有一个 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 取模。

1n,m,l5×106,1k1001 \le n, m, l \le 5 \times 10^6, 1 \le k \le 100

LG P5400

Solution

以下题解修改自ViXbob的题解

显然容斥,从恰好转化为钦定。

f[i]f[i]表示钦定有ii个极大值的方案数,T=n×m×lT = n \times m \times l。则有:

f[i]=(Tg[i])b[i]×h[i]×(Tg[i])!f[i]=\binom{T}{g[i]}b[i]\times h[i]\times (T - g[i])!
  • g[i]g[i]表示和ii个极大值中任意一个至少有一维坐标相同的点的个数。

  • b[i]b[i]表示选出ii个三维坐标都不相同的点的方案数。(ii个极大值)

  • h[i]h[i]表示将g[i]g[i]个数字分配给g[i]g[i]的个点的合法分配的方案数。(合法指的是ii个极大值可以同时存在)

考虑这三个值如何求:

首先,每选择出来一个点(x0,y0,z0)(x_0,y_0,z_0)并钦定它为极大的,则x=x0,y=y0,z=z0x=x_0,y=y_0,z=z_0这三个平面上的点都不能再作为极大的点了,相当于立方体从n×m×ln \times m \times l变成了(n1)×(m1)×(l1)(n-1)\times (m-1)\times (l-1)

因此

g[i]=n×m×l(ni)×(mi)×(li)b[i]=1i!j=0i1(nj)×(mj)×(lj) \begin{aligned} g[i]&=n\times m\times l-(n-i)\times (m-i)\times (l-i) \\ b[i]&=\frac{1}{i!}\prod_{j=0}^{i-1}(n-j)\times(m-j)\times(l-j) \end{aligned}

于是只剩下h[i]h[i]

设需要比第ii个极大的点要小的点上数字集合为SiS_ia[i]a[i]表示第ii个极大的点上的数字。

我们先钦定ii个极大的点上的数字的大小关系,最后再乘i!i!。由数字从小到大的顺序依次考虑这个ii个极大的点。考虑加入第jj个极大的点后方案数的变化:

因为SjS_j和前j1j-1极大的点的集合相交的部分已经确定,且它们一定比a[i]a[i]小(因为aa数组递增),所以只需要选SiS_i和前i1i-1个集合的并集的不交的部分。设其大小为RR,则有:

h[i]=h[i1]×Pg[i]1Rh[i]=h[i-1]\times \text{P}_{g[i] - 1}^{R}

P\text{P}表示排列数

解释一下这个式子,我们就是要在g[i]1g[i] - 1个数中选出RR个用来放SiS_i的这RR个数,要排列。而剩下的数就作为前i1i - 1部分的数,因为已经在h[i1]h[i - 1]考虑过它们的顺序了,所以不需要乘别的东西

g[i]g[i]要减一是因为这所有g[i]g[i]个数中,最大的数已经被a[i]a[i]钦定了

并且有R=g[i]g[i1]1R=g[i]-g[i-1]-1。所以:

h[i]=h[i1]×(g[i]1)!g[i1]!=j=1i(g[j]1)!g[j1]!\begin{aligned}h[i]&=h[i-1] \times \frac{(g[i]-1)!}{g[i-1]!}\\&=\prod_{j=1}^i\frac{(g[j]-1)!}{g[j-1]!}\end{aligned}

然后把所有的东西都代回去(记得乘i!i!),则有:

f[i]=(Tg[i])b[i]×i!j=1i(g[j]1)!g[j1]!×(Tg[i])!=T!(Tg[i])!g[i]!×i!j=1i(g[j]1)!g[j1]!×(Tg[i])!×b[i]\begin{aligned}f[i]&=\binom{T}{g[i]}b[i]\times i!\prod_{j=1}^i\frac{(g[j]-1)!}{g[j-1]!}\times (T - g[i])!\\&=\frac{T!}{(T-g[i])!g[i]!}\times i!\prod_{j=1}^i\frac{(g[j]-1)!}{g[j-1]!}\times (T - g[i])!\times b[i]\end{aligned}

发现我们最后算概率的时候要乘上1T!\frac{1}{T!}并且式子的最后有一个(Tg[i])!(T-g[i])!,则有:

f[i]=i!g[i]!×j=1i(g[j]1)!g[j1]!×b[i]=i!×b[i]j=1i(g[j]1)!g[j]!=i!×b[i]j=1i1g[j]\begin{aligned}f[i]&=\frac{i!}{g[i]!}\times \prod _{j=1}^i\frac{(g[j]-1)!}{g[j-1]!}\times b[i]\\&=i!\times b[i]\prod_{j=1}^i\frac{(g[j]-1)!}{g[j]!}\\&=i!\times b[i]\prod_{j=1}^i\frac{1}{g[j]}\end{aligned}

最后二项式反演一下有:

