BJ模拟 装饰地板【状压dp+特征多项式优化矩阵快速幂】

xiaoxiao2021-02-28  42

题目大意:

给一个 mn m ∗ n 的地板,有 s1 s 1 1×2 1 × 2 的横地砖, s2 s 2 2×1 2 × 1 的竖地砖,问有多少种铺满的方式,对998244353取模。 (m6,n102501s1,s21e9) ( m ≤ 6 , n ≤ 10 2501 , s 1 , s 2 ≤ 1 e 9 )

解题思路:

看到 m m 很小,很容易想到2mn2mn的状压dp,进一步想到 (2m)3log2n ( 2 m ) 3 l o g 2 n 的矩阵快速幂优化dp,但还是不够。 主要问题是如何优化矩阵快速幂复杂度(暴力循环展开就行了)。 考虑 m×m m × m 矩阵的特征多项式 f(λ)=|AλI|=i=0maiλi f ( λ ) = | A − λ I | = ∑ i = 0 m a i λ i ,这是个 m m 次多项式。 根据CayleyHamiltonCayley−Hamilton定理有 f(A)=0 f ( A ) = 0 ,即把矩阵 A A 当作变量带入f(x)f(x)算可以得到一个0矩阵。 所以我们可以把多余的 f(x) f ( x ) 模掉,可以得到 An=Anmodf(A)=i=0ngiAi A n = A n mod f ( A ) = ∑ i = 0 n g i A i ,这样就可以只算 n n 项就可以呢。

那如何求一个矩阵的特征多项式呢?

如果这是一个常系数线性递推矩阵,那么有结论:f(λ)=λni=1nA1iλnif(λ)=λn−∑i=1nA1iλn−i 然而这道题只是个一般的矩阵,所以我们只能插 n+1 n + 1 λ λ 到矩阵 AλI A − λ I 中算出其行列式,再用拉格朗日插值法及多项式乘除法算出特征多项式系数,时间复杂度是 O(n4) O ( n 4 ) 的,可以说是很暴力了。

如果有大佬精通多项式算法,可以把多项式直接带进矩阵高斯消元算,应该是 O(n3) O ( n 3 ) 的。

求出 f(A) f ( A ) 的系数后,就可以用多项式快速幂 O(n2logm) O ( n 2 l o g m ) Ammodf(A) A m mod f ( A ) 的系数,接着求出 Am A m 了。

注意本题n是高精度,要用十进制快速幂。

所以这道题的总复杂度大概是 O(m4m+16m+4mlogn) O ( m 4 m + 16 m + 4 m l o g n ) 的。

其实这道题没有什么思维难度,重点在多项式算法(暴力版),练手是很好的,还可以封装一下模板……。

