前言

关于这个算法的前置知识快速幂和矩阵可以点击链接看我以前的博客

问题

给定\(n \times n\)矩阵\(A\),求\(A^k\)

算法思路

顾名思义，矩阵快速幂就是矩阵乘法 + 快速幂
(这里就不再赘述快速幂的原理，不熟悉的可以去看我以前的博客)
要想实现这个算法，我们首先需要先实现矩阵乘法，设：

\[A = (a_{ij})_{n \times m}, B = (b_{ij})_{m \times p} \]

令

\[AB = C = (c_{ij})_{n \times p} \]

对于

\[i = 1, 2, ..., n, j = 1, 2, ..., p \]

有

\[c_{ij} = \sum_{k = 1}^{m}a_{ik}b_{kj} \]

相比于上述描述，矩阵乘法用代码写出来反而更加通俗易懂：

int a[n][m], b[m][p];
int c[n][p] = {0};

for(int k = 1; i <= m; i++)
    for(int i = 1; i <= n; i++)
        for(int j = 1; j <= p; j++)
            c[i][j] += a[i][k] * b[k][j];

对于快速矩阵幂中遇到的\(A = (a_{ij})_{n \times n}\)的情况则更加简单：\(k, i, j\)均小于等于\(n\)。也就是说：

\[AA = C = (c_{ij})_{n \times n} \]

对于

\[i = 1, 2, ..., n, j = 1, 2, ..., n \]

有

\[c_{ij} = \sum_{k = 1}^{n}a_{ik}b_{kj} \]

现在考虑将矩阵乘法运用到快速幂中，那么我们就要实现两个乘法操作：
\(A\)自乘以及将每部分幂的结果乘起来。
前者用\(A \times A\)即可:

int a[n][n];
int c[n][n] = {0};

for(int k = 1; i <= n; i++)
    for(int i = 1; i <= n; i++)
        for(int j = 1; j <= n; j++)
            c[i][j] += a[i][k] * a[k][j];

后者我们可以使用单位矩阵\(I \times A\)，所以我们还需要定义单位矩阵。
由单位矩阵的定义我们可以写出代码:

int I[n][n] = {0};
for(int i = 1; i <= n; i++) I[i][i] = 1;

然后实现\(I \times A\):

int a[n][n];
int I[n][n] = {0};
int c[n][n] = {0};

for(int i = 1; i <= n; i++) I[i][i] = 1;

for(int k = 1; i <= n; i++)
    for(int i = 1; i <= n; i++)
        for(int j = 1; j <= n; j++)
            c[i][j] += ans[i][k] * a[k][j];

将上述操作添加到快速幂中:在\(n\)(\(n\)指代要求的次数)不等于\(0\)时每次让\(A\)自乘,并且当\(n\)的最后一位二进制是\(1\)时将\(A \times I\);上述操作结束后去掉\(n\)的最后一位二进制位。

代码实现

这里以洛谷P3390 【模板】矩阵快速幂为例

#include <iostream>

#define ll long long

using namespace std;

const int N = 1e9 + 7;

ll a[105][105];
ll ans[105][105] = {0};    //将ans[][]作为保存结果的矩阵
ll n, b;   //b即为次数

void work1() {      //A自乘
    ll c[105][105] = {0};  //用来保存这次运算的结果
    for(int k = 1; k <= n; k++)
        for(int i = 1; i <= n; i++)
            for(int j = 1; j <= n; j++)
                c[i][j] = (c[i][j] + a[i][k] * a[k][j]) % N;
    for(int i = 1; i <= n; i++) 
                for(int j = 1; j <= n; j++) 
                        a[i][j] = c[i][j];
}

void work2() {      //保存结果
    ll c[105][105] = {0};  //用来保存这次运算的结果
    for(int k = 1; k <= n; k++)
        for(int i = 1; i <= n; i++)
            for(int j = 1; j <= n; j++)
                c[i][j] = (c[i][j] + ans[i][k] * a[k][j]) % N;
    for(int i = 1; i <= n; i++) 
        for(int j = 1; j <= n; j++) 
            ans[i][j] = c[i][j];
}

void matrixQuickPow() {
    while(b != 0) {
        if((b & 1) == 1)
            work2();
        work1();
        b >>= 1;
    }
}

int main(void) {
    cin >> n >> b;
    for(int i = 1; i <= n; i++) {
        for(int j = 1; j <= n; j++) {
            cin >> a[i][j];
        }
    }
    for(int i = 1; i <= n; i++) ans[i][i] = 1;  //将ans初始化为单位矩阵
    
    matrixQuickPow();

    for(int i = 1; i <= n; i++) {
        for(int j = 1; j <= n; j++) {
            cout << ans[i][j] % N << " ";
        }
        cout << endl;
    }

    return 0;
}

526互联

矩阵快速幂

前言

问题

算法思路

代码实现