POJ 3233 Matrix Power Series(矩阵快速幂+二分求和)

代码部分参考于HIT《数论及其应用》

#include <cstdio>
using namespace std;
const int maxn=100+5;
typedef struct
{
    int m[maxn][maxn];
}Matrix;
Matrix a,per;
int n,M;

void init()
{
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++){
        scanf("%d",&a.m[i][j]);
        a.m[i][j]%=M;
        per.m[i][j]=(i==j);
    }
}
Matrix add(Matrix a,Matrix b)
{
    Matrix z;
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            z.m[i][j]=(a.m[i][j]+b.m[i][j])%M;
    return z;
}
Matrix Multi(Matrix a,Matrix b)
{
    Matrix z;
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++){
        z.m[i][j]=0;
        for(int k=0;k<n;k++)
            z.m[i][j]+=a.m[i][k]*b.m[k][j];
        z.m[i][j]%=M;
    }
    return z;
}
Matrix power(int k)
{
    Matrix r=per,p=a;
    while(k)
    {
        if(k&1) r=Multi(r,p);
        p=Multi(p,p);
        k>>=1;
    }
    return r;
}
Matrix MatrixSum(int k)///二分求和(等比数列的形式)
{
    if(k==1) return a;
    Matrix tmp,b;
    tmp=MatrixSum(k/2);
    if(k&1)
    {
        b=power(k/2+1);
        tmp=add(tmp,Multi(tmp,b));
        tmp=add(tmp,b);
    }
    else
    {
        b=power(k/2);
        tmp=add(tmp,Multi(tmp,b));
    }
    return tmp;
}
void output(Matrix a)
{
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            printf("%d%c",a.m[i][j],j==n-1?'\n':' ');
}
int main()
{
    int k;
    while(~scanf("%d%d%d",&n,&k,&M))
    {
        init();
        Matrix ans=MatrixSum(k);
        output(ans);
    }
    return 0;
}