相关推荐recommended
求组合数的三种算法
作者:mmseoamin日期:2024-04-30

组合数

  • 一、预处理组合数
  • 二、预处理阶乘
  • 三、卢卡斯定理

    一、预处理组合数

    核心:

    C a b = C a − 1 b + C a − 1 b − 1 C_a^b = C_{a-1}^b + C_{a-1}^{b-1} Cab​=Ca−1b​+Ca−1b−1​

    适用范围: a a a 较小的情况下,如 a ≤ 1 0 3 a \leq 10^3 a≤103。

    算法简析:令 C[n][k] = C n k \text{C[n][k]}=C_n^k C[n][k]=Cnk​,规定 C[0][0] = 1 \text{C[0][0] = 1} C[0][0] = 1,则

    C[n][k] = { 1 , k = = 0 C[n - 1][k] + C[n - 1][k - 1] , 0 < k ≤ n   a n d   n ≥ 1 \begin{split} \text{C[n][k]}=\begin{cases} 1&,k==0\ \text{C[n - 1][k] + C[n - 1][k - 1]}&,01C[n - 1][k] + C[n - 1][k - 1]​,k==0,0

    #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; 
    }
    

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


    二、预处理阶乘

    核心:

    C a b = a ! b ! ( a − b ) ! = a ! ∗ ( b ! ) − 1 ∗ ( ( a − b ) ! ) − 1 C_a^b=\frac{a!}{b!(a-b)!}=a!\ast (b!)^{-1}\ast ((a-b)!)^{-1} Cab​=b!(a−b)!a!​=a!∗(b!)−1∗((a−b)!)−1

    适用范围: a a a 较大的情况下,如 a ≤ 1 0 5 a\leq 10^{5} a≤105。

    算法简析:

    先来看逆元和费马小定理的定义:

    • 1、逆元:对于模 m m m,有两个数 a a a 和 b b b,若 a b ≡ 1 ( mod  m ) ab\equiv1(\text{mod}~m) ab≡1(mod m),即 a a a 和 b b b 的乘积除以 m m m 的余数为1,则 a a a 和 b b b 互为模 m m m 意义下的乘法逆元。记 b = a − 1 ,   a = b − 1 b=a^{-1},~a=b^{-1} b=a−1, a=b−1。
    • 2、费马小定理:若 p p p 是一个素数,且 a a a 不被 p p p 整除,则 a p − 1 ≡ 1 ( mod  p ) a^{p-1}\equiv1(\text{mod}~p) ap−1≡1(mod p),即 a p − 1 a^{p-1} ap−1 除以 p p p 的余数为1。

      由模运算的性质, ≡ \equiv ≡ 两边同乘 a a a 的逆元 a − 1 a^{-1} a−1,得 a p − 2 ≡ a − 1 ( mod  p ) a^{p-2}\equiv a^{-1}(\text{mod}~p) ap−2≡a−1(mod p)。在模 p p p 的意义下, a p − 2 a^{p-2} ap−2 和 a − 1 a^{-1} a−1 等价。求 a − 1 a^{-1} a−1,就转换为求 a p − 2 a^{p-2} ap−2。

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

      • 1、阶乘:由

        a ! = ( a − 1 ) ! ∗ a a!=(a-1)!\ast a a!=(a−1)!∗a

        得,

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

        • 2、阶乘的逆元:由费马小定理,在模 p p p 下,

          ( a ) p − 2 ≡ ( a ) − 1 ( mod  p ) ( a ! ) − 1 = ( ( a − 1 ) ! ) − 1 ∗ a − 1 \begin{split} (a)^{p-2}&\equiv (a)^{-1}(\text{mod}~p) \\ (a!)^{-1}&=((a-1)!)^{-1}\ast a^{-1} \end{split} (a)p−2(a!)−1​≡(a)−1(mod p)=((a−1)!)−1∗a−1​

          得,

          infact[a] = infact[a - 1] ∗ a − 1 a − 1 = pow(a, p - 2) mod p \begin{split} \text{infact[a]}&=\text{infact[a - 1]} \ast a^{-1} \\ a^{-1}&=\text{pow(a, p - 2) mod p} \end{split} infact[a]a−1​=infact[a - 1]∗a−1=pow(a, p - 2) mod p​

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

          #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 = fact[n] * infact[k] * infact[n - k] C_n^k=\text{fact[n] * infact[k] * infact[n - k]} Cnk​=fact[n] * infact[k] * infact[n - k]。注意:计算过程中可能会溢出,要进行模运算。


          三、卢卡斯定理

          核心:卢卡斯定理

          C a b ≡ C a  mod  p b  mod  p   ⋅   C ⌊ a / p ⌋ ⌊ b / p ⌋ ( mod  p ) C_a^b\equiv C_{a~\text{mod}~p}^{b~\text{mod}~p}~·~C_{\lfloor a/p \rfloor}^{\lfloor b/p \rfloor}(\text{mod}~p) Cab​≡Ca mod pb mod p​ ⋅ C⌊a/p⌋⌊b/p⌋​(mod p)

          适用范围: a a a 很大的情况,比如 a ≤ 1 0 18 a \leq 10^{18} a≤1018。

          算法简析:若 a a a 和 b b b 很大,我们可以通过卢卡斯定理缩小 a a a 和 b b b,直至 a ,   b < p a,~b< p a, b

          #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;
          }