快速求组合数的方法

求C(n,m)%mod的方法总结

1.当n,m都很小的时候可以利用杨辉三角直接求。

C(n,m)=C(n-1,m)+C(n-1,m-1)；

2.利用乘法逆元。

乘法逆元：(a/b)%mod=a*(b^(mod-2)) mod为素数。

逆元可以利用扩展欧几里德或欧拉函数求得：

1).扩展欧几里德：b*x+p*y=1 有解，x就是所求

2).费马小定理：b^(p-1)=1(mod p),故b*b^(p-2)=1(mod p),所以x=b^(p-2)

1. n!/(m!*(n-m)! =x%p ,先对算出n！、m！、(n-m)!对p取模的余数，就转换为a/b=x%p;因为p为素数，所以等价于bx+py=a;然后用扩展的欧几里得定理算出 bx’+py’=1的解，x=x’*a，就得到了最终的x的值，即C(m,n)%p得值。

2.逆元 其实如果mod是素数 则b的逆元其实就是b^(mod-2)

也就是 (m!(n-m)!)的逆元为 (m!(n-m)!)p-2 ；

int inv(int a)
{
	//return fpow(a, MOD-2, MOD); 
	return a==1?1:(long long)(MOD-MOD/a)*inv(MOD%a)%MOD; 
}
LL C(LL n,LL m) 
{
	if(m<0)return 0; 
	if(n<m)return 0; 
	if(m>n-m)m=n-m; 
	LL up=1,down=1; 
	for(LL i=0;i<m;i++)
	{
		up=up*(n-i)%MOD; 
		down=down*(i+1)%MOD; 
	}
	return up*inv(down)%MOD; 
} 

3.当n和m比较大，mod是素数且比较小的时候（10^5左右），通过Lucas定理计算

Lucas定理：A、B是非负整数，p是质数。A B写成p进制：A=a[n]a[n-1]…a[0]，B=b[n]b[n-1]…b[0]。 

则组合数C(A,B)与C(a[n],b[n])C(a[n-1],b[n-1])…*C(a[0],b[0]) mod p同余

即：Lucas(n,m,p)=C(n%p,m%p)*Lucas(n/p,m/p,p)

#include<iostream>
//#include<algorithm>
using namespace std;
typedef long long ll;
int quick_power_mod(int a,int b,int m)//pow(a,b)%m
{
	int result = 1;
	int base = a;
	while(b>0)
	{
		if(b&1==1)
		result=(result*base)%m;
		base=(base*base)%m;
		b>>=1;
	}
	return result;
}
//计算组合数取模
ll comp(ll a,ll b,int p)//composite num C(a,b)%p
{
	if(a<b)return 0;
	if(a==b)return 1;
	if(b>a-b)b=a-b;
	int ans=1,ca=1,cb=1;
	for(ll i=0;i<b;i++)
	{
		ca=(ca*(a-i))%p;
		cb=(cb*(b-i))%p;
	}
	ans=(ca*quick_power_mod(cb,p-2,p))%p;
	return ans;
}
ll lucas(ll n,ll m,ll p)
{
	ll ans=1;
	while(n&&m&&ans)
	{
		ans=(ans*comp(n%p,m%p,p))%p;//also can be recusive
		n/=p;
		m/=p;
	}
	return ans;
}
int main()
{
	ll m,n;
	while(cin>>n>>m)
		cout<<lucas(n,m,10007)<<endl;
	return 0;
}

C(n % mod, m % mod) % mod; 如果太大还是利用上面的逆元来处理。

半预处理

由于Lucas定理保证了阶乘的数均小于p，所以可以讲所有的阶乘先预处理，优化C（n,m）

mod的要求：p<10^6，且为素数

有效范围：1<=n,m<=10^9

//半预处理 
const ll MAXN=100000; 
ll fac[MAXN+100]; 
void init(int mod) 
{
	fac[0]=1; 
	for(int i=1;i<mod;i++)
		fac[i]=fac[i-1]*i%mod; 
} 
//半预处理逆元求C（n，m）%mod 
ll C(ll n,ll m) 
{ 
	if (m>n)return 0; 
	return fac[n]*(GetInverse(fac[m]*fac[n-m],mod))%mod; 
}

可以预处理出阶乘和阶乘的逆元。

int C(int n,int m)
{
	return 1ll*fac[n]*inv[m]%p*inv[n-m]%p;
}    
fac[0]=inv[0]=inv[1]=1;
for(int i=1;i<=N;i++)fac[i]=1ll*fac[i-1]*i%p;
for(int i=2;i<=N;i++)inv[i]=1ll*(p-p/i)*inv[p%i]%p;
for(int i=1;i<=N;i++)inv[i]=1ll*inv[i-1]*inv[i]%p;