组合数学初步

组合数学初步

基本计数原理

加法原理:若完成一个事件 An 类方法,第 i 类有 si 种不同的方案,则完成事件 Ai=1nsi 种方案。

乘法原理:若完成一个事件 A 需要 n 步,第 i 步有 ki 种不同的方式,则完成事件 Ai=1nki 种不同的方案。

减法(简单容斥)原理:若完成事件 A 共有 n 次尝试,其中 c 次不能完成,则其中能完成的次数为 nc

排列

n 个人中选 m 个人排成一排,问方案数。

解:

第一个人有 n 种选择,第二个人有 n1 种选择……第 m 个人有 nm+1 种选择。

分步乘法得 i=1mni+1=n!(nm)!,记为 Anm

组合

n 个人中选 m 个人,问方案数。

解:无需排成一排,把重复的去掉即可。

答案为 n!m!(nm)!,记为 Cnm 或者 (nm)

性质:

  • Cnm=Cnnm

  • i=0n=2n

  • k>0,Cnk=Cn1k+Cn1k1

  • Cnm=Cn1m1×nm

  • i=1nCi+mi=Cm+n+1m+11

证明:

i=1nCi+mi=i=1nCi+mm=i=1nCi+mm+Ci+1i+11=Cm+n+1m+11

组合数的一般求法

  • 暴力乘

时间复杂度 O(n)

1.T 组数据,每次求 Cnmmodpp 固定,m,n103

  • 杨辉三角递推

根据归纳恒等式:

k>0,Cnk=Cn1k+Cn1k1

可以递推,预处理 O(n2),单次 O(1)

int C[N][N];
C[0][0]=1;
FOR(i,1,n) {
	C[i][0]=1;
	FOR(j,1,n) C[i][j]=C[i][j-1]+C[i-1][j-1];
}

2.T 组数据,每次求 Cnmmodpp 固定,m,n105

  • 预处理阶乘

在取模时需要求逆元。

预处理 O(n),单次复杂度 O(1)

int C(int n,int m) {
	if(n<m) return 0;
	return fac[n]*inv[m]%Mo*fac[n-m]%Mo;
}
fac[0]=1;
FOR(i,1,n) fac[i]=fac[i-1]*i%Mo;
inv[n]=power(fac[i],Mo-2);
ROF(i,n,1) inv[i-1]=inv[i]*i%Mo;

3.求 Cnmmodpn,m1018,p105,pP

  • 卢卡斯定理

p 为质数,则有:

Cnm=Cn/pm/p×Cnmodpmmodp

证明:

引理:pP,xp+1(x+1)p(modp)

证明:后者用二项式定理展开:

(x+1)p=i=0pCpixi=(i=1p1Cpixi)+1+xp

显然,括号内均为 p 的倍数,可消掉。

证毕。

n=ps+c,m=pt+d

故原命题等价于求证 Cnm=Cst×Ccd

Cnm(1+x)n 中等价于 xm 的系数。

modp 意义下,(1+x)p1+xp(modp)

根据二项式定理,有:

(1+x)n=(1+x)ps+c=(1+x)ps×(1+x)c=(i=0sCsixip)×(i=0cCcixi)

xm=xpt+d=xpt×xd

因为 s<p,所以后面不可能提供 xpt,前面也不可能提供 xd

xm 一定是两式相乘而得。

故原命题成立。

证毕。

const int N=1e5+10;
int t,fac[N];
LL power(LL a,LL b,LL p) {
	LL res=1;
	for(;b;b>>=1) {
		if(b&1) res=(res*a)%p;
		a=a*a%p;
	}
	return res%p;
}
LL C(int n,int m,int p) {
	if(n<m) return 0;
	return fac[n]*power(fac[m],p-2,p)%p*power(fac[n-m],p-2,p)%p;
}
LL lucas(int n,int m,int p) {
	if(n<p&&m<p) return C(n,m,p);
	return C(n%p,m%p,p)*lucas(n/p,m/p,p)%p;
}
int n,m,Mo;
void solve() {
	cin>>n>>m>>Mo;n=n+m;
	FOR(i,1,Mo) fac[i]=fac[i-1]*i%Mo;
	cout<<lucas(n,m,Mo)<<"\n";
}
main() {
	cin>>t;fac[0]=1;
	while(t--) solve();
	return 0;
}

