Description：

高为 $1001$ 宽为 $n$ 的长方形网格，每个位置有 $p$ 的概率为 $0$ ，要求选出的全零矩形与下边界相连，问能选出的最大矩形面积恰好为 $k$ 的概率。

$1\leqslant n\leqslant 10^9,1\leqslant k\leqslant 1000$

Solution：

注意到虽然 $n$ 非常大，但是 $k$ 不是特别大，因此我们可以设计一个复杂度只与 $k$ 有关的辅助 $DP$ 来预处理一些东西，我们可以设计这样一个 $DP$ ，设 $dp[i][j]$ 表示底部宽为 $i$ （高为 $1001$ ），在底下 $j$ 行都是白色，在第 $j+1$ 行有至少一个黑色块，并且所有和底面相连的白色矩形面积都不超过 $k$ 的概率，转移的话可以枚举高度在 $j+1$ 处的第一个矩形，那么这个黑色块的前面高度必须大于 $j$ ，后面至少大于等于 $j$ ： $$ dp[i][j]=[i\times j\leqslant k](1-p)p^{j-1}\sum_{x=1}^i\Bigl((\sum_{j'>j}dp[x-1][j'])\times(\sum_{j'\geqslant j}dp[i-x][j'])\Bigr) $$ 只要把 $dp$ 数组后缀和优化一下，也就是： $$ \begin{align} dp'[i][j]&=\sum_{j'\geqslant j}dp[i][j]\\ dp[i][j]&=[i\times j\leqslant k](1-p)p^j\sum_{x=1}^i(dp'[x-1][j+1]\times dp'[i-x][j]) \end{align} $$ 前面那个概率表示前 $j$ 行都是 $0$ ( $p^j$ )第 $j+1$ 行为 $1$ ( $1-p$ )，后面就是枚举第 $j+1$ 行第一个黑色块的位置。

这个转移是 $O(k^2)$ 的，因为第一维个数为 $k$ ，第二维个数为 $k/i$ ，第三维为 $i$ 。

然后设 $f[n]$ 表示宽为 $n$ 时候的答案，我们枚举最下面最后一段白色的长度，也就是最下一行最后一个黑色块的位置，可以得到转移： $$ f[n]=\sum_{i=0}^kf[n-i-1]\times dp'[i][1]\times (1-p)=\sum_{i=1}^{k+1}f[n-i]\times dp'[i-1][1]\times (1-p) $$ 发现这就是一个常系数齐次线性递推，于是用多项式取模优化即可。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<cctype>

#include<cstring>

using namespace std;

int n,k,x,y;

#define MAXK 2010

int dp[MAXK][MAXK],dp_[MAXK][MAXK];

int p;

#define MOD 998244353

int power(int a,int b)

{

    int res = 1;

    while(b > 0)

    {

        if(b & 1)res = 1ll * res * a % MOD;

        a = 1ll * a * a % MOD;

        b = b >> 1;

    }

    return res;

}

#define MAXN 10000001

int f[MAXN];

int g[MAXK];

int K;

struct matrix

{

	int a[MAXK];

	matrix(){memset(a,0,sizeof(a));}

	int& operator [](int x){return a[x];}

	int operator [](int x)const{return a[x];}

	friend matrix operator * (matrix a,matrix b)

	{

		matrix res;

		for(int i = 0;i < K;++i)

			for(int j = 0;j < K;++j)

				(res[i + j] += 1ll * a[i] * b[j] % MOD) %= MOD;

		for(int i = 2 * K - 2;i >= K;res[i--] = 0)

			for(int j = 1;j <= K;++j)

				(res[i - j] += 1ll * res[i] * g[j] % MOD) %= MOD;

		return res;

	}

};

matrix power(matrix a,int b)

{

	matrix res;res[0] = 1;

	while(b > 0)

	{

		if(b & 1)res = res * a;

		a = a * a;

		b = b >> 1;

	}

	return res;

}

int solve(int k)

{

    memset(dp,0,sizeof(dp));

    memset(dp_,0,sizeof(dp_));

    for(int j = 0;j <= k + 1;++j)dp_[0][j] = 1;

    for(int i = 1;i <= k + 1;++i)

    {

        for(int j = 0;i * j <= k;++j)

        {

            int sum = 0;

            for(int x = 1;x <= i;++x)

            {

                sum = (sum + 1ll * dp_[x - 1][j + 1] * dp_[i - x][j] % MOD) % MOD;

            }

            dp[i][j] = 1ll * (1 - p + MOD) * power(p,j) % MOD * sum % MOD;

            for(int j_ = 0;j_ <= j;++j_)dp_[i][j_] = (dp_[i][j_] + dp[i][j]) % MOD;

        }

    }

    ++k;K = k;

    for(int i = 0;i < k;++i)f[i] = dp_[i][0];

    for(int i = 1;i <= k;++i)g[i] = 1ll * dp_[i - 1][1] % MOD * (1 - p + MOD) % MOD;

    /*for(int x = k;x <= n;++x)

    {

        f[x] = 0;

        for(int i = 1;i <= k;++i)

        {

            f[x] = (f[x] + 1ll * f[x - i] * g[i] % MOD) % MOD;

        }

    }*/

   	matrix m;m[1] = 1;

    m = power(m,n);

    int ans = 0;

    for(int i = 0;i < k;++i)(ans += 1ll * m[i] * f[i] % MOD) %= MOD;

    return ans;

}

int main()

{

    scanf("%d%d%d%d",&n,&k,&x,&y);

    p = 1ll * x * power(y,MOD - 2) % MOD;

    cout << (solve(k) - solve(k - 1) + MOD) % MOD;

    return 0;

}