矩阵快速幂

矩阵快速幂能将 $O(n)$ 的线性递推优化成 $O(\log n)$.

矩阵相当于二维数组.

• 矩阵 $A$ 有 $m$ 行 $n$ 列，称为 $m×n$ 矩阵，简记为 $A_{mn}$.

• 矩阵 $A$ 第 $i$ 行 $j$ 列的元素写作 $a_{ij}$.

$$A=\begin{bmatrix} a_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ \vdots & \vdots & \ddots & \vdots \\ a_{m1}& a_{m2}& \cdots & a_{mn} \end{bmatrix}$$

主对角线上的元素都为 $1$，其余元素为 $0$ 的 $n×n$ 矩阵称为 $n$ 阶单位矩阵，记作 $I_n$ 或 $E_n$.

$$I_n=\begin{bmatrix} 1 & 0 & 0 & \cdots & 0 \\ 0 & 1 & 0 & \cdots & 0 \\ 0 & 0 & 1 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & 1 \end{bmatrix}$$

矩阵的加减法就是将两个矩阵对应位置上的元素相加减.

$$C_{ij}=A_{ij}±B_{ij}$$

定义矩阵 $A×B=C$：

• 必须满足 $A$ 的列数 $=$ $B$ 的行数.

• $A_{mr}×B_{rn}$ 会得到一个 $m×n$ 矩阵.

$$C_{ij}=\sum_{k=1}^nA_{ik}\cdot B_{kj}$$

任意矩阵乘 $I_n$ 都等于它本身.

• 矩阵乘法满足结合律，不满足交换律.

• 只有行数 $=$ 列数的矩阵才能进行乘幂运算.

基于 快速幂 的思想求 $n$ 阶矩阵的 $x$ 次方. 时间复杂度：$O(n^3\log x)$.

矩阵运算可能需要高精度或取模.

const int N = 101; // 矩阵的行数和列数

struct matrix {
    int at[N][N];
    matrix() { memset(at, 0, sizeof at); }
    int* operator[](int key) { return at[key]; }
};

matrix operator *(matrix a, matrix b) { // 矩阵 A × 矩阵 B
    matrix c;
    for (int i = 0; i < N; i ++)
        for (int j = 0; j < N; j ++)
            for (int k = 0; k < N; k ++)
                c[i][j] += a[i][k] * b[k][j];
    return c;
}

matrix operator ^(matrix a, int x) {    // 矩阵 A 的 x 次方
    matrix res;
    for (int i = 0; i < N; i ++)        // 初始化为单位矩阵
        res[i][i] = 1;
    while(x) {
        if(x % 2) res = res * a;
        a = a * a, x /= 2;
    }
    return res;
}

若已知 $f(1),f(2)$，且有矩阵 $A$ 满足

$$A\times \begin{bmatrix} f(n)\\ f(n-1) \end{bmatrix} = \begin{bmatrix} f(n+1)\\ f(n) \end{bmatrix}$$

就可推出 $f(n)$ 的值.

$$(A^{n-2})\times \begin{bmatrix} f(2)\\ f(1) \end{bmatrix} = \begin{bmatrix} f(n)\\ f(n-1) \end{bmatrix}$$

可用矩阵快速幂求 $A^{n-2}$. 时间复杂度：$O(\log n)$.

求斐波那契数列的第 $n$ 项 $f(n)$（$n\geq 10^{12}$）.

$$f(n)= \begin{cases} 1&n=1,2\\ f(n-1)+f(n-2)&n\geq 3 \end{cases}$$

设矩阵 $A=\begin{bmatrix}a_{11}&a_{12}\newline a_{21}&a_{22}\end{bmatrix}$ 满足

$$\begin{bmatrix} a_{11}&a_{12}\\ a_{21}&a_{22} \end{bmatrix} \times \begin{bmatrix} f(n)\\ f(n-1) \end{bmatrix} = \begin{bmatrix} f(n+1)\\ f(n) \end{bmatrix}$$

即

$$\begin{cases} f(n+1)&=&a_{11}\cdot f(n)+a_{12}\cdot f(n-1)\\ f(n)&=&a_{21}\cdot f(n)+a_{22}\cdot f(n-1) \end{cases}$$

解得

$$\begin{cases} a_{11}=1&a_{12}=1\\ a_{21}=1&a_{22}=0 \end{cases} , \quad A= \begin{bmatrix} 1&1\\ 1&0 \end{bmatrix}$$

因此

$$\begin{bmatrix} 1&1\\ 1&0 \end{bmatrix} ^{n-2} \times \begin{bmatrix} f(2)\\ f(1) \end{bmatrix} = \begin{bmatrix} f(n)\\ f(n-1) \end{bmatrix}$$

将 $f(1)=f(2)=1$ 代入即可.

int main() {
    matrix A;
    A[1][1] = 1, A[1][2] = 1;
    A[2][1] = 1, A[2][2] = 0;

    matrix B;
    B[1][1]	= 1;	// f(1) = 1
    B[2][1] = 1;	// f(2) = 1

    cin >> n;
    B = (A ^ (n - 2)) * B;
    cout << B[1][1];
}

求数列的第 $n$ 项 $f(n)$（$n\geq 10^{12}$）.

$$f(n)= \begin{cases} 1&n=1,2\\ f(n-1)+f(n-2)+n-1&n\geq 3 \end{cases}$$

设 $4$ 阶矩阵 $A$ 满足

$$\begin{bmatrix} a_{11}&a_{12}&a_{13}&a_{14}\\ a_{21}&a_{22}&a_{23}&a_{24}\\ a_{31}&a_{32}&a_{33}&a_{34}\\ a_{41}&a_{42}&a_{43}&a_{44} \end{bmatrix} \times \begin{bmatrix} f(n)\\ f(n-1)\\ n\\ 1 \end{bmatrix} = \begin{bmatrix} f(n+1)\\ f(n)\\ n+1\\ 1 \end{bmatrix}$$

解得

$$A= \begin{bmatrix} 1&1&1&0\\ 1&0&0&0\\ 0&0&1&1\\ 0&0&0&1 \end{bmatrix}$$