Link
搞啥不好搞非素数模数。
先枚举港口并将其视为根，那么要算的就是根的儿子数均\(\le k\)的方案数。
考虑枚举根的总儿子数\(x\)，选出这些儿子有\({n-m\choose 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}-{i\choose k}f_{i-k,j-1})\)。
组合数可以用\({n\choose 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\n",mul(ans,C(n,m)));
}
int main(){for(i64 T=read();T;--T) solve();}