【题解】[CTS2019]随机立方体

题目链接：[CTS2019] 随机立方体

题意：

有一个 $n\times m\times l$ 的立方体，立方体中每个格子上都有一个数，如果某个格子上的数比三维坐标至少有一维相同的其他格子上的数都要大的话，我们就称它是极大的。
现在将 $1\sim n\times m\times l$ 这 $n\times m\times l$ 个数等概率随机填入 $n\times m\times l$ 个格子（即任意数字出现在任意格子上的概率均相等），使得每个数恰出现一次，求恰有 $k$ 个极大的数的概率。
$1\le T\le 10, 1\le n\le 5\times 10^6$。

鸽这题不知道鸽了多久，现在来补坑

发现恰好有 $k$ 个极大的不好算，考虑二项式反演之后变成钦定 $k$ 个极大的位置。

设 $f_i$ 是钦定 $i$ 个极大位置的概率，$g_i$ 表示有序地选出 $i$ 个极大位置的方案数。由于显然极大的位置三维坐标都不相同，所以：

$$ g_i=\prod_{j=0}^{i-1} (n-j)(m-j)(l-j) $$

也就是直接给每个位置选三维的坐标。

设 $c_i$ 表示选出 $i$ 个极大位置之后，与这些极大位置有关系的位置的个数。有关系指与至少一个极大位置的至少一维坐标相同。

考虑拿总数减去三维都不和极大位置相同的位置个数，因此：

$$ c_i=nml-(n-i)(m-i)(l-i) $$

设 $h_i$ 表示有序地选出 $i$ 个极大位置之后，有多少种方案给 $c_i$ 个位置分配标号，使得恰好极大位置都是极大的。

考虑从大到小依次填入数字。首先，显然最大的数字一定要填入钦定最大的极大位置，然后便只剩下 $i-1$ 个极大位置。这指引我们将问题转变为 $h_{i-1}$ 的子问题。在 $h_{i-1}$ 中，共 $c_{i-1}$ 个位置有用，$h_i$ 中有 $c_i$ 个，且已经填入了一个最大值，因此还剩下 $(c_i-c_{i-1}-1)$ 个数字可以随便填数字（因为他们和 $h_{i-1}$ 的子问题也没有联系了），因此：

$$ h_i=h_{i-1}\times \dfrac{(c_i-1)!}{c_{i-1}!}=\prod_{j=1}^i \dfrac{(c_j-1)!}{c_{j-1}!} $$

则 $f_i$ 满足式子：

$$ f_i=\dfrac{1}{(nml)!}\binom{nml}{c_i}(nml-c_i)!\times g_i\times h_i $$

表示先选出 $c_i$ 个标号用于极大位置和与极大位置相关的位置，剩下位置任意排列；再选出 $i$ 个极大位置，并在上面合法地填入数字。因为是概率所以要乘 $\dfrac{1}{(nml)!}$

化简一下，可得：

$$ f_i=g_i\prod_{j=1}^i \dfrac{1}{c_j} $$

其中 $\prod_{j=1}^i \dfrac{1}{c_j}$ 可以通过类似阶乘求逆元的方式做到 $\mathcal{O}(n)$ 求解。则答案为：

$$ \text{ans}_k=\sum_{i=k}^{\min(n,m,l)} \binom{i}{k}(-1)^{i-k}f_i $$

时间复杂度 $\mathcal{O}(T\times n)$。

//Code By CXY07
#include<bits/stdc++.h>
using namespace std;

//#define FILE
#define int long long
#define file(FILENAME) freopen(FILENAME".in", "r", stdin), freopen(FILENAME".out", "w", stdout)
#define randint(l, r) (rand() % ((r) - (l) + 1) + (l))
#define LINE() cout << "LINE = " << __LINE__ << endl
#define debug(x) cout << #x << " = " << x << endl
#define abs(x) ((x) < 0 ? (-(x)) : (x))
#define inv(x) qpow((x), mod - 2)
#define lowbit(x) ((x) & (-(x)))
#define ull unsigned long long
#define pii pair<int, int>
#define LL long long
#define mp make_pair
#define pb push_back
#define scd second
#define vec vector
#define fst first
#define endl '\n'

const int MAXN = 5e6 + 10;
const int INF = 2e9;
const double eps = 1e-6;
const double PI = acos(-1);
//const int mod = 1e9 + 7;
const int mod = 998244353;
//const int G = 3;
//const int base = 131;

int T, n, m, l, k, V, lim, Ans;
int fac[MAXN], ifac[MAXN];
int g[MAXN], c[MAXN], ic[MAXN];

template<typename T> inline bool read(T &a) {
    a = 0; char c = getchar(); int f = 1;
    while(c < '0' || c > '9') {if(c == '-') f = -1; c = getchar();}
    while(c >= '0' && c <= '9') {a = a * 10 + (c ^ 48); c = getchar();}
    return a *= f, true;
}

template<typename A, typename ...B>
inline bool read(A &x, B &...y) {return read(x) && read(y...);}

int qpow(int x, int b) {
    int res = 1;
    for(; b; b >>= 1, (x *= x) %= mod) if(b & 1) (res *= x) %= mod;
    return res;
}

int C(int x, int y) {
    if(x < y) return 0;
    return fac[x] * ifac[y] % mod * ifac[x - y] % mod;
}

void init() {
    fac[0] = ifac[0] = 1;
    for(int i = 1; i < MAXN; ++i) fac[i] = fac[i - 1] * i % mod;
    ifac[MAXN - 1] = inv(fac[MAXN - 1]);
    for(int i = MAXN - 2; i; --i) ifac[i] = ifac[i + 1] * (i + 1) % mod;
}

void solve() {
    read(n), read(m), read(l), read(k); 
    V = n * m % mod * l % mod, Ans = 0;
    lim = min(n, min(m, l)); g[0] = ic[0] = 1;
    for(int i = 1; i <= lim; ++i) {
        g[i] = g[i - 1] * (n - i + 1) % mod * (m - i + 1) % mod * (l - i + 1) % mod;
        c[i] = (V - (n - i) * (m - i) % mod * (l - i) % mod + mod) % mod;
        ic[i] = ic[i - 1] * c[i] % mod;
    }
    ic[lim] = inv(ic[lim]);
    for(int i = lim - 1; i; --i) ic[i] = ic[i + 1] * c[i + 1] % mod;
    for(int i = k, opt = 1; i <= lim; ++i) {
        (Ans += C(i, k) * opt % mod * g[i] % mod * ic[i] % mod) %= mod;
        opt = (mod - opt) % mod;
    }
    printf("%lld\n", Ans);
}

signed main () {
#ifdef FILE
    freopen(".in", "r", stdin);
    freopen(".out", "w", stdout);
#endif
    read(T); init();
    while(T--) solve();
    return 0;
}