题目
给定一个 $n imes n$ 的矩阵 $A$ 和正整数 $k$ 和 $m$。求矩阵 $A$ 的幂的和。
$$S = A + A^2 + ... + A^k$$
输出 $S$ 的各个元素对 $M$ 取余后的结果（$1 leq n leq 30, 1 leq k leq 10^9, 1 leq M leq 10^4$）。
分析
数据范围 $n$ 很小，$k$ 很大，不肯能逐一求得。
由于具有等比性质，
设 $S_k = I + A + ... + A^{k-1}$
则有
$$egin{pmatrix} A^k\ S_k end{pmatrix} = egin{pmatrix} A & 0\ I & I end{pmatrix} egin{pmatrix} A^{k-1}\ S_{k-1} end{pmatrix} = egin{pmatrix} A & 0\ I & I end{pmatrix}^kegin{pmatrix} I\ 0 end{pmatrix}$$
对这个新矩阵使用快速幂即可。
代码中的输出，使用了分块矩阵乘法的性质进行了简化。
#include<cstdio> #include<cstring> using namespace std; typedef long long ll; struct matrix { int r, c; int mat[61][61]; matrix(){ memset(mat, 0, sizeof(mat)); } }; int n, k, p; matrix mul(matrix A, matrix B) //矩阵相乘 { matrix ret; ret.r = A.r; ret.c = B.c; for(int i = 0;i < A.r;i++) for(int k = 0;k < A.c;k++) for(int j = 0;j < B.c;j++) { ret.mat[i][j] = (ret.mat[i][j] + A.mat[i][k] * B.mat[k][j]) % p; } return ret; } matrix mpow(matrix A, int n) { matrix ret; ret.r = A.r; ret.c = A.c; for(int i = 0;i < ret.r;i++) ret.mat[i][i] = 1; while(n) { if(n & 1) ret = mul(ret, A); A = mul(A, A); n >>= 1; } return ret; } int main() { scanf("%d%d%d", &n, &k, &p); matrix a, b; a.r = a.c = n; for(int i = 0;i < n;i++) for(int j = 0;j < n;j++) scanf("%d", &a.mat[i][j]); b.r = b.c = 2*n; for(int i = 0;i < n;i++) { for(int j = 0;j < n;j++) b.mat[i][j] = a.mat[i][j]; b.mat[n+i][i] = b.mat[n+i][n+i] = 1; } b = mpow(b, k+1); for(int i = 0;i < n;i++) for(int j = 0;j < n;j++) { int tmp = b.mat[n+i][j] % p; if(i == j) tmp = (tmp + p - 1) % p; printf("%d%c", tmp, j == n-1 ? ' ' : ' '); } }
还有一种经典的分治方法，
例如，
$A+A^2+A^3+A^4 = (A+A^2) + A^2(A + A^2)$，
$A+A^2+A^3+A^4+A^5 = (A+A^2) +A^3 + A^3(A + A^2)$.
因此，分k的奇偶递归一下就可以了。
#include<cstdio> #include<cstring> using namespace std; typedef long long ll; struct matrix { int r, c; int mat[61][61]; matrix(){ memset(mat, 0, sizeof(mat)); } }; int n, k, p; matrix mul(matrix A, matrix B) //矩阵相乘 { matrix ret; ret.r = A.r; ret.c = B.c; for(int i = 0;i < A.r;i++) for(int k = 0;k < A.c;k++) for(int j = 0;j < B.c;j++) { ret.mat[i][j] = (ret.mat[i][j] + A.mat[i][k] * B.mat[k][j]) % p; } return ret; } matrix mpow(matrix A, int n) { matrix ret; ret.r = A.r; ret.c = A.c; for(int i = 0;i < ret.r;i++) ret.mat[i][i] = 1; while(n) { if(n & 1) ret = mul(ret, A); A = mul(A, A); n >>= 1; } return ret; } matrix add(matrix A, matrix B) { matrix ret; ret.r = A.r; ret.c = A.c; for(int i = 0;i < A.r;i++) for(int j = 0;j < A.c;j++) ret.mat[i][j] = (A.mat[i][j] + B.mat[i][j]) % p; return ret; } matrix sum(matrix x, int k) //A+A^2+..+A^k { if(k == 1) return x; matrix s = sum(x, k/2); if(k & 1) { matrix tmp = mpow(x, k/2+1); return add(s, add(tmp, mul(tmp, s))); } else { matrix tmp = mpow(x, k/2); return add(s, mul(tmp, s)); } } int main() { scanf("%d%d%d", &n, &k, &p); matrix a, ans; a.r = a.c = n; for(int i = 0;i < n;i++) for(int j = 0;j < n;j++) scanf("%d", &a.mat[i][j]); ans = sum(a, k); for(int i = 0;i < n;i++) for(int j = 0;j < n;j++) printf("%d%c", ans.mat[i][j], j == n-1 ? ' ' : ' '); }
这个时间复杂度咋算啊？知道的大犇请留言。
参考链接：
1. https://blog.csdn.net/rowanhaoa/article/details/21023599