洛谷P5400 [CTS2019]随机立方体(计数)
题目大意
有一个 (n imes m imes l) 的立方体,立方体中每个格子上都有一个数,如果某个格子上的数比三维坐标至少有一维相同的其他格子上的数都要大的话,我们就称它是极大的。
现在将 (1sim n imes m imes l) 这 (n imes m imes l) 个数等概率随机填入 (n imes m imes l) 个格子(即任意数字出现在任意格子上的概率均相等),使得每个数恰出现一次,求恰有 (k) 个极大的数的概率。答案对 (998244353)(一个质数)取模。
数据范围
对于 (10\%) 的数据,(n,mle 2),(lle 3),(k=1)。
对于 (30\%) 的数据,(n,m,l,kle 12)。
对于 (40\%) 的数据,(n,m,lle 100)。
对于 (50\%) 的数据,(n,m,lle 1000)。
对于 (60\%) 的数据,(n,m,lle 100000),其中有占全部数据 (30\%) 的数据保证 (k=1)。
对于 (80\%) 的数据,(n,m,lle 1000000),其中有占全部数据 (40\%) 的数据保证 (k=1)。
对于 (100\%) 的数据,(1le n,m,lle 5000000),(1le kle 100),(1le Tle 10)。
其中有 (50\%) 的数据保证 (k=1)。
解题思路
一看就是二项式反演,先放个式子
其中 (F(i)) 表示至少有 i 个极大数(设为 A 点)的方案数
(B(i)) 表示与任意极大数至少一维相同的个数,即和极大数有限制关系的点(设为 B 点)个数
(G(i)) 表示选出来 i 个三维都不相同的极大数位置的方案数
(H(i)) 表示填好极大数和和它有限制关系的点的方案数
我们让极大数从小到大填即可,填的过程中,我们会发现当前极大数和之前的 B 点有维相同的话,当前极大数显然大于 B 点,那么剩下的点由于没有限制关系随便填即可
然后要将 (H(i)) 乘上阶乘来还原
因为我们最后要算的是概率,所以我先把 N! 省了
最后一波二项式反演即可
zbk 提供了推 H 的另一种思路,让当前最大点向次大点连边,然后向和之前极大点没有关系的 B 点连边,不难发现是一棵树的形态,求大小排列个数也就是求拓扑序的个数,拓扑序个数有如下公式
这个式子比较好证,不再赘述
#include <queue>
#include <vector>
#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#define MP make_pair
#define ll long long
#define fi first
#define se second
using namespace std;
template <typename T>
void read(T &x) {
x = 0; bool f = 0;
char c = getchar();
for (;!isdigit(c);c=getchar()) if (c=='-') f=1;
for (;isdigit(c);c=getchar()) x=x*10+(c^48);
if (f) x=-x;
}
template<typename F>
inline void write(F x, char ed = '
')
{
static short st[30];short tp=0;
if(x<0) putchar('-'),x=-x;
do st[++tp]=x%10,x/=10; while(x);
while(tp) putchar('0'|st[tp--]);
putchar(ed);
}
template <typename T>
inline void Mx(T &x, T y) { x < y && (x = y); }
template <typename T>
inline void Mn(T &x, T y) { x > y && (x = y); }
const int N = 5000005;
const int P = 998244353;
ll inv[N], jie[N], B[N], m, n, l, Inv, T, k, ans;
ll C(int n, int m) {
return jie[n] * inv[m] % P * inv[n - m] % P;
}
ll F(int k) {
ll t = C(n, k) * C(m, k) % P * C(l, k) % P;
t = t * jie[k] % P * jie[k] % P * jie[k] % P;
t = t * Inv % P * B[k+1] % P;
return t;
}
ll fpw(ll x, ll mi) {
ll res = 1;
for (; mi; mi >>= 1, x = x * x % P)
if (mi & 1) res = res * x % P;
return res;
}
int main() {
jie[0] = jie[1] = inv[0] = inv[1] = 1;
for (int i = 2;i <= 5000000; i++)
jie[i] = jie[i-1] * i % P,
inv[i] = inv[P % i] * (P - P / i) % P;
for (int i = 2;i <= 5000000; i++)
inv[i] = inv[i-1] * inv[i] % P;
for (read(T); T; T--) {
read(n), read(m), read(l), read(k); int lim = min(min(n, m), l);
if (k > lim) { puts("0"); continue; } ll all = n * m % P * l % P;
for (int i = 1;i <= lim; i++)
B[i] = (all + (P - n + i) * (m - i) % P * (l - i)) % P;
B[lim + 1] = 1;
for (int i = lim - 1;i >= 1; i--) B[i] = B[i] * B[i+1] % P;
Inv = fpw(B[1], P - 2), ans = 0;
for (int i = k;i <= lim; i++) {
ll t = F(i) * C(i, k) % P;
if ((i - k) & 1) ans -= t;
else ans += t;
}
write((ans % P + P) % P);
}
return 0;
}