Luogu5400 CTS2019随机立方体（容斥原理）

　　考虑容斥，计算至少有k个极大数的概率。不妨设这k个数对应的格子依次为(k,k,k)……(1,1,1)。那么某一维坐标<=k的格子会对这些格子是否会成为极大数产生影响。先将这样的所有格子和一个数集对应起来，即将答案乘上一个组合数。然后需要考虑的就是这些格子有多少种合法排列顺序。

　　这个排列需要满足的是(i,i,i)之前不能出现某一维坐标为i的格子。可以看做是填完(i,i,i)后，所有三维坐标中最小值为i的格子就可以填了。这样的格子数量容易计算。于是考虑将格子依次塞进排列，显然第一位只能放(k,k,k)，然后所有三维坐标最小值为k的格子被解锁，用一个组合数将他们放在排列中任意位置，再继续放(k-1,k-1,k-1)，以此类推。

　　这样最后化一化得到一些东西，可以发现要计算的是一个数组前缀积的逆元。可以使用经典trick，求出整个数组积的逆元再倒序还原，即可做到线性。

#include<bits/stdc++.h>

using namespace std;

#define ll long long

#define P 998244353

#define N 5000010

char getc(){char c=getchar();while ((c<'A'||c>'Z')&&(c<'a'||c>'z')&&(c<'0'||c>'9')) c=getchar();return c;}

int gcd(int n,int m){return m==0?n:gcd(m,n%m);}

int read()

{

    int x=0,f=1;char c=getchar();

    while (c<'0'||c>'9') {if (c=='-') f=-1;c=getchar();}

    while (c>='0'&&c<='9') x=(x<<1)+(x<<3)+(c^48),c=getchar();

    return x*f;

}

int T,n,m,l,k,fac[N],inv[N],f[N],g[N];

int ksm(int a,int k)

{

    int s=1;

    for (;k;k>>=1,a=1ll*a*a%P) if (k&1) s=1ll*s*a%P;

    return s;

}

int Inv(int a){return ksm(a,P-2);}

void inc(int &x,int y){x+=y;if (x>=P) x-=P;}

int C(int n,int m){if (m>n) return 0;return 1ll*fac[n]*inv[m]%P*inv[n-m]%P;}

int A(int n,int m){if (m>n) return 0;return 1ll*fac[n]*inv[n-m]%P;}

int min(int x,int y,int z){return min(min(x,y),z);}

int main()

{

#ifndef ONLINE_JUDGE

    freopen("a.in","r",stdin);

    freopen("a.out","w",stdout);

    const char LL[]="%I64d\n";

#else

    const char LL[]="%lld\n";

#endif

    T=read();

    fac[0]=1;for (int i=1;i<=N-10;i++) fac[i]=1ll*fac[i-1]*i%P;

    inv[0]=inv[1]=1;for (int i=2;i<=N-10;i++) inv[i]=P-1ll*(P/i)*inv[P%i]%P;

    for (int i=2;i<=N-10;i++) inv[i]=1ll*inv[i]*inv[i-1]%P;

    while (T--)

    {

        n=read(),m=read(),l=read(),k=read();

        int ans=0;

        for (int i=1;i<=min(n,m,l);i++) f[i]=(1ll*(n-i+1)*(m-i+1)%P*(l-i+1)%P-1ll*(n-i)*(m-i)%P*(l-i)%P+P)%P;

        for (int i=1;i<=min(n,m,l);i++) f[i]=(f[i]+f[i-1])%P;g[min(n,m,l)]=1;

        for (int i=1;i<=min(n,m,l);i++) g[min(n,m,l)]=1ll*g[min(n,m,l)]*f[i]%P;

        g[min(n,m,l)]=Inv(g[min(n,m,l)]);

        for (int i=min(n,m,l)-1;i>=1;i--) g[i]=1ll*g[i+1]*f[i+1]%P;

        for (int i=k;i<=min(n,m,l);i++)

        {

            int waytochoosemax=1ll*A(n,i)*A(m,i)%P*A(l,i)%P;

            if (i-k&1) inc(ans,P-1ll*C(i,k)*waytochoosemax%P*g[i]%P);

            else inc(ans,1ll*C(i,k)*waytochoosemax%P*g[i]%P);

        }

        cout<<ans<<endl;

    }

    return 0;

}