矩阵乘法优化递归式

序：

在OI比赛中，很多情况下我们可以能通过打表（找规律）或者某些方式发现一个递归式。

例如：f(n) = f(n - 1)+f(n - 2)，（斐波那契数列）。

通常情况下，我们计算f(n)的时间复杂度就是O(n)（分别计算f(1), f(2) … f(n - 1)）.

但是当n很大又或者还有其他处理的复杂度一叠加便会超时。

如果不学习矩阵乘法优化的话，我们恐怕永远不会想到计算递推式还可以进行优化。

实际上利用矩阵乘法，我们可以将O(n)的计算递归式的复杂度降至O(logn)。

优化递归式的特征

形如f(n) = a1 * f(n - 1) + a2 * f(n - 2) + … + ak * f(n - k)+c (c为常数）

本文讨论的范畴

形如 f(n) = f(n-1) + f(n-2) + .. + f(n-k) 的递归式(1)

形如f(n) = a1 * f(n-1) + a2 * f(n-2) + .. + ak * f(n-k) 的递归式(2)

形如f(n) = a1 * f(n-1) + a2 * f(n-2) + .. + ak * f(n-k) + c 的递归式(3)

实际上理解了最简单的第一个递归式的原理，就很容易理解后面的两种情况。每个前式都是后式的一种特殊情况。

理论基础

首先给出斐波那契数列求第n项的O(logn)的做法，由它引出原理。

已知 :

f(n) = f(n -1) + f(n - 2); – (1)

f(n - 1) = f(n - 1); – (2)

由(1)(2)可以得到这样的一个式子：

[ f(n) ] = [1, 1] * [f(n - 1)]

[f(n - 1)] [1, 0] * [f(n - 2)]

这就核及到矩阵乘法的运算规则。矩阵乘法的计算方式如下所示：

A = X * Y。必须满足row(X) = colom(Y)。

A[i][j] = sum{x[i][k] * x[k][j]}。

row(A) = row(X), colom(A) = colom(Y)..

对于上式，f(n) = f(n - 1) + f(n - 2),f(n - 1) = f(n - 1) + 0，满足斐波那契数列的规则。

我们设右边靠左的式子为A，靠右的为F。

那么计算f(n)，我们只需要计算A^(n - 1) * F。复杂度O(n)。

但是做幂运算，可以通过快速幂将复杂度从O(n)降到O(logn)，因此总复杂度科技降到O(logn)。

通过这个例子我们能发现什么？很显然的便是A与Y（左式与右二式，一下皆简称为A,X,Y）本质上是一个东西，因此通过迭代，直接计算出第n项。

斐波那契数列是一个最简单的例子，它也是(1)的典型例子。

广义斐波那契数列

定义f(n) = a * f(n - 1) + b * f(n - 2).

与之前的区别仅仅在于前面的系数不是1，那么构造出等式只需要照葫芦画瓢即可。

[ f(n) ] = [a, b] [f(n - 1)]

[f(n - 1)] [1, 0] [f(n - 2)],本质没有变化。

那么对于f(n) = a1 * f(n - 1) + …. + ak * f(n - k)，我们只需要构造一个k * k的矩阵。



(*)式便是一般情况的式子，其实只需要使得第一行满足数列公式，其他行满足f(i) = f(i)即可。

因此只需要令f[i][i - 1] = 1(下标从0开始)，其他都是0即可使等式成立。

带常数与系数的递归式

通过之前的经验我们不难看出，矩阵A,Y中的元素一定是每一次计算的时候必需的元素。这次多了一个常数c，因此c也要在A,Y中出现。

对此，我们在需要在最后一行加上c，构造成一个(k + 1) * (k + 1)的矩阵，如下。



观察一下(**)式，第1和k+1行的最后都为1，其他新增的空位全是0，如此便构造完成了。

至于原理，再通过矩阵乘法验证一下不就好了吗？

最后给出代码实现（NOI2012 d1t3）

这道题就是一个裸的(3)式情况。

/*
About: Matrix Fast Exponentiation to calculate Recursion with Coefficient and Constant
From: NOI2012 d1t3
Description:
Input: n, p, c, f(0), Mod, mod;
Recursion: f(n) = (p * f(n-1) + c) % Mod;
Output: f(n) % mod;
Auther: kongse_qi
Date:2017/05/19

Central:
construct a Matrix X[a, 0, 1] and Matrix F[f(1)]. Cal A = X ^ (n-1) * F. The ans is A[0][0] % mod.
[1, 0, 0]             [f(0)]
[0, 0, 1]             [  c ]
[ f(n) ]   [a, 0, 1]   [f(n-1)]
Recursion f(n) = (p * f(n-1) + c) % Mod == [f(n-1)] = [1, 0, 0] * [f(n-2)] (all % mod);
[   c  ]   [0, 0, 1]   [   c  ]
*/

#include <cstdio>
#include <cstring>

#define read(x) scanf("%d", &x)

typedef long long ll;

long long mod;

ll Mul (ll a, ll b)
{
ll ans = 0;
for(; b; b >>= 1, a = (a + a) % mod)
{
if(b & 1) ans = (ans + a) % mod;
}
return ans;
}

struct Matrix{
ll n, m, x[5][5];

Matrix(){}

Matrix(int a, int b):
n(a), m(b){}

Matrix operator * (const Matrix& a)
{
Matrix y(n, a.m);
memset(y.x, 0, sizeof y.x);
for (int i = 0; i < n; i++)
{
for (int k = 0; k < m; k++)
{
if (!x[i][k]) continue;
for (int j = 0; j < a.m; j++)
{
y.x[i][j] = (Mul(x[i][k], a.x[k][j]) + y.x[i][j]) % mod;
}
}
}
return y;
}

Matrix Poww(ll times)
{
Matrix y = *this, z = y;
for(--times; times > 0; times >>= 1, z = z * z)
{
if(times & 1) y = y * z;
}
return y;
}
};

void Create(Matrix& x, Matrix& f, ll a1, ll
be6a
p, ll c)
{
f.n = x.n = x.m = 3;
f.m = 1;
memset(x.x, 0, sizeof x.x);
x.x[0][0] = p;
x.x[0][2] = x.x[1][0] = x.x[2][2] = 1;
f.x[0][0] = (Mul(a1, p) + c) % mod;
f.x[1][0] = a1;
f.x[2][0] = c;
return ;
}

int main()
{
Matrix x, f;
ll a1, n, Mod, c, p;
scanf("%lld%lld%lld%lld%lld%lld", &mod, &p, &c, &a1, &n, &Mod);
Create(x, f, a1, p, c);
Matrix y = x.Poww(n - 1) * f;
printf("%lld\n", y.x[0][0] % Mod);
return 0;
}

至于一些更简单的情况，以下有几道练习。

luogu P1962 斐波那契数列

luogu P1349 广义斐波那契数列

Vijos 1067 Warcraft III 守望者的烦恼

NOI2012 随机数生成器

最后一题值得注意的是，最后几个点乘法过程会炸出longlong，因此需要使用类快速幂的方法快速求和（积），边加边取模。

至此结束。

箜瑟_qi 7:58

2017.05.25