poj3233

题目大意：给定矩阵A，求A + A^2 + A^3 + … + A^k的结果（两个矩阵相加就是对应位置分别相加）。输出的数据mod m。k<=10^9。
    这道题两次二分，相当经典。首先我们知道，A^i可以二分求出。然后我们需要对整个题目的数据规模k进行二分。比如，当k=6时，有：
    A + A^2 + A^3 + A^4 + A^5 + A^6 =(A + A^2 + A^3) + A^3*(A + A^2 + A^3)

对A^3用矩阵快速幂，而(A + A^2 + A^3)可以递归求出，即可一步步推得答案

ps:对角线上为1，其他地方为0的矩阵是单位矩阵（性质：任何矩阵与其相乘均得到本身）

还有，这里矩阵相乘和矩阵相加是不一样的！A,B两矩阵相乘则C(i,j)=∑A[i,k]*B(k,j)；相加则C(i,j)=A(i,j)+B(i,j)

#include<cstdio>
#include<cctype>
#include<iostream>
using namespace std;
int n,m,k;
struct mat{long long m[32][32];}ans,ond,a;

inline void read(int &x){
    char ch=getchar();x=0;
    while(!isdigit(ch))ch=getchar();
    while(isdigit(ch)){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
}

inline void read(long long &x){
    char ch=getchar();x=0;
    while(!isdigit(ch))ch=getchar();
    while(isdigit(ch)){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
}

inline mat mul(mat x,mat y)
{  
    mat c;  
    for(int i=1;i<=n;i++)  
       for(int j=1;j<=n;j++)  
        {  
            c.m[i][j]=0;  
            for(k=1;k<=n;k++)   
                c.m[i][j]=(x.m[i][k]*y.m[k][j]+c.m[i][j])%m;  
        }  
    return c;  
}  

inline mat add(mat x,mat y){
    for(int i=1;i<=n;i++)
       for(int j=1;j<=n;j++)
          x.m[i][j]=(x.m[i][j]+y.m[i][j])%m;
    return x;
}

inline mat powmul(mat x,int k){
    mat b=ond,c=a;
    while(k>0){
        if(k&1)b=mul(b,c);
        c=mul(c,c);
        k>>=1;
    }
    return b;
}

inline mat work(mat now,int k){
    if(k==1)return now;
    if(k&1){
        mat b=work(now,k>>1);
        mat c=powmul(now,(k+1)>>1);
        return add(c,mul(b,add(c,ond)));
    }else{
        mat b=work(now,k>>1);
        return mul(add(powmul(now,k>>1),ond),b);
    }
}

int main(){
    read(n);read(k);read(m);
    for(int i=1;i<=n;i++){ond.m[i][i]=1;for(int j=1;j<=n;j++)read(a.m[i][j]);}
    ans=work(a,k);
    for(int i=1;i<=n;i++)
       for(int j=1;j<=n;j++)
          if(j!=n)printf("%lld ",ans.m[i][j]);else printf("%lld\n",ans.m[i][j]);
}