Lucas 与拓展 Lucas 求组合数

前往 算竞题解 分类 || 前往 算竞模板 分类

简介

$Lucas$ 定理用于处理带模数排列数的计算，$O(n)$ 预处理，$O(\log{n})$ 询问。

$Lucas$ 定理要求模数 $p$ 为质数，而扩展 $Lucas$ 算法不要求，后者时间复杂度略高于前者，为 $O(\log{p}\log{m})$。

模板题链接

下方模板自取，相关证明见下方链接相关题解。

前置：各种求逆元、逆元求组合数、中国剩余定理、$Wilson$ 定理。

洛谷 P3807 【模板】Lucas 定理

洛谷 P4720 【模板】exLucas

Lucas 定理

// Lucas 定理求排列组合
// 时间: O(log_MOD(n)) + 预处理: O(MOD)
// 要求 MOD 是质数且较小, MAX 大于 MOD
const int MAX = 1000005;
const int MOD = 999983;

// 预处理逆元
int inv[MAX];
void load_inv(int n) {
    inv[1] = 1;
    for (int i = 2; i <= n; i++)
        inv[i] = (MOD - MOD / i) * inv[MOD % i] % MOD;
}

// 预处理阶乘
int fac[MAX];
void load_fac(int n) {
    fac[0] = 1;
    for (int i = 1; i <= n; i++)
        fac[i] = fac[i - 1] * i % MOD;
}

// 预处理阶乘的逆元
int ifc[MAX];
void load_ifc(int n) {
    ifc[0] = 1;
    for (int i = 1; i <= n; i++)
        ifc[i] = ifc[i - 1] * inv[i] % MOD;
}

// 初始化绑定
void init(int n) {
    load_inv(n);
    load_fac(n);
    load_ifc(n);  // 必须放在最后
}

// 询问排列数与组合数
inline int A(int m, int n) { return 0 <= n and n <= m ? fac[m] * ifc[m - n] % MOD : 0; }
inline int C(int m, int n) { return 0 <= n and n <= m ? fac[m] * ifc[n] % MOD * ifc[m - n] % MOD : 0; }

// Lucas 定理. 时间: O(log_MOD(n)) + 预处理: O(MOD)
// 要求模是质数
int Lucas(int m, int n) {
    if (n < 0 or n > m) 
        return 0;
    int res = 1;
    while (n) {
        res = res * C(m % MOD, n % MOD) % MOD;
        m /= MOD;
        n /= MOD;
    }
    return res;
}

扩展 Lucas 算法

// exLucas 求排列组合
// 时间: O(log(MOD)log(m)) + 预处理: O(MOD)
// 要求 MOD 较小, MAX 大于 MOD
const int MAX = 1000005;
const int MOD = 999983;

// 分解质因数(无筛优化)
pair<int, int> divisors[128];
int prime_divisors(int x) {
    int cnt = 0;
    for (int i = 2; i <= x / i; i++)
        if (x % i == 0) {
            divisors[cnt] = {i, 0};
            while (x % i == 0) {
                x /= i;
                divisors[cnt].second++;
            }
            cnt++;
        }
    if (x > 1)
        divisors[cnt++] = {x, 1};
    return cnt;
}

// 预处理阶乘 及 分解质因数
int fac[MAX];
int cnt_divisors;
void init_exLucas() {
    fac[0] = 1;
    for (int i = 1; i <= MOD; i++)
        fac[i] = fac[i - 1] * i % MOD;
    cnt_divisors = prime_divisors(MOD);
}

// 快速幂
int quick_pow(int base, int exponent) {
    __int128_t res = 1;
    while (exponent) {
        if (exponent & 1)
            res *= base;
        base *= base;
        exponent >>= 1;
    }
    return res;
}
int quick_pow(int base, int exponent, int modulo) {
    int res = 1 % modulo;
    base %= modulo;
    while (exponent) {
        if (exponent & 1)
            res = res * base % modulo;
        base = base * base % modulo;
        exponent >>= 1;
    }
    return res;
}

// 扩展欧几里得
int _exgcd(int a, int b, int &x, int &y) {
    int d = a;
    if (b != 0) {
        d = _exgcd(b, a % b, y, x);
        y -= (a / b) * x;
    } else {
        x = 1;
        y = 0;
    }
    return d;
}
inline bool exgcd(int a, int b, int c, int &x, int &y) {
    int _x, _y, g;
    g = _exgcd(a, b, _x, _y);
    if (c % g) return false;
    int p = b / g;
    x = (c / g * _x % p + p) % p;
    y = (c - a * x) / b;
    return true;
}

// 线性同余方程
inline int linear_congruence_theorem(int a, int c, int m) {
    int x, y;
    bool flag = exgcd(a, m, c, x, y);
    return flag ? x : -1;
}

// 扩欧求逆元
inline int inv(int primal, int modulo) { return linear_congruence_theorem(primal, 1, modulo); }

// 中国剩余定理(改)
int CRT(int rest[], int modulo[], int n) {
    int res = 0, M = MOD;  // 本例中 M 直接取 MOD 即可
    for (int i = 0; i < n; i++) {
        int mi = modulo[i], Mi = M / mi;
        res = (res + rest[i] * Mi * inv(Mi % mi, mi)) % M;
    }
    return res;
}

// 扩展 Lucas 定理. 时间: O(log(MOD)log(m)) + 预处理: O(MOD)
int exLucas(int m, int n) {
    if (n < 0 or n > m)
        return 0;
    static int rest[128], modulo[128];
    for (int i = 0; i < cnt_divisors; i++) {
        int p = divisors[i].first, M = quick_pow(p, divisors[i].second);
        int cnt = 0, R = 1;
        for (int i = m; i; i /= p) {
            cnt += i / p;
            R = R * fac[i % p] % M;
        }
        for (int i = n; i; i /= p) {
            cnt -= i / p;
            R = R * inv(fac[i % p], M) % M;
        }
        for (int i = m - n; i; i /= p) {
            cnt -= i / p;
            R = R * inv(fac[i % p], M) % M;
        }
        R = R * quick_pow(fac[p], cnt, M) % M;
        rest[i] = R;
        modulo[i] = M;
    }
    return CRT(rest, modulo, cnt_divisors);  // 时间: O(log(MOD))
}