bzoj4944: [Noi2017]泳池（dp+常系数齐次线性递推）

传送门
首先转化问题，变成计算能选择的最大的游泳场的面积不超过$K$的概率是多少然后差分即可。

设$f_i$表示强制第$i$个格子是危险的，前$i-1$个格子的面积不超过$K$的概率。

那么我们最后可以强制第$n+1$个格子是危险的计算出前$n$个格子不超过$K$的概率，答案就是$\frac{f_{n+1}}{1-q}$

考虑怎么递推这个$f$，由于最多连续有$K$个格子是连续的，因此不难想到这个$f$数组有形如$f_n={\sum }_{i=1}^{k+1}coe{f}_i{f}_{k-i}$的递推式，如果我们能推出$coef$数组与$f_{0,1,2,...,k}$，就能在$O(k^3\log n)/O(k^2\log n)/O(k\log k\log n)$的时间内推出$f_{n+1}$。

其中$O(k^3\log n)$的矩阵快速幂只能拿到部分分，而$O(k^2\log n)/O(k\log k\log n)$的多项式取模可以拿到$100pts$的好成绩。

不会用多项式取模优化递推的可以看我这篇讲常系数齐次线性递推的博客。

现在假装我们已经会递推了
于是问题变成考虑$coef$怎么求。

考虑到我们枚举$f$转移的意义相当于是枚举了长度为$1001$，宽度为$n-i$的矩形中选出来挨着下边界的矩形的最大面积不超过$K$的概率。

那么定义状态$g_{i,j}$表示对于一个宽度为$i$的矩形，最矮的限制点的高度为$j$使得从中选出来挨着下边界的矩形的最大面积不超过$K$的概率。

如何转移？

我们从左到右枚举第一个出现的最矮的限制，然后就转化成了子问题。

即$g_{i,j}={\sum }_{p=1}^i({\sum }_{k<p}g(i-1,k))({\sum }_{k\le p}g(n-i,k))$
发现可以后缀和优化转移，最后我们的$coe{f}_i=g_{i,1}\ast p$

然后再用多项式取模优化递推即可。

代码：

