POJ3233：Matrix Power Series

对n<=30（其实可以100）大小的矩阵A求A^1+A^2+……+A^K，K<=1e9，A中的数%m。

从K的二进制位入手。K分解二进制，比如10110，令F[i]=A^1+A^2+……+A^(2^i)，那么答案就是F[10000]*A^110+F[100]*A^10+F[10]+A^0。也就是说如果知道F就可以得答案。

F亦可递推，F[i]=F[i-1]*(A^(2^i-1)+A^0)。完美！什么log方，都是假的！

#include<stdio.h>
#include<string.h>
#include<algorithm>
#include<math.h>
//#include<iostream>
using namespace std;

int n,K,mod;
#define maxn 111
#define LL long long
typedef LL mat[maxn][maxn];
mat f,a,last,tmp,t,ans,E;
void init(mat &a)
{
    memset(a,0,sizeof(a));
    for (int i=1;i<=n;i++) a[i][i]=1;
}
void add(mat a,mat b,mat &ans)
{
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
            ans[i][j]=(a[i][j]+b[i][j])%mod;
}
void mul(mat a,mat b,mat &ans)
{
    mat t;
    memset(t,0,sizeof(t));
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
            for (int k=1;k<=n;k++)
                t[i][j]=(t[i][j]+a[i][k]*b[k][j]%mod)%mod;
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
            ans[i][j]=t[i][j];
}
int main()
{
    scanf("%d%d%d",&n,&K,&mod);
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
            scanf("%lld",&a[i][j]),tmp[i][j]=f[i][j]=a[i][j];
    init(last);init(E);
    memset(ans,0,sizeof(ans));
    while (K)
    {
        if (K&1)
        {
            mul(f,last,t);
            add(ans,t,ans);
            mul(last,tmp,last);
        }
        add(tmp,E,t);
        mul(f,t,f);
        mul(tmp,tmp,tmp);
        K>>=1;
    }
    for (int i=1;i<=n;i++)
    {
        for (int j=1;j<=n;j++)
            printf("%lld ",ans[i][j]);
        puts("");
    }
    return 0;
}