luogu5400 CTS2019 随机立方体

题意

$n\times m\times l$个格子的立方体，随机填入$1$到$n\times m\times l$的数。
一个格子如果比所有与他有最少一维相同的格子上面的数都大，那么它是极大的。
问恰有$k$个极大格子的概率，$mod998244353$，$T$组数据

数据范围

$n,m,l\le 5\times 10^6,k\le 100,T\le 10,12s$

解法

我这里有一个思考起来比较简单的做法。
其实一个随机排列，可以把它变成随机一个无限长的格子序列，每个格子等概率的随机，按照第一次出现的顺序依次填数字。这样是等价的。
随机格子就好办一些，可以看成三个序列，分别在值域$[1,n],[1,m],[1,l]$随机，记作$A,B,C$
容易发现对于极大的格子，设其为$(a,b,c)$，设$i$是最小的位置满足$A_i=a,B_i=b,C_i=c$。$i$就是$a$在$A$中，$b$在$B$中，$c$在$C$中第一次出现的位置。
不容斥还是没法算的，现在钦定$k$个位置，$A,B,C$在这些位置上都有新出现的数字。考虑所有钦定$k$个位置的方案，首先可以用$A_n^k{A}_m^k{A}_l^k(\frac{1}{nml}{)}^k$来固定这些位置上填的数。然后钦定两两位置之间的间隔，由于是无限序列，这个是可以随便定的，从$0$到$+\infty$，容易发现这相当于一个等比数列，第$i$个位置和第$i+1$个位置之间的位置，填的数满足条件的概率是$\frac{(n-k+i)(m-k+i)(l-k+i)}{nml}$这个就是公比，把这些等比数列求和的值相乘就得到了钦定$k$个位置所有情况的概率。
最后是容斥系数，需要注意的是要把一定是极大格子的第一个位置去掉，可以得到$(-1{)}^{(i-k)}(\genfrac{}{}{0ex}{}{i-1}{k-1})$的系数。
预处理阶乘和阶乘逆元，能够$O(1)$得到排列数和组合数。求等比数列的时候需要求$O(n)$个数的逆元，这个用luoguP5431的方法可以做到$O(n)$。
时间复杂度$O(Tn)$，基本上用不了$1.5s$就过了。。。

代码

#include<iostream>
#include<cstring>
#include<cassert>
#include<cmath>
#include<map>
#include<set>
#include<queue>
#include<stack>
#include<vector>
#include<time.h>
#include<bitset>
#include<cstdio>
#include<algorithm>
using namespace std;
#define REP(i,x,y) for(ll i=x;i<=y;i++)
#define rep(i,n) REP(i,1,n)
#define rep0(i,n) REP(i,0,n-1)
#define repG(i,x) for(ll i=pos[x];~i;i=e[i].next)
#define ll long long
#define db double
const ll N=5e6+507;
const ll INF=1e9+7;
const ll mod=998244353;
ll Pow(ll x,ll y){
	ll ans=1,now=x;
	while(y){
		if(y&1)ans=ans*now%mod;
		now=now*now%mod;
		y>>=1;
	}
	return ans;
}
ll T,m,n,l,k;
ll fac[N],inv[N],p[N],q[N],Fc[N];
ll C(ll x,ll y){
	if(y<0||x<y)return 0;
	return fac[x]*inv[y]%mod*inv[x-y]%mod;
}
ll A(ll x,ll y){
	if(y<0||x<y)return 0;
	return fac[x]*inv[x-y]%mod;
}
int main(){
	scanf("%lld",&T);
	ll Max=N-3;
	fac[0]=1;
	rep(i,Max)fac[i]=fac[i-1]*i%mod;
	inv[Max]=Pow(fac[Max],mod-2);
	for(ll i=Max;i;i--)inv[i-1]=inv[i]*i%mod;
	while(T--){
		scanf("%lld%lld%lld%lld",&n,&m,&l,&k);
		ll o=min(n,min(m,l));
		if(k>o){
			puts("0");
			continue;
		}
		Fc[0]=p[1]=1;
		ll nml=n*m%mod*l%mod;
		ll inml=Pow(nml,mod-2);
		rep(i,o-1){
			q[i]=(n-i)*(m-i)%mod*(l-i)%mod*inml%mod;
			q[i]=(1-q[i]+mod)%mod;
			q[i]=q[i]*nml%mod; 
			Fc[i]=Fc[i-1]*q[i]%mod;
		}
		p[o]=Pow(Fc[o-1],mod-2);
		for(ll i=o-1;i;i--)p[i]=p[i+1]*q[i]%mod;
		ll ans=0;
		REP(i,k,o){
			p[i]=p[i]*inml%mod*A(n,i)%mod*A(m,i)%mod*A(l,i)%mod;
			if((i-k)&1)ans=(ans-p[i]*C(i-1,k-1))%mod;
			else ans=(ans+p[i]*C(i-1,k-1))%mod;
		}
		ans=(ans%mod+mod)%mod;
		printf("%lld\n",ans);
	}
	return 0;
}