#include <iostream>
#include <cstring>
#include <cstdio>
using namespace std;
#define maxn 35
int n, k, m;//A是n*n阶的矩阵,k是要求的sum = A^1 + A^2 +.....A^k; m 是被模的元素
struct Mat {
int val[maxn][maxn];
void unit() { //单位矩阵
for (int i = 0; i < maxn; i++)
val[i][i] = 1;
}
void zero() { //零矩阵
memset(val, 0, sizeof(val));
}
} x;
Mat operator *(const Mat &a, const Mat &b) { //矩阵乘法
Mat tmp;
tmp.zero();
for (int k = 1; k <= n; k++) {
for (int i = 1; i <= n; i++)
if (a.val[i][k])
for (int j = 1; j <= n; j++) {
tmp.val[i][j] += a.val[i][k] * b.val[k][j];
if (tmp.val[i][j] >= m)
tmp.val[i][j] %= m;
}
}
return tmp;
}
Mat operator ^(Mat x, int n) { //矩阵快速幂
Mat tmp;
tmp.zero();
tmp.unit();
while (n) {
if (n & 1)
tmp = tmp * x;
x = x * x;
n >>= 1;
}
return tmp;
}
Mat operator +(const Mat &a, const Mat &b) { //矩阵加法
Mat tmp;
for (int i = 1; i <= n; i++)
for (int j = 1; j <= n; j++)
tmp.val[i][j] = (a.val[i][j] + b.val[i][j]) % m;
return tmp;
}
Mat sum(int k) {//求sum S = A + A2 + A3 + … + Ak.
if (k == 1)
return x;
else {
Mat tmp = sum(k >> 1);
if (k & 1) {
Mat tmp2 = x ^ ((k >> 1) + 1);
return tmp + tmp2 + tmp * tmp2;
} else {
Mat tmp2 = x ^ (k >> 1);
return tmp + tmp * tmp2;
}
}
}
int main(){
while(scanf("%d%d%d",&n,&k,&m)!=EOF){
int i,j;
for(i = 1 ; i <= n ; ++i){
for(j = 1 ; j <= n ; ++j){
scanf("%d",&x.val[i][j]);
x.val[i][j] %= m;
}
}
Mat ans = sum(k);
for(i = 1 ; i <= n ; ++i){
for(j = 1 ; j < n ; ++j){
printf("%d ",ans.val[i][j]);
}
printf("%d\n",ans.val[i][n]);
}
}
return 0;
}
(Relax 矩阵快速幂 1.2)POJ 3233 Matrix Power Series(用矩阵加法+矩阵快速幂来求sum= A + A2 + A3 + … + Ak)
最新推荐文章于 2022-05-08 15:24:05 发布