#include<bits/stdc++.h> #define ll long long using namespace std; const int N=150,M=2505,mod=998244353; int m,s1,s2,mx,l1,l2,ans; char s[M],t[M]; int l[M],r[M]; int Pow(ll x,int y) { ll res=1; for(;y;y>>=1,x=x*x%mod) if(y&1)res=res*x%mod; return res; } struct Poly { ll a[N];int deg; Poly(){memset(a,0,sizeof(a));} inline friend Poly operator + (const Poly &A,const Poly &B) { Poly res;res.deg=max(A.deg,B.deg); for(int i=0;i<=res.deg;i++)res.a[i]=(A.a[i]+B.a[i])%mod; return res; } inline friend Poly operator * (const Poly &A,const int &b) { Poly res=A; for(int i=0;i<=res.deg;i++)res.a[i]=A.a[i]*b%mod; return res; } inline friend Poly operator * (const Poly &A,const Poly &B) { Poly res;res.deg=A.deg+B.deg; for(int i=0;i<=A.deg;i++) for(int j=0;j<=B.deg;j++) res.a[i+j]=(res.a[i+j]+A.a[i]*B.a[j])%mod; return res; } inline friend Poly operator % (const Poly &A,const Poly &B) { Poly res=A;ll inv=Pow(B.a[B.deg],mod-2),t; for(int i=res.deg;i>=B.deg;i--) { t=res.a[i]*inv%mod; for(int j=0;j<=B.deg;j++) res.a[i-B.deg+j]=(res.a[i-B.deg+j]-t*B.a[j])%mod; } res.deg=B.deg; return res; } inline friend Poly operator / (Poly A,const Poly &B) { Poly res;res.deg=A.deg-B.deg; ll inv=Pow(B.a[B.deg],mod-2),t; for(int i=A.deg;i>=B.deg;i--) { res.a[i-B.deg]=A.a[i]*inv%mod; for(int j=0;j<=B.deg;j++) A.a[i-B.deg+j]=(A.a[i-B.deg+j]-res.a[i-B.deg]*B.a[j])%mod; } return res; } inline friend Poly Pow(Poly A,int *b,int len,const Poly &C) { Poly res,t;res.a[0]=1,res.deg=0; for(int i=0;i<len;i++) { t=A; for(int j=1;j<10;j++) { if(j==b[i])res=res*A%C; A=A*t%C; } } return res; } }F; struct matrix { ll a[N][N]; matrix(){memset(a,0,sizeof(a));} inline friend matrix operator * (const matrix &A,const matrix &B) { matrix res; for(int i=0;i<=mx;i++) for(int k=0;k<=mx;k++) for(int j=0;j<=mx;j++) res.a[i][j]=(res.a[i][j]+A.a[i][k]*B.a[k][j])%mod; return res; } inline friend matrix operator * (const matrix &A,const int &b) { matrix res; for(int i=0;i<=mx;i++) for(int j=0;j<=mx;j++) res.a[i][j]=A.a[i][j]*b%mod; return res; } inline friend matrix operator + (const matrix &A,const matrix &B) { matrix res; for(int i=0;i<=mx;i++) for(int j=0;j<=mx;j++) res.a[i][j]=(A.a[i][j]+B.a[i][j])%mod; return res; } inline friend matrix operator - (const matrix &A,const matrix &B) { matrix res; for(int i=0;i<=mx;i++) for(int j=0;j<=mx;j++) res.a[i][j]=(A.a[i][j]-B.a[i][j])%mod; return res; } inline friend int det(matrix A) { ll res=1; for(int i=0;i<=mx;i++) { if(!A.a[i][i]) { int bz=0; for(int j=i+1;j<=mx;j++) if(A.a[j][i]) { bz=1; for(int k=i;k<=mx;k++) swap(A.a[j][k],A.a[i][k]); break; } if(!bz)return 0; } res=res*A.a[i][i]%mod; int inv=Pow(A.a[i][i],mod-2); for(int j=i+1;j<=mx;j++) { int b=A.a[j][i]*inv%mod; for(int k=i;k<=mx;k++) A.a[j][k]=(A.a[j][k]-A.a[i][k]*b)%mod; } } return res; } }A,I; void pre() { l1=strlen(s); for(int i=0;i<l1;i++)l[i]=s[i]-'0'; reverse(l,l+l1); l2=strlen(t); for(int i=0;i<l2;i++)r[i]=t[i]-'0'; reverse(r,r+l2); r[0]++; for(int i=0;i<l2;i++) if(r[i]>=10)r[i+1]+=r[i]/10,r[i]%=10; else break; if(r[l2])l2++; } void Init() { mx=1<<m; A.a[0][mx]=A.a[mx][mx]=1; for(int i=0;i<mx;i++) for(int j=0;j<mx;j++) { int k;ll res=1; for(k=1;k<=m;) { int a=(i&(1<<k-1)),b=(j&(1<<k-1)),c=(i&(1<<k)),d=(j&(1<<k)); if(k==m)c=d=1; if(a&&!b)k++; else if(!a&&b)res=res*s1%mod,k++; else if(!a&&!b&&!c&&!d)res=res*s2%mod,k+=2; else break; } if(k>m)A.a[i][j]=res; } for(int i=0;i<=mx;i++)I.a[i][i]=1; } Poly get_Poly(matrix A) { static ll x[N],y[N]; for(int i=1;i<=mx+2;i++)x[i]=i,y[i]=det(A-I*i); Poly F,G;G.a[0]=1,G.deg=0,F.deg=mx+1; for(int i=1;i<=mx+2;i++) { Poly t;t.deg=1,t.a[1]=1,t.a[0]=-x[i]; G=G*t; } for(int i=1;i<=mx+2;i++) { ll a=1; for(int j=1;j<=mx+2;j++)if(j!=i) a=a*(x[i]-x[j])%mod; a=Pow(a,mod-2),a=a*y[i]%mod; Poly t;t.deg=1,t.a[1]=1,t.a[0]=-x[i]; t=G/t;F=F+t*a; } return F; } int calc(matrix A,int *b,int len) { Poly G;G.a[1]=1,G.deg=1; G=Pow(G,b,len,F); matrix t=A,P=I;memset(A.a,0,sizeof(A.a)); for(int i=0;i<=mx;i++) A=A+P*G.a[i],P=P*t; matrix res;res.a[0][0]=1;res=res*A; return (res.a[0][mx]+mod)%mod; } int main() { //freopen("lx.in","r",stdin); scanf("%s%s",s,t);pre(); scanf("%d%d%d",&m,&s1,&s2); Init(); F=get_Poly(A); ans=(calc(A,r,l2)-calc(A,l,l1)+mod)%mod; printf("%d\n",ans); return 0; }
转载请注明原文地址: https://www.6miu.com/read-2622707.html

最新回复(0)