求组合数的四种方法
求组合数
公式法
- 公式
\(C_a^b=C_{a-1}^{b-1}+C_{a-1}^{b}\) - 公式理解
在b个要抽取的物品中任选一个x物品，第一次抽取情况无非两种，抽到了x和没有抽到x，如果抽到了x那么就只需要在剩下的a-1个中再抽取b-1个即可，如果没有抽到x，就在剩下的a-1个中抽取b个，两次相加即为抽取的总种数。 - 代码思路及实现
我们只需要将数据范围n内的所有组合数都预处理出来，就可以实现求任一个组合数。时间复杂度\(o(n^2)\),本质是一个从1到n的等差数列求和。
给定n组数，每组包含一对a,b,求每一组\(C_a^b\)
#include<iostream>
using namespace std;
const int N=2010;
long long c[N][N];
const int mod=1e9+7;
int main()
{
int n;
cin>>n;
for(int i=0;i<=2000;i++)
{
for(int j=0;j<=i;j++)
{
if(j==0)
{
c[i][j]=1;
}
else
{
c[i][j]=(c[i-1][j]+c[i-1][j-1])%mod;
}
}
}
while(n--)
{
int a,b;
scanf("%d %d",&a,&b);
printf("%ld\n",c[a][b]);
}
return 0;
}
- 代码细节
注意当b=0时的边界判定
定义法
- 算法思路
由组合数定义\(C_a^b=\frac{a(a-1)(a-2)\dots(a-b+1)}{b!}=\frac{a!}{b!(a-b)!}\),这种结果一般会很大所以会要求模上一个数，所以我们算出a!后就不去算除法而是用逆元来算。 - 代码实现
当 \(0\le b\le a\le{10}^5\)如果 用公式法时间复杂度\(o(n^2)\)会超时，而定义法时间复杂度为\(o(nlogn)\)不会超时。
#include<iostream>
using namespace std;
const int mod=1e9+7;
const int N=1e5+10;
typedef long long LL;
LL fant[N],infant[N];//fant[i]代表i的阶乘，infant[i]代表i的阶乘的逆元
LL quick_mod(int a,int b,int p)//这道题mod是一个质数所以可以用快速幂来求
{
LL res=1;
while(b)
{
if(b&1)
{
res=res*a%p;
}
b>>=1;
a=(LL)a*a%p;
}
return res;
}
int main()
{
int n;
cin>>n;
fant[0]=infant[0]=1;
for(int i=1;i<N;i++)
{
fant[i]=(LL)fant[i-1]*i%mod;
infant[i]=(LL)infant[i-1]*quick_mod(i,mod-2,mod)%mod;//第29行
}
while(n--)
{
int a,b;
scanf("%d %d",&a,&b);
printf("%lld\n",(LL)fant[a]*infant[b]%mod*infant[a-b]%mod);
}
return 0;
}
- 代码细节
第29行求逆元的过程本质上就是\(\frac{a}{bc}\equiv a*b^{-1}*c^{-1}(mod\enspace p)\),即除以两个数同余先乘上一个数的逆元再乘上后一个数的逆元。
卢卡斯定理
- 定理
\(C_a^b\equiv C_{a\space mod\space p}^{b\space mod\space p}*C_{a/p}^{b/p}(mod\enspace p)\) - 数学推导
要用到一个公式\((1+x)^p\equiv 1+x^p(mod\enspace p)\),这个公式的证明把左边展开即可：\((1+x)^p=C_p^0+C_p^1x+C_P^2x^2+\dots+C_p^px^p\)因为从\(C_p^1x\)到\(c_p^{p-1}x^{p-1}\)都是p的倍数[1]模上一个p之后都为0，所以\((1+x)^p\equiv 1+x^p(mod\enspace p)\)
将a和b用p进制表示为：
\(a=a_kp^k+a_{k-1}p^{k-1}+a_{k-2}p^{k-2}+\dots+a_0p^0\)
\(b=b_kp^k+b_{k-1}p^{k-1}+b_{k-2}p^{k-2}+\dots+b_0p^0\)
对于\((1+x)^a\)有\((1+x)^a=(1+x)^{a_kp^k+a_{k-1}p^{k-1}+a_{k-2}p^{k-2}+\dots+a_0p^0}=(1+x)^{a_kp^k}(1+x)^{a_{k-1}p^{k-1}}\dots(1+x)^{a_0p_0}=(1+x^{p^k})^{a_k}(1+x^{p^{k-1}})^{a_{k-1}}\dots(1+x^{p^0})^{a_0}\)，由等式左边可以看出\(x^b\)的系数为\(C_a^b\),而\(x^b=x^{b_kp^k+b_{k-1}p^{k-1}+b_{k-2}p^{k-2}+\dots+b_0p^0}=(1+x)^{b_kp^k}(1+x)^{b_{k-1}p^{k-1}}\dots(1+x)^{b_0p_0}\),而从等式右边可以看出\(x^{b_kp^k}\)的系数为\(C_{a_k}^{b_k}\)以此类推得\(x^b\)的系数\(C_a^b\)\(=C_{a_k}^{b_k}C_{a_{k-1}}^{b_{k-1}}\dots C_{a_0}^{b_0}(mod\enspace p)\)，因为a,b是用p进制表示的则有：\(a_0=a\space mod\space p,b_0=b\space mod\space p,\\ \lfloor a/p \rfloor\)实际意义为a在p进制下向右移一位,\(\lfloor b/p \rfloor\)实际意义为b在p进制下向右移一位。那么对\(C_{\lfloor a/p \rfloor}^{\lfloor b/p \rfloor}\)再进行一次上面的操作可以得到\(C_{\lfloor a/p \rfloor}^{\lfloor b/p \rfloor}=C_{a_k}^{b_k}C_{a_{k-1}}^{b_{k-1}}\dots C_{a_1}^{b_1}(mod\enspace p)\)
综上有：\(C_a^b=C_{\lfloor a/p \rfloor}^{\lfloor b/p \rfloor}*C_{a\space mod\space p}^{a\space mod\space p}(mod\enspace p)\)
-代码实现
卢卡斯算法时间复杂度对于\(0\le b\le a\le {10}^{18},0\le p \le{10}^5\)为\(o(klog_pn*n*log_2p)\)这里的数据组数k一般比较小
#include<iostream>
#include<algorithm>
using namespace std;
typedef long long LL;
int n;
int qmi(int a,int b,int p)
{
int res=1;
while(b)
{
if(b&1)
{
res=(LL)res*a%p;
}
b>>=1;
a=(LL)a*a%p;
}
return res;
}
int C(LL a,LL b,int p)
{
int res=1;
for(int i=1,j=a;i<=b;i++,j--)
{
res=(LL)res*j%p;
res=(LL)res*qmi(i,p-2,p)%p;
}
return res;
}
int lucas(LL a,LL b,int p)
{
if(b>a)
{
return 0;
}
if(a<p)
{
return C(a,b,p);
}
return (LL)C(a%p,b%p,p)*lucas(a/p,b/p,p)%p;//a%p和b%p一定小于p所以直接用C来算组合数即可
}
int main()
{
cin>>n;
while(n--)
{
LL a,b;
int p;
scanf("%lld %lld %d",&a,&b,&p);
printf("%d\n",lucas(a,b,p));
}
return 0;
}
- 代码细节
当a,b都小于p时候才用C函数来算，其余都继续递归直到小于p。
大精度计算
- 针对类型
可计算a,b较大且结果不模上一个数的组合数 - 算法思路
\(C_a^b=\frac{a!}{b!(a-b)!}\)对其分解质因数可得\(C_a^b=\frac{a!}{b!(a-b)!}=p_1^{\alpha_1}p_2^{\alpha_2}\dots p_n^{\alpha_n}(p_1p_2\dots p_n均为质数)\)对于任一个\(p_i\)它的次数都等于分子\(a!\)所含\(p_i\)的个数减去分母\(b!*(a-b)!\)所含的\(p_i\)的个数，求一个阶乘的某个质因子的个数有公式:\(cnt(a!)=\lfloor \frac{a}{p} \rfloor+\lfloor \frac{a}{p^2} \rfloor+\dots\)[2]一直加到分母比a大时停止(因为分母比a大时向下取整为0加了和没加一样),则结果为\(C_a^b=\prod \limits_{i=1}^{n}p_i^{\alpha_i}\) - 实现步骤
先筛出2-a的所有质数，然后算出每个质数对应的次数，再算出\(\prod \limits_{i=1}^{n}p_i^{\alpha_i}\)即可 - 代码实现
输入a,b求出\(C_a^b(0\le b\le a\le 5000)\)
#include<iostream>
#include<algorithm>
#include<vector>
using namespace std;
typedef long long LL;
const int N=5010;
int primes[N];
int sum[N];
int cnt;
bool str[N];
int get(int x,int p)
{
int res=0;
while(x)
{
res+=x/p;
x/=p;
}
return res;
}
vector<int> mul(vector<int> a,int b)
{
vector<int> c;
int t=0;
for(int i=0;i<a.size();i++)
{
t+=a[i]*b;
c.push_back(t%10);
t/=10;
}
while(t)//这里得写成while,当b大于10的时候t在这里就可能大于10从而需要多次进位
{
c.push_back(t%10);
t/=10;
}
}
void get_primes(int n)
return c;
{
for(int i=2;i<=n;i++)
{
if(!str[i])
{
primes[cnt++]=i;
}
for(int j=0;i<=n/primes[j];j++)
{
str[i*primes[j]]=true;
if(i%primes[j]==0)
{
break;
}
}
}
}
int main()
{
int a,b;
cin>>a>>b;
get_primes(a);
for(int i=0;i<cnt;i++)
{
int p=primes[i];
sum[p]=get(a,p)-get(b,p)-get(a-b,p);
}
vector<int> res;
res.push_back(1);
for(int i=0;i<cnt;i++)
{
for(int j=0;j<sum[primes[i]];j++)
{
res=mul(res,primes[i]);
}
}
for(int i=res.size()-1;i>=0;i--)
{
printf("%d",res[i]);
}
return 0;
}
-代码细节
注意求大精度数相乘的模板
当p为质数时，有公式：\(C_p^b\equiv 0(mod\enspace p)\)证明：\(C_p^b=\frac{p!}{b!(p-b)!}=\frac{p(p-1)!}{b(b-1)!(p-b)!}=\frac{p}{b}*C_{p-1}^{b-1}\)由于p是质数，则b与p互质，b在模p条件下存在逆元，则有：\(C_p^b=p*infact(b)*C_{p-1}^{b-1}(mod\enspace p)\)显然是一个P的倍数则模上p后为0. ↩︎
\(\lfloor \frac{a}{p} \rfloor\)就表示a!中的a项乘积中是p的倍数的个数即a项乘积中能提出的p的个数，同理\(\lfloor \frac{a}{p^2} \rfloor\)就表示a!中的a项乘积中是\(p^2\)的倍数的个数即a项乘积中能提出的\(p^2\)的个数，以此类推就可以得到a!中含p的个数，这里你可能会有疑问，不应该加上二倍\(p^2\)的个数，三倍\(p^3\)的个数...才是所含p的个数吗，其实加上一倍就可以了，因为比如\(p^2\)的倍数的个数在算p的倍数的时候就已经被算过了，毕竟\(p^2\)的k倍数一定是\(p\)的k*p倍数。 ↩︎
浙公网安备 33010602011771号