一、内容
https://blog.csdn.net/qq_41280600/article/details/103073544
二、思路
- 我们规定f(A, k) 这个函数为求A + A^2^ + A^3^ ..... A^k^ 的值
- 当k为偶数时:
S = A + A^2^ + A^3^ ..... A^k^ = (A + A^1^....A^k/2^) + (A^k/2+1^+.....A^k^) = (A + A^1^....A^k/2^) +A^k/2^(A + A^1^....A^k/2^) = (1 + A^k/2^) * (A + A^1^....A^k/2^) = (1 + A^k/2^) * f(A, k / 2)。 - 同理当k为奇数时:
S = (A + A^1^....A^(k+1)/2^) +A^(k+1)/2^(A + A^1^....A^(k+1)/2-1^) = (1 + A^(k+1)/2^) * (A + A^1^....A^(k-1)/2^) = (1 + A^(k+1)/2^) * f(A, (k-1)/ 2) - 最后利用矩阵快速幂求解即可。
三、代码
#include <cstdio>
#include <cstring>
typedef long long ll;
using namespace std;
const int N = 35;
int n, k, mod;
struct Node {
int m[N][N];
};
Node ONE; //单位矩阵
Node mul(Node a, Node b) {
Node ans;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
ans.m[i][j] = 0; //注意初始化。
for (int k = 0; k < n; k++) {
ans.m[i][j] = (ans.m[i][j] + (ll)a.m[i][k] * b.m[k][j]) % mod;
}
}
}
return ans;
}
Node add(Node a, Node b) {
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
a.m[i][j] = (a.m[i][j] + (ll)b.m[i][j]) % mod;
}
}
return a;
}
Node qpow(Node a, int b) {
Node ans = ONE;
while (b) {
if (b & 1) ans = mul(ans, a);
a = mul(a, a);
b >>= 1;
}
return ans;
}
Node f(Node A, int k) {
if (k <= 1) return A;
if (k & 1) {
//奇数情况下
Node B = qpow(A, (k + 1) / 2);
Node C = f(A, (k - 1) / 2);
return add(mul(add(ONE, B), C), B);
} else {
Node B = f(A, k / 2);
return mul(add(qpow(A, k / 2), ONE), B);
}
}
int main() {
scanf("%d%d%d", &n, &k, &mod);
Node rec;
for (int i = 0; i < n; i++) {
ONE.m[i][i] = 1;
for (int j = 0; j < n; j++) {
scanf("%d", &rec.m[i][j]);
}
}
Node ans = f(rec, k);
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
printf("%d ", ans.m[i][j]);
}
puts("");
}
return 0;
}