【BZOJ4944】【UOJ316】【NOI2017】泳池

【题目链接】

• BZOJ
• UOJ

【思路要点】

• 考虑对于已知网格，如何计算最大的安全区域的面积。
• 我们会选用笛卡尔树。
• 因此，我们考虑枚举区间最小值的位置，进行笛卡尔树DP。
• 设\(dp_{i,j,0}(i*j≤k)\)表示当前考虑的部分宽度为\(j\)，且靠近沙滩的\(i*j\)个格子被确认是安全的情况下最大安全面积恰好为\(k\)的概率。
• 设\(dp_{i,j,1}(i*j≤k)\)表示当前考虑的部分宽度为\(j\)，且靠近沙滩的\(i*j\)个格子被确认是安全的情况下最大安全面积不超过\(k\)的概率。
• 枚举第\(i+1\)行第一次出现不安全的格子的位置\(pos\)，或者在\(i+1\)层未出现不安全的格子。
• 有$$dp_{i,j,0}=q^j*dp_{i+1,j,0}+\sum_{pos=1}^{j}(1-q)*q^{pos-1}*(dp_{i+1,pos-1,0}*dp_{i,j-pos,1}+dp_{i+1,pos-1,1}*dp_{i,j-pos,0}+dp_{i+1,pos-1,0}*dp_{i,j-pos,0})$$ $$dp_{i,j,1}=q^j*dp_{i+1,j,1}+\sum_{pos=1}^{j}(1-q)*q^{pos-1}*dp_{i+1,pos-1,1}*dp_{i,j-pos,1}$$
• 注意到\(i*j≤k\)，因此这个DP除了\(i=0\)时有\(O(N)\)个状态，其余的总状态数为\(\sum_{i=1}^{k}\frac{k}{i}\)。
• 因此其时间复杂度为\(O(N^2+\sum_{i=1}^{k}(\frac{k}{i})^2)=O(N^2+K^2)\)，可以得到70分。
• 考虑对于\(dp_{0,i,0}\)运行BM算法，我们发现这个数列存在一个长度不超过\(2*K\)的常系数线性递推式。
• 用Cayley-Hamilton定理优化递推即可。
• 时间复杂度\(O(K^2LogN)\)。

【代码】

