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

数学 算法 组合数学

---

主页地址 ：月光森林

---

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

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

另外要说明的一点是，这种方法仅适用于常系数线性递推关系，即递推关系

$$h_n = a_1h_{n-1} + a_2h_{n-2} + \cdots + a_kh_{n-k} + b$$

中，$a_1, a_2, \cdots, a_k$ 以及$b$ 均为常数(可以为0)的递推关系，且$h_n, h_{n-1}, h_{n-2}, \cdots, h_{n - k}$ 的指数均为1。

这里我们先考虑递推关系为齐次递推关系的情况，即忽略常数项b(令 $b = 0$ ), 由线性代数的知识可以知道，在递推关系

$$h_n = a_1h_{n-1} + a_2h_{n-2} + \cdots + a_kh_{n-k}$$

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

$$\begin{bmatrix} h_{s+k} \\ h_{s+k-1} \\h_{s+k-2} \\ \vdots \\ h_{s+2} \\ h_{s + 1} \\ \end{bmatrix} = \ A \begin{bmatrix} h_{s+k-1} \\ h_{s+k-2} \\ h_{s+k-3} \\ \vdots \\ h_{s+1} \\h_s \\ \end{bmatrix}$$

其中， A矩阵为：

$$\begin{bmatrix} a_1 & a_2 & a_ 3 & \cdots & a_{k-1} & a_k \\ 1 & 0 & 0 & \cdots & 0 & 0 \\ 0 & 1 & 0 & \cdots & 0 & 0 \\ 0 & 0 & 1 & \cdots & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 & 0 \\ \end{bmatrix}$$

可以看出， $A$ 矩阵是一个$k + 1$ 阶方阵，且第一行由递推关系等号右边的k个元素$h_{n-1}, h_{n-2}, h_{n-3}, \cdots, h_{n-k}$ 系数构成。从第2行到第k行为行矩阵$e_i$ ——第 i 列为1、其余元素均为 0 的行矩阵，并且 $i \in [1, k-1]$。同样，这个$(k-1) \times k$ 矩阵也可以看作是一个k阶单位矩阵去掉最后一行(即第 k 行)得到的。

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

$$\begin{equation} \begin{aligned} \begin{bmatrix} h_{s+k} \\ h_{s+k-1} \\h_{s+k-2} \\ \vdots \\ h_{s+2}\\ h_{s + 1} \\ \end{bmatrix} = \ A & \begin{bmatrix} h_{s+k-1} \\ h_{s+k-2} \\ h_{s+k-3} \\ \vdots \\ h_{s+1} \\h_s \\ \end{bmatrix} \\ = \ A^2 & \begin{bmatrix} h_{s+k-2} \\ h_{s+k-3} \\ h_{s+k-4} \\ \vdots \\h_s \\ h_{s-1} \\ \end{bmatrix} \\ \vdots \\ = \ A^s & \begin{bmatrix} h_k \\ h_{k-1} \\ h_{k-2} \\ \vdots \\h_2 \\ h_{1} \\ \end{bmatrix} \\ \end{aligned} \end{equation}$$

于是可以得到等式

$$\begin{bmatrix} h_n \\ h_{n-1} \\h_{n-2} \\ \vdots \\ h_{n-k+2}\\ h_{n-k+1} \\ \end{bmatrix} = \ A^{n-k} \begin{bmatrix} h_k \\ h_{k-1} \\ h_{k-2} \\ \vdots \\h_2 \\ h_{1} \\ \end{bmatrix}$$

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

$$h_n = a_1h_{n-1} + a_2h_{n-2} + \cdots + a_kh_{n-k} + b$$

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

$$\begin{bmatrix} h_{s+k} \\ h_{s+k-1} \\h_{s+k-2} \\ \vdots \\ h_{s+2} \\ h_{s + 1} \\ b \\ \end{bmatrix} = \ A \begin{bmatrix} h_{s+k-1} \\ h_{s+k-2} \\ h_{s+k-3} \\ \vdots \\ h_{s+1} \\h_s \\ b \\ \end{bmatrix}$$

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

$$\begin{bmatrix} a_1 & a_2 & a_ 3 & \cdots & a_{k-1} & a_k & 1 \\ 1 & 0 & 0 & \cdots & 0 & 0 & 0\\ 0 & 1 & 0 & \cdots & 0 & 0 & 0\\ 0 & 0 & 1 & \cdots & 0 & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots\\ 0 & 0 & 0 & \cdots & 1 & 0 & 0 \\ 0 & 0 & 0 & \cdots & 0 & 0 & 1 \\ \end{bmatrix}$$

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

通过迭代，就可以得到：

$$\begin{bmatrix} h_n \\ h_{n-1} \\h_{n-2} \\ \vdots \\ h_{n-k+2}\\ h_{n-k+1} \\ b \\ \end{bmatrix} = \ A^{n-k} \begin{bmatrix} h_k \\ h_{k-1} \\ h_{k-2} \\ \vdots \\h_2 \\ h_{1} \\ b \\ \end{bmatrix}$$

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

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

$$\left\{ \begin{aligned} h_n =& \ 2h_{n-1} - h_{n-2} + 3 \\ h_1 = & \ 1, h_2 = 1 \end{aligned} \right.$$

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

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

$$\begin{bmatrix} 2 & -1 & 1 \\ 1 & 0 & 0 \\ 0 & 0 & 1 \\ \end{bmatrix}$$

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

#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即可，万变不离其宗。

以上です。