矩阵快速幂在常系数线性递推关系中的应用

先引入一下，知乎上有一个问题 关于斐波拉契数列的一个低级问题 。题主询问了关于求解斐波拉契数列第n项对10007取模的结果。而这个n，可以达到106甚至109 。

解法已经在排名第一的回答中给出了，主要思路就是快速幂和矩阵乘法的结合律，亦即矩阵快速幂。具体方法这里也就不再给出。但可以依托此思想，拓展出在O(logn) 的时间下计算一个递推关系的第n项。

另外要说明的一点是，这种方法仅适用于常系数线性递推关系，即递推关系hn=a1hn−1+a2hn−2+⋯+akhn−k+b

中，a1,a2,⋯,ak 以及b 均为常数(可以为0)的递推关系，且hn,hn−1,hn−2,⋯,hn−k 的指数均为1。

这里我们先考虑递推关系为齐次递推关系的情况，即忽略常数项b(令 b=0 ), 由线性代数的知识可以知道，在递推关系hn=a1hn−1+a2hn−2+⋯+akhn−k

中，如果令 n=s+k, 则原递推关系可以转换为矩阵问题：

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢hs+khs+k−1hs+k−2⋮hs+2hs+1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥= A⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢hs+k−1hs+k−2hs+k−3⋮hs+1hs⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥

其中， A矩阵为：

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢a1100⋮0a2010⋮0a3001⋮0⋯⋯⋯⋯⋱⋯ak−1000⋮1ak0000⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

可以看出， A 矩阵是一个k+1 阶方阵，且第一行由递推关系等号右边的k个元素hn−1,hn−2,hn−3,⋯,hn−k 系数构成。从第2行到第k行为行矩阵ei ——第 i 列为1、其余元素均为 0 的行矩阵，并且 i∈[1,k−1]。同样，这个(k−1)×k 矩阵也可以看作是一个k阶单位矩阵去掉最后一行(即第 k 行)得到的。

接下来就要运用到矩阵乘法的结合律 来迭代上述矩阵等式：

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢hs+khs+k−1hs+k−2⋮hs+2hs+1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥= A= A2⋮= As⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢hs+k−1hs+k−2hs+k−3⋮hs+1hs⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢hs+k−2hs+k−3hs+k−4⋮hshs−1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢hkhk−1hk−2⋮h2h1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥

于是可以得到等式

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢hnhn−1hn−2⋮hn−k+2hn−k+1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥= An−k⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢hkhk−1hk−2⋮h2h1⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥

现在来考虑非齐次情况的递推关系，即：

hn=a1hn−1+a2hn−2+⋯+akhn−k+b

中b≠0 。实际上读者们已经可以推出，只需要将原来的初始矩阵改变一下即可，为了方便，仍然令n=s+k, 递推关系可以转化为矩阵问题：

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢hs+khs+k−1hs+k−2⋮hs+2hs+1b⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥= A⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢hs+k−1hs+k−2hs+k−3⋮hs+1hsb⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

此时A矩阵变成了 k+1 阶方阵：

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢a1100⋮00a2010⋮00a3001⋮00⋯⋯⋯⋯⋱⋯⋯ak−1000⋮10ak000⋮001000⋮01⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

这里的第一行仍然是由递推关系等号右边的k个元素hn−1,hn−2,hn−3,⋯,hn−k 的系数以及一个1构成。从第2行到第k+1行为行矩阵ei ——第 i 列为1、其余元素均为 0 的行矩阵，并且 i∈[1,k−1]∪{k+1}。类似地，这个k×(k+1) 矩阵也可以看作是一个(k+1)阶单位矩阵去掉倒数第二行(即第 k 行)得到的。

通过迭代，就可以得到：

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢hnhn−1hn−2⋮hn−k+2hn−k+1b⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥= An−k⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢hkhk−1hk−2⋮h2h1b⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

这样，便可以通过矩阵快速幂算法 快速得到递推关系hn 的值。做一次矩阵运算的复杂度是O(k3)， 快速幂的复杂度为O(logn) 。由此可以得到求解 hn 的时间复杂度是 O(k3logn)。

例 已知一个常系数线性递推关系

{hn=h1= 2hn−1−hn−2+3 1,h2=1

求它的第n项，n的范围是[0, 1e9]。因为数值可能过大，所以输出得数对10007取模即可。

由上面的推理我们可以得到这个常系数线性递推关系的A 矩阵为

⎡⎣⎢210−100101⎤⎦⎥

由此可以得到求解hn 的代码：

#include <iostream>
#include <cstdio>
#include <cstring>

using namespace std;

const maxn(5), mod(10007);
struct matrix{
int c, r;
int mat[maxn][maxn];

matrix(){}

matrix(int _r, int _c)
:r(_r), c(_c) {memset(mat, 0, sizeof(mat));}

matrix operator * (const matrix &t)const
{
matrix ans(r, t.c);
if(c != t.r) return ans;
for(int i = 0; i < r; i++)
for(int j = 0; j < t.c; j++)
{
int &temp = ans.mat[i][j];
for(int k = 0; k < c; k++)
temp += mat[i][k] * t.mat[k][j];
}
return ans;
}
}A, h;

void init()
{
int _ma[maxn][maxn] =
{
2, -1, 1, 0, 0,
1, 0, 0, 0, 0,
0, 0, 1, 0, 0,
0, 0, 0, 0, 0,
0, 0, 0, 0, 0
};
memcpy(A.mat, _ma, sizeof(_ma));
A.c = A.r = 3;
h.mat[0][0] = h.mat[1][0] = 1;
h.mat[2][0] = 3;
h.r = 3, h.c = 1;
}

matrix pow(matrix t, int n)
{
matrix ans(3, 3);
for(int i = 0; i < maxn; i++)
ans.mat[i][i] = 1;
while(n)
{
if(n & 1) ans = ans * t;
t = t * t;
for(int i = 0; i < t.r; i++)
for(int k = 0;  k < t.c; k++)
t.mat[i][k] %= mod;
n >>= 1;
}
return ans;
}

int main()
{
int n, ca = 1;
while(~scanf("%d", &n) && n)
{
if(n == 1 || n == 2)
{
printf("1\n");
continue;
}
init();
matrix ans = pow(A, n - 2);
ans = ans * h;
printf("Case %d: %d\n",ca++, ans.mat[0][0]);
}
return 0;
}

至于矩阵快速幂的编写方法就和上面给出的例子中类似，只需要改变matrix结构体的r和c即可，万变不离其宗。

以上です。