[LOJ#3119][Luogu5400][CTS2019]随机立方体(容斥+DP)

https://www.cnblogs.com/cjyyb/p/10900993.html

#include<cstdio>
#include<algorithm>
#define rep(i,l,r) for (int i=(l); i<=(r); i++)
using namespace std;

const int N=5000010,mod=998244353;
int n,m,l,d,V,M,k,T,ans,v[N],w[N],s[N],is[N],fac[N],inv[N];

int C(int n,int m){ return n<m ? 0 : 1ll*fac[n]*inv[m]%mod*inv[n-m]%mod; }

int ksm(int a,int b){
    int res=1;
    for (; b; a=1ll*a*a%mod,b>>=1)
        if (b & 1) res=1ll*res*a%mod;
    return res;
}

void init(int n){
    fac[0]=1; rep(i,1,n) fac[i]=1ll*fac[i-1]*i%mod;
    inv[n]=ksm(fac[n],mod-2); for (int i=n; i; i--) inv[i-1]=1ll*inv[i]*i%mod;
}

int main(){
    freopen("cube.in","r",stdin);
    freopen("cube.out","w",stdout);
    init(N-10);
    for (scanf("%d",&T); T--; ){
        scanf("%d%d%d%d",&n,&m,&l,&k); V=1ll*n*m%mod*l%mod; M=min(min(n,m),l); ans=0; d=1;
        rep(i,1,M) v[i]=(V-1ll*(n-i)*(m-i)%mod*(l-i)%mod+mod)%mod;
        rep(i,1,M) w[i]=1ll*C(n,i)*C(m,i)%mod*C(l,i)%mod*fac[i]%mod*fac[i]%mod*fac[i]%mod;
        s[0]=1; rep(i,1,M) s[i]=1ll*s[i-1]*v[i]%mod;
        is[M]=ksm(s[M],mod-2); for (int i=M; i; i--) is[i-1]=1ll*is[i]*v[i]%mod;
        rep(i,k,M) ans=(ans+1ll*d*C(i,k)%mod*w[i]%mod*is[i])%mod,d=mod-d;
        printf("%d
",ans);
    }
    return 0;
}