[BZOJ]4162: shlw loves matrix II

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())<'0'||c>'9');
    for(x=c-'0';(c=getchar())>='0'&&c<='9';)x=x*10+c-'0';
    return x;
}
#define ML 10000
#define MK 50
#define MN 100
#define MOD 1000000007
char s[ML+5];
int k,c[MK+5][MK+5],x[MK+5],ans[MK+5][MK+5];
inline int mo1(int x){return x<MOD?x:x-MOD;}
inline int mo2(int x){return x<0?x+MOD:x;}
int inv(int x)
{
    int r=1,y=MOD-2;
    for(;y;y>>=1,x=1LL*x*x%MOD)if(y&1)r=1LL*r*x%MOD;
    return r;
}
struct poly
{
    int n,a[MN+5];
    poly(int n=0,int a0=0,int a1=0):n(n){memset(a,0,sizeof(a));a[0]=a0;a[1]=a1;}
    friend poly operator+(poly a,const poly&b)
    {
        a.n=max(a.n,b.n);
        for(int i=0;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=0;i<=c.n;c.a[i++]=x%MOD)for(x=j=0;j<=i&&j<=a.n;++j)
            x+=1LL*a.a[j]*b.a[i-j],x>8e18?x%=MOD:0;
        return c;
    }
    friend poly operator*(poly a,int b)
    {
        for(int i=0;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=1;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=0;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-1;
        return a;
    }
}f,a,p(0,1);
struct mat
{
    int z[MK+5][MK+5];
    mat(){memset(z,0,sizeof(z));}
    friend mat operator*(const mat&a,const mat&b)
    {
        mat c;int i,j,k;ll x;
        for(i=0;i<=MK;++i)for(j=0;j<=MK;c.z[i][j++]=x%MOD)
            for(x=k=0;k<=MK;++k)x+=1LL*a.z[i][k]*b.z[k][j],x>8e18?x%=MOD:0;
        return c;
    }
}m,n;
int cal(int x)
{
    int i,j,l,ans=1;
    for(i=1;i<=k;++i)for(j=1;j<=k;++j)c[i][j]=m.z[i][j];
    for(i=1;i<=k;++i)c[i][i]=mo2(c[i][i]-x);
    for(i=1;i<=k;++i)
    {
        for(j=i;j<=k;++j)if(c[j][i])break;
        if(j>k)return 0;
        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+1);k=read();
    for(i=1;i<=k;++i)for(j=1;j<=k;++j)m.z[i][j]=read();
    for(i=0;i<=k;++i)x[i]=cal(i);
    for(i=0;i<=k;++i)p=p*poly(1,mo2(-i),1);
    for(i=0;i<=k;++i)
    {
        for(v=x[i],j=0;j<=k;++j)if(i!=j)v=1LL*v*inv(mo2(i-j))%MOD;
        f=f+p/poly(1,mo2(-i),1)*v;
    }
    a=poly(0,1);p=poly(1,0,1);
    for(i=strlen(s+1);i;--i,p=p*p%f)if(s[i]>'0')a=a*p%f;
    for(i=0;i<=k;++i)n.z[i][i]=1;
    for(l=0;l<k;++l,n=n*m)for(i=1;i<=k;++i)for(j=1;j<=k;++j)
        ans[i][j]=(ans[i][j]+1LL*a.a[l]*n.z[i][j])%MOD;
    for(i=1;i<=k;printf("%d
",ans[i++][j]))for(j=1;j<k;++j)printf("%d ",ans[i][j]);
}