【POJ3233】Matrix Power Series

题面

给出矩阵A，求S = A + A² + A³ + … + A^k.

分析

矩阵的乘方是可以通过快速幂很快的推出来，主要是相加的问题

但是别忘了，虽然没有等比求和，但是矩阵是满足结合律的

因此求出了A + A² + A³ + … + A^k/2后只需要乘A^k/2，就可以得到A^(k/2)+1 + A^(k/2)+2 + A^(k/2)+3 + … + A^k

类似地，一直二分下去。（上面忽略了k的奇偶，写的时候要注意）

代码

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<iostream>
using namespace std;
struct email
{
    int x[33][33];
}a,o;
int n,m,k;
email mul(email a,email b)
{
    email c;
    memset(c.x,0,sizeof(c.x));
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            for(int k=1;k<=n;k++)
                c.x[i][j]=(a.x[i][k]*b.x[k][j]+c.x[i][j])%m;
    return c;
}
email add(email a,email b)
{
    int i,j;
    for(i=1;i<=n;i++)
        for(j=1;j<=n;j++)
            a.x[i][j]+=b.x[i][j],a.x[i][j]%=m; 
    return a;
}

email ksm(email a,int k)
{
    if(k==1)return a;
      email b=ksm(a,k/2);
    b=mul(b,b);
    if(k&1)b=mul(b,a);
    return b;
}
email divide(email a,int k)
{
    if(k==1)return a;
    if(k%2==0)
    {
        email b=divide(a,k/2);
        return mul(add(ksm(a,k/2),o),b);
    }
    else
    {
        email b=divide(a,(k)/2);
        return add(mul(add(ksm(a,k/2),o),b),ksm(a,k));
    }
}
int main()
{
    scanf("%d%d%d",&n,&k,&m);
    for(int i=1;i<=n;i++)o.x[i][i]=1;
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n;j++)
            scanf("%d",&a.x[i][j]);
    email ans=divide(a,k);
    for(int i=1;i<=n;i++)
    {
        for(int j=1;j<=n;j++)
            printf("%d ",ans.x[i][j]);
        printf("
");
    }
    return 0;
}