论组合数：胎教级教学

1前言

本文将介绍求组合数的方式
主要方法有：动态规划，公式，卢卡斯定理，高精度+线性筛
由于前置知识过多，建议先阅读我数学专题的所有文章
我还会填我之前费马小定理的坑
同时会介绍二项式定理
本文力求详细，所以篇幅较长
建议配合acwing基础课使用

2问题

现有n个物品，从中取出x个
无序则为C(x,n)（组合），有序则为A(x,n)（排列）
由于组合数的应用较广，我们只讨论组合数（其实排列和组合是很像的）

3动态规划

我们根据定义的，对于一个组合C(a,b)
只有a<=b才有意义，而且a等于0时C(a,b)为1，这个很好证
我们可以表示出状态dp[x][y] = C(y,x)即从x个数取y个数的组合
那么，状态转移方程怎么写？
根据容斥原理，如图

由图得状态转移方程dp[i][j] = dp[i-1][j]+dp[i-1][j-1]
有人可能会疑惑，这样递推，不是把dp[i-1][j]中不含j的算了两遍吗？
其实不然，因为有i-1，默认选取j，所以数量是正确的
适宜解决多次询问，数据较小的问题
复杂度O(n^2)，询问O(1）

附代码

#include<bits/stdc++.h>
using namespace std;
int dp[2222][2222];
const int MOD = 1e9+7;
void init(){
	for(int i = 0;i<=2010;i++){
		for(int j = 0;j<=i;j++){
			if(j==0){
				dp[i][j] = 1;
			}else{
				dp[i][j] = (dp[i-1][j-1]+dp[i-1][j])%MOD;
			}
		}
	}
}
int main(){
	int T;
	cin>>T;
	init();
	while(T--){
		int x,y;
		cin>>x>>y;
		cout<<dp[x][y]<<endl;
	}
	return 0;
}

4公式

直接给公式
C(a,b) = b!/((b-a)!a!)
这个是怎么推的呢，要证就需要A（排列）来帮助了
首先，n个数全排列（即所有顺序）有n!种
每一项都在剩下的数取一个
n*(n-1)....*2*1 = n!
然后考虑A(x,y)，x<y，不可以取完
则有y*(y-1)...*(y-x+1)
化简为y!/(y-x)!，现有x个数，有序变无序，直接除以全排列，简单粗暴，推出公式
那么问题来了，怎么算
首先，我们需要用到阶乘，可以预处理
然后，要模一个大的质数（具体看题，一般为1e9+7）,这下除法还要逆元
我们发现模的这个数是质数，这真好办，直接快速幂，可保证有解（参考费马小定理）
费马小定理
依旧给出结论
若存在整数 a ，p 且g c d ( a ， p ) = 1，即二者互质则
a^(p−1) ≡ 1 (modp)
证明
还记得小学乘法表吗，可以发现3,7,9在表内的最后一位均不同，把表扩大，发现这是周期性变化
而且周期就是10
这又是为什么
原来3和10互质，最小公倍数为30，30的倍数mod10 = 0
30的倍数-3mod10就等于7
以此类推
但是证出这个有什么用呢？
因为3^0 mod 10 = 1
所以3^1 mod 10 = 3
继续推可得3^x mod 10 = (x*3)mod10
只要x*3mod10是周期性变化的，3^x也会是周期性变化的，周期还相同，都为10
根据上面的乘法表，x*3mod10就是周期性的
至此，a = 3,p = 10的特殊值证完了，那么怎么保证所有a^(p-1)modp都等于1呢
很简单，任何数的0次幂都等于1！
周期为p-1，则a^(p-1)modp=1
我们提出一个a，留下a^(p-2)，直接成了a的逆元
作者为了保障尽可能通俗易懂，篇幅较长，建议多次阅读证明部分
适宜解决组数和数据范围较均衡的问题
O(n log n)
又附代码

#include<bits/stdc++.h>
using namespace std;
const long long MOD = 1e9+7;
long long f[114514];
long long c[114514]; 
long long quick_pow(long long a,long long b){
    if(b==0){
    	return 1;	
	}
    if(b&1){
    	return a*quick_pow(a,b-1)%MOD;
	}
    long long x = quick_pow(a,b/2);
    return x*x%MOD;
}
void init(){
	f[0] = 1;
	c[0] = 1;
	for(long long i = 1;i<=100010;i++){
		f[i] = f[i-1]*i%MOD;
		c[i] = quick_pow(f[i],MOD-2); 
	}
}
int main(){
	init();
	long long T;
	cin>>T;
	while(T--){
		long long a,b;
		cin>>a>>b;
		long long h1 = c[a-b]%MOD;
		long long h2 = c[b]%MOD;
		long long x = (((f[a]*h1)%MOD)*h2)%MOD;
		cout<<x<<endl; 
	}
	return 0;
}

5卢卡斯定理

还是先给结论
C(a,b) ≡ C(a%p,b%p)*C(a/p,b/p)（mod p)
其中保证p为质数，不然要用exlucas
证明

