meaning of the title

n × m × l n\times m\times l n × m × l lattice cubes, filled randomly 1 1 1 to n × m × l n\times m\times l n × m × The number of l.
If a lattice is larger than the number on all the lattices with the same dimension at least, it is maximal.
Q QIAOYOU k k The probability of k maximal lattices, m o d   998244353 mod\ 998244353 mod 998244353， T T Group T Data

Data range

n , m , l ≤ 5 × 1 0 6 , k ≤ 100 , T ≤ 10 , 12 s n,m,l\leq 5\times 10^6,k\leq 100,T\leq 10,12s n,m,l≤5×106,k≤100,T≤10,12s

solution

I have a simple way to think about it.
In fact, a random arrangement can be turned into a random sequence of infinite lattices. Each lattice is random with equal probability, and numbers are filled in order of the first occurrence. This is equivalent.
Random lattice is easier to do. It can be regarded as three sequences in the value range [ 1 , n ] , [ 1 , m ] , [ 1 , l ] [1,n],[1,m],[1,l] [1,n],[1,m],[1,l] random, recorded as A , B , C A,B,C A,B,C
It is easy to find that for a very large lattice, set it to ( a , b , c ) (a,b,c) (a,b,c), set i i i is the minimum position satisfied A i = a , B i = b , C i = c A_i=a,B_i=b,C_i=c Ai​=a,Bi​=b,Ci​=c. i i i is a a a in A A In A, b b b in B B In B, c c c in C C The first position in C.
It's not to be denounced or calculated. Now it's decided by the emperor k k k locations, A , B , C A,B,C A. B and C have new numbers in these positions. Consider all imperial appointments k k The k-location scheme can be used first A n k A m k A l k ( 1 n m l ) k A_n^kA_m^kA_l^k(\frac{1}{nml})^k Ang ， Amk ， Alk (nml1) k to fix the numbers filled in these positions. Then the interval between the two positions is determined. Because it is an infinite sequence, this can be determined at will 0 0 0 to + ∞ +\infty +∞, it is easy to find that this is equivalent to an equal ratio sequence, the second i i i position and i + 1 i+1 For the position between i+1 positions, the probability that the filled number meets the conditions is ( n − k + i ) ( m − k + i ) ( l − k + i ) n m l \frac{(n-k+i)(m-k+i)(l-k+i)}{nml} nml(n − k+i)(m − k+i)(l − k+i) is the common ratio. By multiplying the sum values of these equal ratio sequences, we get the imperial ratio k k Probability of all cases at k locations.
Finally, the inclusion exclusion coefficient. It should be noted that the first position of the maximum lattice must be removed to obtain ( − 1 ) ( i − k ) ( i − 1 k − 1 ) (-1)^{(i-k)}\tbinom{i-1}{k-1} (− 1)(i − k)(k − 1i − 1).
Preprocessing factorial and factorial inverse elements can O ( 1 ) O(1) O(1) obtains the permutation number and combination number. It is required to calculate the proportional sequence O ( n ) O(n) The inverse of the number of O(n) can be achieved by using the method of luoguP5431 O ( n ) O(n) O(n).
Time complexity O ( T n ) O(Tn) O(Tn), basically useless 1.5 s 1.5s After 1.5s...

code

#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;
}