## $C_a^b$的多种场景下的求法 ###[一、$AcWing$ $885$. 求组合数 $I$](https://www.acwing.com/problem/content/887/) **理论依据**:$\large C_a^b=C_{a-1}^b+C_{a-1}^{b-1}$ 适合场景:$\large a<=2000,b<=2000$ **不严谨的证明**: 有$a$个苹果,现在需要选出$b$个苹果。一共有多少种选法呢? **答**:走到第一个苹果面前,面临两个选择: - 选择它 - 放弃它 如果选择它,将要面对$a-1$个苹果中选$b-1$个苹果的问题 如果放弃它,将要面对$a-1$个苹果中选$b$个苹果的问题 根据加法原理,两种选择方法加在一起,就是方法总数,也就是上面的递推式。 **实现代码** ```cpp {.line-numbers} #include using namespace std; const int N = 2010; const int p = 1e9 + 7; int c[N][N]; void init() { /* 功能:第0列初始化为1,现实意义:C(n,0)=1 理解:从n个元素中,每次取出0个元素,也就是每次都不取出元素,有几种取法? 很明显只有1种取法,就是什么都不取这样一种情况。类似:C(n,0)=C(n,n)=1 这是递归的base case,如果没有这个1打底,以后再怎么加法原理都是白费 */ for (int i = 1; i < N; i++) c[i][0] = 1, c[i][i] = 1; // base case for (int i = 1; i < N; i++) for (int j = 1; j < i; j++) // 递归式最重要!这是根据数学公式得到的递推式 // 从递推式由右向左看过去,需要由小到大填充a,b,并且,m>n c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % p; } int n; int main() { init(); cin >> n; while (n--) { int a, b; cin >> a >> b; printf("%d\n", c[a][b]); } return 0; } ``` ### [二、$AcWing$ $886$. 求组合数 $II$](https://www.acwing.com/problem/content/888/) 上一个办法是把$C_a^b$的值预处理出来了。用的是$c[N][N]$,$N$最大$2010$. 本题的$a$,$b$都是上限$10^5$,如果按上题来,就是$c[10^5][10^5]$, 直接报 $Memory\ Limit\ Exceeded$. 所以不能直接递推求出所有解。 适合场景:$a<=1e5,b<=1e5$ 理论依据:$\large C_a^b=\frac{a!}{(a-b)! * b!}$ 举个栗子: $C_3^2=\frac{3!}{2!\times 1!}=\frac{3 \times 2 \times 1}{2 \times 1}=3$ - 本题由于怕数据太大,要求结果 $mod\ (1e9+7)$,这个 **$MOD=1e9+7$是质数(一般取模的都是质数)**,可以直接使用 **费马小定理+快速幂** 求逆元 - 现在要求计算出 $a!$, ${b!}^{-1}$, ${(a-b)!}^{-1}$,使用两个数组来递推: $fact[i] = i!\ \% \ MOD$ $infact[i] = (i!)^{-1}\ \% \ MOD$ 所以: $$\large C_a^b=\frac{a!}{(a-b)! * b!} = (fact[a] * infact[a-b] * infact[b] )\ \% \ MOD$$ **实现代码** ```cpp {.line-numbers} #include using namespace std; typedef long long LL; const int N = 100010; const int p = 1e9 + 7; int fact[N]; // 阶乘结果数组 int infact[N]; // 阶乘逆元结果数组 // 快速幂 int qmi(int a, int k) { int res = 1; while (k) { if (k & 1) res = (LL)res * a % p; a = (LL)a * a % p; k >>= 1; } return res; } int main() { // base case fact[0] = infact[0] = 1; // 0!=1 (0!)^-1=1 for (int i = 1; i < N; i++) { fact[i] = (LL)fact[i - 1] * i % p; // 阶乘递推 // ①用qmi(i, mod - 2, mod)求出i在% mod场景下的逆元 // ②得到的逆元再乘以infact[i-1]累乘到infact中 // ③注意取模 // ④递归式最重要!这是根据数学公式得到的递推式 infact[i] = (LL)infact[i - 1] * qmi(i, p - 2) % p; } int n; cin >> n; while (n--) { int a, b; cin >> a >> b; printf("%d\n", (LL)fact[a] * infact[b] % p * infact[a - b] % p); } return 0; } ``` > **注意**:因为$AcWing\ 886$的数据范围比$AcWing\ 885$更大,所以$AcWing\ 886$的代码在$AcWing\ 885$中可以直接$AC$,但反过来不行。 ###[$AcWing$ $887$. 求组合数 $III$](https://www.acwing.com/problem/content/889/) #### 1.$Lucas$公式 $$\LARGE C_a^b \equiv C_{a\%p}^{b\%p} * C_{a \div p}^{b \div p} (mod \ p)$$ 适用场景:数据范围:$b<=a<=1e18$,$p<=1e5$,$p \in prime$ 预处理出 $1...p$ 的阶乘$fact[i]$和阶乘的逆元$infacf[i]$,用卢卡斯定理进行回答。 注意:这个$p$是小于$1e5$的,看到要求模$1e9+7$绕行啊!不是干哪个的啊!适合$a,b$很大,$p$很小的场景啊!!! > **$Lucas$定理的理解** 实现的时候,只需要对这个东西继续递归调用$Lucas$定理即可。也就是,对于$C_{a÷p}^{b÷p}$继续调用,最后就会出现一个连乘。 于是我们发现,这个东西又除又模的,很像进制拆解。 没错。 也就是说,这个$Lucas$定理就相当于把$n,m$变成了$p$进制数,然后对$p$进制下的每一位都计算组合数,最后乘起来。 #### 2.如何求解$C_a^b$ $C_a^b=\frac{a!}{(a-b)! \times b!}=\frac{a \times (a-1) \times (a-2) \times ...\times(a-b+1) \times (a-b) \times ... \times 1}{(a-b) \times (a-b-1) \times ... \times 1 \times b!}= \frac{a \times (a-1) \times (a-2) \times ...(a-b+1)}{b!}=\frac{a \times (a-1) \times (a-2) \times ...(a-b+1)}{b \times (b-1) \times (b-2) \times ... \times 2 \times 1}$ 根据这个结论: - 对于$a$来讲,需要变量$j$从$a$一直遍历到$a-b+1$ - 对于$b$来讲,需要变量$i$从$1$一起遍历到$b$ 而且$a$到$a-b+1$其实就是$b$次,比如$10$遍历到$8$,就是$10,9,8$三次,即$b=3$,也就是$10-3+1=8$。 而$i$也是遍历$b$次,真有太有意思了,它们两个可以在一个$b$次的循环中一起变动 (代码能省则省啊~)! ```cpp {.line-numbers} int C(int a, int b, int p) { if (a < b) return 0; int down = 1, up = 1; for (int i = a, j = 1; j <= b; i--, j++) { up = (LL)up * i % p; down = (LL)down * j % p; } return (LL)up * qmi(down, p - 2, p) % p; } ``` **实现代码** ```cpp {.line-numbers} #include using namespace std; typedef long long LL; // 快速幂 int qmi(int a, int k, int p) { int res = 1; while (k) { if (k & 1) res = (LL)res * a % p; a = (LL)a * a % p; k >>= 1; } return res; } // 功能:组合数模板 int C(int a, int b, int p) { if (a < b) return 0; int down = 1, up = 1; for (int i = a, j = 1; j <= b; i--, j++) { up = (LL)up * i % p; down = (LL)down * j % p; } return (LL)up * qmi(down, p - 2, p) % p; } // Lucas公式 int lucas(LL a, LL b, int p) { if (a < p && b < p) return C(a, b, p); return (LL)C(a % p, b % p, p) * lucas(a / p, b / p, p) % p; // 递归套用公式 } int n, p; int main() { cin >> n; while (n--) { // n组询问 LL a, b; cin >> a >> b >> p; cout << lucas(a, b, p) << endl; // 利用lucas公式计算组合数 } return 0; } ``` ### [四、$AcWing$ $888$. 求组合数 $IV$](https://www.acwing.com/problem/content/890/) **不让取模怎么办?** #### 1、与前三道题的区别 不再是数据上限的问题了,而是一直不$mod$,有多大都保留,这无疑可能会爆$long\ long$,需要使用高精度。 #### 2、公式一【从定义出发的组合数公式】: $\large C_a^b=\frac{a \times (a-1) \times ... \times(a-b+1)}{b\times (b-1) \times ... \times 1}=\frac{a \times (a-1) \times ... \times (a-b+1) \times (a-b)!}{ b \times (b-1) \times ... \times 1 \times (a-b)!} =\frac{a!}{b! \times (a-b)!}$ #### 3、是不是我们利用高精度+上面的组合数公式直接算就行了? 不是的,因为$yxc$(大雪菜老师)说,这样的效率太低,不可以,需要再想一个好办法。 #### 4、公式二 【算术基本定理】 算术基本定理:$\large C_a^b=p_1^{\alpha1} \times p_2^{\alpha2} \times p_3^{\alpha3} ... \times p_k^{\alpha k}$ 其中$p_1$,$p_2$...是质数,$\alpha 1$,$\alpha 2$,...是指质数$p_1$,$p_2$,...的个数。 如果能分解质因数成功的话,那么就可以通过 **高精度乘法** 解决掉这个问题。 (1) 那么$p_1$和$p_2$这些东东怎么求? 思路:因为$a,b$都是在$[1,5000]$之间的数,所以可以通过筛质数的方法,提前获取到一个质数数组,然后逐个去看,是不是含有这个质数。就能知道有哪些$p_1,p_2,...$了。 (2) $\alpha 1$,$\alpha 2,...,\alpha k$又该怎么求呢? #### 5、公式三【求质数$p$在$a!$中出现的次数】 $$\large cnt=\lfloor \frac{a}{p} \rfloor + \lfloor \frac{a}{p^2} \rfloor + \lfloor \frac{a}{p^3} \rfloor + ...+ \lfloor \frac{a}{p^k} \rfloor$$ [参考资料](https://blog.csdn.net/spidy_harker/article/details/88414504) **计算$n$的阶乘中质数$p$的个数**: ```cpp {.line-numbers} // 获得n!中包含质数p的个数,本质就是转化为p进制嘛 int get(int n, int p) { int res = 0; while (n) { res += n / p; n /= p; } return res; } ``` **举个栗子**: $ 5!=1 \times 2 \times 3 \times 4 \times 5$,求$2$的个数, 则$1 \sim 5$之间含$2$的个数就是 $\lfloor \frac{5}{2} \rfloor$,就是$2$个。 则$1 \sim 5$之间含$2^2$的个数就是 $\lfloor \frac{5}{2^2} \rfloor$,就是$1$个。 $2^3$就大于$5$,就不再继续。那么一共的个数就是$2+1=3$个。表示在 $5!$这个数字中,存在$3$个$2$,就是$2$中有$1$个$2$,$4$中有两个$2$。 #### 6、算法流程 (1)、筛素数,把$1-5000$之间的素数筛出来 (2)、计算$C_a^b$中每个已求得素数的个数 (3)、利用高精度,计算$C_a^b$$=p_1^{\alpha1} \times p_2^{\alpha2} \times p_3^{\alpha3} ... \times p_k^{\alpha k}$的值 **实现代码** ```cpp {.line-numbers} #include using namespace std; const int N = 5010; int primes[N], cnt; bool st[N]; // 欧拉筛 void get_primes(int n) { for (int i = 2; i <= n; i++) { if (!st[i]) primes[cnt++] = i; for (int j = 0; primes[j] <= n / i; j++) { st[primes[j] * i] = true; if (i % primes[j] == 0) break; } } } // 高精乘低精 void mul(int a[], int &al, int b) { int t = 0; for (int i = 1; i <= al; i++) { t += a[i] * b; a[i] = t % 10; t /= 10; } while (t) { a[++al] = t % 10; t /= 10; } } /** * 功能:n的阶乘中包含的质因子p的个数 * @param n n的阶乘 * @param p 质因子p * @return 有多少个 */ int get(int n, int p) { int cnt = 0; while (n) { // p^1的个数,p^2的个数,p^3的个数... cnt += n / p; n /= p; } return cnt; } // C(a,b)的结果,高精度保存到c数组,同时,返回c数组的长度len void C(int a, int b, int c[], int &cl) { // 乘法的基数是1 c[1] = 1, cl = 1; // 由于高精度数组中只有一位,是1,所以长度也是1 for (int i = 0; i < cnt; i++) { // 枚举区间内所有质数 int p = primes[i]; /* C(a,b)=a!/(b! * (a-b)!) a!中有多少个质数因子p 减去(a-b)!的多少个质数因子p, 再减去b!的质数因子p的个数,就是总个数 s记录了p这个质数因子出现的次数 */ int s = get(a, p) - get(b, p) - get(a - b, p); while (s--) mul(c, cl, p); // 不断的乘p,结果保存到数组c中。len将带回c的有效长度 } } int a, b; int c[N], cl; int main() { cin >> a >> b; // 筛质数 get_primes(N - 1); // 计算组合数 C(a, b, c, cl); // 输出高精度结果 for (int i = cl; i >= 1; i--) printf("%d", c[i]); return 0; } ``` #### 8、这段代码如何理解? ``` s = get(a, p) - get(a - b, p) - get(b, p); ``` 看公式一: $C_a^b=\frac{a \times (a-1) \times ... \times(a-b+1)}{b\times (b-1) \times ... \times 1}$ $=\frac{a \times (a-1) \times ... \times (a-b+1) \times (a-b)!}{ b \times (b-1) \times ... \times 1 \times (a-b)!}$ $=\frac{a!}{b! \times (a-b)!}$ $a!$中质数$p$的个数,减去$b!$中质数$p$的个数,再减去$(a-b)!$中质数$p$的个数,就是公因子消掉的意思。