poj 3233 Matrix Power Series---矩阵快速幂
要求矩阵A的k次幂，矩阵快速幂加上二分求和
其中，矩阵相乘二分：A^2k=A^k*A^k,
A^(2k+1)=A^k*A^k*A.
求和二分：A+A^2+A...+A^(2k+1)= A+A^2+...+A^k+A^(k+1)+A^(k+1)*(A+A^2+...+A^k).
A+A^2+...+A^2k = A+A^2+...+A^k+A^k*(A+A^2+...+A^k).
#include <cstdlib>#include <cstring>#include <cstdio>#include <iostream>using namespace std;int N,k,m;struct matrix{ int a[35][35];}origin,res;matrix multiply(matrix x,matrix y){ matrix temp; memset(temp.a,0,sizeof(temp.a)); for(int i=0;i<N;i++) { for(int j=0;j<N;j++) { for(int k=0;k<N;k++) { temp.a[i][j]+=(x.a[i][k]*y.a[k][j]); } temp.a[i][j]%=m; } } return temp;}/*matrix calc(int n){ while(n) { if(n&1) { res=multiply(res,origin); n--; } else { n>>=1; origin=multiply(origin,origin); } } return multiply(res,origin);}*/matrix calc(int exp){ matrix p,q; p=origin; q=res; while (exp!=1) { if (exp&1) { exp--; q=multiply(p,q); } else { exp>>=1; p=multiply(p,p); } } return multiply(p,q);}matrix add(matrix x,matrix y){ matrix tmp; int i,j; for(i=0;i<N;i++) for(j=0;j<N;j++) { tmp.a[i][j]=x.a[i][j]+y.a[i][j]; tmp.a[i][j]%=m; } return tmp;}void init(){ int i,j; memset(res.a,0,sizeof res.a); for(i=0;i<N;i++) for(j=0;j<N;j++) { scanf("%d",&origin.a[i][j]); origin.a[i][j]%=m; } for(i=0;i<N;i++) { res.a[i][i]=1; }}matrix sum(int n){ matrix tmp,now; if(n==1) return origin; tmp=sum(n/2); if(n&1) { now=calc(n/2+1); tmp=add(tmp,multiply(tmp,now)); tmp=add(tmp,now); } else { now=calc(n/2); tmp=add(tmp,multiply(tmp,now)); } return tmp;}int main(){ int i,j; while(~scanf("%d%d%d",&N,&k,&m)) { init(); matrix tmp=sum(k); for(i=0;i<N;i++) { printf("%d",tmp.a[i][0]); for(j=1;j<N;j++) { printf(" %d",tmp.a[i][j]); } putchar('\n'); } } return 0;}