自动搬运

查看原文

来自洛谷，原作者为

	Salamander **

搬运于2025-08-24 21:53:07，当前版本为作者最后更新于2018-05-23 16:43:34，作者可能在搬运后再次修改，您可在原文处查看最新版

自动搬运只会搬运当前题目点赞数最高的题解，您可前往洛谷题解查看更多

以下是正文

首先，恰好等于$k$的概率可以通过差分转化为小于等于$k$的概率减去小于等于$k-1$的概率。

考虑怎么求小于等于$k$的概率。设$f_n$为底部宽为$n$时的答案，那么我们可以枚举最后连续的一段白色的长度，即$\displaystyle f_n=\sum_{i=0}^{k}f_{n-i-1}\cdot\sum_{j=2}^{\lfloor{k/i}\rfloor+1}dp_{i,j}$。其中$dp_{i,j}$表示一个宽为$i$的矩形，高度为$1$至$j-1$都全部是白色，高度为$j$的一行开始出现障碍，并且包含的所有贴住底部的白色矩形(包括宽不一定为$i$和高超过$j$的那些)大小都不超过$k$的概率。记$g_{i,j}=\sum_{l\geq j}dp_{i,l}$。

如果我们可以求出$dp$和$g$数组，那么转移方程就变成了齐次常系数线性递推的形式：$\displaystyle f_n=\sum_{i=0}^{k}f_{n-i-1}\cdot g_{i,2}$。注意到$k$比较大，传统的矩阵乘法为$O(k^3\log n)$无法通过，可以用多项式取模的方法来做。

先考虑怎么求出$dp$和$g$。枚举高度为$j$处的第一个障碍，

$\displaystyle dp_{i,j}=[i\cdot(j-1)\leq k](1-p)p^{j-1}\sum_{l=1}^{i}\left(\sum_{q>j}dp_{l-1,q}\right)\left(\sum_{q\geq j}dp_{i-l,q}\right)$

前缀和优化一下，就可以做到$O(k^2\log k)$了，因为合法的状态中要求$i\cdot(j-1)\leq k$，所以一共有$\sum_{i=1}^{k}\lfloor{k/i}\rfloor=O(k\log k)$个非零状态，每次转移为$O(k)$，已经可以通过此题。实际上还可以用FFT优化到$k\log^2k$。

$~$

下面是用多项式取模优化齐次常系数线性递推的方法：

对于递推$f_n=\sum_{i=1}^ka_kf_{n-k}$

设矩阵乘法中的转移矩阵为$A$，设初始状态矩阵为$B[f_k,f_{k-1},...,f_1]$，我们现在要求$BA^n$，即求$A^n$。

考虑$A$的特征多项式$\lambda(x)=det(xI-A)$，$I$为单位矩阵。由于转移矩阵的特殊性，在此处我们可以手动求一下特征多项式，$\lambda(x)=x^k-\sum_{i=1}^ka_ix^{k-i}$，是一个$k$次多项式，并且满足$\lambda(A)=0$。

如果我们把$A^n$表示成若干倍的$\lambda(A)$，即$A^n=\lambda(A)g(A)+r(A)$。由于$\lambda(A)=0$，所以$A^n=r(A)$。而$r(A)=A^n\mod\lambda(A)$，我们只需要求出$x^n$模特征多项式后的多项式，然后将$A$带入即可求出$A^n$。

$~$

多项式取模可以使用Picks的方法，用FFT做到$k\log k$，加上多项式快速幂复杂度就是$O(k\log k\log n)$，边乘边取模。当然在这里没有必要，直接暴力取模即可。

$~$

如何暴力取模?假设要求$f(x)\mod g(x)$，设$f(x)=g(x)h(x)+r(x)$，我们就是要求$r(x)$，并且要求$r(x)$的次数小于$g(x)$的次数。不妨将$g(x)$看成$0$，那么此时$f(x)$就等于$r(x)$。设$g(x)=\sum_{i=0}^kb_ix^i=0$，那么我们可以得出$\displaystyle x^k=-\sum_{i=0}^{k-1}\frac{b_i}{b_k}x^i$，我们只需要把$f(x)$中系数过大的一部分迭代，即可将系数一步步压到比$k$小。

两个$k$次多项式相乘，超出的部分系数是$O(k)$的，所以暴力取模的复杂度为$O(k^2)$，套上快速幂就是$O(k^2\log n)$，可以通过此题。

$~$

然后我们考虑怎么统计答案。我们已经用一个关于$A$的$k-1$次多项式$A^n=R(A)=\sum_{i=0}^{k-1}r_iA^i$表示出了$A^n$。

