POJ 3233 二分求等比数列 矩阵快速幂

#include <cstdio>
#include <iostream>
#include <cstring>
#define sf scanf
#define pf printf
using namespace std;
const int maxn = 35;
int n;
__int64 mod;
struct Mat{
__int64 C[maxn][maxn];
Mat(){memset(C,0,sizeof C);for(int i = 0;i < n;++i) C[i][i] = 1;}
void out(){
for(int i = 0;i < n;++i){
pf("%I64d",C[i][0]);
for(int j = 1;j < n;++j) pf(" %I64d",C[i][j]);
pf("\n");
}
}
};
Mat operator*(const Mat& A,const Mat& B){
Mat ret;
for(int i = 0;i < n;++i){
for(int j = 0;j < n;++j){
ret.C[i][j] = 0;
for(int k = 0;k < n;++k){
ret.C[i][j] += A.C[i][k]*B.C[k][j];
ret.C[i][j] %= mod;
}
}
}
return ret;
}
Mat operator+(const Mat& A,const Mat& B){
Mat ret;
for(int i = 0;i < n;++i)
for(int j = 0;j < n;++j)
{ret.C[i][j] = A.C[i][j] + B.C[i][j];ret.C[i][j] %= mod;}
return ret;
}
Mat operator^(Mat A,int k){
Mat ret;
while(k){
if(k & 1) ret = ret * A;
A = A * A;
k = k >> 1;
}
return ret;
}

Mat getSum(const Mat A,int k){
if(k == 1){
return A;
}
else{
Mat ret = getSum(A,k / 2);
if(k & 1){
Mat tmp = A ^ (k / 2 + 1);
return ret + (ret * tmp) + tmp;
}
else{
return ret + ret * (A ^ (k / 2));
}
}
}

int main(){
int k;
while( ~sf("%d%d%I64d",&n,&k,&mod)){
Mat st;
for(int i = 0;i < n;++i)
for(int j = 0;j < n;++j) {sf("%I64d",&st.C[i][j]);st.C[i][j] %= mod;}
getSum(st,k).out();
}
return 0;

}