求组合数的三种算法

一、预处理组合数

核心：

$$C_a^b=C_{a-1}^b+C_{a-1}^{b-1}$$

适用范围：$a$ 较小的情况下，如 $a\le 10^3$。
算法简析：令 $\text{C[n][k]}=C_n^k$，规定 $\text{C[0][0] = 1}$，则

$$\begin{array}{r}\text{C[n][k]}=\{\begin{array}{ll}1& ,k==0\\ \text{C[n - 1][k] + C[n - 1][k - 1]}& ,0<k\le nandn\ge 1\end{array}\end{array}$$

#define MAX 4000
#define MOD 6662333

int C[MAX][MAX], n;

void solve(void)
{
	for (int i = 0; i <= n; 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; 
}

预处理后，直接访问 $\text{C[n][k]}$，就能得到 $C_n^k$ 的值。

---

二、预处理阶乘

核心：

$$C_a^b=\frac{a!}{b!(a-b)!}=a!\ast (b!{)}^{-1}\ast ((a-b)!{)}^{-1}$$

适用范围：$a$ 较大的情况下，如 $a\le 10^5$。
算法简析：
先来看逆元和费马小定理的定义：

• 1、逆元：对于模 $m$，有两个数 $a$ 和 $b$，若 $ab\equiv 1(\text{mod}m)$，即 $a$ 和 $b$ 的乘积除以 $m$ 的余数为1，则 $a$ 和 $b$ 互为模 $m$ 意义下的乘法逆元。记 $b=a^{-1},a=b^{-1}$。
• 2、费马小定理：若 $p$ 是一个素数，且 $a$ 不被 $p$ 整除，则 $a^{p-1}\equiv 1(\text{mod}p)$，即 $a^{p-1}$ 除以 $p$ 的余数为1。
  由模运算的性质，$\equiv$ 两边同乘 $a$ 的逆元 $a^{-1}$，得 $a^{p-2}\equiv a^{-1}(\text{mod}p)$。在模 $p$ 的意义下，$a^{p-2}$ 和 $a^{-1}$ 等价。求 $a^{-1}$，就转换为求 $a^{p-2}$。

现在，我们的重点是求阶乘和阶乘的逆元。我们用两个数组 fact[] 和 infact[] 分别表示阶乘(fact[a] 为 $a!$)和阶乘的逆元(infact[a] 表示 $(a!{)}^{-1}$)。规定 $\text{fact[0] = infact[0] = 1}$。

• 1、阶乘：由

$$a!=(a-1)!\ast a$$

得，

$$\text{fact[a] = fact[a - 1] * a}$$

• 2、阶乘的逆元：由费马小定理，在模 $p$ 下，

$$\begin{array}{rl}(a{)}^{p-2}& \equiv (a{)}^{-1}(\text{mod}p)\\ (a!{)}^{-1}& =((a-1)!{)}^{-1}\ast a^{-1}\end{array}$$

得，

$$\begin{array}{rl}\text{infact[a]}& =\text{infact[a - 1]}\ast a^{-1}\\ a^{-1}& =\text{pow(a, p - 2) mod p}\end{array}$$

注：计算逆元时，可以通过快速幂来提高算法效率。

#define MAX 4000
#define MOD 6662333

int fact[MAX], infact[MAX], n;

typedef long long ll;

int qum(int x, int n, int mod)
{
	ll ret = 1;
	while (n > 0)
	{
		if (n & 1)    ret = ret * x % MOD;
		x = (ll)x * x % MOD;
		n >>= 1;
	}
	return ret;	
} 

void solve(void)
{
	fact[0] = infact[0] = 1;
	for (int i = 1; i <= n; i++)
	{
		fact[i] = (ll)fact[i - 1] * i % MOD;
		infact[i] = (ll)infact[i - 1] * qum(i, MOD - 2, MOD) % MOD;
	}
}

预处理后，$C_n^k=\text{fact[n] * infact[k] * infact[n - k]}$。注意：计算过程中可能会溢出，要进行模运算。

---

三、卢卡斯定理

核心：卢卡斯定理

$$C_a^b\equiv C_{a\text{mod}p}^{b\text{mod}p}\cdot C_{\lfloor a/p\rfloor }^{\lfloor b/p\rfloor }(\text{mod}p)$$

适用范围：$a$ 很大的情况，比如 $a\le 10^{18}$。
算法简析：若 $a$ 和 $b$ 很大，我们可以通过卢卡斯定理缩小 $a$ 和 $b$，直至 $a,b<p$ ($p$ 一般是较小的素数)。这时，在使用前两种方法求解。

#define MAX 4000
#define MOD 6662333

int fact[MAX], n;

typedef long long ll;

int qum(int x, int n, int mod)
{
	ll ret = 1;
	while (n > 0)
	{
		if (n & 1)    ret = ret * x % MOD;
		x = (ll)x * x % MOD;
		n >>= 1;
	}
	return ret;	
} 

void init(void)
{
	fact[0] = 1;
	for (int i = 1; i <= n; i++)
	{
		fact[i] = (ll)fact[i - 1] * i % MOD;
	}
}

int C(int a, int b, int mod)
{
	if (a < b)    return 0;
	return (ll)fact[a] * qum(fact[b], mod - 2, mod) % mod * qum(fact[a - b], mod - 2, mod) % mod;
}

int lucas(int a, int b, int mod)
{
	if (a < mod && b < mod)    return C(a, b, mod);
	return (ll)C(a % mod, b % mod, mod) * lucas(a / mod, b / mod, mod) % mod;
}

---

完