根据矩阵乘法的分配律，$B\cdot R(A)=\sum_{i=0}^{k-1}r_i\cdot B\cdot A^i$。我们要的结果是矩阵的第一个元素，而$B\cdot A^i$的第一个元素就是$f_{k+i}$，所以我们要先求出$f_k,f_{k+1},...,f_{2k}$，然后直接扫一遍统计答案即可。

$~$

复杂度$O(k^2(\log k+\log n))$，其实可以做到$O(k\log k(\log k+\log n))$。

我代码中的dp部分和上面写的转移不太一样。

#include<bits/stdc++.h>

#define For(i,_beg,_end) for(int i=(_beg),i##end=(_end);i<=i##end;++i)
#define Rep(i,_beg,_end) for(int i=(_beg),i##end=(_end);i>=i##end;--i)

template<typename T>T Max(const T &x,const T &y){return x<y?y:x;}
template<typename T>T Min(const T &x,const T &y){return x<y?x:y;}
template<typename T>int chkmax(T &x,const T &y){return x<y?(x=y,1):0;}
template<typename T>int chkmin(T &x,const T &y){return x>y?(x=y,1):0;}
template<typename T>void read(T &x){
	T f=1;char ch=getchar();
	for(;ch<'0'||ch>'9';ch=getchar())if(ch=='-')f=-1;
	for(x=0;ch>='0'&&ch<='9';ch=getchar())x=x*10+ch-'0';
	x*=f;
}

typedef long long LL;
const int N=1010,mod=998244353;
int n,m;
LL a,b,p,q,dp[N][N],g[N][N],pw[N];
LL A[N],f[N<<1];

LL power(LL,LL);
LL Solve(int);

int main(){
	read(n);read(m);read(a);read(b);
	p=a*power(b,mod-2)%mod;q=(mod+1-p)%mod;
	pw[0]=1;
	For(i,1,m) pw[i]=pw[i-1]*p%mod;
	
	printf("%lld\n",(Solve(m)-Solve(m-1)+mod)%mod);
	return 0;
}

LL Solve(int k){
	memset(dp,0,sizeof dp);
	memset(g,0,sizeof g);
	For(i,1,k+2) g[0][i]=dp[0][i]=1;
	For(i,1,k) Rep(j,k/i+1,2){
		For(l,1,i) dp[i][j]=(dp[i][j]+g[l-1][j+1]*g[i-l][j]%mod*pw[l-1]%mod*q)%mod;
		g[i][j]=(g[i][j+1]*pw[i]+dp[i][j])%mod;
	}
	memset(A,0,sizeof A);
	For(i,0,k) A[i+1]=q*g[i][2]%mod*pw[i]%mod;
	memset(f,0,sizeof f);
	f[0]=1;
	For(i,1,k){
		f[i]=g[i][2]*pw[i]%mod;
		For(j,1,i) f[i]=(f[i]+A[j]*f[i-j])%mod;
	}
	k++;
	For(i,k,k<<1) For(j,1,k) f[i]=(f[i]+A[j]*f[i-j])%mod;
	if(n<=k)return f[n];
	int y=n-k,len=1,L=0;
	LL res[N<<2],tmp[N<<2],x[N<<2];
	memset(res,0,sizeof res);
	memset(x,0,sizeof x);
	res[0]=1;x[1]=1;
	for(;y;y>>=1){
		if(y&1){
			For(i,0,L+len) tmp[i]=0;
			For(i,0,L) For(j,0,len) tmp[i+j]=(tmp[i+j]+res[i]*x[j])%mod;
			L+=len;
			Rep(i,L,k) For(j,1,k) tmp[i-j]=(tmp[i-j]+tmp[i]*A[j])%mod;
			chkmin(L,k-1);
			For(i,0,L) res[i]=tmp[i];
		}
		For(i,0,len+len) tmp[i]=0;
		For(i,0,len) For(j,0,len) tmp[i+j]=(tmp[i+j]+x[i]*x[j])%mod;
		len<<=1;
		Rep(i,len,k) For(j,1,k) tmp[i-j]=(tmp[i-j]+tmp[i]*A[j])%mod;
		chkmin(len,k-1);
		For(i,0,len) x[i]=tmp[i];
	}
	LL ans=0;
	For(i,0,k-1) ans=(ans+res[i]*f[i+k])%mod;
	return ans;
}
LL power(LL x,LL y){
	LL res=1;
	for(;y;y>>=1,x=x*x%mod) if(y&1) res=res*x%mod;
	return res;
}