原博客地址

预处理结果

问题类型

n组查询：$1 \leq n \leq 10^4$
每次查询$C_a^b$ : $1 \leq b \leq a \leq 2000$
如果按照公式计算，最坏时间复杂度是O(n * b)，当n取$10^4$，a 和 b取 $10^5$ 时，总运算次数为 $10^4 * 2 * 10^5 = 2 * 10^9$

算法原理及流程

根据公式 $C_a^b = c_{a - 1}^b + C_{a - 1}^{b - 1}$ 进行递推

代码实现

#include <iostream>

using namespace std;

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

int n;
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; // c_i^0都等于0
            else c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]) % mod;
}
int main()
{
    init();
    cin >> n;
    while (n --)
    {
        int a, b;
        cin >> a >> b;
        cout << c[a][b] << endl;
    }
    return 0;
}

预处理中间值

问题类型

n组查询：$1 \leq n \leq 10^4$
每次查询$C_a^b$ : $1 \leq b \leq a \leq 10^5$
之所以不能预处理结果，是因为a和b的范围过大，用数组下标已经存不下了

算法原理及流程

预处理阶乘和阶乘的逆元。由于根据公式计算涉及到除法，所以取模需要用到逆元

代码实现

解释一下代码中计算阶乘的逆元时的两种写法等价的原因：$(k!)^{-1} = (k!)^{mod - 2} = ((k - 1)! * k)^{mod - 2} = ((k - 1)!)^{mod - 2} * (k)^{mod - 2} = ((k - 1)!)^{-1} * (k)^{-1}$

#include <iostream>

using namespace std;

typedef long long LL;

const int N = 1e5 + 10, mod = 1e9 + 7;

int n;
LL fact[N], infact[N]; // fact[i]:i的阶层 infact[i]:i的阶乘的逆元

int qpow(int a, int b, int mod)
{
    int res = 1;
    while (b)
    {
        if (b & 1) res = (LL)res * a % mod;
        a = (LL)a * a % mod;
        b >>= 1;
    }
    return res;
}
void init()
{
    fact[0] = infact[0] = 1;
    for (int i = 1; i < N; ++ i) 
    {
        fact[i] = (LL)fact[i - 1] * i % mod;
        infact[i] = qpow(fact[i], mod - 2, mod); // 因为mod是一个质数，所以可以使用费马小定理计算逆元
        // 或者 infact[i] = (LL)infact[i - 1] * qpow(i, mod - 2, mod); 原因见上
    }
}
int main()
{
    init();
    cin >> n;
    while (n --)
    {
        int a, b;
        cin >> a >> b;
        cout << ((LL)fact[a] * infact[b] % mod * infact[a - b] % mod) << endl;
        // 两个大int可能会爆long long，所以需要先%一下
    }
    return 0;
}

卢卡斯优化

问题类型

n组查询：$1 \leq n \leq 20$
每次查询$C_a^b$ : $1 \leq b \leq a \leq 10^{18}$
mod p: $1 \leq p \leq 10^5$
查询次数不多；该类型具有一个很明显的特征就是 $C_a^b$中的a和b都是大于模数p的，所以可以采用卢卡斯定理进行优化

算法原理及流程

卢卡斯(lucas)定理：$C_a^b \equiv C_{a / p}^{b / p} * C_{a \% p}^{b \% p} \pmod {p}$

首先调用卢卡斯定理进行优化，然后直接使用公式进行计算，除法取模需要求逆元

代码实现

#include <iostream>

using namespace std;

typedef long long LL;

int n;

int qpow(int a, int b, int p)
{
    int res = 1;
    while (b)
    {
        if (b & 1) res = (LL)res * a % p;
        a = (LL)a * a % p;
        b >>= 1;
    }
    return res;
}
int c(int a, int b, int p) // 调用组合数时候就一定保证a和b是小于p的了，所以a和b定为int
{
    if (b > a) return 0;

    // 按照公式计算
    int res = 1;
    for (int i = a - b + 1, j = 1; j <= b; ++ i, ++ j)
    {
        res = (LL)res * i % p;
        // 保证b是小于p的，因为p质数，所以1到b的所有数都一定与p互质，所以计算逆元可以采用费马小定理
        res = (LL)res * qpow(j, p - 2, p) % p; // 理解原理需要用到“预处理中间值-代码实现”中提到的等式
    }

    return res;
}
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; // %p保证运算后的结果小于p，无法再使用lucas优化了，所以直接调用组合式函数即可，但/运算则无法保证
}
int main()
{
    cin >> n;
    while (n --)
    {
        int p;
        LL a, b;
        cin >> a >> b >> p;

        cout << lucas(a, b, p) << endl;
    }
    return 0;
}

质因数分解

问题类型

计算结果不能取模，结果可能很大

算法原理及流程

原理
很明显的是，结果很大需要用高精度，如果单纯通过公式进行计算，$c_a^b = \frac{(a \ - \ b \ + \ 1) \ * \ … \ * \ a}{1 \ * … \ * \ b}$，需要用到高精度乘法和除法，整起来就有些复杂了，所以考虑将$C_a^b$ 转化为 $ p_1^{k_1} * p_2^{k_2} * \ … \ * p_n^{k_n} $，，从而将运算转变为单纯的乘法。
如果对公式 $c_a^b = \frac{(a \ - \ b \ + \ 1) \ * \ … \ * \ a}{1 \ * … \ * \ b}$ 中的每个数进行一次质因数分解，复杂度就是 $O(n \sqrt[2]{n})$，估计也没什么大问题，但是代码写起来还是有点复杂的，所以这里利用一个很神奇的计算方法：

a! 中 质因子 p 的个数为 $\lfloor \frac{a}{p} \rfloor \ + \ \lfloor \frac{a}{p^2} \rfloor \ + \ … \ + \ \lfloor \frac{a}{p^n} \rfloor$.

不妨举个例子来理解这个方法，计算6！中质因子2的个数: $6! = 1 * 2 * 3 * 4 * 5 * 6$,其中2，4，6是2的倍数，4是$2^2$的倍数，而2，6对质因子2个数的贡献值为1，4的贡献值为2。4中有两个2，在$\lfloor \frac{a}{p} \rfloor$时考虑到的是它其中一个2，$\lfloor \frac{a}{p^2} \rfloor$考虑的则是另一个2。

流程
$C_a^b = \frac{a!}{b! \ * \ (a - b)!}$,我们首先预处理出所有质因子，然后分别计算出所有质因子的个数，最后调用高精度乘法计算得出最后结果

代码实现

#include <iostream>
#include <vector>

using namespace std;

const int N = 5010;

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

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

void get_primes(int n) // 打素数表
{
    st[0] = st[1] = true;
    for (int i = 2; i <= n; ++ i)
    {
        if (!st[i]) primes[cnt ++]  = i;
        for (int j = 0; primes[j] <= n / i; ++ j)
        {
            st[i * primes[j]] = true;
            if (i % primes[j] == 0) break;
        }
    }
}
int get(int a, int p) // 返回a!中质因子p的个数
{
    int res = 0;
    while (a)
    {
        res += a / p;
        a /= p;
    }
    return res;
}

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(b, p) - get(a - 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; -- i) cout << res[i];
    cout << endl;

    return 0;
}