题目链接

\(S=k\)可以拆成\(S\le k\)减去\(S\le 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)\)，注意\(i\le 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_1\lambda^{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('\n')
#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\n",(solve(K,n)-solve(K-1,n)+mod)%mod);
    return 0;
}