预处理 O(p),单次 O(logn)

  • Cnmn,m103,不取模

高精度除法直接秒了!

发现无论除的是什么 妖魔鬼怪,最后总会是一个非负整数。

所以考虑对阶乘分解质因数,将除法算式转化为质数幂的指数相减的形式。

最后只需要做高精度乘法了。

单次复杂度 O(n)

(本代码以洛谷P1771为例)

const int N=1e5+10;
int n,x,v[N],phi[N];
vector<int>prime;
LL power(LL a,LL b,LL p) {
	LL res=1%p;
	for(;b;b>>=1) {
		if(b&1) res=(res*a)%p;
		a=(a*a)%p;
	}
	return res;
}
void Euler() {
	phi[1]=1;
	FOR(i,2,100000) {
		if(!v[i]) phi[i]=i-1,v[i]=i,prime.pb(i);
		for(int p:prime) {
			if(i*p>100000||v[i]<p) break;
			v[i*p]=p;
			phi[i*p]=phi[i]*(i%p?p-1:p); 
		}
	}
}
int sum[N];
int Get(int x,int p) {
	int cnt=0;
	while(x) {
		cnt+=x/p;
		x/=p;
	}
	return cnt;
}
vector<int>Mul(vector<int>x,int y) {
	LL tol=0;vector<int>res;
	for(auto c:x) {
		tol+=(LL)y*c;
		res.pb(tol%10);
		tol/=10;
	}
	while(tol) res.pb(tol%10),tol/=10;
	return res;
}
main() {
	cin>>n>>x;
	x=power(x,x,1000);
	Euler();
	for(auto p:prime) sum[p]=Get(x-1,p)-Get(n-1,p)-Get(x-n,p);
	vector<int>ans;ans.pb(1);
	for(auto p:prime) while(sum[p]--) ans=Mul(ans,p);
	reverse(ans.begin(),ans.end());
	for(auto x:ans) write(x);putchar('\n');
	return 0;
}
  • Cnmmodp 的值,其中 n,m109p 不一定为质数,但保证 p 的所有质因子幂次均为 1

p 分解质因数,分别求出原式对 p 的各个质因子取模的结果。

最后用 中国剩余定理 解同余方程组即可。

(这里以洛谷 P2840 为例)

如果你不会中国剩余定理

const int Mo=999911659,N=4e4+10;
int p[4]={2,3,4679,35617},M[4],t[4],a[4],inv[N],fac[N];
int n,g;
int power(int x,int y,int p) {
	int res=1;
	for(;y;y>>=1) {
		if(y&1) res=res*x%p;
		x=x*x%p;
	}
	return res;
}
int C(int n,int m,int p) {
	if(n<m) return 0;
	return fac[n]*inv[m]%p*inv[n-m]%p;
}
int Lucas(int n,int m,int p) {
	if(n<p&&m<p) return C(n,m,p);
	return C(n%p,m%p,p)*Lucas(n/p,m/p,p)%p;
}
vector<int>c; 
void Init(int id) {
	int P=p[id];
	fac[0]=1;
	FOR(i,1,P-1) fac[i]=fac[i-1]*i%P;
	inv[P-1]=power(fac[P-1],P-2,P);
	ROF(i,P-1,1) inv[i-1]=inv[i]*i%P;
	for(int x:c) (a[id]+=Lucas(n,x,P))%=P;
} 
int Crt() {
	int m=Mo-1;
	FOR(i,0,3) M[i]=m/p[i],t[i]=power(M[i],p[i]-2,p[i]);
	int ans=0;
	FOR(i,0,3) (ans+=a[i]*M[i]%m*t[i]%m)%=m;
	return ans;
}
int solve() {
	for(int i=1;i*i<=n;++i) {
		if(n%i) continue;
		c.pb(i);
		if(i==n/i) continue;
		c.pb(n/i);
	}
	FOR(i,0,3) Init(i);
	return Crt();
}
main() {
	cin>>n>>g;
	if(g==Mo) {
		puts("0");
		return 0;
	}
	cout<<power(g,solve(),Mo)<<"\n";
	return 0;
}
  • Cnmmodp 的值,其中 p 不一定为质数,n,m105

