[bzoj4162]shlw loves matrix II

来自FallDream的博客，未经允许，请勿转载，谢谢

---

给定矩阵k*k的矩阵M，请计算 M^n，并将其中每一个元素对 1000000007 取模输出。 k<=50 n<=2^10000

考虑求出矩阵的特征多项式，这点我们可以通过带入$lambda=x0..xk$，求出矩阵的行列式，然后通过插值求出多项式。

然后搬出一个很厉害的定理:f(A)=0 其中f(x)是特征多项式，A是矩阵，所以我们可以把所求的$x^{n}$对f(x)取膜，从而让答案变成一个k-1次多项式。然后我们把原矩阵带进去就行了。

插值的复杂度是$O(n^{4})$,然后后面那部分的复杂度是$k^{2}logn$ 

然后做了这道题好像懂了怎么在O(klogklogn)内做完"常系数线性递推"，貌似还要nlogn的多项式取余 真麻烦TAT

#include<iostream>
#include<cstdio>
#include<cstring>
#define ll long long
#define mod 1000000007
using namespace std;
inline int read()
{
    int x = 0; char ch = getchar();
    while(ch < '0' || ch > '9')ch = getchar();
    while(ch >= '0' && ch <= '9'){x = x * 10 + ch - '0';ch = getchar();}
    return x;
}

char st[10001];
int y[51],k;

int pow(int x,int k)
{
    int sum=1;
    for(;k;k>>=1,x=1LL*x*x%mod)
        if(k&1) sum=1LL*sum*x%mod;
    return sum;
}
struct Matrix
{
    int s[51][51];
    Matrix operator*(Matrix b)
    {
        Matrix c;memset(c.s,0,sizeof(c.s));
        for(int i=1;i<=k;i++)
            for(int l=1;l<=k;l++) if(s[i][l])
                for(int j=1;j<=k;j++)
                    c.s[i][j]=(c.s[i][j]+1LL*s[i][l]*b.s[l][j])%mod;
        return c;
    }
    int calc()
    {
        int sum=1,j;
        for(int i=1;i<=k;i++)
        {
            for(j=i;j<=k;j++)
                if(s[j][i])
                {
                    if(j!=i)
                    {
                        sum=mod-sum;
                        for(int l=1;l<=k;l++)
                            swap(s[j][l],s[i][l]);
                    }
                    break;
                }
            if(j==k+1) return 0;
            for(int j=i+1;j<=k;j++)
            {
                int inv=1LL*pow(s[i][i],mod-2)*s[j][i]%mod;
                for(int l=i;l<=k;l++)
                    s[j][l]=(1LL*s[i][l]*inv%mod-s[j][l]+mod)%mod;
            }
        }
        for(int i=1;i<=k;i++) sum=1LL*sum*s[i][i]%mod;
        return sum;
    }
}a,b,ans;

struct poly
{
    int s[101];
    poly(int x=0){memset(s,0,sizeof(s));s[0]=x;}
    poly operator^(int x)
    {
        poly c;
        for(int i=k<<1;i;i--)
            c.s[i]=(s[i-1]+1LL*x*s[i])%mod;
        c.s[0]=(1LL*s[0]*x)%mod;
        return c;
    }
    poly operator*(int x)
    {
        poly c(0);
        for(int i=0;i<=k<<1;i++)
            c.s[i]=1LL*s[i]*x%mod;
        return c;
    }
    poly operator+(poly b)
    {
        poly c(0);
        for(int i=0;i<=k<<1;i++)
            c.s[i]=(s[i]+b.s[i])%mod;
        return c;
    }
    poly operator*(poly b)
    {
        poly c(0);
        for(int i=0;i<=k;i++)
            for(int j=0;j<=k;j++)
                c.s[i+j]=(c.s[i+j]+1LL*s[i]*b.s[j])%mod;
        return c;
    }
    friend poly operator%(poly a,poly b)
    {
        for(int i=k;i>=0;i--)
        {
            int t=1LL*a.s[i+k]*pow(b.s[k],mod-2)%mod;
            for(int j=0;j<=k;j++)
                a.s[i+j]=(a.s[i+j]-1LL*b.s[j]*t%mod+mod)%mod;
        }
        return a;
    }
}F(0),Ans(1);

int main()
{
    scanf("%s",st+1);k=read();
    for(register int i=1;i<=k;i++)
        for(register int j=1;j<=k;j++)
            b.s[i][j]=a.s[i][j]=read();
    for(register int t=0;t<=k;t++)
    {
        memcpy(b.s,a.s,sizeof(b.s));
        for(int i=1;i<=k;i++)
            b.s[i][i]=(b.s[i][i]-t+mod)%mod;
        y[t]=b.calc();
    }
    for(register int t=0;t<=k;t++)
    {
        poly tmp(1);
        for(register int i=0;i<=k;i++) if(i!=t)
            tmp=tmp^(mod-i),tmp=tmp*pow((t-i+mod)%mod,mod-2);
        tmp=tmp*y[t];F=F+tmp;
    }
    poly x(0);x.s[1]=1;
    for(int i=strlen(st+1);i;i--)
    {
        if(st[i]=='1') Ans=Ans*x%F;
        x=x*x%F;
    }
    memset(b.s,0,sizeof(b.s));
    for(int i=1;i<=k;i++) b.s[i][i]=1;
    for(int t=0;t<k;t++,b=a*b)
        for(int i=1;i<=k;i++)
            for(int j=1;j<=k;j++)
                ans.s[i][j]=(ans.s[i][j]+1LL*Ans.s[t]*b.s[i][j])%mod;
    for(register int i=1;i<=k;i++)
        for(register int j=1;j<=k;j++)
            printf("%d%c",ans.s[i][j],j!=k?' ':'
');
    return 0;
}