步骤1
我们可以发现，a可以表示为p进制数:ak*p^k + ak-1*p^k-1......+a0*p^0
举个例子，13表示为二进制:1*8+1*4+0*2+1*1
同理可得b = bk*p^k + bk-1*p^k-1......+b0*p^0
这个等会会用到

步骤2
二项式定理，对于(x+y)^n，展开之后每项x,y的指数和为n
那么，当y的指数为k时，x的指数为（n-k），合并同类项后系数为C(k,n)，根据组合数定义即可证

步骤3
C(a,b)(0<a<b)
则C(a,b)mod b = 0
这个特别好证，用公式展开，分子都含有b了 

步骤4
（1+x）^p
可用二项式定理展开
1无论几次方都是1，不用看，所以可得原式等于
C(k,p)x^k(0<=k<=p)求和
由步骤3得
在modp的情况下，所有(0<k<p)可消
剩下情况组合都等于1，也可消成1作为系数
则原式等于1+x^p
整理一下(1+x)^p = 1+x^p
可以直接开括号

步骤5

铺垫得够多了，开始证明
由二项式定理得C(b,a)为(1+x)^a中x^b的系数
初中数学学过同底数指数幂相乘，底数不变指数相加，结合步骤一，我们可以发现
（1+x）^a等价于(1+x)^a0((1+x）^p)^a1......((1+x)^(p^k))^ak     (我们发现展开后有k位，k取决于log(p,a）即以p为底a的对数，可算时间复杂度)
相当于是展开了
我们单看每项
((1+x)^(p^k))^ak
我们直接逆用步骤4
(1+x)^(p^k*ak)
这不是二项式定理的一般形式吗
直接二项式定理
C(bk/p,ak/p)为(1+x)^(p^k*ak)中x^(p^k*bk)的系数
把p消掉则有C(bk,ak)为(1+x)^(p^k*ak)中x^(p^k*bk)的系数
单独解决了每一项，我们乘起来，运用步骤1把x^(p^k*bk)合并
可得原式即(1+x)^a0((1+x）^p)^a1......((1+x)^(p^k))^ak可化为
C(b0,a0)C(b1,a1)......C(bk,ak)x^b
上面证过（1+x）^a等价于C(b0,a0)C(b1,a1)......C(bk,ak)x^b（等量代换）
x^b的系数既可以表示为C(b,a)
又可以表示为C(b0,a0)C(b1,a1)......C(bk,ak)
则两者相等
以为对于每个ak,bk都需要除p取模（参考10进制的分割每一位）
可得递归式C(a,b) ≡ C(a%p,b%p)*C(a/p,b/p)（mod p)
证毕
回顾一下，这个证明过程真的长，给作者CPU干烧了
建议lucas模板变成蓝，CRT变成绿，难度换一下
lucas证明最妙的一步是把a分解
然后产生两个等价的式子，不断套用公式，整出来的组合数
再逆用公式化简，竟然搞出来了一堆组合数
分解以后快了不少，适用于p小的情况
复杂度（p log p log n)（log n 以p为底，log p 以2为底）
所以可近似看作(p log n)（此处就都以2为底）
代码要靠快速幂求逆元+递归写出
码量不大，建议理解之后背诵
叒附代码

#include<bits/stdc++.h>
using namespace std;
int p;
long long qmi(int a,int k){
	int res = 1;
	while(k){
		if(k&1){
			res = (long long)res*a%p;
		}
		a = (long long)a*a%p;
		k>>=1; 
	}
	return res;
}
long long C(int a,int b){
	int res = 1;
	for(int i = 1,j = a;i<=b;i++,j--){
		res = (long long)res*j%p;
		res = (long long)res*qmi(i,p-2)%p;
	}
	return res;
}
long long lucas(long long a,long long b){
	if(a<p&&b<p){
		return C(a,b);
	}else{
		return (long long)(C(a%p,b%p)*lucas(a/p,b/p))%p;
	}
}
int main(){
	int n;
	cin>>n;
	while(n--){
		long long a,b;
		cin>>a>>b>>p;
		cout<<lucas(a+b,b)<<endl;
	} 
	return 0;
}

6高精度+线性筛

不让模的问题真的恶心死人了
直接从定义出发
C(b,a) = (a*(a-1)...*(a-b+1))/(b*(b-1)...*1)
a,b都是低精度，int可存
所以要实现高精乘和高精除
高精除？不可能写！
因为组合数必定为整数，直接分解质因数，这样就剩高精乘int
组合数必为整数的证明
1定义，无需说
2我自己研发的奇妙证法
情况1:b>=(a/2)
此时2b是分子因数
情况2:b<(a/2)
此时b是分子的因数
证毕
那么，怎么分解呢
首先，线性筛筛质数
然后，用a依次除以p^k(1<=k<=log(p,a)),加和即可
适用问题：不能模的问题
无须考虑复杂度，没TLE之前数组先炸了，会导致无解，这种情况不存在
我知道大家需要什么
叕附代码

#include<bits/stdc++.h>
using namespace std;
const int N = 5555;
int prime[N],cnt,st[N],sum[N];
void get_prime(int n){
    for(int i=2;i<=n;i++){
        if(!st[i]){
        	prime[cnt++] = i;
		}
        for(int j=0;prime[j]<=n/i;j++){
            st[prime[j]*i] = true;
            if(i%prime[j]==0) break;
        }
    }
}
int aaa(int n,int p){
    int res = 0;
    while(n){
        res+=n/p;
        n/=p;
    }
    return res;
}
vector<int> mul(vector<int>&A,int b){
    vector<int> C;
    int t = 0;
    for (int i = 0;i<A.size()||t;i++){
        if (i < A.size()) t += A[i] * b;
        C.push_back(t % 10);
        t /= 10;
    }
    while(C.size()>1&&C.back()==0){
    	C.pop_back();
	}
    return C;
}
int main(){
    int a,b;
    cin>>a>>b;
    get_prime(a);
    for(int i=0;i<cnt;i++){
    	sum[i] = aaa(a,prime[i])-aaa(a-b,prime[i])-aaa(b,prime[i]);
	}
    vector<int> ans;
    ans.push_back(1);
    for(int i=0;i<cnt;i++){
        for(int j=0;j<sum[i];j++){
            ans = mul(ans,prime[i]);
        }
    }
    for(int i=ans.size()-1;i>=0;i--){
    	cout<<ans[i];
	}
}

7后记
本文一共提及了四种求组合数的方式
有的代码难写，有的思路难想
不同的方法没有优劣，适合不同的题目
还有exlucas没有完成，需要excrt作为前置知识，又双叒叕挖坑了，会补的
本文作者是个蒟蒻，欢迎各位大牛指点