https://www.cnblogs.com/cjyyb/p/10900993.html
1 #include<cstdio> 2 #include<algorithm> 3 #define rep(i,l,r) for (int i=(l); i<=(r); i++) 4 using namespace std; 5 6 const int N=5000010,mod=998244353; 7 int n,m,l,d,V,M,k,T,ans,v[N],w[N],s[N],is[N],fac[N],inv[N]; 8 9 int C(int n,int m){ return n<m ? 0 : 1ll*fac[n]*inv[m]%mod*inv[n-m]%mod; } 10 11 int ksm(int a,int b){ 12 int res=1; 13 for (; b; a=1ll*a*a%mod,b>>=1) 14 if (b & 1) res=1ll*res*a%mod; 15 return res; 16 } 17 18 void init(int n){ 19 fac[0]=1; rep(i,1,n) fac[i]=1ll*fac[i-1]*i%mod; 20 inv[n]=ksm(fac[n],mod-2); for (int i=n; i; i--) inv[i-1]=1ll*inv[i]*i%mod; 21 } 22 23 int main(){ 24 freopen("cube.in","r",stdin); 25 freopen("cube.out","w",stdout); 26 init(N-10); 27 for (scanf("%d",&T); T--; ){ 28 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; 29 rep(i,1,M) v[i]=(V-1ll*(n-i)*(m-i)%mod*(l-i)%mod+mod)%mod; 30 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; 31 s[0]=1; rep(i,1,M) s[i]=1ll*s[i-1]*v[i]%mod; 32 is[M]=ksm(s[M],mod-2); for (int i=M; i; i--) is[i-1]=1ll*is[i]*v[i]%mod; 33 rep(i,k,M) ans=(ans+1ll*d*C(i,k)%mod*w[i]%mod*is[i])%mod,d=mod-d; 34 printf("%d ",ans); 35 } 36 return 0; 37 }