从fibonacci数列开始说起

给定一个nn，求fibonacci数列的第nn项模一个数PP，即F_n \mod PFn​modP。

题目：洛谷1962

较朴素的方法：递推求数

int fibonacci(int n) {
    static int f[N];
    f[1] = 1;
    for (int i = 2; i <= n; ++i)
        f[i] = f[i - 1] + f[i - 2] % P;
    return f[n];
}

这样，我们可以发现其时间复杂度与空间复杂度都是O(n)O(n)，当然，我们可以发现每个数都只需要它之前的两个数，所以我们可以通过这一性质进行空间压缩，将空间复杂度降到O(1)O(1)。这个效率确实不错，但是当n的取值范围提升到10^{12}1012这个数量级时便无能为力了。

事实上，有一个方法可以以O(\log n)O(logn)的效率将F_nFn​算出。

将数列递推式转化为线性变换

我们可以发现，将空间压缩后的每个时刻状态，可以当做一个矢量：

A\_n = \left( \begin{array}{cc} F\_{n-1} & F\_n \end{array} \right)A_n=(F_n−1​F_n​)

在转化成下一个状态的过程中：

A\_{n+1} = \left( \begin{array}{cc} F\_{n} & F\_{n+1} \end{array} \right) = \left( \begin{array}{cc} A\_{n2} & A\_{n1}+A\_{n2} \end{array} \right)A_n+1=(F_n​F_n+1​)=(A_n2​A_n1+A_n2​)

我们便可以将其改写成一个矩阵变换，设

M = \left( \begin{array}{cc} 0 & 1 \\ 1 & 1 \end{array} \right)M=(01​11​)

则有

A_{n+1} = A_n MAn+1​=An​M

根据归纳法，我们可以得知

A_n=A_1 M^{n-1}An​=A1​Mn−1

进而有

F\_n=(A\_1 M^{n-1})\_2=(M^{n-1})\_{2,2}F_n=(A_1Mn−1)_2=(Mn−1)_2,2

那我们要问了，矩阵乘法的复杂度是O(r^3)O(r3)的，这里rr为2即乘方运算大致是8n8n次加法，而原本的求法是nn次加法，这岂不是常数反而更大了？实际上，矩阵乘方和数字的普通乘方运算，都可以通过快速幂来进行加速。那么便可以通过O(\log n)O(logn)的复杂度完成乘方运算。接下来，我将给出证明。

代码样例：

#include <cstdio>
#include <cstring>

using namespace std;

typedef long long ll;

struct matrix {
    int a[2][2];

    matrix()
    { memset(a, 0, sizeof(a)); }

    matrix(const matrix& x)
    { memcpy(a, x.a, sizeof(a)); }

    matrix operator*(const matrix& x) const;
};

const int P = 1e9+7;

int main() {
    matrix ans, m;
    ll n;
    scanf("%lld", &n);
    --n;
    ans.a[0][0] = ans.a[1][1] = 1;
    m.a[0][1] = m.a[1][0] = m.a[1][1] = 1;
    while (n) {
        if (n & 1)
            ans = ans * m;
        m = m * m;
        n >>= 1;
    }
    printf("%d\n", ans.a[1][1]);
    return 0;
}

matrix matrix::operator*(const matrix& x) const {
    matrix ret;
    for (int i = 0; i < 2; ++i)
        for (int j = 0; j < 2; ++j)
            for (int k = 0; k < 2; ++k)
                ret.a[i][j] = (((ll) a[i][k] * x.a[k][j]) + ret.a[i][j]) % P;
    return ret;
}

矩阵乘法的基本定理

矩阵乘法的结合律

设有矩阵A,B,CA,B,C分别的大小为n \times m, m \times p, p \times qn×m,m×p,p×q。求证(AB)C=A(BC)(AB)C=A(BC)，进一步为(AB)C_{i,j}=A(BC)_{i,j}(AB)Ci,j​=A(BC)i,j​

根据矩阵乘法的定义，有

AB\_{i,j} = \sum\_{k=1}^m A\_{i,k} \times B\_{k,j}AB_i,j=∑_k=1mA_i,k×B_k,j

