首页 > 代码库 > [BZOJ]4162: shlw loves matrix II
[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
5 9 5
5 4 0
8 8 8
Sample Output
110 121 65
45 61 25
144 168 104
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\n",ans[i++][j]))for(j=1;j<k;++j)printf("%d ",ans[i][j]); }
[BZOJ]4162: shlw loves matrix II
声明:以上内容来自用户投稿及互联网公开渠道收集整理发布,本网站不拥有所有权,未作人工编辑处理,也不承担相关法律责任,若内容有误或涉及侵权可进行投诉: 投诉/举报 工作人员会在5个工作日内联系你,一经查实,本站将立刻删除涉嫌侵权内容。