Link
搞啥不好搞非素数模数。
先枚举港口并将其视为根，那么要算的就是根的儿子数均(le k)的方案数。
考虑枚举根的总儿子数(x)，选出这些儿子有({n-mchoose x})种方案，把(x)个儿子挂在(m)个根上的方案数设为(f_{x,m})，然后把根删掉，剩下的就是(n-m)个点构成(x)棵有根树的方案数，这个可以Generalized Caylay算。
(f_{x,m})等价于(x)个球染(m)种色，每种颜色的个数(le k)的方案数，转移考虑补集容斥，(f_{i,j}=j(f_{i-1,j}-{ichoose k}f_{i-k,j-1}))。
组合数可以用({nchoose m}=frac{n^{underline m}}m)算，然后模数不是质数所以代码多了1K。
#include<cstdio>
#include<cstring>
using i64=long long;
const int N=101;
i64 read(){i64 x;scanf("%lld",&x);return x;}
i64 P,tot,f[N][N*N],inv[N*N],fac[N],cnt[N],s1[N],s2[N],c[N*N][N];
i64 inc(i64 a,i64 b){return a+=b,a>=P? a-P:a;}
i64 dec(i64 a,i64 b){return a-=b,a<0? a+P:a;}
i64 mul(i64 a,i64 b){return 1ull*a*b%P;}
i64 pow(i64 a,i64 k){i64 r=1;for(;k;k>>=1,a=mul(a,a))if(k&1)r=mul(a,r);return r;}
void exgcd(i64 a,i64 b,i64&x,i64&y){!b? (x=1,y=0):(exgcd(b,a%b,y,x),y-=a/b*x);}
i64 cal(i64 x,i64 d,i64*s){for(int i=1;i<=tot;++i) while(!(x%fac[i])) s[i]=inc(s[i],d),x/=fac[i];return x;}
i64 product(i64*s){i64 r=1;for(int i=1;i<=tot;++i) r=mul(r,pow(fac[i],s[i]));return r;}
i64 C(int n,int m)
{
i64 a=1,b=1,x,y;memset(s1,0,sizeof s1);
for(int i=n-m+1;i<=n;++i) a=mul(a,cal(i,1,s1));
for(int i=1;i<=m;++i) b=mul(b,cal(i,P-1,s1));
return exgcd(b,P,x,y),mul(mul(a,inc(x,P)),product(s1));
}
void FAC(i64 t)
{
memset(cnt+1,0,sizeof cnt),tot=0;
for(i64 i=2;i*i<=t;++i) if(!(t%i)) for(fac[++tot]=i;!(t%i);t/=i,++cnt[tot]);
if(t^1) fac[++tot]=t,cnt[tot]=1;
}
void solve()
{
int n=read(),m=read(),k=read();i64 ans=0;
FAC(P=read()),memset(f,0,sizeof f),memset(s2,0,sizeof s2);
for(int i=0,j;i<=10000;++i) for(j=c[i][0]=1;j<=i&&j<=100;++j) c[i][j]=inc(c[i-1][j-1],c[i-1][j]);
for(int i=1;i<=m;++i)
{
f[i][0]=1;
for(int j=1;j<=m*k;++j)
{
f[i][j]=mul(i,f[i][j-1]);
if(j>k) f[i][j]=dec(f[i][j],mul(mul(i,c[j-1][k]),f[i-1][j-k-1]));
}
}
for(i64 i=0,a=1,b=1,x,y;i<=m*k;++i)
if(n>i+m) ans=inc(ans,mul(mul(mul(mul(f[m][i],pow(n-m,n-i-m-1)),i),a),product(s2))),a=mul(a,cal(n-m-i,1,s2)),b=cal(i+1,P-1,s2),exgcd(b,P,x,y),b=inc(x,P),a=mul(a,b);
else if(n==i+m) ans=inc(ans,f[m][i]);
printf("%lld
",mul(ans,C(n,m)));
}
int main(){for(i64 T=read();T;--T) solve();}