Description

T次询问，每次询问给出n,m,k,P，问n个点的所有有标号生成森林中，连通块数为m的方案中，从每棵树中选择一个度数<=k的点的方案数对P取模的方案数
n<=10^9,m,k<=100

Solution

我们先选关键点，再数生成树
显然任意m个点的方案数是相等的，我们只需要指定m个关键点，答案乘上$\binom{n}{m}$(mn​)
考虑枚举这m个点的度数，设总度数为D，剩下的部分为，n-m个点，D棵有标号有根树的生成森林个数
考虑如何做n个点m棵有标号有根树，因为是有根树，新建一个点n+1向所有根连边，这个等价于n+1个点，点n+1的度数为m的树的个数
用prufer序计数，相当于选择m-1个位置填n+1，剩下的位置填1~n，方案数为$\binom{n-1}{m-1}n^{n-m}$(m−1n−1​)nn−m
接下来问题相当于，将n个物品分为m组，每组数量<=k
设F[n][m]表示答案，考虑第n个物品在哪一组，如果加入后不合法相当于这一组本来就有k个，枚举是哪k个剩余的相当于是F[n-1-k][m-1]
所以我们有$F[n][m]=F[n-1][m]*m-F[n-1-k][m-1]*\binom{n-1}{k}*m$F[n][m]=F[n−1][m]∗m−F[n−1−k][m−1]∗(kn−1​)∗m
最后将两部分合并，有一个$\binom{n-1}{m-1}$(m−1n−1​)由于P不一定为质数所以直接暴力分解质因数
单次询问复杂度为$O(m^2k+mk^2+mk\log^2 P)$O(m2k+mk2+mklog2P)

Code

#include <cstdio>
#include <cstring>
#include <algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define fd(i,a,b) for(int i=a;i>=b;i--)
using namespace std;

typedef long long ll;

const int N=1e4+5;

int ty,n,m,k,Mo,lim,cnt[N],p[N];
ll f[105][N],C[N][105];

ll pwr(ll x,int y) {
	ll z=1;
	for(;y;y>>=1,x=x*x%Mo)
		if (y&1) z=z*x%Mo;
	return z;
}

ll calc() {
	ll ret=1;
	fo(i,1,p[0]) ret=ret*pwr(p[i],cnt[i])%Mo;
	return ret;
}

ll ins(int x,int f) {
	if (!x) return 1;
	fo(i,1,p[0]) while (!(x%p[i])) x/=p[i],cnt[i]+=f;
	return x;
}

void exgcd(ll a,ll b,ll &x,ll &y) {
	if (!b) {x=1;y=0;return;}
	ll xx,yy;
	exgcd(b,a%b,xx,yy);
	x=yy;y=xx-a/b*yy;
}

ll inv(ll x) {
	if (x==1) return 1;
	ll a,b;
	exgcd(x,Mo,a,b);
	return (a%Mo+Mo)%Mo;
}

ll Comb(int n,int m) {
	if (m>n-m) m=n-m;
	fo(i,1,p[0]) cnt[i]=0;ll ret=1;
	fo(i,n-m+1,n) ret=ret*ins(i,1)%Mo;
	fo(i,1,m) ret=ret*inv(ins(i,-1))%Mo;
	ret=ret*calc()%Mo;
	return ret;
}

void inc(ll &x,ll y) {x=x+y>=Mo?x+y-Mo:x+y;}
void dec(ll &x,ll y) {x=x-y<0?x-y+Mo:x-y;}

int main() {
	freopen("islands.in","r",stdin);
	freopen("islands.out","w",stdout);
	for(scanf("%d",&ty);ty;ty--) {
		scanf("%d%d%d%d",&n,&m,&k,&Mo);
		if (Mo==1) {puts("0");continue;}
		if (n==m) {puts("1");continue;}
		p[0]=0;int x=Mo;
		for(int i=2;i*i<=x;i++) {
			if (!(x%i)) p[++p[0]]=i;
			while (!(x%i)) x/=i;
		}
		if (x>1) p[++p[0]]=x;
		fo(i,0,m*k) {
			C[i][0]=1;
			fo(j,1,min(i,k)) C[i][j]=(C[i-1][j]+C[i-1][j-1])%Mo;
		}
		fo(i,0,m) fo(j,0,i*k) f[i][j]=0;f[0][0]=1;
		fo(i,1,m) {
			f[i][0]=1;
			fo(j,1,i*k) {
				inc(f[i][j],f[i][j-1]*i%Mo);
				if (j>=k+1) dec(f[i][j],f[i-1][j-k-1]*C[j-1][k]%Mo*i%Mo);
			}
		}
		fo(i,1,p[0]) cnt[i]=0;
		ll ans=0,C=1;
		fo(i,1,m*k) {
			if (i>n-m) break;
			(ans+=f[m][i]%Mo*pwr(n-m,n-m-i)%Mo*calc()%Mo*C%Mo)%=Mo;
			C=C*ins(n-m-i,1)%Mo;
			C=C*inv(ins(i,-1))%Mo;
		}
		printf("%lld\n",ans*Comb(n,m)%Mo);
	}
	return 0;
}