矩阵快速幂，矩阵加法，矩阵乘法

转载自大神脑钊文：http://blog.csdn.net/wzw1376124061/article/details/52492665

#include<iostream>  
#include<cstdio>  
#include<cstdlib>  
#include<cmath>  
#include<cstring>  
#include<algorithm>  
using namespace std;  
int n,k;  
int mod; //n*n阶的矩阵;k是sum=A^1+A^2+.....+A^k;m是被模的元素    
struct Matrix{  
    int matrix[50][50];  
    Matrix(int a=0){  
        memset(matrix,0,sizeof(matrix)); //矩阵清零   
        if(a==1) //单位矩阵  
            for(int i=0;i<50;i++)  
                matrix[i][i]=1;  
    }  
}m;  
Matrix operator * (Matrix a,Matrix b){ //矩阵乘法  
    Matrix res;  
    memset(res.matrix,0,sizeof(res.matrix));  
    for(int i=1;i<=n;i++)  
        for(int j=1;j<=n;j++)  
            for(int k=1;k<=n;k++)  
                res.matrix[i][j]=(res.matrix[i][j]+(a.matrix[i][k]*b.matrix[k][j])%mod)%mod;  
    return res;  
}  
Matrix operator + (Matrix a,Matrix b){ //矩阵加法  
    Matrix res;  
    for(int i=1;i<=n;i++)  
        for(int j=1;j<=n;j++)  
            res.matrix[i][j]=(a.matrix[i][j]+b.matrix[i][j])%mod;   
    return res;  
}  
Matrix operator ^ (Matrix a,int k){ //矩阵快速幂  
    bool flag=false;  
    Matrix ans=1;  
    while(k){  
        if(k&1){  
            if(flag)ans=ans*a;  
            else ans=a;  
            flag=true;  
        }  
        a=a*a;  
        k>>=1;  
    }  
    return ans;  
}  
Matrix sum(int k){ //求sum( S(k) )= A + A2 + A3 + …+ Ak  
    if(k==1)return m; //递归回去   
    else{  
        Matrix tmp=sum(k>>1); //sum( S(k/2) )  
        if(k&1){ //k为奇数   
            Matrix tmp2=m^((k>>1)+1);   
            return tmp+tmp2+tmp*tmp2;  
        }  
        else{  
            Matrix tmp2=m^(k>>1);  
            return tmp+tmp*tmp2; //sum( S(k/2) )*(tmp2+1)  
        }  
    }  
}  
int main(){  
    int i,j;  
    Matrix ans;  
    scanf("%d%d%d",&n,&k,&mod);  
    for(i=1;i<=n;i++)  
        for(j=1;j<=n;j++){  
            scanf("%d",&m.matrix[i][j]);  
            m.matrix[i][j]%=mod;  
        }  
    ans=sum(k);  
    for(i=1;i<=n;i++,puts(""))  
        for(j=1;j<=n;j++)  
            printf("%d ",ans.matrix[i][j]);  
    return 0;  
}