组合数学初步
组合数学初步
基本计数原理
加法原理:若完成一个事件 有 类方法,第 类有 种不同的方案,则完成事件 有 种方案。
乘法原理:若完成一个事件 需要 步,第 步有 种不同的方式,则完成事件 有 种不同的方案。
减法(简单容斥)原理:若完成事件 共有 次尝试,其中 次不能完成,则其中能完成的次数为 。
排列
从 个人中选 个人排成一排,问方案数。
解:
第一个人有 种选择,第二个人有 种选择……第 个人有 种选择。
分步乘法得 ,记为 。
组合
从 个人中选 个人,问方案数。
解:无需排成一排,把重复的去掉即可。
答案为 ,记为 或者 。
性质:
-
。
-
。
-
。
-
。
-
证明:
组合数的一般求法
- 暴力乘
时间复杂度 。
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. 组数据,每次求 , 固定,。
- 预处理阶乘
在取模时需要求逆元。
预处理 ,单次复杂度 。
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.求 ,。
- 卢卡斯定理
若 为质数,则有:
证明:
引理:
证明:后者用二项式定理展开:
显然,括号内均为 的倍数,可消掉。
证毕。
设
故原命题等价于求证 。
在 中等价于 的系数。
在 意义下,。
根据二项式定理,有:
又 。
因为 ,所以后面不可能提供 ,前面也不可能提供 。
故 一定是两式相乘而得。
故原命题成立。
证毕。
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;
}
预处理 ,单次 。
- 求 ,,不取模
高精度除法直接秒了!
发现无论除的是什么 妖魔鬼怪,最后总会是一个非负整数。
所以考虑对阶乘分解质因数,将除法算式转化为质数幂的指数相减的形式。
最后只需要做高精度乘法了。
单次复杂度 。
(本代码以洛谷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;
}
- 求 的值,其中 , 不一定为质数,但保证 的所有质因子幂次均为 。
将 分解质因数,分别求出原式对 的各个质因子取模的结果。
最后用 中国剩余定理 解同余方程组即可。
(这里以洛谷 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;
}
- 求 的值,其中 不一定为质数,。
跟高精度的思想差不多,乘的时候顺带取模就行了,代码不放了。
- 求 的值,其中
这个刺激了
类似上述思想,将 分解质因数,设 的唯一分解为 。
求出同余方程组
问题转换成如何求出 的值。
这个式子等价于:
这玩意没办法用逆元,考虑转化。
由于 仅有一个本质不同的质因子 ,考虑在 与 除掉 的若干次幂,使得分母与模数互质,这样就可以用扩展欧几里得求逆元了。
但是,剩下被除掉的部分还是不和模数互质,故还需要在分子上除掉若干次幂消掉分子的被除掉的贡献。
这样式子就变成了
问题转换成求 。
注意到上面的 最后一定要被除掉,并发现上述 还可能有 因子存在,故设递推式:
这样求解目标就是 ,递归边界为 。
问题转换成求 后面带的那一串乘法算式。
注意到后面的所有乘数都不是 的倍数,所以 其实在 意义下是等价的,都是 。
所以只求小于 部分的乘积,中间部分直接对小于 的部分做快速幂,剩下的 的部分暴力乘即可。
求一次 复杂度是 的,这个 是以 为底的,所以比一般的 要小。
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
}
问题来了,怎么求 ?
其实就是求 中有多少 这个质因子。
事实上,在 的过程中,已经顺便求出来了。
的这个 ,将递归过程中的 加起来就行了。
int g(int x,int p) {//求 x! 中有多少 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);//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;//乘积。
}
最后,使用中国剩余定理合并方程组,即可求出答案,总时间复杂度 。
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 没有任何关系。
组合数经典问题
直接求很不好做。
考虑上面的减法原理,用总数减去不会越狱的方案数。
总数显然是 。
而不会越狱当且仅当相邻两个监狱的信仰不同。
所以第一个可以有 种选择,第二个可有 种选择(不能选第一个人选的),第三个有 种选择(不能选第二个人的)……
所以答案为 。需要注意减法取模的情况。
例2.错排问题
设 表示当排列长度为 时该问题答案。
考虑递推,假设 放在了第 个位置,则有 种选择。
则此时第 个信封有两种选择:
-
放在第 个位置,则剩余 个再次错排,有 种方案。
-
不放在第 个位置,则剩余 个错排,有 种方案。
上述先放 ,再放 ,分两步,用乘法原理。第二步分两类,用加法原理。
则 。
先选几个固定后,剩下的错排即可。
很小,使得我们可以杨辉三角递推组合数。
注意到 固定,所以递推时可以对 取模,一个组合数能够整除 当且仅当该数字 后为 。
最后做二维前缀和即可做到每次 查询。
例5.方程的解
给你一个不定方程 ,求正整数解个数。
首先, 可以通过快速幂得到。
考虑这样一个问题:
给你 个相同的球与 个不同的盒子,每个盒子至少放一个,求将这 个球全部放入 个盒子的方案数。
将这 个球排成一排,最左侧和最右侧分别放上一个板子,如下:
这样问题转换成放入 个板子(不算两边的板子),两个板子之间至少有一个球,求方案数。
如下即为一个 的合法方案:
最终每个盒子放入的球就是两个板子之间的球的数量。
注意到球之间共有 个空隙,我们只要在选择 个空隙,对每个空隙分别插入一个板子就一定能够构造出一个合法方案。
问题转换成从 个东西中选出 个东西,求方案数。
答案很显然了,!
回到上面的问题,我们抽象一下,将 个 排成一排,有 个数,初始为 ,每个数要选择一些 加到自己上面,不能不选,求方案数。
这不就是放球问题吗?
这样问题就完美解决,答案即为 。
不过,该题没有模数,需要写高精。
扩展:给定不定方程 ,求非负整数解的个数。
直接求显然不好做。
考虑做一些数学变换。
上面的 是换元的 。
发现 的取值与 一一对应。
所以我们求出 的解数就可以求出 的解数。
而因为 ,所以这题做完了。
例6.[USACO09FEB] Bulls And Cows S
枚举放多少头公牛,假设放 头。
则母牛剩下 头。
问题变成放 个板子隔成,左、右留的位置无限制,但板子之间必须空 个以上。
写成不定方程的形式就是
考虑像换元一样,将左边除左右两边的所有项减去 ,原式可化为:
这样就是不定方程非负整数解板子了,放置公牛 个的答案即:
最终答案即:
按原式计算即可。
例7.车的放置
先考虑一个 行 列的普通矩形放 个怎么做。
首先肯定是在 行中选 列,方案数为 。
其次,选第一个所占得行时有 种方案,第二个有 种。
所以方案数有 种。
接着将原图形分成两个矩形,这里设分成 与 的矩形。
先放第一个,枚举放多少个(设为 个),方案有 种。
再放第二个,设第二个放了 个,但是放第二个时由于第一个放 个后第二个的行只有 个了。
所以方案数是 。
乘起来就是答案。
一个 的网格上有 个格点,为了叙述方便,以下记 。
先不考虑不在同一条直线上的限制,就是在所有点中任选 个,方案显然有 种。
考虑在此基础上减掉不合法的方案。
显然,不合法方案分成横着、竖着、斜着的。
对于横着,一共 行,每一行都有 种方案,乘起来即可。
竖着的同理。
斜着的有两种情况:斜率正与斜率负。
显然他们是对称的,即方案数相同。
只需要求出斜率为正的方案,最终乘 即可。
注意到我们可以枚举两个不在同一水平线且不在同一竖直线上的点,这样它们中间的点的个数(不包括两端)就是答案。
因为任选中间的一个点与两端的两个点总是可以在同一水平线上。
引理:在格点上,若一个点坐标为 ,另一个点坐标为 ,则它们连线上格点个数(不包括两端)为:
证明:
为叙述方便,以下把绝对值去掉。
对于斜边上任意一点 ,满足 。
来一张图。
记 。
则:
若 ,则 。
则一定存在一个 ,使得 。
则 。
经过一堆不等式变换后(或者按原不等式取值范围)得到 。
最终得到 。
即 最多有 种取值,即最多有 种合法的 。
掐头去尾即可证明原命题。
证毕。
直接枚举斜边的端点,计算即可。
事实上,直接枚举是 的,需要进一步优化。
实际上,我们固定形状后,格子图中有多少个这类形状的三角形已经确定了。
假设底、高分别为 ,则数量显然为 。
证明留给读者思考。
复杂度
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;
}
多重集排列数
多重集是一种广义的集合。
普通的集合是不可重的,即集合中每个元素都不相同。
多重集是一种可重集合,可以表示成这么一种形式:
上述可重集 就是描述了一个集合内有 个 、 个 …… 个 。
所谓排列数,就是求 有多少个本质不同的排列。
记 。
则上述问题的答案为:
简证:
不考虑数字相同的情况,共有 种。
每一个相同的数字内部总共会有 种,需要除掉。
证毕。
将 与 之间连有向边,容易发现一定会形成若干个有向环。
若交换完成,则一定是 个自环。
而最优解一定是每次交换都能够让至少一个数字复位,在图上的表现就是让一个大环断成两个小环。
设 表示让一个长为 的大环分解成 个自环的方案数。
再设 表示将一个长为 的环分成长为 和长为 的环的方案数。
注意到若 ,,否则 。
因为第一种可能会有重复情况,画图即可理解,这里不再赘述。
有转移:
这样就完了。……吗?
不然,考虑上述乘法原理,我们先把 环消掉再做 环,显然漏了情况。
其实可以先将 环分成两个,再将 环分成两个……
交替着做也是一种方案!肿么办?
发现一个长为 的环操作次数必然是 。
则上面的操作总数一共是 ,要从中选 次操作给第一个环,剩下的给第二个环,有多少方案?
这不就组合数吗!
或者也可以按多重集排列数理解, 个 , 个 有多少排列方案。
最终答案是一样的。
所以真正的方程应为:
最后 dfs 找环贡献加起来就行。
但是这样时间复杂度是 的,无法通过。
事实上, 有个通项公式:
这样用快速幂就可以 求出答案了。
等等……为啥是 ?
设第 个环的长度为 ,一共有 个环。
则有 。
所以千万不要预处理 数组!
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;
}
最后的问题, 的通项是怎么来的?
观察可得、oeis大法
证明:
引理:有标号无根树的方案数为 (Cayley 公式)。
证明:根据 Prüfer序列 的性质,任何长为 的值域为 的序列都可以唯一确定成一棵有标号树。
根据乘法原理,方案数显然为 。
证毕。
考虑一个有 个点的环,显然需要 次交换就可以分成 个自环。
假设有一张图 ,开始时没有边。
每次若交换 ,则将 之间连一条无向边。
显然,这个图不会出现环,否则一定是断开的环之间重连,一定不是最优解。
而又因为该图有 条边,所以操作完成后 必定是一棵树。
考虑将操作倒过来,不断让两个环之间交换,最后合成一个大环,这样就转化为在树上删边。
若将树的每条边给一个顺序,按照这个顺序进行删边并合并环,则最后得到的环一定是不变的。
根据引理,边没有顺序的有标号树有 棵,加上边的顺序后就是 。
这样就把边有标号的环和边有标号的有标号树建立了一一对应关系。
又因为合并顺序不同的环在不考虑合并顺序的差异后均为原环,本质上是一样的,有重复。
所以要去重,要在原来基础上除掉边重复的方案即 。
所以方案数就是 。
证毕。
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】凌霞软件回馈社区,博客园 & 1Panel & Halo 联合会员上线
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】博客园社区专享云产品让利特惠,阿里云新客6.5折上折
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 在鹅厂做java开发是什么体验
· 百万级群聊的设计实践
· WPF到Web的无缝过渡:英雄联盟客户端的OpenSilver迁移实战
· 永远不要相信用户的输入:从 SQL 注入攻防看输入验证的重要性
· 浏览器原生「磁吸」效果!Anchor Positioning 锚点定位神器解析