bzoj 4162: shlw loves matrix II 拉格朗日插值法+矩阵乘法

题意

给定矩阵 M，请计算 M^n，并将其中每一个元素对 1000000007 取模输出。
对于 100% 数据，满足 n <= 2^10000;k <= 50; 0 <= Mij < 10^9 +7

分析

我们可以带入k+1个不同的$\lambda$算出所有的$|A-\lambda I|$，然后用拉格朗日插值法把M的特征多项式$f(x)$插出来。
因为根据Cayley-Hamilton定理有$f(M)=0$，所以$x^n=x^n\bmod f(x)$，用快速幂算出来模意义下的值后把A带进去算即可。

代码

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
using namespace std;

typedef long long LL;

const int N=55;
const int MOD=1000000007;

int n,a[N][N],b[N*2],c[N][N],d[N][N],r[N*2],tmp[N*2],f[N],g[N],mat[N][N],ans[N][N];
pair<int,int> pts[N];
char str[10005];

int read()
{
    int x=0,f=1;char ch=getchar();
    while (ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=getchar();}
    while (ch>='0'&&ch<='9'){x=x*10+ch-'0';ch=getchar();}
    return x*f;
}

int ksm(int x,int y)
{
    int ans=1;
    while (y)
    {
        if (y&1) ans=(LL)ans*x%MOD;
        x=(LL)x*x%MOD;y>>=1;
    }
    return ans;
}

int det()
{
    int ans=1;
    for (int i=1;i<=n;i++)
    {
        int l=i;
        for (int j=i;j<=n;j++) if (mat[j][i]) {l=j;break;}
        if (l!=i)
        {
            for (int j=i;j<=n;j++) swap(mat[i][j],mat[l][j]);
            ans=-ans;
        }
        for (int j=i+1;j<=n;j++)
            if (mat[j][i])
            {
                int w=(LL)mat[j][i]*ksm(mat[i][i],MOD-2)%MOD;
                for (int k=i;k<=n;k++) (mat[j][k]-=(LL)mat[i][k]*w%MOD)%=MOD;
            }
        ans=(LL)ans*mat[i][i]%MOD;
    }
    ans+=ans<0?MOD:0;
    return ans;
}

void mul(int *a,int *b,int *c)
{
    for (int i=0;i<=n*2-2;i++) tmp[i]=0;
    for (int i=0;i<n;i++)
        for (int j=0;j<n;j++)
            (tmp[i+j]+=(LL)a[i]*b[j]%MOD)%=MOD;
    int ny=ksm(f[n],MOD-2);
    for (int i=n*2-2;i>=n;i--)
        for (int j=n-1;j>=0;j--)
            (tmp[i-n+j]-=(LL)f[j]*tmp[i]%MOD*ny%MOD)%=MOD;
    for (int i=0;i<n;i++) c[i]=tmp[i];
}

int main()
{
    scanf("%s",str);n=read();
    for (int i=1;i<=n;i++)
        for (int j=1;j<=n;j++)
            a[i][j]=read();
    for (int i=0;i<=n;i++)
    {
        memcpy(mat,a,sizeof(a));
        for (int j=1;j<=n;j++) mat[j][j]-=i;
        pts[i]=make_pair(i,det());
    }
    for (int i=0;i<=n;i++)
    {
        for (int j=0;j<=n;j++) g[j]=0;
        g[0]=pts[i].second;
        for (int j=0;j<=n;j++) if (i!=j) g[0]=(LL)g[0]*ksm(pts[j].first-pts[i].first,MOD-2)%MOD;
        for (int j=0;j<=n;j++)
            if (i!=j)
            {
                for (int k=n;k>=1;k--) g[k]=((LL)g[k]*pts[j].first%MOD-g[k-1])%MOD;
                g[0]=(LL)g[0]*pts[j].first%MOD;
            }
        for (int j=0;j<=n;j++) (f[j]+=g[j])%=MOD;
    }
    r[0]=1;b[1]=1;
    int len=strlen(str);
    for (int i=len-1;i>=0;i--)
    {
        if (str[i]=='1') mul(r,b,r);
        mul(b,b,b);
    }
    for (int i=1;i<=n;i++) c[i][i]=1;
    for (int i=0;i<n;i++)
    {
        for (int j=1;j<=n;j++)
            for (int k=1;k<=n;k++)
                (ans[j][k]+=(LL)c[j][k]*r[i]%MOD)%=MOD;
        memset(d,0,sizeof(d));
        for (int j=1;j<=n;j++)
            for (int k=1;k<=n;k++)
                for (int l=1;l<=n;l++)
                    (d[j][k]+=(LL)c[j][l]*a[l][k]%MOD)%=MOD;
        memcpy(c,d,sizeof(d));
    }
    for (int i=1;i<=n;i++)
    {
        for (int j=1;j<n;j++) printf("%d ",(ans[i][j]+MOD)%MOD);
        printf("%d\n",(ans[i][n]+MOD)%MOD);
    }
    return 0;
}