poj 3233 Matrix Power Series 矩阵快速幂 二分

计算 s(k)= a+a^2+a^3+.....a^k

a^k 用矩阵的快速幂

因为k比较大，不能直接循环相加，因吃应用二分的方法

s(k)=s(k/2)+s(k/2)*a^k/2   (k为奇，偶略有不同)

#include<iostream>
#include<cstdio>
using namespace std;
const int maxn=40;
int n,m,k;
typedef struct
{
    int m[maxn][maxn];
} Matrax;
Matrax a,per;
Matrax multi(Matrax a,Matrax b,int m) //矩阵乘法，(a*b)%m
{
    Matrax 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]+=a.m[i][k]*b.m[k][j];
            c.m[i][j]%=m;
        }
    }
    return c;
}
Matrax power(Matrax a,int k,int m) //矩阵快速幂，(a^k)%m
{
    Matrax p,ans;
    p=a;
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++) ans.m[i][j]=(i==j);//初始化ans为单位矩阵
    while(k)
    {
        if(k&1)
        {
            ans=multi(ans,p,m);
            k--;
        }
        k/=2;
        p=multi(p,p,m);
    }
    return ans;
}
Matrax addm(Matrax a,Matrax b,int m) //矩阵加法 (a+b)%m
{
    Matrax 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])%m;
    return c;
}
void cc(Matrax a )
{
    for(int i=0;i<n;i++)
    {
        cout<<a.m[i][0];
        for(int j=1;j<n;j++) cout<<" "<<a.m[i][j];
        cout<<endl;
    }

}
Matrax erfen(Matrax a,int k,int m)
{
    if(k==1) return a;
    Matrax temp,ans,p;
    temp=erfen(a,k/2,m);
    if(k%2)
    {
        p=power(a,k/2+1,m);
        ans=addm(temp,multi(temp,p,m),m);
        ans=addm(ans,p,m);
    }
    else
    {
        p=power(a,k/2,m);
        ans=addm(temp,multi(temp,p,m),m);
        //cc(ans);
    }
    return ans;
}
int main()
{
    while(cin>>n>>k>>m)
    {
        Matrax a;
        for(int i=0;i<n;i++)
        {
            for(int j=0;j<n;j++) 
            {
                cin>>a.m[i][j];
                a.m[i][j]%=m;
            }
        }
        a=erfen(a,k,m);
        for(int i=0;i<n;i++)
        {
            cout<<a.m[i][0];
            for(int j=1;j<n;j++) cout<<" "<<a.m[i][j];
            cout<<endl;
        }
    }
}