【CTS2019】随机立方体（容斥）

【CTS2019】随机立方体（容斥）

题面

LOJ
洛谷

题解

做这道题目的时候不难想到容斥的方面。
那么我们考虑怎么计算至少有(k)个极大值的方案数。
我们首先可以把(k)个极大值的位置给确定出来，方案数是(displaystyle {nchoose k}{mchoose k}{lchoose k}(k!)^3)，乘上(k!)是为了确定之间的顺序关系，即我们先确定(xyz)三维，然后把这三维要一一对应到点才行。假设这个值是(w[k])。
剩下要填的是两个部分，一个是剩下的((n-k)(m-k)(l-k))个没有什么影响的位置，以及和极大值有至少一维坐标相交的点。
所以方案数大概可以写成({nmlchoose nml-(n-k)(m-k)(l-k)}w[k]((n-k)(m-k)(l-k))!*h[k])的样子。
其中(h[k])是分配极大值所在的(k*k*k)的这个小立方体上的数字的方案数。
为了方便(V=nml,v[k]=V-(n-k)(m-k)(l-k))。
所以我们只需要考虑怎么分配这个挖掉(k)层的合法方案数。
首先每次确定掉一个极大值之后，我们就可以把它所在的这三个面直接丢掉，变成小一圈的一个立方体，而对于极大值而言，因为它要比所有同层的数都要大，所以我们从大往小，从内往外考虑填数。
首先最值的位置一定是三个面的交点，并且最值一定是当前所有剩余可填的数中的最大值。所以只需要确定剩下的数的数就行了，这个时候还有(v[i]-1)个数可以选，要选(v[i]-v[i-1]-1)个数，所以填进去的方案数就是(frac{(v[i]-1)!}{v[i-1]!})。
所以(h[k]=h[k-1]*frac{(v[k]-1)!}{v[k-1]!})。
所以(displaystyle h[k]=prod_{i=1}^{k} frac{(v[i]-1)!}{v[i-1]!})
然后把答案式掏出来，是：

[egin{aligned} & displaystyle {Vchoose v[k]}w[k](V-v[k])!h[k]\ &=frac{V!}{v[k]!}w[k]h[k]\ &=V!frac{1}{v[k]!}w[k]h[k]\ &=V!frac{1}{v[k]!}w[k]prod_{i=1}^k(v[i]-1)!prod_{i=0}^{k-1}frac{1}{v[i]!}\ &=V!w[k]prod_{i=1}^k frac{1}{v[i]} end{aligned}]

然后题目要求的是概率，所以和(V!)就没有关系了。
那么就只有后半部分。
这样子就很容易计算了。

#include<iostream>
#include<cstdio>
using namespace std;
#define MAX 5001000
#define MOD 998244353
inline int read()
{
	int x=0;bool t=false;char ch=getchar();
	while((ch<'0'||ch>'9')&&ch!='-')ch=getchar();
	if(ch=='-')t=true,ch=getchar();
	while(ch<='9'&&ch>='0')x=x*10+ch-48,ch=getchar();
	return t?-x:x;
}
int jc[MAX],jv[MAX],inv[MAX];
int n,m,l,V,M,k,ans;
int v[MAX],w[MAX],s[MAX],invs[MAX];
int fpow(int a,int b){int s=1;while(b){if(b&1)s=1ll*s*a%MOD;a=1ll*a*a%MOD;b>>=1;}return s;}
int C(int n,int m){return 1ll*jc[n]*jv[m]%MOD*jv[n-m]%MOD;}
int main()
{
	jc[0]=jv[0]=inv[0]=inv[1]=1;
	for(int i=2;i<MAX;++i)inv[i]=1ll*inv[MOD%i]*(MOD-MOD/i)%MOD;
	for(int i=1;i<MAX;++i)jc[i]=1ll*jc[i-1]*i%MOD;
	for(int i=1;i<MAX;++i)jv[i]=1ll*jv[i-1]*inv[i]%MOD;
	int T=read();
	while(T--)
	{
		n=read();m=read();l=read();k=read();V=1ll*n*m%MOD*l%MOD;M=min(min(n,m),l);ans=0;
		for(int i=1;i<=M;++i)v[i]=(V-1ll*(n-i)*(m-i)%MOD*(l-i)%MOD+MOD)%MOD;
		for(int i=1;i<=M;++i)w[i]=1ll*C(n,i)*C(m,i)%MOD*C(l,i)%MOD*jc[i]%MOD*jc[i]%MOD*jc[i]%MOD;
		s[0]=1;for(int i=1;i<=M;++i)s[i]=1ll*s[i-1]*v[i]%MOD;
		invs[M]=fpow(s[M],MOD-2);for(int i=M-1;i;--i)invs[i]=1ll*invs[i+1]*v[i+1]%MOD;
		for(int i=k,d=1;i<=M;++i,d=MOD-d)ans=(ans+1ll*d*C(i,k)%MOD*w[i]%MOD*invs[i])%MOD;
		printf("%d
",ans);
	}
	return 0;
}