【數學】求組合數的4種方法

組合數公式：（圖來自百度百科）

---

1.迭代法（預處理）求組合數

適用於\(C_a^b\)中\(a\) 和\(b\)不是很大的情況,一般\(1 \leq a,b \leq 10^4\)
所以可以直接預處理出來\(C_a^b\)，用的時候直接查表即可。

#include <iostream>

using namespace std;

const int N = 2010, mod = 1e9 + 7;

int c[N][N];

void init()
{
    for(int i = 0; i < N; i ++ )
        for(int j = 0; j <= i; j ++ )
            if(!j) c[i][j] = 1;
            else c[i][j] = (c[i - 1][j - 1] + c[i - 1][j]) % mod;
}

int main()
{
    init();
    int n;
    scanf("%d", &n);
    while(n --)
    {
        int a,b;
        scanf("%d%d",&a,&b);
        printf("%d\n", c[a][b]);
    }
    return 0;
}

---

2.利用乘法逆元求組合數

\(C_n^m = \frac{n!}{(n-m)m!}\)，此時\(1\leq m,n \leq10^5\)

對乘法逆元不熟悉的可以看這里
將組合數公式轉化為除法形式：n! 表示為fact[n],(n-m)!表示為infact[n - m]，m!表示為infact[m]
所以組合數公式可以寫成:\(C_n^m\) = fact[n] * infact[m] * infact[n - m]
infact表示逆元數組，模數是質數的情況下可以用費馬小定理+快速冪求逆元

費馬小定理兩邊同除\(a\)得：
所以\(a\)模\(p\)的逆元就是\(a^{p-2}\)，這個可以用快速冪求。

#include <iostream>

using namespace std;

typedef long long LL;

const int N = 100010, mod = 1e9 + 7;
int fact[N],infact[N];

LL qmi(LL a,int k,int p)
{
    LL res = 1;
    while(k)
    {
        if(k & 1) res = res * a % p;
        k >>= 1;
        a = a * a % p;
    }
    return res;
}

int main()
{
    //預處理階乘和階乘的逆元
    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] * qmi(i, mod - 2, mod) % mod; //這里是關鍵，把組合數的公式轉化為乘法形式
    }

    int n;
    scanf("%d",&n);
    while(n --)
    {
        int a,b;
        scanf("%d%d",&a,&b);
        printf("%lld\n", (LL)fact[a] * infact[a - b] % mod * infact[b] % mod);
        //因為3個1e9級別的數相乘會爆longlong，所以乘了兩個后要mod一下1e9+7，不影響結果
    }
    return 0;
}

---

3.Lucas定理求組合數

此時\(1\leq a,b \leq 10^{18}\)
Lucas定理：\(C_a^b \equiv C_{a \% p}^{b \% p} \cdot C_{a / p}^{b / p} (mod \ p)\)
中間組合數按定義算即可：\(C_a^b = \frac{a!}{(a-b)!\ b!}\)

#include <iostream>

using namespace std;

typedef long long LL;

LL qmi(LL a,int k, int p)
{
    LL res = 1;
    while(k)
    {
        if(k & 1) res = res * a % p;
        k >>= 1;
        a = a * a % p;
    }
    return res;
}

int C(int a,int b,int p)
{
    if(a < b) return 0;

    LL res = 1;
    for(int i = 1, j = a; i <= b; i ++ , j -- )
    {
        res = res * j % p;
        res = res * qmi(i, p - 2, p) % p;
    }
    return res;
}

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

int main()
{
    int n;
    cin >> n;
    while(n --)
    {
        LL a,b;
        int p;
        cin >> a >> b >> p;
        cout << lucas(a, b, p) << endl;
    }
    return 0;
}

---

4.高精度求組合數

這類題與其他三類題不同的地方在於，題目沒有讓求\(C_a^b\)模某個數，而是直接讓你求出\(C_a^b\)的值
\(1 \leq a,b \leq 5000\)，這就需要用到高精度了。

求組合數的話還是從定義出發，\(C_a^b = \frac{a!}{(a-b)!\ b!}\)
由算術基本定理，
同理，組合數公式中分子和分母都是用階乘表示，階乘當然也可以分解為有限個質數的乘積
這樣我們可以把分子分母都分解質因數，然后分子和分母上的質因數還可以消掉，顯然分子是一定大於分母的(組合數算下來總不能小於1吧，當然a<b的情況我們會特判掉)，所以分母中的質因子一定會被消掉，這樣這個組合數公式就變成一個由若干個質因子乘積表示的式子了，所以我們只需要寫一個高精度乘法就可以了。

分解質因數：
這里要求\(n!\)中質因子\(p\)出現的次數，公式如下：
質因子p的次數為：\(\lfloor \frac{n}{p} \rfloor +\lfloor \frac{n}{p^2} \rfloor + \ldots+\lfloor \frac{n}{p^n} \rfloor\)
上式中，第一項的含義為\(1\sim n\)中\(p\)的倍數的個數，第二項為\(1\sim n\)中\(p^2\)的倍數，依次類推...
當然，這個式子中並不一定n項都存在，大部分情況是只存在個別幾項。

#include <iostream>
#include <vector>

using namespace std;

const int N = 5010;

int primes[N],cnt;
bool st[N];
int sum[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;
        }
    }
}

int get(int n, int p)//分解質因數，求n的階乘中有多少個質因子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(); i ++ )
    {
        t += a[i] * b;
        c.push_back(t % 10);
        t /= 10;
    }
    while(t)
    {
        c.push_back(t % 10);
        t /= 10;
    }
    return c;
}

int main()
{
    int a,b;
    cin >> a >> b;
    get_primes(a);

    for(int i = 0; i < cnt; i ++ )
    {
        int p = primes[i];
        sum[i] = get(a, p) - get(a - b, p) - get(b, p);
    }

    vector<int> res;
    res.push_back(1);

    for(int i = 0; i < cnt; i ++ )
        for(int j = 0; j < sum[i]; j ++ )
            res = mul(res, primes[i]);


    for(int i = res.size() - 1; i >= 0; i -- ) printf("%d", res[i]);
    puts("");

    return 0;
}

---

階乘分解

另外再補充一道階乘分解，也需要用到分解質因數

#include <iostream>
#include <algorithm>
#include <cstring>
#include <cstdio>

using namespace std;

const int N = 1e6 + 10;

int primes[N], cnt;
bool st[N];

void init(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;
        }
    }
}

int main()
{
    int n;
    scanf("%d", &n);
    init(n);

    for(int i = 0; i < cnt; i ++ )
    {
        int p = primes[i];
        int s = 0;
        for(int j = n; j; j /= p) s += j / p; //s就是每個質數的指數冪
        printf("%d %d\n", p, s);
    }

    return 0;
}