【矩阵快速幂+二分】Matrix Power Series POJ - 3233

#include <cstdio>
#include <cstring>
#include <algorithm>

using namespace std;

const int N = 40;

typedef struct Matrax{
int m

;
}matrax;

matrax e, per;
int n, m, mod;

void read();/*读取矩阵+初始化单位矩阵*/
matrax Add(matrax a, matrax b);/*矩阵加法*/
matrax multi(matrax a, matrax b);/*矩阵乘法*/
matrax power(int k);/*矩阵快速幂*/
matrax matrax_sum(int k);/*等比矩阵前n项和*/
void Pri(matrax ans);/*输出矩阵*/

int main(){
while(~scanf("%d %d %d", &n, &m, &mod)){
read();
matrax ans = matrax_sum(m);
Pri(ans);
}
return 0;
}
void read(){/*读取矩阵+初始化单位矩阵*/
for(int i = 0; i < n; i++){
for(int j = 0; j < n; j++){
scanf("%d", &e.m[i][j]);
e.m[i][j] %= mod;
per.m[i][j] = (i == j);
}
}
}
matrax Add(matrax a, matrax b){/*矩阵加法*/
matrax c;
for(int i = 0; i < n; i++){
for(int j = 0; j < n; j++){
c.m[i][j] = (a.m[i][j] + b.m[i][j])%mod;
}
}
return c;
}
matrax multi(matrax a, matrax b){/*矩阵乘法*/
matrax c;
for(int i = 0; i < n; i++){
for(int j = 0; j < n; j++){
c.m[i][j] = 0;
for(int k = 0; k < n; k++){
c.m[i][j] += (a.m[i][k] * b.m[k][j])%mod;
}
c.m[i][j] %= mod;
}
}
return c;
}
matrax power(int k){/*矩阵快速幂*/
matrax p, ans = per;
p = e;
while(k){
if(k & 1)
ans = multi(ans, p);
p = multi(p, p);
k >>= 1;
}
return ans;
}
matrax matrax_sum(int k){/*等比矩阵前n项和*/
if(k == 1)
return e;
matrax tmp, b;
tmp = matrax_sum(k>>1);
if(k & 1){
b = power((k>>1)+1);
tmp = Add(tmp, multi(b, tmp));
tmp = Add(tmp, b);
}
else {
b = power(k>>1);
tmp = Add(tmp, multi(b, tmp));
}
return tmp;
}
void Pri(matrax ans){/*输出矩阵*/
for(int i = 0; i < n; i++){
for(int j = 0; j < n; j++){
printf("%d%c", ans.m[i][j], j == n-1? '\n': ' ');
}
}
}