[BZOJ]4162: shlw loves matrix II

时间：2022-05-03 16:24:59

Time Limit: 30 Sec  Memory Limit: 128 MB

Description

　　给定矩阵 M，请计算 M^n，并将其中每一个元素对 1000000007 取模输出。

Input

　　第 1 行包含两个整数 n,k，其中 n 使用二进制表示，可能含有前导零;
　　余下 k 行描述了一个 k * k 的矩阵 M。

Output

　　输出题目描述中要求的矩阵，格式同输入。

Sample Input

　　010 3
　　5 9 5
　　5 4 0
　　8 8 8

Sample Output

　　110 121 65
　　45 61 25
　　144 168 104

HINT

　　对于 100% 数据，满足 n <= 2^10000;k <= 50; 0 <= Mij < 10^9 +7

Solution

　　矩阵乘法特征多项式优化，矩阵M的特征多项式是$f(\lambda)=|M-\lambda I|$，随便带入k+1个$\lambda$，高斯消元求出行列式的值，然后插值就能求出M的特征多项式（这里的总复杂度为$O(k^4)$，主要在计算行列式上面）。根据Cayley-Hamilton定理，$f(M)=0$，也就是说对同一个式子减去若干个$f(M)$后值不变，我们即可用M的k-1次多项式表示一个M的次幂，两个多项式相乘后对$f(M)$取模，然后就可以快速幂了，多项式取模和乘法用$O(k^2)$暴力就够了，总复杂度$O(k^4+k^2logn)$。

Code

#include<cstdio>

#include<cstring>

#include<algorithm>

using namespace std;

#define ll long long

inline int read()

{

    int x;char c;

    while((c=getchar())<''||c>'');

    for(x=c-'';(c=getchar())>=''&&c<='';)x=x*+c-'';

    return x;

}

#define ML 10000

#define MK 50

#define MN 100

#define MOD 1000000007

char s[ML+];

int k,c[MK+][MK+],x[MK+],ans[MK+][MK+];

inline int mo1(int x){return x<MOD?x:x-MOD;}

inline int mo2(int x){return x<?x+MOD:x;}

int inv(int x)

{

    int r=,y=MOD-;

    for(;y;y>>=,x=1LL*x*x%MOD)if(y&)r=1LL*r*x%MOD;

    return r;

}

struct poly

{

    int n,a[MN+];

    poly(int n=,int a0=,int a1=):n(n){memset(a,,sizeof(a));a[]=a0;a[]=a1;}

    friend poly operator+(poly a,const poly&b)

    {

        a.n=max(a.n,b.n);

        for(int i=;i<=a.n;++i)a.a[i]=mo1(a.a[i]+b.a[i]);

        return a;

    }

    friend poly operator*(const poly&a,const poly&b)

    {

        poly c(a.n+b.n);int i,j;ll x;

        for(i=;i<=c.n;c.a[i++]=x%MOD)for(x=j=;j<=i&&j<=a.n;++j)

            x+=1LL*a.a[j]*b.a[i-j],x>8e18?x%=MOD:;

        return c;

    }

    friend poly operator*(poly a,int b)

    {

        for(int i=;i<=a.n;++i)a.a[i]=1LL*a.a[i]*b%MOD;

        return a;

    }

    friend poly operator/(poly a,const poly&b)

    {

        poly c(a.n-b.n);int i,j;

        for(i=a.n;i>=b.n;--i)

        {

            c.a[i-b.n]=1LL*a.a[i]*inv(b.a[b.n])%MOD;

            for(j=;j<=b.n;++j)a.a[i-j]=mo2((a.a[i-j]-1LL*b.a[b.n-j]*c.a[i-b.n])%MOD);

        }

        return c;

    }

    friend poly operator%(poly a,const poly&b)

    {

        int i,j,x;

        for(i=a.n;i>=b.n;--i)

        {

            x=1LL*a.a[i]*inv(b.a[b.n])%MOD;

            for(j=;j<=b.n;++j)a.a[i-j]=mo2((a.a[i-j]-1LL*b.a[b.n-j]*x)%MOD);

        }

        a.n=b.n-;

        return a;

    }

}f,a,p(,);

struct mat

{

    int z[MK+][MK+];

    mat(){memset(z,,sizeof(z));}

    friend mat operator*(const mat&a,const mat&b)

    {

        mat c;int i,j,k;ll x;

        for(i=;i<=MK;++i)for(j=;j<=MK;c.z[i][j++]=x%MOD)

            for(x=k=;k<=MK;++k)x+=1LL*a.z[i][k]*b.z[k][j],x>8e18?x%=MOD:;

        return c;

    }

}m,n;

int cal(int x)

{

    int i,j,l,ans=;

    for(i=;i<=k;++i)for(j=;j<=k;++j)c[i][j]=m.z[i][j];

    for(i=;i<=k;++i)c[i][i]=mo2(c[i][i]-x);

    for(i=;i<=k;++i)

    {

        for(j=i;j<=k;++j)if(c[j][i])break;

        if(j>k)return ;

        if(j>i)for(ans=MOD-ans,l=i;l<=k;++l)swap(c[i][l],c[j][l]);

        ans=1LL*ans*c[i][i]%MOD;

        for(j=i;++j<=k;)

        {

            x=1LL*c[j][i]*inv(c[i][i])%MOD;

            for(l=i;l<=k;++l)c[j][l]=mo2((c[j][l]-1LL*c[i][l]*x)%MOD);

        }

    }

    return ans;

}

int main()

{

    int i,j,l,v;

    scanf("%s",s+);k=read();

    for(i=;i<=k;++i)for(j=;j<=k;++j)m.z[i][j]=read();

    for(i=;i<=k;++i)x[i]=cal(i);

    for(i=;i<=k;++i)p=p*poly(,mo2(-i),);

    for(i=;i<=k;++i)

    {

        for(v=x[i],j=;j<=k;++j)if(i!=j)v=1LL*v*inv(mo2(i-j))%MOD;

        f=f+p/poly(,mo2(-i),)*v;

    }

    a=poly(,);p=poly(,,);

    for(i=strlen(s+);i;--i,p=p*p%f)if(s[i]>'')a=a*p%f;

    for(i=;i<=k;++i)n.z[i][i]=;

    for(l=;l<k;++l,n=n*m)for(i=;i<=k;++i)for(j=;j<=k;++j)

        ans[i][j]=(ans[i][j]+1LL*a.a[l]*n.z[i][j])%MOD;

    for(i=;i<=k;printf("%d\n",ans[i++][j]))for(j=;j<k;++j)printf("%d ",ans[i][j]);

}

相关文章