POJ3233:Matrix Power Series(矩阵快速幂+递推式)

传送门

题意

给出n,m,k,求

\[\sum_{i=1}^kA^i \]

A是矩阵

分析

我们首先会想到等比公式，然后得到这样一个式子：

\[\frac{A^{k+1}-E}{A-E} \]

发现要用矩阵除法，可以用求矩阵逆来做，现在我们换一种做法，我们发现有这样一个性质：

\[\left[ \begin{matrix} A & E \\ 0 & E \\ \end{matrix} \right]^n= \left[ \begin{matrix} A^{n} & \sum_{i=0}^{n-1}A^i \\ 0 & E \\ \end{matrix} \right] \]

那么我们将原先矩阵扩大成四倍，对矩阵求k+1次幂，然后取右上角减去一个单位阵即可

trick

代码

#include<cstdio>
#include<cstring>
#include<iostream>
using namespace std;
struct Matrix
{
    int matrix[80][80];
}ans,ret;
int t,n,k,mod;
Matrix multi(Matrix a,Matrix b)
{
    Matrix tmp;
    for(int i=1;i<=n;++i)for(int j=1;j<=n;++j)
    {
        tmp.matrix[i][j]=0;
        for(int k=1;k<=n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
    }
    for(int i=1;i<=n;++i)for(int j=1;j<=n;++j)
    {
        for(int k=n+1;k<=2*n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
    }
     for(int i=1;i<=n;++i)for(int j=n+1;j<=2*n;++j)
    {
        tmp.matrix[i][j]=0;
        for(int k=1;k<=n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
    }
    for(int i=1;i<=n;++i)for(int j=n+1;j<=2*n;++j)
    {
        for(int k=n+1;k<=2*n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
    }
    for(int i=1+n;i<=2*n;++i)for(int j=1;j<=n;++j)
    {
        tmp.matrix[i][j]=0;
        for(int k=1;k<=n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
    }
    for(int i=1+n;i<=2*n;++i)for(int j=1;j<=n;++j)
    {
        for(int k=n+1;k<=2*n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
    }
    for(int i=n+1;i<=2*n;++i)for(int j=n+1;j<=2*n;++j)
    {
        tmp.matrix[i][j]=0;
        for(int k=1;k<=n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
    }
    for(int i=1+n;i<=2*n;++i)for(int j=n+1;j<=2*n;++j)
    {
        for(int k=n+1;k<=2*n;++k) (tmp.matrix[i][j]+=a.matrix[i][k]*b.matrix[k][j])%=mod;
    }
    return tmp;
}
void fast_mod(int p)
{
    memset(ans.matrix,0,sizeof(ans.matrix));
    for(int i=1;i<=2*n;++i) ans.matrix[i][i]=1;
    for(;p;p>>=1,ret=multi(ret,ret)) if(p&1) ans=multi(ans,ret);
}


int main()
{
    //freopen("data.in","w",stdout);
    while(scanf("%d %d %d",&n,&k,&mod)!=EOF)
    {
        for(int i=1;i<=n;++i)for(int j=1;j<=n;++j) scanf("%d",&ret.matrix[i][j]);
        for(int i=1;i<=n;++i)for(int j=n+1;j<=2*n;++j) if((j-i)==n) ret.matrix[i][j]=1;else ret.matrix[i][j]=0;
        for(int i=n+1;i<=2*n;++i)for(int j=1+n;j<=2*n;++j) if(i==j) ret.matrix[i][j]=1;else ret.matrix[i][j]=0;
        for(int i=n+1;i<=2*n;++i)for(int j=1;j<=n;++j)   ret.matrix[i][j]=0;
        //for(int i=1;i<=2*n;++i)for(int j=1;j<=2*n;++j) printf("%d%c",ret.matrix[i][j],j==2*n?'\n':' ');
        fast_mod(k+1);
        //for(int i=1;i<=2*n;++i)for(int j=1;j<=2*n;++j) printf("%d%c",ans.matrix[i][j],j==2*n?'\n':' ');
        for(int i=1;i<=n;++i)for(int j=n+1;j<=2*n;++j) if((j-i)==n) (ans.matrix[i][j]+=(mod-1))%=mod;
        for(int i=1;i<=n;++i)for(int j=1+n;j<=2*n;++j) printf("%d%c",ans.matrix[i][j],j==2*n?'\n':' ');
    }
    return 0;
}