uoj316【NOI2017】泳池

题目链接
(S=k)可以拆成(Sle k)减去(Sle k-1)。用((i,j))表示第i行第j列。
设(g(i,j))表示前i行前j列都安全其他未知满足条件的概率，(h(i,j))表示前i行前j列是安全的但是((i+1,j))是危险的，其他未知，满足条件的概率。当(i*j>k)时两个数组的值都是0。
(g(i,j))的转移可以在(i+1)行枚举最右的危险格子来转移：(g(i,j)=sum_{k=0}^jh(i,k)*g(i+1,j-k))。(h(i,j))同理枚举(i+1)行除了((i+1,j))以外的最右的危险格子来转移：(h(i,j)=sum_{k=0}^{j-1}h(i,k)*g(i+1,j-k-1)*(1-q)*q^i)，这里dp的复杂度是(sum_{i}(frac{k}{i})^2=O(k^2))的。

然后令(f(i))表示前i列的最大矩形(le k)的概率，转移时枚举第一行最长的连续安全区的长度(为(j-1))：(f(i)=sum_{j=1}^{k+1}f(i-j)*g(1,j-1)*(1-q))，注意(ile k)时还要从(g(1,i))转移，因为可能它就是第一个连续的安全区。

这是常系数线性递推，k只有1000，可以在(O(k^2logn))的时间内完成，考虑(f_i=sum_{j=1}^kf_{i-j}*a_j)这样一个递推公式，有一个定理是说设这个矩阵的特征多项式为(g(lambda))，(对于一般的矩阵(g(lambda)=|lambda I-A|)，而这个矩阵(g(lambda)=lambda^k-a_1lambda^{k-1}-...-a_{k-1}lambda-a_k))，那么(g(A)=0)。因此(A^k=a_1A^{k-1}+...+a_k)也就是说(A^k)可以用(A^{k-1}...A^0)线性表示出来，然后假如我们可以用那k个矩阵线性表示出(A^i)，那么(A^i)乘上(A)后再把多出来的(A^k)项用那个定理拆掉，所以对于任意的(A^i)都可以用(A^{k-1}...A^0)线性表示出来。
我们预处理(A^k)~(A^{2k-2})的线性表示，两个矩阵相乘相当于两个k-1次多项式相乘，乘出来大于k-1次方的项就拆掉。

最后，设(A^n=b_0A_0+b_1A_1+...+b_{k-1}A_{k-1})，由于我们要求(A^n*X)的某一项的值，即(b_0A_0X+b_1A_1X+...+b_{k-1}A_{k-1}X)的那一项的值，(k^2)暴力算出(A_0)到(A_{k-1})的值即可。
复杂度(O(k^2logn))，其实两部分都可以进一步优化。

#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<iostream>
#include<vector>
#include<algorithm>
#include<cmath>
#define P puts("lala")
#define cp cerr<<"lala"<<endl
#define ln putchar('
')
#define pb push_back
#define fi first
#define se second
using namespace std;
inline int read()
{
    char ch=getchar();int g=1,re=0;
    while(ch<'0'||ch>'9') {if(ch=='-')g=-1;ch=getchar();}
    while(ch<='9'&&ch>='0') re=(re<<1)+(re<<3)+(ch^48),ch=getchar();
    return re*g;
}
typedef long long ll;
typedef pair<int,int> pii;

const int mod=998244353;
const int N=1050;
ll qpow(ll a,int n)
{
	ll ans=1;
	for(;n;n>>=1,a=a*a%mod) if(n&1) ans=ans*a%mod;
	return ans;
}
int m[N<<1][N],C[N<<1];

void mul(int *a,int *b,int k)
{
	for(int i=0;i<=(k-1<<1);++i) C[i]=0;
	for(int i=0;i<k;++i) for(int j=0;j<k;++j) C[i+j]=(C[i+j]+1ll*a[i]*b[j])%mod;
	for(int i=k;i<=(k-1<<1);++i) if(C[i])
		for(int j=0;j<k;++j) C[j]=(C[j]+1ll*m[i][j]*C[i])%mod;
	for(int i=0;i<k;++i) a[i]=C[i];
}

int f[N<<1],g[N][N],h[N][N],a[N];
int q,pwq[N],A[N<<1],S[N<<1];

int solve(int K,int n)
{
	if(!K) return qpow((1-q+mod)%mod,n);
	memset(g,0,sizeof(g)); memset(h,0,sizeof(h));
	memset(f,0,sizeof(f)); memset(a,0,sizeof(a));
	memset(m,0,sizeof(m)); memset(A,0,sizeof(A));
	memset(S,0,sizeof(S));
	g[K][1]=1ll*pwq[K]*(1-q+mod)%mod;
	for(int i=1;i<=K;++i) g[i][0]=1,h[i][0]=1;
	for(int i=K-1;i>=1;--i)
	{
		for(int j=1;j*i<=K;++j)
		{
			for(int k=0;k<j;++k)
				h[i][j]=(h[i][j]+1ll*h[i][k]*g[i+1][j-k-1]%mod*pwq[i]%mod*(1-q+mod))%mod;
			for(int k=0;k<=j;++k)
				g[i][j]=(g[i][j]+1ll*g[i+1][j-k]*h[i][k])%mod;
		}
	}
	K++;
	for(int i=1;i<=K;++i) a[i]=1ll*g[1][i-1]*(1-q+mod)%mod;

	f[0]=1;
	for(int i=1;i<K;++i) 
	{
		for(int j=1;j<=i;++j) f[i]=(f[i]+1ll*f[i-j]*a[j])%mod;
		f[i]=(f[i]+g[1][i])%mod;//!!!
	}
	//for(int i=K;i<=n;++i) for(int j=1;j<=K;++j) f[i]=(f[i]+1ll*f[i-j]*a[j])%mod;
	//return f[n];

	for(int i=0;i<K;++i) m[K][i]=a[K-i];
	for(int i=K+1;i<=(K-1<<1);++i)
	{
		for(int j=1;j<=K;++j) m[i][j]=m[i-1][j-1];
		for(int j=0;j<K;++j) m[i][j]=(m[i][j]+1ll*m[i][K]*a[K-j])%mod;
		m[i][K]=0;
	}
	S[0]=1; A[1]=1;
	for(;n;n>>=1,mul(A,A,K)) if(n&1) mul(S,A,K);
	int fn=0;
	for(int i=0;i<K;++i) fn=(fn+1ll*f[i]*S[i])%mod;
	return fn;
}

int K,n;
int main()
{
#ifndef ONLINE_JUDGE
	freopen("pool.in","r",stdin);freopen("pool.out","w",stdout);
#endif
	n=read(); K=read(); q=1ll*read()*qpow(read(),mod-2)%mod;
	pwq[0]=1;
	for(int i=1;i<=K;++i) pwq[i]=1ll*pwq[i-1]*q%mod;
	printf("%d
",(solve(K,n)-solve(K-1,n)+mod)%mod);
	return 0;
}