B.矩阵快速幂
等比数列二分求和

#include <iostream>
#include<cstdio>
#include<cstring>
using namespace std;
int n,k,mod;
const int MAXN=35;
struct Mat{
int n,m;
int mat[MAXN][MAXN];
Mat(){
memset(mat,0,sizeof(mat));
n=MAXN,m=MAXN;
}
Mat(int x,int y){
memset(mat,0,sizeof(mat));
n=x,m=y;
};
Mat operator*(Mat b){
Mat c(n,b.m);
for(int i=0;i<n;i++)
for(int j=0;j<b.m;j++){
for(int k=0;k<m;k++){
c.mat[i][j]+=(mat[i][k]*b.mat[k][j])%mod;
c.mat[i][j]%=mod;
}
}
return c;
}
Mat operator+(Mat b){
Mat c(n,m);
for(int i=0;i<n;i++)
for(int j=0;j<m;j++)
c.mat[i][j]=(mat[i][j]+b.mat[i][j])%mod;
return c;
}
void in(){
for(int i=0;i<n;i++)
for(int j=0;j<m;j++)
scanf("%d",&mat[i][j]);
}
void out(){
for(int i=0;i<n;i++){
for(int j=0;j<m;j++)
printf("%d%c",mat[i][j],j==m-1?'\n':' ');
}
}
}E;
Mat qpow(Mat a,int b){Mat ret=E; while(b){ if(b&1) ret=ret*a; a=a*a; b>>=1;} return ret;}
Mat cal(Mat a,int k){//二分求等比数列和--a+a^2+a^3+...+a^k
if(k==1) return a;
else if(k&1){
Mat cur=qpow(a,k/2+1);
return cal(a,k/2)*(E+cur)+cur;
}
else return cal(a,k/2)*(E+qpow(a,k/2));
}
int main()
{
scanf("%d%d%d",&n,&k,&mod);
E.n=n,E.m=n;
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
E.mat[i][j]=(i==j?1:0);
Mat now(n,n);
now.in();
now=cal(now,k);
now.out();
return 0;
}
#include <iostream>
#include<cstdi
