POJ 3233 Matrix Power Series （矩阵快速幂+二分求解）

#include <iostream>
#include <cstdio>
#include <string.h>

using namespace std;
const int maxn=31;
int mod;
int n,k,m;
struct Matrix{
    int m[maxn][maxn];
    void init(){
        memset(m,0,sizeof(m));
    }
    void initE(){
        memset(m,0,sizeof(m));
        for(int i=0;i<n;i++)
            m[i][i]=1;
    }
}A;
//重载+运算符
Matrix operator+(Matrix a,Matrix b){
    Matrix 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;
}
//重载*运算符
Matrix operator*(Matrix a,Matrix b){
    Matrix 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]=(c.m[i][j]+a.m[i][k]*b.m[k][j]%mod)%mod;
            }
        }
    }
    return c;
}
//矩阵快速幂
Matrix MquickPow(Matrix A,int b){
    Matrix ret;
    ret.initE();
    while(b){
        if(b&1)
            ret=ret*A;
        A=A*A;
        b=b>>1;
    }
    return ret;
}
Matrix p;
Matrix dfs(Matrix A,int k){
    if(k==1){
        p=A;
        return A;
    }
    Matrix ret,ans;
    ret=dfs(A,k/2);
    //Matrix p=MquickPow(A,k/2);如果用快速幂计算p=A^(k/2)，则要1500ms，而直接在递归的时候同时计算p，则只要188ms。
    if(k&1){
        //return ret+ret*p+p*p*A;
        ans=ret+ret*p+p*p*A;
        p=p*p*A;
    }
    else{
        //return ret+ret*p;
        ans=ret+ret*p;
        p=p*p;
    }
    return ans;
}
int main()
{
    scanf("%d%d%d",&n,&k,&m);
    mod=m;
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            scanf("%d",&A.m[i][j]);
        }
    }
    Matrix ans;
    ans=dfs(A,k);
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++){
            printf("%d ",ans.m[i][j]);
        }
        printf("\n");
    }
    return 0;
}