「CTS 2019」随机立方体【容斥原理 + 概率与期望】

Time Limit: 12 Sec Memory Limit: 512 MB

Description

有一个 n×m×ln \times m \times l 的立方体,每个格子上都有一个数,如果某个格子上的数比三维坐标至少有一维相同的其他格子上的数都要大的话,我们就称它是极大的。

现在将 1n×m×l1 \sim n \times m \times ln×m×ln \times m \times l 个数等概率随机填入 n×m×ln \times m \times l 个格子,使得每个数均出现一次,求恰有 kk 个极大的数的概率。答案对 998244353998244353 取模。

Input

输入包含多组数据。第一行包含一个整数 TT,表示数组组数。

接下来 TT 行,每行包含四个整数 n,m,l,kn, m, l, k,表示一次询问。

Output

对于每次询问,输出一行一个整数,表示答案对 998244353998244353 取模的余数。

Sample Input

1
2
3
4
5
6
5
1 1 1 1
2 2 2 1
7 8 9 3
123 456 789 1
1000 1000 1000 10

Sample Output

1
2
3
4
5
1
142606337
736950806
246172965
189652652

Constraints

对于 30%30\% 的数据,满足 n,m,l,k12n, m, l, k \leqslant 12

对于 80%80\% 的数据,满足 n,m,l106n, m, l \leqslant 10^6

对于 100%100\% 的数据,满足 n,m,l5×106,n, m, l \leqslant 5 \times 10^6, k100,k \leqslant 100, T10T \leqslant 10

Solution

显然位置 (x,y,z)(x, y, z) 控制的集合是:

{(a,b,c)a=x or b=y or c=z}\big\{ (a, b, c) \mid a = x \text{ or } b = y \text{ or } c = z \big\}

恰好 kk 个位置较难计算,考虑用容斥计算强制 kk 个位置是极大值的概率。

设这些数为 a1,a2,,aka_1, a_2, \cdots, a_k (a1>a2>>ak)(a_1 > a_2 > \cdots > a_k),那么需要满足 aiaia_i \geqslant a_i 控制范围的所有数。但是 a1,a2,,aka_1, a_2, \cdots, a_k 控制范围有交,概率不独立,不能直接相乘。

可以将这个条件改为 aiai,ai+1,,aka_i \geqslant a_i, a_{i + 1}, \cdots, a_k 控制范围的并的所有数。这样处理之后概率就是独立的了。控制范围的并可以直接容斥计算。

f(i)f(i) 表示强制 ii 个位置是极大值的概率,g(i)g(i) 表示恰好 ii 个位置是极大值。由 二项式反演 可得:

f(i)=j=in(ji)g(j) g(i)=j=in(1)ji(ji)f(j)\begin{aligned} & f(i) = \sum_{j = i}^n \binom{j}{i} g(j) \\ \Rightarrow\ & g(i) = \sum_{j = i}^n (-1)^{j - i} \binom{j}{i} f(j) \end{aligned}

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
#include <bits/stdc++.h>
using namespace std;

const int maxn = 5000010, P = 998244353;
int T, n, m, l, k, v[maxn], fact[maxn], inv[maxn], f[maxn], g[maxn];

int A(int n, int m) { return 1LL * fact[n] * inv[n - m] % P; }
int C(int n, int m) { return 1LL * fact[n] * inv[m] % P * inv[n - m] % P; }

int qp(int x, int y) {
int z = 1;
for (; y; y >>= 1, x = 1LL * x * x % P) {
if (y & 1) z = 1LL * z * x % P;
}
return z;
}

int main() {
scanf("%d", &T);
for (int i = fact[0] = inv[0] = 1; i <= 5000000; i++) {
fact[i] = 1LL * i * fact[i - 1] % P;
}
inv[5000000] = qp(fact[5000000], P - 2);
for (int i = 4999999; i; i--) {
inv[i] = 1LL * (i + 1) * inv[i + 1] % P;
}
while (T--) {
scanf("%d %d %d %d", &n, &m, &l, &k);
if (m > l) swap(m, l);
if (n > m) swap(n, m);
int x = (1LL * n * m + 1LL * n * l + 1LL * m * l) % P, ans = 0;
for (int i = g[n] = 1; i <= n; i++) {
v[i] = (1LL * i * i % P * i % P + 1LL * i * x % P) % P;
v[i] = (v[i] - 1LL * i * i % P * (n + m + l) % P) % P;
g[n] = 1LL * g[n] * v[i] % P;
}
g[n] = qp(g[n], P - 2);
for (int i = n - 1; i; i--) {
g[i] = 1LL * v[i + 1] * g[i + 1] % P;
}
for (int i = k; i <= n; i++) {
f[i] = 1LL * A(n, i) * A(m, i) % P * A(l, i) % P * g[i] % P;
}
for (int i = k; i <= n; i++) {
if ((i - k) & 1) ans = (ans - 1LL * C(i, k) * f[i] % P) % P;
else ans = (ans + 1LL * C(i, k) * f[i] % P) % P;
}
printf("%d\n", (ans + P) % P);
}
return 0;
}