二项式反演好题。
首先看到“恰好 \(k\) 个极大值点”，我们可以套路地想到二项式反演，具体来说我们记 \(f_i\) 为钦定 \(i\) 个点为极大值点的方案数，那么
\[ans=\dfrac{1}{(nml)!}\sum\limits_{i=k}^{\min(n,m,l)}f_i(-1)^{i-k}\dbinom{i}{k} \]考虑怎么求 \(f_i\)，首先我们肯定要选出 \(i\) 个极大的位置。我们假设 \(g_i\) 为选出 \(i\) 个极大的位置的方案数，那么显然 \(g_i\) 就是把每一个位置选择的方案数都乘起来。但这样会算重，具体来说，对于每一个合法的选出 \(i\) 个极大的位置的方案数，我们这样算相当于给全部 \(i\) 个极大点都上了标号，因此一种合法的方案会被重复计算 \(i!\) 次，总方案数还需除以 \(i!\)，即
\[g_i=\dfrac{1}{i!}\prod\limits_{j=0}^{i-1}(n-j)(m-j)(l-j) \]接下来考虑怎么给与这 \(i\) 个极大值点在同一行/列/高度的位置填上数。首先我们设 \(c_i\) 表示与 \(i\) 个极大点在同一行/列/高度的格子数量。显然 \(c_i\) 可以通过总格子数量减去与 \(i\) 个极大值点都不在同一行/列/高度的格子数量，即 \(c_i=nml-(n-i)(m-l)(l-i)\)。再设 \(h_i\) 表示为给与这 \(i\) 个极大值点的任意一个极大值点在同一行/列/高度的位置填上 \(1\sim 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)!\)，因此我们可以得到：
\[h_i=h_{i-1}·\dbinom{c_i-1}{c_i-c_{i-1}-1}·(c_i-c_{i-1}-1)! \]即
\[h_i=h_{i-1}·\dfrac{(c_i-1)!}{c_{i-1}!} \]递推一下有：
\[h_i=\prod\limits_{j=1}^i\dfrac{(c_j-1)!}{c_{j-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\)，因此
\[f_i=\dbinom{nml}{c_i}i!g_ih_i(nml-c_i)! \]展开来可以得到：
\[f_i=\dfrac{(nml)!}{c_i!(nml-c_i)!}·\dfrac{1}{i!}\prod\limits_{j=0}^{i-1}(n-j)(m-j)(l-j)·\prod\limits_{j=1}^i\dfrac{(c_j-1)!}{c_{j-1}!}·i!·(nml-c_i)! \]发现有一堆东西可以怼掉，\(i!\) 和 \(\dfrac{1}{i!}\) 怼掉了，\((nml-c_i)!\) 和 \(\dfrac{1}{(nml-c_i)!}\) 怼掉了，我们还可以发现，前面的 \(\dfrac{1}{c_i!}\) 和后面的 \(\prod\limits_{j=1}^i\dfrac{(c_j-1)!}{c_{j-1}!}\) 拼起来变成 \(\dfrac{\prod\limits_{j=1}^i(c_j-1)!}{\prod\limits_{j=1}^ic_j!}\)，显然这东西等于 \(\prod\limits_{j=1}^i\dfrac{1}{c_j}\)，于是
\[f_i=(nml)!·\prod\limits_{j=0}^{i-1}(n-j)(m-j)(l-j)·\prod\limits_{j=1}^i\dfrac{1}{c_j} \]带到最一开始的答案的式子中
\[ans=\dfrac{1}{(nml)!}\sum\limits_{i=k}^{\min(n,m,l)}(nml)!·(\prod\limits_{j=0}^{i-1}(n-j)(m-j)(l-j)·\prod\limits_{j=1}^i\dfrac{1}{c_j})(-1)^{i-k}\dbinom{i}{k} \]\((nml)!\) 与 \(\dfrac{1}{(nml)!}\) 又怼掉了，剩下的式子中不含我们组合数学中不能直接算的东西（上下底数都很大的组合数、超过 \(10^7\) 的阶乘等），因此预处理出 \(\prod\limits_{j=0}^{i-1}(n-j)(m-j)(l-j)\)，以及 \(\prod\limits_{j=1}^i\dfrac{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\n",res);
}
int main(){
init_fac(MAXN);
int qu;scanf("%d",&qu);while(qu--) solve();
return 0;
}