G3 矩阵快速幂 • EH's Blog

题目情境#

题目链接：https://www.luogu.com.cn/problem/P3390 ↗

题目描述

给定 $n\times n$ 的矩阵 $A$ ，求 $A^k$ 。

输入格式

第一行两个整数 $n,k$ 。
接下来 $n$ 行，每行 $n$ 个整数，第 $i$ 行的第 $j$ 的数表示 $A_{i,j}$ 。

输出格式

输出 $A^k$

共 $n$ 行，每行 $n$ 个数，第 $i$ 行第 $j$ 个数表示 $(A^k)_{i,j}$ ，每个元素对 $10^9+7$ 取模。

输入 #1

2 1
1 1
1 1

plaintext

输出 #1

1 1
1 1

plaintext

输入 #2

plaintext

输出 #2

121824 149688 177552
275886 338985 402084
429948 528282 626616

plaintext

对于 $100\%$ 的数据， $1\le n \le 100$ ， $0 \le k \le 10^{12}$ ， $|A_{i,j}| \le 1000$ 。

算法解析：#

矩阵快速幂：把标量快速幂推广到矩阵。初始为单位矩阵，转移矩阵自乘平方，指数 $k$ 按二进制拆分，当前位为 1 则乘上对应矩阵。矩阵乘法三重循环 $O(n^3)$ ，总复杂度 $O(n^3\log k)$ ，模 $10^9+7$ 。

Python代码实现#

import sys
MOD = 10**9 + 7
def mat_mult(X, Y, n):  # 矩阵乘法：返回 X * Y
    Z = [[0] * n for _ in range(n)]
    for i in range(n):
        for k in range(n):
            if X[i][k]:  # 跳过零元，加快运算
                xik = X[i][k]
                for j in range(n):
                    Z[i][j] = (Z[i][j] + xik * Y[k][j]) % MOD
    return Z

def mat_pow(A, exp, n):  # 快速幂：A^exp
    # 初始化为单位矩阵
    res = [[1 if i == j else 0 for j in range(n)] for i in range(n)]
    base = A
    while exp:
        if exp & 1:
            res = mat_mult(res, base, n)
        base = mat_mult(base, base, n)
        exp >>= 1
    return res
n, k = map(int, input().split())
# 读取矩阵，转换为 0-index 列表
A = []
for i in range(n):
    A.append(list(map(int, input().split())))
# 计算 A^k 模 MOD
result = mat_pow(A, k, n)
# 输出结果矩阵
out = []
for row in result:
    out.append(' '.join(str(v) for v in row))
sys.stdout.write('\n'.join(out))

python

C++代码实现#

#include <iostream>
#include <cstring>
#include <algorithm>
using namespace std;

typedef long long LL;
const int mod=1000000007;
struct matrix{
    LL c[101][101];
    matrix(){memset(c, 0, sizeof c);}
} A, res;
LL n, k;

matrix operator*(matrix &x, matrix &y){ //矩阵乘法
    matrix t; //临时矩阵
    for(int i=1; i<=n; i++)
      for(int j=1; j<=n; j++)
        for(int k=1; k<=n; k++)
          t.c[i][j]=(t.c[i][j]+x.c[i][k]*y.c[k][j])%mod;
    return t;
}
void quickpow(LL k){ //快速幂
    for(int i=1; i<=n; i++) res.c[i][i]=1; //单位矩阵
    while(k){
        if(k & 1) res = res*A;
        A = A*A;
        k >>= 1;
    }
}
int main(){
    scanf("%d%lld",&n,&k);
    for(int i=1; i<=n; i++)
        for(int j=1; j<=n; j++)
            scanf("%d",&A.c[i][j]);
    quickpow(k);
    for(int i=1; i<=n; i++){
        for(int j=1; j<=n; j++)
            printf("%d ",res.c[i][j]);
        puts("");
    }
    return 0;
}