Address

洛谷P3824
BZOJ4944
UOJ#316
LOJ#2304

Solution…

一、差分 容斥

要限制最大值恰好为一个定值往往是不好做的。
所以考虑容 (cha) 斥 (fen) ，把询问 $=K$ 拆成 $\le K$ 和 $\le K-1$ ，这样就转化成了出现过的每一个安全子矩形面积都不超过 $K$ （或不超过 $K-1$ ）。
为了方便，下面我们只讨论每一个安全子矩形面积都不超过 $K$ 的求法。
（求 $\le K$ 和求 $\le K-1$ 的方法大致相同，但是要特殊处理 $\le 0$ 和 $\le-1$ 的情况，否则在 UOJ 上会被 hack ）
$\le 0$ 的答案为：
$(1-q)^n$
即最下面 $n$ 个格子全部不安全。
$\le -1$ 的答案就不用解释了，为 $0$ 。

二、 $g[i][j]$ 和 $s[i][j]$

状态 $g[i][j]$ 的定义为底边长为 $i$ ，高为 $1001$ ，最下面恰好 $j$ 行全部安全并且这个底边长为 $i$ 高为 $1001$ 的矩形内不存在任意安全子矩形面积大于 $K$ 的概率。
$s[i][j]$ 表示底边长为 $i$ ，高为 $1001$ ，最下面至少 $j$ 行全部安全并且不存在任意安全子矩形面积大于 $K$ 的概率。相当于是 $g[i][j]$ 的后缀和。
设 $a[k]$ 为第 $k$ 列的最下面有多少个安全的格子。
那么显然第 $l$ 列到第 $r$ 列内的最大安全子矩形为 $r-l+1$ 乘上 $a$ 在区间 $[l,r]$ 内的最小值。
这涉及到所有区间的最小值，故可以采用笛卡尔树的思想，找到最小值并向最小值的两边分治处理。
对于 $g[i][j]$ ，我们先枚举 $k\in[0,i-1]$ 表示 $a$ 中最左的最小值在 $k+1$ 位置。显然 $a[k+1]$ 应该为 $j$ （ $g$ 的定义为最小值恰好为 $j$ ）。
则 $a[1\dots k]$ 必须严格大于 $a[k+1]$ 即 $j$ （因为 $a[k+1]$ 是最左的最小值），
即 $s[k][j+1]$ 。
$a[k+2\dots i]$ 必须大于等于 $a[k+1]=j$ ，即 $s[i-1-k][j]$ 。
$a[k]=j$ 的概率就是下面 $j$ 个格子全部安全而第 $j+1$ 个格子不安全的概率。
即 $q^j(1-q)$ 。
得出转移：
$g[i][j]=\sum_{k=0}^{i-1}q^j\times(1-q)\times s[k][j+1]\times s[i-1-k][j]$
$s$ 是 $g$ 的后缀和：
$s[i][j]=s[i][j+1]+g[i][j]$
边界：
$s[0][0...K+1]=1$
复杂度：看上去是 $O(K^3)$ ，但是任何一个 $g[i][j]\ne 0$ 都必须满足 $i\times j\le K$ ，所以 $j$ 的上界只有 $\lfloor\frac Ki\rfloor$ 。
所以复杂度：
$K\sum_{i=1}^K\lfloor\frac Ki\rfloor=O(K^2\log K)$

三、 $f[i]$

$f[i]$ 的定义是底边长为 $i$ 高为 $1001$ ，任意安全子矩形面积都不超过 $K$ 的概率。
分两种情况：
（1）第 $i$ 列最下面的格子不安全。
$f[i]+=(1-q)\times f[i-1]$
（2）第 $[i-j+1,i]$ 列最下面的格子安全且产生的最大安全子矩形面积不超过 $K$ ，而第 $i-j$ 列的格子不安全。
$f[i]+=(1-q)\times s[j][1]\times f[i-j-1]$
得出转移：
$f[i]=\sum_{j=1}^{K+1}(1-q)\times s[j-1][1]\times f[i-j]$
边界 $f[0]=1$ 。
当然，如果 $i-j<0$ 就不要转移。
这是一个常系数线性递推式。
我们可以先暴力计算出 $f[0\dots K]$ ，然后利用快速幂对特征多项式取模得到答案。
关于常系数线性递推可以右转：蒟蒻的 blog
复杂度 $O(K^2\log n)$ 。

