Description：
(1<=n<=1e9,1<=m,k<=100)
模数不是质数。
题解：
先选m个点，最后答案乘上(C_{n}^m)。
不妨枚举m个点的度数和D，那么我们需要解决两个问题：
- 一共m个有标号盒子，D个有标号小球放到盒子里，且每个盒子的球数不超过k的方案数。
- n-m个有标号点的D棵有根树的森林划分
Task1：
事实上这个东西可以直接NTT卷起来，效率应该是最高的，但是因为模数不是质数，所以不行。
设(f[i][j])表示i个盒子，j个小球的方案数。
不难得到一个容斥的转移：
(f[i][j]=f[i][j-1]*i-f[i-1[j-(k+1)]*i*C_{j-1}^k)
组合数直接杨辉三角预处理。
复杂度：(O(m^2k))
Task2：
利用扩展Cayley公式：
n个点，m棵树，且1-m的点在不同的树里的方案数：
拓展prufer序列的定义，现在是取出森林中最大的叶子，输出与它相邻的点，删掉它，直到剩下1..m
发现序列长度是n-m，且前n-1-m个位置可以填1..n，最后一个只能填1..m，所以：
(F(n,m)=m*n^{n-1-m})
那直接乘上一个(C_n^m)来把1..m换成其它根即是我们要求的。
或者说直接推导：
新建一个虚点n+1，让所有的根连向n+1，那么就可以做树上purfer（取编号最小的叶子），直到剩下n+1一个点。
由于n+1的度数是m，所以在prufer序列中出现m-1次。
因此(=C_{n-1}^{m-1}*n^{n-m})
和上面是等价的。
Task3：
(Ans=C(n,m)*sum_{i=0}^{mk}f[m][i]*C_{n-m-1}^{i-1}*(n-m)^{n-m-i})
模数不是质数，所以组合数需要分解质因数来算。
Code：
#include<bits/stdc++.h>
#define fo(i, x, y) for(int i = x, B = y; i <= B; i ++)
#define ff(i, x, y) for(int i = x, B = y; i < B; i ++)
#define fd(i, x, y) for(int i = x, B = y; i >= B; i --)
#define ll long long
#define pp printf
#define hh pp("
")
using namespace std;
int T, n, m, k, mo;
ll ksm(ll x, ll y) {
ll s = 1;
for(; y; y /= 2, x = x * x % mo)
if(y & 1) s = s * x % mo;
return s;
}
int u[105], v[105], u0;
int pmo;
void fen(int x) {
pmo = x;
u0 = 0;
for(int i = 2; i * i <= x; i ++) if(x % i == 0) {
u[++ u0] = i; v[u0] = 0;
for(; x % i == 0; x /= i) v[u0] ++;
}
if(x > 1) u[++ u0] = x, v[u0] = 0;
fo(i, 1, u0) pmo = pmo / u[i] * (u[i] - 1);
}
ll c[10005][105];
ll inv(int x) { return ksm(x, pmo - 1);}
struct nod {
int v[11];
};
nod operator * (nod a, nod b) {
a.v[0] = (ll) a.v[0] * b.v[0] % mo;
fo(i, 1, u0) a.v[i] += b.v[i];
return a;
}
nod operator / (nod a, nod b) {
a.v[0] = (ll) a.v[0] * inv(b.v[0]) % mo;
fo(i, 1, u0) a.v[i] -= b.v[i];
return a;
}
nod p[10005], q[10005];
void gg(int x, nod &p) {
if(!x) {
fo(i, 1, u0) p.v[i] = 0;
p.v[0] = 1;
return;
}
fo(i, 1, u0) {
p.v[i] = 0;
for(; x % u[i] == 0; x /= u[i]) p.v[i] ++;
}
p.v[0] = x;
}
void build(int n) {
gg(0, p[0]);
gg(0, q[0]);
fo(i, 1, min(10000, n)) {
gg(i, p[i]);
p[i] = p[i - 1] * p[i];
gg(n - i + 1, q[i]);
q[i] = q[i - 1] * q[i];
}
}
ll C(int x) {
nod w = q[x] / p[x];
ll s = w.v[0];
fo(i, 1, u0) s = s * ksm(u[i], w.v[i]) % mo;
return s;
}
ll f[105][10005];
int main() {
freopen("islands.in", "r", stdin);
freopen("islands.out", "w", stdout);
scanf("%d", &T);
fo(ii, 1, T) {
scanf("%d %d %d %d", &n, &m, &k, &mo);
if(n - m == 0) {
pp("%d
", 1 % mo);
continue;
}
fen(mo);
fo(i, 0, 10000) {
c[i][0] = 1;
fo(j, 1, min(i, 100)) c[i][j] = (c[i - 1][j - 1] + c[i - 1][j]) % mo;
}
build(n - m);
memset(f, 0, sizeof f);
fo(i, 1, m) {
f[i][0] = 1;
fo(j, 1, i * k) {
f[i][j] = f[i][j - 1] * i;
if(j >= k + 1) f[i][j] -= f[i - 1][j - (k + 1)] * c[j - 1][k] % mo * i;
f[i][j] = (f[i][j] % mo + mo) % mo;
}
}
ll ans = 0;
fo(i, 0, m * k) {
if(n - m - 1 - i >= 0) ans += f[m][i] * C(i) % mo * i % mo * ksm(n - m, n - m - 1 - i) % mo;
if(n - m - 1 - i == -1) ans += f[m][i] * C(i) % mo;
}
build(n);
ans = ans % mo * C(m) % mo;
pp("%lld
", ans);
}
}