跟高精度的思想差不多,乘的时候顺带取模就行了,代码不放了。

  • Cnmmodq 的值,其中 n,m1018,q106

洛谷p4720 扩展lucas定理

这个刺激了

类似上述思想,将 q 分解质因数,设 q 的唯一分解为 p1c1p2c2p3c3poco

求出同余方程组

{xCnm(modp1c1)xCnm(modp2c2)xCnm(modp3c3)xCnm(modpoco)

问题转换成如何求出 Cnmmodpc 的值。

这个式子等价于:

n!m!(nm)!modpc

这玩意没办法用逆元,考虑转化。

由于 pc 仅有一个本质不同的质因子 p,考虑在 m!(nm)! 除掉 p 的若干次幂,使得分母与模数互质,这样就可以用扩展欧几里得求逆元了。

但是,剩下被除掉的部分还是不和模数互质,故还需要在分子上除掉若干次幂消掉分子的被除掉的贡献。

这样式子就变成了

n!pxm!py(nm)!pz×pxyzmodpc

问题转换成求 n!px

n!=n×(n1)×(n2)××1=(p×2p×3p)×(1×2×3)=pnp×np!×(1×2×3)

注意到上面的 pnp 最后一定要被除掉,并发现上述 np! 还可能有 p 因子存在,故设递推式:

f(n)=f(np)×(1×2×3)

这样求解目标就是 f(n),递归边界为 f(0)=1

问题转换成求 f(n) 后面带的那一串乘法算式。

注意到后面的所有乘数都不是 p 的倍数,所以 l<pc,l,pc+l,2pc+l(npc1)×pc+l 其实在 modpc 意义下是等价的,都是 l

所以只求小于 pc 部分的乘积,中间部分直接对小于 pc 的部分做快速幂,剩下的 npc×pcn 的部分暴力乘即可。

求一次 f(n) 复杂度是 O(pclogn) 的,这个 log 是以 p 为底的,所以比一般的 log 要小。

int f(int n,int p,int pk) {
	if(!n) return 1;//边界
	int st=1,ed=1;
	FOR(i,1,pk) if(i%p) (st*=i)%=pk;//暴力求开始部分
	st=power(st,n/pk,pk);//快速幂求中间贡献
	FOR(i,pk*(n/pk),n) if(i%p) (ed*=(i%pk))%=pk;//余项
	return f(n/p,p,pk)*st%pk*ed%pk;//递归求解n/p
}

问题来了,怎么求 pxyz?

其实就是求 n! 中有多少 p 这个质因子。

事实上,在 f(n) 的过程中,已经顺便求出来了。

pnp×np!×(1×2×3)

的这个 pnp,将递归过程中的 np 加起来就行了。

int g(int x,int p) {//求 x! 中有多少 p
	int cnt=0;
	while(x) {
		cnt+=x/p;
		x/=p;
	}
	return cnt;
} 

求出 f(n) 本身与 f(m),f(nm)modpc 意义下的逆元,乘起来后对 pc 取模就是 Cnmmodpc

int inv(int a,int p) {
	int x,y;
	exgcd(a,p,x,y);
	return (x+p)%p;
}
int C(int n2,int m2,int p,int pk) {
	int nc=f(n2,p,pk),mc=inv(f(m2,p,pk),pk),nmc=inv(f(n2-m2,p,pk),pk);//n!,1/m!,1/(n-m)!
	int pa=power(p,g(n2,p)-g(m2,p)-g(n2-m2,p),pk);//p^(x-y-z)
	return nc*mc%pk*nmc%pk*pa%pk;//乘积。
}

最后,使用中国剩余定理合并方程组,即可求出答案,总时间复杂度 O(qlogn)

int exLucas() {
	int nw=p;
	for(int i=2;i*i<=p;++i) {
		if(nw%i) continue;
		int cy=0; 
		while(nw%i==0) ++cy,nw/=i;
		vec.pb(mkp(i,cy));
	}
	if(nw>1) vec.pb(mkp(nw,1));
	for(auto o:vec) {
		int x=o.fr,y=o.se,pk=1;
		FOR(i,1,y) pk*=x;
		crt.pb(mkp(C(n,m,x,pk),pk)); 
	}
	int m=1,ans=0;
	for(auto o:crt) m*=o.se;
	for(auto o:crt) {
		int x=o.se,M=m/x;
		int l,r;exgcd(M,x,l,r);
		l=(l+m)%m;
		(ans+=l%m*M%m*o.fr%m)%=m;
	}
	return ans;
}

全部代码如下:

const int N=1e6+10;
int n,m,p,M[N];
void exgcd(int a,int b,int &x,int &y) {
	if(!b) return x=1,y=0,void();
	exgcd(b,a%b,x,y);
	int t=x;x=y;y=t-(a/b)*y;
}
int mul(int a,int b,int p) {
	int res=0;a%=p,b%=p;
	for(;b;b>>=1) {
		if(b&1) res=(res+a)%p;
		a=a*2LL%p;
	}
	return res;
}
int power(int a,int b,int p) {
	int res=1%p;a%=p;
	for(;b;b>>=1) {
		if(b&1) res=mul(res,a,p)%p;
		a=mul(a,a,p)%p;
	}
	return res;
}
int f(int n,int p,int pk) {
	if(!n) return 1;
	int st=1,ed=1;
	FOR(i,1,pk) if(i%p) (st*=i)%=pk;
	st=power(st,n/pk,pk);
	FOR(i,pk*(n/pk),n) if(i%p) (ed*=(i%pk))%=pk;
	return f(n/p,p,pk)*st%pk*ed%pk;
}
int g(int x,int p) {
	int cnt=0;
	while(x) {
		cnt+=x/p;
		x/=p;
	}
	return cnt;
} 
int inv(int a,int p) {
	int x,y;
	exgcd(a,p,x,y);
	return (x+p)%p;
}
int C(int n2,int m2,int p,int pk) {
	int nc=f(n2,p,pk),mc=inv(f(m2,p,pk),pk),nmc=inv(f(n2-m2,p,pk),pk);
	int pa=power(p,g(n2,p)-g(m2,p)-g(n2-m2,p),pk);
	return nc*mc%pk*nmc%pk*pa%pk;
}
vector<pair<int,int> >vec,crt;
int exLucas() {
	int nw=p;
	for(int i=2;i*i<=p;++i) {
		if(nw%i) continue;
		int cy=0; 
		while(nw%i==0) ++cy,nw/=i;
		vec.pb(mkp(i,cy));
	}
	if(nw>1) vec.pb(mkp(nw,1));
	for(auto o:vec) {
		int x=o.fr,y=o.se,pk=1;
		FOR(i,1,y) pk*=x;
		crt.pb(mkp(C(n,m,x,pk),pk)); 
	}
	int m=1,ans=0;
	for(auto o:crt) m*=o.se;
	for(auto o:crt) {
		int x=o.se,M=m/x;
		int l,r;exgcd(M,x,l,r);
		l=(l+m)%m;
		(ans+=l%m*M%m*o.fr%m)%=m;
	}
	return ans;
}
main() {
	cin>>n>>m>>p;
	cout<<exLucas()<<"\n";
	return 0;
}

讲个笑话,扩展 lucas 和 lucas 没有任何关系。

组合数经典问题

例1.[HNOI2008] 越狱

直接求很不好做。

考虑上面的减法原理,用总数减去不会越狱的方案数。

总数显然是 mn

而不会越狱当且仅当相邻两个监狱的信仰不同。

所以第一个可以有 m 种选择,第二个可有 m1 种选择(不能选第一个人选的),第三个有 m1 种选择(不能选第二个人的)……

所以答案为 mnm×(m1)n1。需要注意减法取模的情况。

例2.错排问题

f(n) 表示当排列长度为 n 时该问题答案。

考虑递推,假设 n 放在了第 k 个位置,则有 n1 种选择。

则此时第 k 个信封有两种选择:

  • 放在第 n 个位置,则剩余 n2 个再次错排,有 f(n2) 种方案。

  • 不放在第 n 个位置,则剩余 n1 个错排,有 f(n1) 种方案。

上述先放 n,再放 k,分两步,用乘法原理。第二步分两类,用加法原理。

f(n)=(n1)×(f(n1)+f(n2))

例3.[SDOI2016] 排列计数

先选几个固定后,剩下的错排即可。

例4.[NOIP2016 提高组] 组合数问题

n,m 很小,使得我们可以杨辉三角递推组合数。

注意到 k 固定,所以递推时可以对 k 取模,一个组合数能够整除 k 当且仅当该数字 modk 后为 0

最后做二维前缀和即可做到每次 O(1) 查询。

例5.方程的解

给你一个不定方程 i=1kxi=g(n),求正整数解个数。

首先,g(n) 可以通过快速幂得到。

考虑这样一个问题:

给你 n 个相同的球与 k 个不同的盒子,每个盒子至少放一个,求将这 n 个球全部放入 k 个盒子的方案数。

将这 n 个球排成一排,最左侧和最右侧分别放上一个板子,如下:

|ooooooooooooooo|

这样问题转换成放入 k1 个板子(不算两边的板子),两个板子之间至少有一个球,求方案数。

如下即为一个 k=3 的合法方案:

|ooo|ooooooo|ooooo|

最终每个盒子放入的球就是两个板子之间的球的数量。

注意到球之间共有 n1 个空隙,我们只要在选择 k1 个空隙,对每个空隙分别插入一个板子就一定能够构造出一个合法方案。

问题转换成从 n1 个东西中选出 k1 个东西,求方案数。

答案很显然了,Cn1k1

回到上面的问题,我们抽象一下,将 g(x)1 排成一排,有 k 个数,初始为 0,每个数要选择一些 1 加到自己上面,不能不选,求方案数。

这不就是放球问题吗?

这样问题就完美解决,答案即为 Cg(x)1k1

不过,该题没有模数,需要写高精。

扩展:给定不定方程 i=1kxi=n,求非负整数解的个数。

直接求显然不好做。

考虑做一些数学变换。

i=1kxi=ni=1k(xi+1)=n+ki=1kyi=n+k

上面的 yi 是换元的 xi+1

发现 yi 的取值与 xi 一一对应。

所以我们求出 yi 的解数就可以求出 xi 的解数。

而因为 yi1,所以这题做完了。

Cn+k1k1

例6.[USACO09FEB] Bulls And Cows S

枚举放多少头公牛,假设放 p 头。

则母牛剩下 np 头。

问题变成放 p 个板子隔成,左、右留的位置无限制,但板子之间必须空 k 个以上。

写成不定方程的形式就是

i=1p+1xi=np,x1,xp+10,i[2,p]Z,xik

考虑像换元一样,将左边除左右两边的所有项减去 k,原式可化为:

i=1p+1xi=npk×(p1),i[1,p+1]Z,xi0

这样就是不定方程非负整数解板子了,放置公牛 p 个的答案即:

Cnk×(p1)+1p

最终答案即:

p=0nCnk×(p1)+1p

按原式计算即可。

例7.车的放置

先考虑一个 nm 列的普通矩形放 k 个怎么做。

首先肯定是在 m 行中选 k 列,方案数为 Cmk

其次,选第一个所占得行时有 n 种方案,第二个有 n1 种。

所以方案数有 Cmk×Ank 种。

接着将原图形分成两个矩形,这里设分成 a×b(a+c)×d 的矩形。

先放第一个,枚举放多少个(设为 k 个),方案有 Cak×Abk 种。

再放第二个,设第二个放了 z 个,但是放第二个时由于第一个放 k 个后第二个的行只有 a+ck 个了。

所以方案数是 Ca+ckz×Adz

乘起来就是答案。

例8.[CQOI2014] 数三角形

一个 N×M 的网格上有 (N+1)×(M+1) 个格点,为了叙述方便,以下记 n=N+1,m=M+1

先不考虑不在同一条直线上的限制,就是在所有点中任选 3 个,方案显然有 Cnm3 种。

考虑在此基础上减掉不合法的方案。

显然,不合法方案分成横着、竖着、斜着的。

对于横着,一共 n 行,每一行都有 Cm3 种方案,乘起来即可。

竖着的同理。

斜着的有两种情况:斜率正与斜率负。

显然他们是对称的,即方案数相同。

只需要求出斜率为正的方案,最终乘 2 即可。

注意到我们可以枚举两个不在同一水平线且不在同一竖直线上的点,这样它们中间的点的个数(不包括两端)就是答案。

因为任选中间的一个点与两端的两个点总是可以在同一水平线上。

引理:在格点上,若一个点坐标为 (a,b),另一个点坐标为 (c,d),则它们连线上格点个数(不包括两端)为

gcd(|ac|,|bd|)1

证明:

为叙述方便,以下把绝对值去掉。

对于斜边上任意一点 (x,y),满足 ybxa=dbca,x[a,c],y[b,d]

来一张图。

g=gcd(db,ca)

则:

y=dbca×(xa)+b=dbgcag×(xa)+b=dbg×xacag+b

yN,则 cagxa

则一定存在一个 q,使得 cag×q=xa

x=cag×q+a

经过一堆不等式变换后(或者按原不等式取值范围)得到 cag×q+a[a,c]

最终得到 q[0,g]Z

q 最多有 g+1 种取值,即最多有 g+1 种合法的 x

掐头去尾即可证明原命题。

证毕。

直接枚举斜边的端点,计算即可。

事实上,直接枚举是 O(n2m2logn) 的,需要进一步优化。

实际上,我们固定形状后,格子图中有多少个这类形状的三角形已经确定了。

假设底、高分别为 a,b,则数量显然为 (na+1)×(mb+1)

证明留给读者思考。

复杂度 O(nmlogn)

int ans,n,m,s[1010][1010],a[1010][1010];
int gcd(int x,int y) {
	return y?gcd(y,x%y):x;
}
int calc(int x) {
	return x*(x-1)*(x-2)/6;
} 
main() {
	cin>>n>>m;
	LL c=(n+1)*(m+1); 
	ans=calc(c)-(m+1)*calc(n+1)-(n+1)*calc(m+1);
	FOR(i,1,n+1) FOR(j,1,m+1) ans-=2LL*(gcd(i,j)-1)*(n-i+1)*(m-j+1);
	cout<<ans<<"\n";
	return 0;
}

多重集排列数

多重集是一种广义的集合。

普通的集合是不可重的,即集合中每个元素都不相同。

多重集是一种可重集合,可以表示成这么一种形式:

A={n1×a1,n2×a2,,nk×ak}

上述可重集 A 就是描述了一个集合内有 n1a1n2a2…… nkak

所谓排列数,就是求 A 有多少个本质不同的排列。

n=i=1kni

则上述问题的答案为:

n!i=1kni!

简证:

不考虑数字相同的情况,共有 n! 种。

每一个相同的数字内部总共会有 ni! 种,需要除掉。

证毕。

例9.Counting swap

ipi 之间连有向边,容易发现一定会形成若干个有向环。

若交换完成,则一定是 n 个自环。

而最优解一定是每次交换都能够让至少一个数字复位,在图上的表现就是让一个大环断成两个小环。

fi 表示让一个长为 i 的大环分解成 i 个自环的方案数。

再设 g(x,y) 表示将一个长为 x+y 的环分成长为 x 和长为 y 的环的方案数。

注意到若 x=yg(x,y)=x,否则 g(x,y)=x+y

因为第一种可能会有重复情况,画图即可理解,这里不再赘述。

有转移:

fix=1i2fx×fy×g(x,y)

这样就完了。……吗?

不然,考虑上述乘法原理,我们先把 x 环消掉再做 y 环,显然漏了情况。

其实可以先将 x 环分成两个,再将 y 环分成两个……

交替着做也是一种方案!肿么办?

发现一个长为 x 的环操作次数必然是 x1

则上面的操作总数一共是 x+y2,要从中选 x1 次操作给第一个环,剩下的给第二个环,有多少方案?

这不就组合数吗!

或者也可以按多重集排列数理解,x11y10 有多少排列方案。

最终答案是一样的。

所以真正的方程应为:

fix=1i2fx×fy×g(x,y)×Cx+y2x1

最后 dfs 找环贡献加起来就行。

但是这样时间复杂度是 O(n2) 的,无法通过。

事实上,f 有个通项公式:

fn=nn2

这样用快速幂就可以 O(n) 求出答案了。

等等……为啥是 O(n)

设第 i 个环的长度为 li,一共有 k 个环。

则有 i=1kli=ni=1klogli

所以千万不要预处理 f 数组!

const int N=1e5+10,Mo=1e9+9;
int n,cnt,col[N],ct[N];
LL inv[N],fac[N];
bool vis[N];
vector<int>e[N];
void dfs(int x) {
	col[x]=cnt;
	for(int y:e[x]) {
		if(col[y]) continue;
		dfs(y);
	}
} 
LL power(LL a,LL b) {
	LL res=1;
	for(;b;b>>=1) {
		if(b&1) res=res*a%Mo;
		a=a*a%Mo;
	}
	return res;
}
void solve() {
	cin>>n;cnt=0;
	FOR(i,1,n) {
		int x=read();
		e[x].pb(i);
	}
	FOR(i,1,n) if(!col[i]) ++cnt,dfs(i);
	FOR(i,1,n) ct[col[i]]++;
	LL ans=fac[n-cnt];
	FOR(i,1,cnt) {
		(ans*=inv[ct[i]-1])%=Mo;
		if(ct[i]<2) continue;
		(ans*=power(ct[i],ct[i]-2))%=Mo;
	}
	cout<<ans<<"\n";
	FOR(i,1,n) e[i].clear(),ct[i]=col[i]=0;
}
main() {
	fac[0]=1;
	FOR(i,1,N-10) fac[i]=fac[i-1]*i%Mo;
	inv[N-10]=power(fac[N-10],Mo-2);
	ROF(i,N-10,1) inv[i-1]=inv[i]*i%Mo;
	int t=read();
	while(t--) solve();
	return 0;
}

最后的问题,fi 的通项是怎么来的?

观察可得、oeis大法

证明:

引理:有标号无根树的方案数为 nn2(Cayley 公式)。

证明:根据 Prüfer序列 的性质,任何长为 n2 的值域为 [1,n] 的序列都可以唯一确定成一棵有标号树。

根据乘法原理,方案数显然为 nn2

证毕。

考虑一个有 n 个点的环,显然需要 n1 次交换就可以分成 n 个自环。

假设有一张图 G,开始时没有边。

每次若交换 px,py,则将 x,y 之间连一条无向边。

显然,这个图不会出现环,否则一定是断开的环之间重连,一定不是最优解。

而又因为该图有 n1 条边,所以操作完成后 G 必定是一棵树。

考虑将操作倒过来,不断让两个环之间交换,最后合成一个大环,这样就转化为在树上删边。

若将树的每条边给一个顺序,按照这个顺序进行删边并合并环,则最后得到的环一定是不变的。

根据引理,边没有顺序的有标号树有 nn2 棵,加上边的顺序后就是 nn2×(n1)!

这样就把边有标号的环和边有标号的有标号树建立了一一对应关系。

又因为合并顺序不同的环在不考虑合并顺序的差异后均为原环,本质上是一样的,有重复。

所以要去重,要在原来基础上除掉边重复的方案即 (n1)!

所以方案数就是 nn2

证毕。

posted @   cannotmath  阅读(120)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 在鹅厂做java开发是什么体验
· 百万级群聊的设计实践
· WPF到Web的无缝过渡:英雄联盟客户端的OpenSilver迁移实战
· 永远不要相信用户的输入:从 SQL 注入攻防看输入验证的重要性
· 浏览器原生「磁吸」效果!Anchor Positioning 锚点定位神器解析
点击右上角即可分享
微信分享提示