四、Other

注意到 $g$ 的转移是卷积形式。
如果利用多项式求逆进行 $g$ 的转移，并且利用 NTT 进行多项式取模，
那么复杂度可以优化到 $O(K\log K\log n)$ 。

Code

#include <cmath>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#define For(i, a, b) for (i = a; i <= b; i++)
#define Rof(i, a, b) for (i = a; i >= b; i--)
using namespace std;

inline int read()
{
	int res = 0; bool bo = 0; char c;
	while (((c = getchar()) < '0' || c > '9') && c != '-');
	if (c == '-') bo = 1; else res = c - 48;
	while ((c = getchar()) >= '0' && c <= '9')
		res = (res << 3) + (res << 1) + (c - 48);
	return bo ? ~res + 1 : res;
}

const int N = 1007, M = N << 1, ZZQ = 998244353;
int n, K, p, q[N], f[N][N], g[N], h[M], a[M], tmp[M], s[N][N], st[N];

int qpow(int a, int b)
{
	int res = 1;
	while (b)
	{
		if (b & 1) res = 1ll * res * a % ZZQ;
		a = 1ll * a * a % ZZQ;
		b >>= 1;
	}
	return res;
}

void prod(int *a, int *b, int K)
{
	int i, j;
	memset(tmp, 0, sizeof(tmp));
	For (i, 0, K) For (j, 0, K)
		tmp[i + j] = (tmp[i + j] + 1ll * a[i] * b[j] % ZZQ) % ZZQ;
}

void xmod(int *a, int n, int K)
{
	int i, j;
	Rof (i, n, K + 1)
	{
		int tmp = a[i];
		For (j, i - K - 1, i)
			a[j] = (a[j] - 1ll * tmp * g[j - i + K + 1] % ZZQ + ZZQ) % ZZQ;
	}
}

int jiejuediao(int K)
{
	if (K == -1) return 0;
	if (K == 0) return qpow((1 - p + ZZQ) % ZZQ, n);
	int i, j, k;
	q[0] = 1;
	For (i, 1, K) q[i] = 1ll * q[i - 1] * p % ZZQ;
	For (i, 0, K + 1) For (j, 0, K + 1) f[i][j] = s[i][j] = 0;
	For (j, 0, K + 1) s[0][j] = 1;
	For (i, 1, K) Rof (j, K / i, 0)
	{
		For (k, 0, i - 1) f[i][j] = (f[i][j] +
			1ll * s[k][j + 1] * s[i - 1 - k][j] % ZZQ
			* q[j] % ZZQ * (1 - p + ZZQ) % ZZQ) % ZZQ;
		s[i][j] = (s[i][j + 1] + f[i][j]) % ZZQ;
	}
	g[K + 1] = 1;
	For (i, 0, K)
		g[i] = (ZZQ - 1ll * s[K - i][1] * (1 - p + ZZQ) % ZZQ) % ZZQ;
	memset(a, 0, sizeof(a)); memset(h, 0, sizeof(h));
	a[h[0] = 1] = 1;
	int tm = n;
	while (tm)
	{
		if (tm & 1)
		{
			prod(h, a, K);
			For (i, 0, K << 1) h[i] = tmp[i];
			xmod(h, K << 1, K);
		}
		prod(a, a, K);
		For (i, 0, K << 1) a[i] = tmp[i];
		xmod(a, K << 1, K);
		tm >>= 1;
	}
	int res = 0;
	st[0] = 1;
	For (i, 1, K)
	{
		st[i] = s[i][1];
		For (j, 0, i - 1)
			st[i] = (st[i] + 1ll * st[j] * (1 - p + ZZQ) % ZZQ
			* s[i - j - 1][1] % ZZQ) % ZZQ;
	}
	For (i, 0, K) res = (res + 1ll * h[i] * st[i] % ZZQ) % ZZQ;
	return res;
}

int main()
{
	n = read(); K = read();
	int x = read(), y = read();
	p = 1ll * x * qpow(y, ZZQ - 2) % ZZQ;
	cout << (jiejuediao(K) - jiejuediao(K - 1) + ZZQ) % ZZQ << endl;
	return 0;
}