二项式反演好题。
首先看到“恰好 (k) 个极大值点”，我们可以套路地想到二项式反演，具体来说我们记 (f_i) 为钦定 (i) 个点为极大值点的方案数，那么
考虑怎么求 (f_i)，首先我们肯定要选出 (i) 个极大的位置。我们假设 (g_i) 为选出 (i) 个极大的位置的方案数，那么显然 (g_i) 就是把每一个位置选择的方案数都乘起来。但这样会算重，具体来说，对于每一个合法的选出 (i) 个极大的位置的方案数，我们这样算相当于给全部 (i) 个极大点都上了标号，因此一种合法的方案会被重复计算 (i!) 次，总方案数还需除以 (i!)，即
接下来考虑怎么给与这 (i) 个极大值点在同一行/列/高度的位置填上数。首先我们设 (c_i) 表示与 (i) 个极大点在同一行/列/高度的格子数量。显然 (c_i) 可以通过总格子数量减去与 (i) 个极大值点都不在同一行/列/高度的格子数量，即 (c_i=nml-(n-i)(m-l)(l-i))。再设 (h_i) 表示为给与这 (i) 个极大值点的任意一个极大值点在同一行/列/高度的位置填上 (1sim c_i) 的数的方案数。我们考虑填上 (c_i) 的那个极值点，显然如果去掉那个极值点，那问题就转化为 (i-1) 的情况，方案数自然就是 (h_{i-1})。而增加这个极值点则会多出 (c_i-c_{i-1}) 个与极大值点在同一行/列/高度的位置，由于第 (i) 个极值点已经填好了数，因此我们还需选出 (c_i-c_{i-1}-1) 个数，方案数为 (dbinom{c_i-1}{c_i-c_{i-1}-1})，填好这 (c_i-c_{i-1}-1) 个数后还可以将它们随意排列，方案数就是 ((c_i-c_{i-1}-1)!)，因此我们可以得到：
即
递推一下有：
最后考虑怎样求 (f_i)，首先我们钦定 (i) 个极大值的位置，方案数 (g_i)，我们还要选出 (c_i) 个数安排给与 (i) 个极大值点在同一行/列/高度的点们，方案数 (dbinom{nml}{c_i})，由于 (i) 个极值点顺序可以调换，因此还要乘上 (i!)；剩余 (nml-c_i) 个点可以随便填，方案数 ((nml-c_i)!)，最后我们还要将这 (c_i) 个数填到对应的位置上去，方案数 (h_i)，因此
展开来可以得到：
发现有一堆东西可以怼掉，(i!) 和 (dfrac{1}{i!}) 怼掉了，((nml-c_i)!) 和 (dfrac{1}{(nml-c_i)!}) 怼掉了，我们还可以发现，前面的 (dfrac{1}{c_i!}) 和后面的 (prodlimits_{j=1}^idfrac{(c_j-1)!}{c_{j-1}!}) 拼起来变成 (dfrac{prodlimits_{j=1}^i(c_j-1)!}{prodlimits_{j=1}^ic_j!})，显然这东西等于 (prodlimits_{j=1}^idfrac{1}{c_j})，于是
带到最一开始的答案的式子中
((nml)!) 与 (dfrac{1}{(nml)!}) 又怼掉了，剩下的式子中不含我们组合数学中不能直接算的东西（上下底数都很大的组合数、超过 (10^7) 的阶乘等），因此预处理出 (prodlimits_{j=0}^{i-1}(n-j)(m-j)(l-j))，以及 (prodlimits_{j=1}^idfrac{1}{c_j})，然后直接算上式的值是没问题的。有一个注意点，就是直接对每个 (c_i) 算一波逆元复杂度会多个 (log)，然后你就会获得 80pts 的好成绩。不过考虑借鉴求阶乘及其逆元的套路，我们先递推求出 (c_j) 的前缀积，然后对最后一项求一波逆元，然后再从后往前递推即可求出 (dfrac{1}{c_j}) 的前缀积，这大概也算是组合数学中一个优化的小套路了吧（
复杂度线性。
const int MAXN=5e6;
const int MOD=998244353;
int qpow(int x,int e){
int ret=1;
for(;e;e>>=1,x=1ll*x*x%MOD) if(e&1) ret=1ll*ret*x%MOD;
return ret;
}
int fac[MAXN+5],ifac[MAXN+5];
void init_fac(int n){
for(int i=(fac[0]=ifac[0]=ifac[1]=1)+1;i<=n;i++) ifac[i]=1ll*ifac[MOD%i]*(MOD-MOD/i)%MOD;
for(int i=1;i<=n;i++) fac[i]=1ll*fac[i-1]*i%MOD,ifac[i]=1ll*ifac[i-1]*ifac[i]%MOD;
}
int binom(int x,int y){return 1ll*fac[x]*ifac[y]%MOD*ifac[x-y]%MOD;}
int n,m,l,k,g[MAXN+5],pre[MAXN+5],b[MAXN+5],inv_pre[MAXN+5],f[MAXN+5];
int calc(int x){return 1ll*(n-x)*(m-x)%MOD*(l-x)%MOD;}
void solve(){
scanf("%d%d%d%d",&n,&m,&l,&k);b[0]=pre[0]=1;
int lim=min(min(n,m),l),sz=1ll*n*m%MOD*l%MOD,res=0;
for(int i=1;i<=lim;i++) g[i]=(sz-calc(i)+MOD)%MOD;
for(int i=1;i<=lim;i++) pre[i]=1ll*pre[i-1]*g[i]%MOD;
inv_pre[lim]=qpow(pre[lim],MOD-2);
for(int i=lim-1;~i;i--) inv_pre[i]=1ll*inv_pre[i+1]*g[i+1]%MOD;
for(int i=1;i<=lim;i++) b[i]=1ll*b[i-1]*calc(i-1)%MOD;
for(int i=1;i<=lim;i++) f[i]=1ll*b[i]*inv_pre[i]%MOD;
for(int i=k;i<=lim;i++){
int mul=1ll*binom(i,k)*f[i]%MOD;
if((i-k)&1) res=(res-mul+MOD)%MOD;
else (res+=mul)%=MOD;
} printf("%d
",res);
}
int main(){
init_fac(MAXN);
int qu;scanf("%d",&qu);while(qu--) solve();
return 0;
}