POJ 3233 Matrix Power Series 矩阵快速幂

if(n%2==0)s(n)=s(n/2)+s(n/2)*A^(n/2)

else s(n)=s((n-1)/2)+s((n-1)/2)*A^((n-1)/2) + A^n。

//#pragma comment(linker, "/STACK:1024000000,1024000000")
#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<algorithm>
#include<iostream>
#include<sstream>
#include<cmath>
#include<climits>
#include<string>
#include<map>
#include<queue>
#include<vector>
#include<stack>
#include<set>
using namespace std;
typedef long long ll;
typedef pair<int,int> pii;
#define pb(a) push_back(a)
#define INF 0x1f1f1f1f
#define lson idx<<1,l,mid
#define rson idx<<1|1,mid+1,r
#define PI  3.1415926535898
template<class T> T min(const T& a,const T& b,const T& c) {
    return min(min(a,b),min(a,c));
}
template<class T> T max(const T& a,const T& b,const T& c) {
    return max(max(a,b),max(a,c));
}
void debug() {
#ifdef ONLINE_JUDGE
#else

    freopen("d:\\in.txt","r",stdin);
   // freopen("d:\\out1.txt","w",stdout);
#endif
}
int getch() {
    int ch;
    while((ch=getchar())!=EOF) {
        if(ch!=' '&&ch!='\n')return ch;
    }
    return EOF;
}

int n,m;
struct Matrix
{
    int da[30][30];
    Matrix operator * (const Matrix &ans)
    {
        Matrix res;
        for(int i=0;i<n;i++)
        {
            for(int j=0;j<n;j++)
            {
                res.da[i][j]=0;
                for(int k=0;k<n;k++)
                    res.da[i][j]=(res.da[i][j]+da[i][k]*ans.da[k][j])%m;
            }
        }
        return res;
    }
    Matrix operator +(const Matrix &ans)
    {
        Matrix res;
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
                res.da[i][j]=(da[i][j]+ans.da[i][j])%m;
        return res;
    }
};
Matrix base;
Matrix f(int n)
{
    if(n==1)return base;
    if(n%2==0)
    {
        Matrix res=f(n/2);
        return res*res;
    }else
    {
        Matrix res=f(n/2);
        return  res*res*base;
    }
}
Matrix s(int k)
{
    if(k==1)
        return base;
    if(k%2==0)
    {
        Matrix res;
        res=s(k/2);
        return res+res*f(k/2);
    }else
    {
        Matrix res;
        res=s(k/2);
        return res+res*f(k/2)+f(k);
    }
}
int main()
{
    int k;
    while(scanf("%d%d%d",&n,&k,&m)!=EOF)
    {
        for(int i=0;i<n;i++)
            for(int j=0;j<n;j++)
            {
                scanf("%d",&base.da[i][j]);
                base.da[i][j]%=m;
            }
        Matrix x=f(3);
        Matrix res=s(k);
        for(int i=0;i<n;i++)
        {
            for(int j=0;j<n;j++)
                printf("%d%c",res.da[i][j],j+1==n?'\n':' ');
        }
    }
    return 0;
}