LOJ3119 随机正方体

神仙的二项式反演，中间推导有点麻烦。

题意：随机在三维 $n\times m \times l$ 的立方体中，每个格子随机放入 $1\sim n\times m\times l$ 使得每个数恰好出现一次。若三维坐标中至少有一个和当前坐标相同的坐标所填的数中，没有比当前数大的数，则称当前坐标是一个极大值。求恰好有 $k$ 个极大值的概率。$T(T\leq 10)$ 组数据，$n, m, l\leq 5\times 10 ^ 6, k\leq 100$。

首先声明，$k\leq 100$ 似乎没起到作用……

看到恰好，果断二项式反演，下面 $f(k)$ 表示钦定 $k$ 个极大值的方案数。显然 $k$ 不能超过 $\min\{n, m, l\}$，令其为 $lim$。则答案为：

$$ans = \sum_{i = k} ^ {lim} (-1) ^ {i - k} \binom ik f(i)$$

下面考虑如何计算 $f(k)$。

首先我们需要三维都选出 $k$ 个坐标，贡献就是 $n ^ {\underline k}m ^ {\underline k}l ^ {\underline k}$（注意是有序的）。与极大值无关的地方随便填，设为 $a_k = (n - k)(m - k)(l - k)$，为 $(nml) ^ {\underline {a_k}}$。

然后现在顺序这些都不需要考虑，我们现在只需要填入 $k$ 个最大值。

如果我们先填入 $k$ 个极大值中最大的那一个，填的时候我们将前 $k - 1$ 个所覆盖的（指至少有一维坐标相同）的先不管。这启示我们可以递推下去，假设这个函数为 $g(k)$。

我们填的时候，有影响的位置就应该是 $a_{k - 1} - a_k$，然后我们需要在 $nml - a_k - 1$ 个数（没有这个极大值，因为他一定选最大的那一个）中随意选出 $a_{k - 1} - a_k - 1$ 个数，然后除了极大值的位置固定以外，其他的又随意填，贡献就是：

$$\binom{nml - a_k - 1}{a_{k - 1} - a_k - 1}(a_{k - 1} - a_k - 1)!$$

通过这个，我们很容易从 $g(k - 1)$ 转移到 $g(k)$，而转移系数就是上面那个式子。

接下来就是把上面的式子一通合并（防止重名，题目中的 $k$ 更换为 $st$）：

$$\begin{aligned} ans =& \dfrac 1{nml!}\sum_{k = st} ^ {lim} (-1) ^ {k - st} \binom k{st} n ^ {\underline k}m ^ {\underline k} l ^ {\underline k} nml ^ {\underline{a_k}} g(k)\\ =& \dfrac 1{nml!}\sum_{k = st} ^ {lim} (-1) ^ {k - st} \binom k{st} n ^ {\underline k}m ^ {\underline k} l ^ {\underline k} nml ^ {\underline{a_k}} \prod_{i = 1} ^ k \dfrac{(nml - a_k - 1)!} {(nml - a_{k - 1})!} \\ =& \sum_{k = st} ^ {lim} (-1) ^ {k - st} \binom k{st} n ^ {\underline k} m ^ {\underline k} l ^ {\underline k} \prod_{i = 1} ^ k \dfrac{1}{(nml - a_k)!} \end{aligned}$$

这样我们线性预处理前缀逆元即可 $O(n)$ 计算。

void work()
{
	std::cin >> n >> m >> l >> k, lim = std::min({n, m, l});
	if (k > lim) return void(puts("0"));
	for (int i = 1; i <= lim; ++ i)
		adj(r[i] = (LL) n * m % Mod * l % Mod - (LL) (n - i) * (m - i) % Mod * (l - i) % Mod);
	int tot = 1;
	for (int i = 1; i <= lim; ++ i) tot = (LL) tot * r[i] % Mod;
	tot = qpow(tot), pre[lim] = tot;
	for (int i = lim - 1; i; -- i) pre[i] = (LL) pre[i + 1] * r[i + 1] % Mod;
	int res = 0;
	for (int i = k, op = 1; i <= lim; ++ i, op = Mod - op)
		res = (res + (LL) op * C(n, i) % Mod * C(m, i) % Mod * C(l, i)
			% Mod * qpow(fact[i], 3) % Mod * pre[i] % Mod * C(i, k)) % Mod;
	printf("%d\n", res);
}