ans=i=kmin(n,m,l)(1)ik(ik)f[i]=i=kmin(n,m,l)(1)ik(ik)i!×b[i]j=1i1g[j] \begin{aligned} ans &= \sum_{i=k}^{\min(n,m,l)}(-1)^{i-k}\binom{i}{k}f[i]\\ &= \sum_{i=k}^{\min(n,m,l)}(-1)^{i-k}\binom{i}{k}i!\times b[i]\prod_{j=1}^i\frac{1}{g[j]} \end{aligned}

最后需要线性预处理g[j]g[j]前缀积的逆元(求g[j]g[j]的逆元也可以),是个很经典的trick,类似线性预处理阶乘的逆,处理出gg数组的前缀积,再用一次快速幂求最后一个数的逆元,再依次乘回来即可

复杂度O(min(n,m,l))O(\min(n,m,l))

Code

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
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
#include <iostream>
#include <ctime>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <vector>
#include <climits>
#include <queue>
#include <set>
#include <cmath>
#include <map>
#include <bitset>
#include <fstream>
#include <tr1/unordered_map>
#include <assert.h>

#define x first
#define y second
#define y0 Y0
#define y1 Y1
#define mp make_pair
#define pb push_back
#define DEBUG(x) cout << #x << " = " << x << endl;

using namespace std;

typedef long long LL;
typedef pair <int, int> pii;

template <typename T> inline int Chkmax (T &a, T b) { return a < b ? a = b, 1 : 0; }
template <typename T> inline int Chkmin (T &a, T b) { return a > b ? a = b, 1 : 0; }
template <typename T> inline T read ()
{
T sum = 0, fl = 1; char ch = getchar();
for (; !isdigit(ch); ch = getchar()) if (ch == '-') fl = -1;
for (; isdigit(ch); ch = getchar()) sum = (sum << 3) + (sum << 1) + (ch - '0');
return sum * fl;
}

inline void proc_status ()
{
ifstream t ("/proc/self/status");
cerr << string (istreambuf_iterator <char> (t), istreambuf_iterator <char> ()) << endl;
}

const int MAXN = 5e6;
const int MOD = 998244353;

int N, M, L, K;

namespace MATH
{
int fac[MAXN + 5], ifac[MAXN + 5], inv[MAXN + 5];

inline void ADD (int &a, int b) { if ((a += b) >= MOD) a -= MOD; }

inline int Pow (int a, int b)
{
int ans = 1;
for (int i = b; i; i >>= 1, a = (LL) a * a % MOD) if (i & 1) ans = (LL) ans * a % MOD;
return ans;
}

inline void init (int n = MAXN + 1)
{
fac[0] = 1; for (int i = 1; i <= n; ++i) fac[i] = (LL) fac[i - 1] * i % MOD;
ifac[n] = Pow (fac[n], MOD - 2); for (int i = n - 1; i >= 0; --i) ifac[i] = (LL) ifac[i + 1] * (i + 1) % MOD;
for (int i = 0; i < n; ++i) inv[i + 1] = (LL) fac[i] * ifac[i + 1] % MOD;
}

inline int C (int n, int m) { if (n < m || n < 0 || m < 0) return 0; return (LL) fac[n] * ifac[m] % MOD * ifac[n - m] % MOD; }
}

using namespace MATH;

inline void get_inv (int *A, int *IA, int n)
{
int sum = 1;
for (int i = 1; i <= n; ++i) sum = (LL) sum * A[i] % MOD;
IA[n] = Pow (sum, MOD - 2); for (int i = n - 1; i >= 0; --i) IA[i] = (LL) IA[i + 1] * A[i + 1] % MOD;
}

inline void Solve ()
{
int lim = min (N, min (M, L));

static int G[MAXN + 5], IG[MAXN + 5];

for (int i = 1; i <= lim; ++i) G[i] = ((LL) N * M % MOD * L % MOD - (LL) (N - i) * (M - i) % MOD * (L - i) % MOD + MOD) % MOD;
get_inv (G, IG, lim);

static int B[MAXN + 5];

B[0] = 1;
for (int i = 1; i <= lim; ++i) B[i] = (LL) B[i - 1] * (N - i + 1) % MOD * (M - i + 1) % MOD * (L - i + 1) % MOD * inv[i] % MOD;

int ans = 0;
for (int i = K; i <= lim; ++i)
{
int coef = (LL) C (i, K) * fac[i] % MOD * B[i] % MOD * IG[i] % MOD;
// cout << G[i] << ' ' << IG[i] << endl;
if ((i - K) & 1) ADD (ans, MOD - coef);
else ADD (ans, coef);
}

cout << ans << endl;
}

inline void Input ()
{
N = read<int>(), M = read<int>(), L = read<int>(), K = read<int>();
}

int main ()
{

#ifdef hk_cnyali
freopen("A.in", "r", stdin);
freopen("A.out", "w", stdout);
#endif

int Te = read<int>();

MATH :: init ();
while (Te--)
{
Input ();
Solve ();
}

return 0;
}