#include<bits/stdc++.h>
#define ri register int
using namespace std;
const int rlen=1<<18|1;
inline char gc(){
	static char buf[rlen],*ib,*ob;
	(ib==ob)&&(ob=(ib=buf)+fread(buf,1,rlen,stdin));
	return ib==ob?-1:*ib++;
}
inline int read(){
	int ans=0;
	char ch=gc();
	bool f=1;
	while(!isdigit(ch))f^=ch=='-',ch=gc();
	while(isdigit(ch))ans=((ans<<2)+ans<<1)+(ch^48),ch=gc();
	return f?ans:-ans;
}
typedef vector<int> poly;
typedef long long ll;
const int mod=998244353;
inline int add(int a,int b){return (a+=b)>=mod?a-mod:a;}
inline int dec(int a,int b){return (a-=b)<0?a+mod:a;}
inline int mul(int a,int b){return (ll)a*b%mod;}
inline void Add(int&a,int b){(a+=b)>=mod?a-=mod:a;}
inline void Dec(int&a,int b){(a-=b)<0?a+=mod:a;}
inline void Mul(int&a,int b){a=(ll)a*b%mod;}
inline int ksm(int a,int p){int ret=1;for(;p;p>>=1,Mul(a,a))if(p&1)Mul(ret,a);return ret;}
int w[23],invv[23],lim,tim;
vector<int>rev[23];
inline void init_w(){
	w[22]=ksm(3,(mod-1)>>23);
	for(ri i=21;~i;--i)w[i]=mul(w[i+1],w[i+1]);
	invv[0]=1;
	for(ri i=1,iv=mod+1>>1;i<23;++i)invv[i]=mul(invv[i-1],iv);
}
inline void init(const int&up){
	lim=1,tim=0;
	while(lim<up)lim<<=1,++tim;
	if(rev[tim].size())return;
	rev[tim].resize(lim);
	for(ri i=0;i<lim;++i)rev[tim][i]=(rev[tim][i>>1]>>1)|((i&1)<<(tim-1));
}
inline void ntt(poly&a,const int&type){
	for(ri i=0;i<lim;++i)if(i<rev[tim][i])swap(a[i],a[rev[tim][i]]);
	for(ri i=1,t=0,a0,a1;i<lim;i<<=1,++t)for(ri j=0,len=i<<1;j<lim;j+=len)
	for(ri k=0,mt=1;k<i;++k,Mul(mt,w[t]))a0=a[j+k],a1=mul(a[j+k+i],mt),a[j+k]=add(a0,a1),a[j+k+i]=dec(a0,a1);
	if(~type)return;
	reverse(++a.begin(),a.end());
	for(ri i=0;i<lim;++i)Mul(a[i],invv[tim]);
}
inline poly operator-(poly a,poly b){
	int n=a.size(),m=b.size();
	poly c(max(n,m));
	for(ri i=0;i<n;++i)Add(c[i],a[i]);
	for(ri i=0;i<m;++i)Dec(c[i],b[i]);
	return c;
}
inline poly operator*(poly a,poly b){
	int n=a.size(),m=b.size(),t=n+m-1;
	if(t<=128){
		poly c(t);
		for(ri i=0;i<n;++i)for(ri j=0;j<m;++j)Add(c[i+j],mul(a[i],b[j]));
		return c;
	}
	init(t);
	a.resize(lim),ntt(a,1);
	b.resize(lim),ntt(b,1);
	for(ri i=0;i<lim;++i)Mul(a[i],b[i]);
	return ntt(a,-1),a.resize(t),a;
}
inline poly poly_inv(poly a,int k){
	poly b(1,ksm(a[0],mod-2)),c;
	for(ri i=4,up=k<<2;i<up;i<<=1){
		init(i);
		c=a,c.resize(i>>1);
		b.resize(lim),ntt(b,1);
		c.resize(lim),ntt(c,1);
		for(ri j=0;j<lim;++j)Mul(b[j],dec(2,mul(b[j],c[j])));
		ntt(b,-1),b.resize(i>>1);
	}
	return b.resize(k),b;
}
inline poly operator/(poly a,poly b){
	int n=a.size(),m=b.size(),t=n-m+1;
	reverse(a.begin(),a.end());
	reverse(b.begin(),b.end()),b=poly_inv(b,t);
	return a=a*b,a.resize(t),reverse(a.begin(),a.end()),a;
}
inline poly operator%(poly a,poly b){if(a.size()<b.size())return a;poly c=a-a/b*b;return c.resize(b.size()),c;}
const int N=1005;
int a[N],b[N],ss[N][N];
inline int calc(int*coe,int*a,int k,int n){
	if(n+1<k)return a[n+1];
	poly md(k+1),f(2),g(1,1);
	f[1]=md[k]=1;
	for(ri i=1;i<=k;++i)md[k-i]=coe[i]?mod-coe[i]:0;
	for(++n;n;n>>=1,f=f*f%md)if(n&1)g=g*f%md;
	int ret=0;
	for(ri i=0;i<k;++i)Add(ret,mul(g[i],a[i]));
	return ret;
}
int n,k,q,p,pw[N];
inline int solve(int n,int K){
	memset(b,0,sizeof(b));
	memset(ss,0,sizeof(ss));
	for(ri i=1;i<=K+1;++i)ss[0][i]=1;
	for(ri j=K;j;--j)for(ri t=0,i=1,up=K/j;i<=up;++i,t=0){
		for(ri k=1;k<=i;++k)Add(t,mul(ss[k-1][j+1],ss[i-k][j]));
		ss[i][j]=add(ss[i][j+1],mul(t,mul(pw[j],p)));
	}
	++K;
	a[1]=p;
	for(ri i=1;i<K;++i)a[i+1]=mul(p,ss[i][1]);
	b[0]=1;
	for(ri i=1;i<K;++i)for(ri j=0;j<i;++j)Add(b[i],mul(b[j],a[i-j]));
	return calc(a,b,K,n);
}
int main(){
	init_w();
	n=read(),k=read(),q=read(),Mul(q,ksm(read(),mod-2)),p=dec(1,q);
	pw[0]=1;
	for(ri i=1;i<=k;++i)pw[i]=mul(pw[i-1],q);
	int a=solve(n,k),b=solve(n,k-1);
	cout<<mul(ksm(p,mod-2),dec(a,b));
	return 0;
}