故

$$ ```cpp \begin{eqnarray} (AB)C_{i,j} & = & \sum_{l=1}^p AB_{i,l} \times C_{l,j} \\ & = & \sum_{l=1}^p (\sum_{k=1}^m A_{i,k} \times B_{k,l}) \times C_{l,j} \\ & = & \sum_{l=1}^p \sum_{k=1}^m A_{i,k} \times B_{k,l} \times C_{l,j} \\ & = & \sum_{k=1}^m \sum_{l=1}^p A_{i,k} \times B_{k,l} \times C_{l,j} \\ & = & \sum_{k=1}^m (\sum_{l=1}^p B_{k,l} \times C_{l,j}) \times A_{i,k} \\ & = & \sum_{k=1}^m A_{i,k} \times BC_{k,j} \\ & = & A(BC)_{i,j} \end{eqnarray} ``` $$

证毕。

矩阵快速幂的正确性

设方阵AA(行数为rr)和自然数nn，可以将nn分解为一个mm位的二进制数

\overline{b_mb_{m-1}...b_1b_0}bm​bm−1​...b1​b0​​，则

$$ \begin{eqnarray} A^n ```cpp & = & A^{\sum_{i=0}^m b_i 2^i} \\ & = & \prod_{i=0}^m A^{b_i 2^i} \\ & = & \prod_{b_i = 1} A^{2^i} \end{eqnarray} ``` $$

从第一步到第二步的变换是不平凡的，它只有基于结合律才成立。因此，我们可以发现通过反复平方法可以以O(r^3 \log n)O(r3logn)的效率求出A^nAn。

重定义运算符

矩阵乘法中的运算涉及两种操作，即乘法和加法，我们在对结合律进行证明的过程中可以发现，我们依赖了：两种运算符各自符合交换律、结合律，乘法对加法符合分配率。形式化描述，即有运算\oplus⊕加法与\otimes⊗乘法，只要它们符合以下性质

$$ ```cpp \begin{eqnarray} x \oplus y & = & y \oplus x \\ (x \oplus y) \oplus z & = & x \oplus (y \oplus z) \\ x \otimes y & = & y \otimes x \\ (x \otimes y) \otimes z & = & x \otimes (y \otimes z) \\ x \otimes (y \oplus z) & = & (x \otimes y) \oplus (x \otimes z) \end{eqnarray} ``` $$

那么将该运算替换掉原本矩阵乘法的加法和乘法，交换律依然成立，快速幂也仍然可以使用。比如说，我们知道带余加法和带余乘法仍然满足该性质，所以我们也可以置\times_p \rightarrow \otimes, +_p \rightarrow \oplus×p​→⊗,+p​→⊕。

当然，如果还想要让这种重定义的矩阵可逆等，就需要该矩阵的元素所属集合与两种运算构成域(field)，也就是还要有逆元和单位元。

矩阵快速幂解决最短路问题

我们尝试置+ \rightarrow \otimes, \min \rightarrow \oplus+→⊗,min→⊕后，可以发现这仍然符合矩阵运算所需要的性质，而这时

AB\_{i,j} = \bigoplus\_{k=1}^n A\_{i,k} \otimes B\_{k,j} = \min\_{k=1}^n (A\_{i,k} + B\_{k,j})AB_i,j=⨁_k=1nA_i,k⊗B_k,j=min_k=1n(A_i,k+B_k,j)

（这个minmin的特别用法是我自己编的，嘻嘻）

我们惊奇地发现，这就等同于一次边(i,j)(i,j)的松弛操作！于是我们对图的邻接矩阵AdjAdj计算n-1n−1次幂，就可以得出两两边的最短路了。其中如果(u,v)(u,v)无边，则置Adj_{u,v} \leftarrow +\inftyAdju,v​←+∞。计算Adj^{n-1}Adjn−1的时间复杂度为O(n^3 \log n)O(n3logn)，虽然逊色于floyd算法的O(n^3)O(n3)，但是这也是一个非常启发式的方法。