从fibonacci数列开始说起

给定一个nn,求fibonacci数列的第nn项模一个数PP,即F_n \mod PFnmodP。

题目:洛谷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_n1F_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_nF_n+1)=(A_n2A_n1+A_n2)

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

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

则有

A_{n+1} = A_n MAn+1=AnM

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

A_n=A_1 M^{n-1}An=A1Mn1

进而有

F\_n=(A\_1 M^{n-1})\_2=(M^{n-1})\_{2,2}F_n=(A_1Mn1)_2=(Mn1)_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}bmbm1...b1b0,则

$$ \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,kB_k,j=min_k=1n(A_i,k+B_k,j)

(这个minmin的特别用法是我自己编的,嘻嘻)

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