组合数计算

\[{n \choose m}={n! \over (n-m)!m!} \]

1 .定义法求组合数

Code

LL C(LL a,LL b,LL MOD)
{
	if(b>a) return 0;
	LL res=1;
	for(LL i=1,j=a;i<=b;i++,j--) {
		res=res*j%MOD;
		res=res*quickpow(i,MOD-2,MOD)%MOD;
	}
	return res%MOD;
}

时间复杂度：\(O(m)\)

2 .递推法求组合数

Code

LL c[N][N];
void init()
{
    c[0][0]=1;
    for(int i=1;i<=2000+3;i++)
        for(int j=0;j<=i;j++) {
            if(j) c[i][j]=(c[i][j]+c[i-1][j-1])%MOD;
            c[i][j]=(c[i][j]+c[i-1][j])%MOD; 
        }
}

时间复杂度:\(O(nm)-O(1)\).

3 .阶乘逆元法求组合数

Code

LL power(LL x,LL k,LL MOD)
{
    LL res=1; x%=MOD;
    while(k) {
        if(k&1) res=res*x%MOD;
        x=x*x%MOD; k>>=1;
    }
    return res%MOD;
}

inline LL C(LL n, LL m)
{
    return fact[n]*power(fact[n-m]*fact[m],MOD-2,MOD)%MOD;
}

预处理部分

	fact[0]=1;
	for(LL i=1;i<N;i++) 
		fact[i]=fact[i-1]*i%MOD;

时间复杂度：\(O(n)-O(\log n)\).
如果预处理时同时处理逆元的话，
可以做到 \(O(n \log n)-O(1)\).

4 .lucas定理求组合数

只适用于模质数的情况
设预处理时间为 \(f(p)\)，单次求组合数的时间为 \(g(p)\)
时间复杂度为 \(\mathcal O(f(p)+g(p)\log _{p}m)\)

故一般要求模数要小。

Code

typedef long long LL;

const int N=2e5+5;

LL power(LL x,LL k,LL MOD)
{
	LL res=1; x%=MOD;
	while(k) {
		if(k&1) res=res*x%MOD;
		x=x*x%MOD; k>>=1;
	}
	return res%MOD;
}

LL fact[N],infact[N];

LL C(int n,int m,LL p) 
{ 
	if(n<m) return 0; // 注意这个特判。
	else return fact[n]*infact[n-m]%p*infact[m]%p; 
}
LL lucas(int n,int m,LL p)
{
	if(n<p && m<p) return C(n,m,p);
	else return C(n%p,m%p,p)*lucas(n/p,m/p,p)%p;
}

预处理部分

	fact[0]=infact[0]=1;
	for(int i=1;i<=p;i++) {
		fact[i]=fact[i-1]*i%p;
		infact[i]=power(fact[i],p-2,p);
	}

y 总原来那个复杂度好像不大对。
目前这个是 \(\mathcal O(p \log p + \log_pm)\)
如果只递推阶乘，现算逆元的话可以 \(\mathcal O(p + (\log_pm)\cdot(\log_2 p))\)

Code:

LL power(LL x,LL k,LL MOD)
{
	LL res=1; x%=MOD;
	while(k) {
		if(k&1) res=res*x%MOD;
		x=x*x%MOD; k>>=1;
	}
	return res%MOD;
}

LL fact[N];

LL C(int n,int m,LL p) 
{ 
	if(n<m) return 0;
	else return fact[n]*power(fact[n-m],p-2,p)%p*power(fact[m],p-2,p)%p; 
}
LL lucas(int n,int m,LL p)
{
	if(n<p && m<p) return C(n,m,p);
	else return C(n%p,m%p,p)*lucas(n/p,m/p,p)%p;
}

fact[0]=1;
for(int i=1;i<=p;i++) fact[i]=fact[i-1]*i%p;

5. 扩展卢卡斯定理求组合数 exlucas

P4720 【模板】扩展卢卡斯
题解

要求模的数较小，不过可以不是质数。

\(n,m\le 10^{18},p \le 10^6\)

Code：

typedef long long LL;

LL muler(LL x,LL k,LL MOD)
{
	LL res=0; x=(x%MOD+MOD)%MOD; k=(k%MOD+MOD)%MOD; 
	while(k) {
		if(k&1) res=(res+x)%MOD;
		x=(x+x)%MOD; k>>=1;
	}
	return res%MOD;
}

