题面
思路
首先，这道题是可以暴力min-max反演+NTT做出来的......但是这个不美观，我来讲一个做起来舒服一点的做法
一个非常basic的idea：我们发现在一只鸽子吃饱以后再喂给它的玉米都是“无效”的，并且我们如此认为，那么有效的玉米数量是确定的：$nk$
吃饱序列和投喂序列
那么，我们考虑一个序列$r_i$，表示第$i$次喂完玉米之前，有多少只鸽子是吃饱的，我们称之为吃饱序列
注意到本题中每只鸽子互不相同，因此我们再确定一个“有效喂鸽子操作”的序列，我们称之为投喂序列
特别注意：吃饱序列的构造虽然部分依赖于投喂序列，但是投喂序列的出现概率是完全依赖于吃饱序列的
显然，对于一组操作序列和吃饱序列，我们可以得到这组序列出现的总概率：$Prob=\prod_{i=1}^{nk}P_{r_i}$
其中$P_i$表示吃饱了$i$个的情况下，下一个投喂选到我们的目标鸽子的概率
那么我们现在实际上把投喂序列变成了可以由吃饱序列求出来，而吃饱序列的下一项又反过来由投喂序列确定，这么一个情况
我们如果要考虑总贡献，我们发现还需要考虑最终成功完成一次有效投喂（注意因为前面算的是概率，这里只要随便投喂一个没吃饱的鸽子就可以了）的期望时间
这个时间$T_i=E_{r_i}$，其中$E_i=\frac{n}{n-i}$，这里表示在$i$个鸽子吃饱的前提下的有效投喂期望时间
那么，我们可以得到我们在确定了一个吃饱序列和对应的投喂序列时最终答案的表达式：$Ans=\prod_{i=1}{nk} E_{r_i}$
上式的两个部分分别代表每一个投喂序列出现的概率，以及这个吃饱序列的期望完成时间
转化为DP
这个东西不好处理，因为我们没有办法直接知道每次成功投喂以后会不会使吃饱序列的下一项+1（也就是有一只鸽子吃饱了）
注意到贡献都只和$r_i$有关系，而和目前没吃饱的鸽子吃掉的玉米的分配没有关系！
所以我们大可以随意分配这些没吃饱的鸽子吃掉的玉米，下文中简称为白玉米
那么我们可以基于上面的表达式得到一个$dp$的做法：
设$f[i][j]$表示投喂了$i$次，有$j$个鸽子吃饱了的总贡献，$g[i][j]$则表示上述情况出现的概率（也就是只考虑表达式中含$P_{r_i}$的部分）
那么我们把表达式转化一下：
$f[i][j]=\sum_{\lbrace r\rbrace} \prod_{x=1}{i} E_{r_y}$
$f[i][j]=\sum_{\lbrace r\rbrace} (P_j\prod_{x=1}{i-1} E_{r_y})$
$f[i][j]=P_j\ast(\sum_{\lbrace r\rbrace} \prod_{x=1}{i-1} E_{r_y})+P_j\ast E_j\ast(\sum_{\lbrace r\rbrace} \prod_{x=1}^{i-1}P_{r_x})$
$f[i][j]=P_jf[i-1][j]+P_jE_jg[i-1][j]$
这样我们就完成了没有新鸽子吃饱的情况下的$f[i][j]$的转移
那么对于$g[i][j]$的转移，很显然是$g[i][j]=P_jg[i-1][j]$，不再赘述
对于新加入的玉米使得一只鸽子吃饱的情况，我们需要对目前存在的白玉米进行染色，此时染色的方案数显然为$\binom{i-jk}{k-1}$
所以对于从$f[i][j]$到$f[i+1][j+1]$的转移，只需要在上面的转移的基础上乘上上述组合系数即可
若仍有疑问，可以参见代码的实现
Code
#include<iostream>
#include<cstdio>
#include<cstring>
#include<algorithm>
#include<cassert>
#define MOD 998244353
#define ll long long
using namespace std;
inline int read(){
int re=0,flag=1;char ch=getchar();
while(!isdigit(ch)){
if(ch=='-') flag=-1;
ch=getchar();
}
while(isdigit(ch)) re=(re<<1)+(re<<3)+ch-'0',ch=getchar();
return re*flag;
}
inline int add(int a,int b){
a+=b;
if(a>=MOD) a-=MOD;
return a;
}
inline void addd(int &a,int b){
a+=b;
if(a>=MOD) a-=MOD;
}
inline int qpow(int a,int b){
int re=1;
while(b){
if(b&1) re=1ll*re*a%MOD;
a=1ll*a*a%MOD;b>>=1;
}
return re;
}
int n,m,fac[1000010],finv[1000010],inv[1000010];
inline void init(){
int i,len=1000000;
fac[0]=fac[1]=finv[0]=finv[1]=inv[1]=1;
for(i=2;i<=len;i++) inv[i]=1ll*(MOD-MOD/i)*inv[MOD%i]%MOD;
for(i=2;i<=len;i++) fac[i]=1ll*fac[i-1]*i%MOD;
finv[len]=qpow(fac[len],MOD-2);
for(i=len;i>2;i--) finv[i-1]=1ll*finv[i]*i%MOD;
}
inline int C(int x,int y){
if(x<0||y<0||x<y) return 0;
return 1ll*fac[x]*finv[y]%MOD*finv[x-y]%MOD;
}
int f[100010][110],g[100010][110],p[100010],e[100010];
int main(){
n=read();m=read();int i,j,tf,tg,tt;
init();
for(i=0;i<=n;i++) p[i]=inv[n-i],e[i]=1ll*n*inv[n-i]%MOD;
f[0][0]=0;g[0][0]=1;
for(i=0;i<n*m;i++){
for(j=0;j*m<=i;j++){
tg=1ll*g[i][j]*p[j]%MOD;
tf=add(1ll*f[i][j]*p[j]%MOD,1ll*p[j]*e[j]%MOD*g[i][j]%MOD);
tt=C(i-j*m,m-1);
addd(f[i+1][j],tf);
addd(g[i+1][j],tg);
addd(f[i+1][j+1],1ll*tf*tt%MOD);
addd(g[i+1][j+1],1ll*tg*tt%MOD);
}
}
// for(i=0;i<=n*m;i++) for(j=0;j*m<=i;j++) cout<<i<<' '<<j<<' '<<f[i][j]<<' '<<g[i][j]<<'\n';
cout<<(1ll*fac[n]*f[n*m][n]%MOD)<<'\n';
}