題目鏈接
dp部分是問的別人的。\(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\)

然後令\(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|\)

最後，設\(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^2+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;
}

uoj316【NOI2017】泳池