#include<bits/stdc++.h>
using namespace std;
const int MAXN = 4505;
const int P = 998244353;
template <typename T> void chkmax(T &x, T y) {x = max(x, y); }
template <typename T> void chkmin(T &x, T y) {x = min(x, y); } 
template <typename T> void read(T &x) {
	x = 0; int f = 1;
	char c = getchar();
	for (; !isdigit(c); c = getchar()) if (c == '-') f = -f;
	for (; isdigit(c); c = getchar()) x = x * 10 + c - '0';
	x *= f;
}
template <typename T> void write(T x) {
	if (x < 0) x = -x, putchar('-');
	if (x > 9) write(x / 10);
	putchar(x % 10 + '0');
}
template <typename T> void writeln(T x) {
	write(x);
	puts("");
}
namespace CayleyHamilton {
	const int MAXN = 5005;
	int n, k;
	int a[MAXN], now[MAXN];
	int res[MAXN], h[MAXN];
	void times(int *x, int *y) {
		static int tmp[MAXN];
		memset(tmp, 0, sizeof(tmp));
		for (int i = 0; i <= k - 1; i++)
		for (int j = 0; j <= k - 1; j++)
			tmp[i + j] = (tmp[i + j] + 1ll * x[i] * y[j]) % P;
		for (int i = 2 * k - 2; i >= k; i--) {
			for (int j = 1; j <= k; j++)
				tmp[i - j] = (tmp[i - j] + 1ll * tmp[i] * a[j]) % P;
			tmp[i] = 0;
		}
		for (int i = 0; i <= k - 1; i++)
			x[i] = tmp[i];
	}
	void init(vector <int> &res, int *val) {
		k = res.size();
		for (int i = 1; i <= k; i++)
			a[i] = res[i - 1];
		for (int i = 1; i <= k * 2; i++)
			h[i] = val[i];
	}
	int query(int n) {
		if (n <= 2 * k) return h[n];
		memset(res, 0, sizeof(res));
		memset(now, 0, sizeof(now));
		n -= k; res[0] = now[1] = 1;
		while (n != 0) {
			if (n & 1) times(res, now);
			n >>= 1; times(now, now);
		}
		int ans = 0;
		for (int i = 0; i <= k - 1; i++)
			ans = (ans + 1ll * h[k + i] * res[i]) % P;
		return ans;
	}
}
namespace BerlekampMassey {
	const int MAXN = 5005;
	int n, val[MAXN], cnt, fail[MAXN], delta[MAXN];
	vector <int> ans[MAXN];
	int power(int x, int y) {
		if (y == 0) return 1;
		int tmp = power(x, y / 2);
		if (y % 2 == 0) return 1ll * tmp * tmp % P;
		else return 1ll * tmp * tmp % P * x % P;
	}
	void work() {
		ans[cnt = 0].clear();
		for (int i = 1; i <= n; i++) {
			int now = val[i];
			for (unsigned j = 0; j < ans[cnt].size(); j++)
				now = (now - 1ll * val[i - j - 1] * ans[cnt][j] % P + P) % P;
			delta[i] = now; if (now == 0) continue;
			fail[cnt] = i;
			if (cnt == 0) {
				ans[++cnt].clear();
				ans[cnt].resize(i);
				continue;
			}
			ans[++cnt].clear();
			ans[cnt].resize(i - fail[cnt - 2] - 1);
			int mul = 1ll * now * power(delta[fail[cnt - 2]], P - 2) % P;
			ans[cnt].push_back(mul);
			for (unsigned j = 0; j < ans[cnt - 2].size(); j++)
				ans[cnt].push_back(1ll * ans[cnt - 2][j] * (P - mul) % P);
			if (ans[cnt].size() < ans[cnt - 1].size()) ans[cnt].resize(ans[cnt - 1].size());
			for (unsigned j = 0; j < ans[cnt - 1].size(); j++)
				ans[cnt][j] = (ans[cnt][j] + ans[cnt - 1][j]) % P;
		}
	}
}
int n, k, p, q;
int dp[MAXN][MAXN][2];
int power(int x, int y) {
	if (y == 0) return 1;
	int tmp = power(x, y / 2);
	if (y % 2 == 0) return 1ll * tmp * tmp % P;
	else return 1ll * tmp * tmp % P * x % P;
}
int main() {
	freopen("UOJ316.in", "r", stdin);
	freopen("UOJ316.out", "w", stdout);
	read(n), read(k), read(p), read(q);
	p = 1ll * p * power(q, P - 2) % P;
	for (int i = 0; i <= k + 1; i++)
		dp[i][0][1] = 1;
	for (int i = k; i >= 0; i--)
	for (int j = 1, s = i; j <= n && j < MAXN && s <= k; j++, s += i) {
		int now = 1;
		for (int pos = 1; pos <= j; pos++) {
			int tmp = now * (P + 1ll - p) % P;
			if (s == k) {
				dp[i][j][0] = (1ll * dp[i + 1][pos - 1][1] * tmp % P * dp[i][j - pos][1] + dp[i][j][0]) % P;
				dp[i][j][1] = (1ll * dp[i + 1][pos - 1][1] * tmp % P * dp[i][j - pos][1] + dp[i][j][1]) % P;
			} else {
				dp[i][j][0] = (1ll * dp[i + 1][pos - 1][0] * tmp % P * dp[i][j - pos][1] + dp[i][j][0]) % P;
				dp[i][j][0] = (1ll * dp[i + 1][pos - 1][1] * tmp % P * dp[i][j - pos][0] + dp[i][j][0]) % P;
				dp[i][j][0] = ((P - 1ll) * dp[i + 1][pos - 1][0] % P * tmp % P * dp[i][j - pos][0] + dp[i][j][0]) % P;
				dp[i][j][1] = (1ll * dp[i + 1][pos - 1][1] * tmp % P * dp[i][j - pos][1] + dp[i][j][1]) % P;
			}
			now = 1ll * now * p % P;
		}
		dp[i][j][0] = (1ll * dp[i + 1][j][0] * now + dp[i][j][0]) % P;
		dp[i][j][1] = (1ll * dp[i + 1][j][1] * now + dp[i][j][1]) % P;
	}
	if (n < MAXN) {
		writeln(dp[0][n][0]);
		return 0;
	}
	BerlekampMassey :: n = MAXN - 1;
	for (int i = 1; i < MAXN; i++)
		BerlekampMassey :: val[i] = dp[0][i][0];
	BerlekampMassey :: work();
	CayleyHamilton :: init(BerlekampMassey :: ans[BerlekampMassey :: cnt], BerlekampMassey :: val);
	writeln(CayleyHamilton :: query(n));
	return 0;
}