LL power(LL x,LL k,LL MOD)
{
	LL res=1; x%=MOD;
	while(k) {
		if(k&1) res=muler(res,x,MOD);
		x=muler(x,x,MOD); k>>=1;
	}
	return res%MOD;
}

LL exgcd(LL a,LL b,LL& x,LL& y)
{
	if(b==0) {
		x=1; y=0;
		return a;	
	}
	LL z=exgcd(b,a%b,y,x);
	y-=a/b*x;
	return z;
}

LL inv(LL x,LL p) 
{
	LL y,z; exgcd(x,p,y,z);
	return (y%p+p)%p;  
}

LL excrt(int n,LL b[],LL a[])
{
	LL m=a[1],ans=b[1];
	for(int i=2;i<=n;i++) {
		LL y,z,d=exgcd(m,a[i],y,z);
		if((b[i]-ans)%d!=0) return -1;
		y=muler(y,(b[i]-ans)/d,a[i]/d);
		
		ans+=y*m;
		m=a[i]/d*m;
		ans=(ans%m+m)%m;
	}
	return ans;
}

LL divide_p(LL n,LL p,LL pk)
{
	if(!n) return 1;
	LL res=1;
	for(LL i=1;i<=pk;i++) 
		if(i%p) res=res*i%pk;
	res=power(res,n/pk,pk); 
	for(LL i=n%pk;i>=1;i--) 
		if(i%p) res=res*i%pk;
	return res*divide_p(n/p,p,pk)%pk;
}

LL Getp(LL n,LL p) 
{
	LL res=0;
	while(n) res+=n/p,n/=p;
	return res;
}

LL C(LL n,LL m,LL p,LL pk) 
{
	return divide_p(n,p,pk)*inv(divide_p(n-m,p,pk),pk)%pk
				*inv(divide_p(m,p,pk),pk)%pk
				*power(p,Getp(n,p)-Getp(n-m,p)-Getp(m,p),pk)%pk;
}

LL exlucas(LL n,LL m,LL p)
{
	static LL A[1024],B[1024];
	
	if(m>n) return 0;
	int tot=0;
	LL t=p;
	for(LL i=2;i*i<=t;i++) {
		if(t%i==0) {
			A[++tot]=1;
			for(;t%i==0;t/=i) A[tot]*=i;
			B[tot]=C(n,m,i,A[tot]);
		}
	}
	if(t>1) A[++tot]=t,B[tot]=C(n,m,t,t);
	return excrt(tot,B,A);
}

6.阶乘分解求组合数

适用于 没有模数 或者 模数是合数并且 \(n,m\) 较小的情况。

\(n,m \le 10^6,p \le 10^{18}\)

高精度 Code:

vector<int> mul(vector<int> a,int b)
{
	vector<int> c;
	int t=0,i;
	for(i=0;i<(int)a.size();i++) {
		t=t+a[i]*b;
		c.push_back(t%10);
		t/=10;
	}
	while(t>0) {
		c.push_back(t%10);
		t/=10;
	}
	return c;
}

const int N=10010;

int p[N],tot=0;
bool tag[N];

void primes(int n)
{
	int i,j;
	for(i=2;i<=n;i++) {
		if(!tag[i]) p[++tot]=i;
		for(j=1;j<=tot && p[j]*i<=n;j++) {
			tag[i*p[j]]=true;
			if(i%p[j]==0) break;
		}
	}
	return;
}

int get(int n,int p)
{
	int cnt=0;
   while(n>=p) cnt+=n/p,n/=p;
	return cnt;
}

vector<int> C(int n,int m)
{
	int cnt;
	vector<int> c; c.push_back(1);
	for(int j=1;j<=tot && p[j]<=n;j++) {
		cnt=get(n,p[j])-get(m,p[j])-get(n-m,p[j]);
		while(cnt--) c=mul(c,p[j]);
	}
	return c;
}

时间复杂度: \(O(\log ^{2}C(n,m)+n)\).
可以近似认为是 \(O(n^2)\).
时间复杂度分析：

• 质因数分解复杂度 \(O(n)\)
• \(C(n,m)\) 共有 \(\lg C(n,m)\) 位，因此每次乘法mul需要 \(O(\lg C(n,m))\)
• \(C(n,m)\) 最多有 \(\log_{2}C(n,m)\) 个质因子，因此需要mul \(\log_{2}C(n,m)\) 次

综上，时间复杂度: \(O(\log^{2}C(n,m)+n)\).

如果模的是合数的话吧高精度去掉就可以了，时间复杂度 \(O(n)-O({n \over \ln n})\)

materials

• blog
• simple problems