卢卡斯Lucas & 扩展卢卡斯exLucas

jiazhaopeng · 2020-02-05 15:14:29 · 个人记录

卢卡斯定理

适用条件：

只算数量较少的C(n, m), 模数P较小（约1e5），且P为质数。

递推式：C(n, m) = C(n/P, m/P) * C(n%P, m%P) % P;

当 n < m 的时候， C_n^m = 0！

当 n == 0 或者 m == 0 的时候就可以停止了。

Code:

//阶乘已经预处理好，逆元用快速幂
long long Lucas(long long n, long long m) {//n > m 
    if (m == 0 || n == 0)   return 1;
    long long res = Lucas(n / p, m / p) * C(n % p, m % p) % p;
    return res;
}

扩展卢卡斯

适用范围：

只求少量组合数，m, n巨大，但模数 M = \prod p^q ，且 p^q 不大。

复杂度：O(\log_Pn * (p^q))

扩展卢卡斯实际上还可以搞任何阶乘除阶乘的类似问题。比如：P2183 [国家集训队]礼物。主要思路是对每个 p_q 单独考虑，对因子 p 和其它因子单独考虑。提出来 p 的倍数，剩下的利用暴力前缀积搞，p 的倍数部分递归搞。

Code:

//P2183礼物
int realP;
int n, m, w[10];
int p[44], c[44], ptot, rp[44];
inline void Div(int x) {
    for (int i = 2; i * i <= x; ++i) {
        if (x % i == 0) {
            p[++ptot] = i; rp[ptot] = 1;
            while (x % i == 0)  x /= i, ++c[ptot], rp[ptot] *= i;
        }
    }
    if (x > 1)  p[++ptot] = x, c[ptot] = 1, rp[ptot] = x;
}

inline int quickpow(int x, int k, int P)

int exgcd(int a, int b, int &x, int &y)

int calc(int a, int b) {//ax + by = 1 (gcd(a,b)=1), return x
    int x, y; exgcd(a, b, x, y);
    return (x % b + b) % b;
}

int nw;
int yu, ct;
int Prod[101000];
int fakeans[44];
void sol(int n) {
    if (!n) return ;
    yu = yu * quickpow(Prod[rp[nw]], n / rp[nw], rp[nw]) % rp[nw] * Prod[n % rp[nw]] % rp[nw];
    ct += n / p[nw];
    sol(n / p[nw]);
}

inline void merge() {
    int ans = 0;
    for (int i = 1; i <= ptot; ++i) {
        int tmp = calc(realP / rp[i], rp[i]);
        ans = (ans + fakeans[i] * (realP / rp[i]) % realP * tmp) % realP;
    }
    printf("%lld\n", (ans % realP + realP) % realP);
}

signed main() {//calculate (n!)/(w_1! w_2!...)
    Div(realP);
    for (int i = 1; i <= ptot; ++i) {
        nw = i;
        yu = 1, ct = 0;
        Prod[0] = 1;
        for (int j = 1; j <= rp[i]; ++j)
            if (j % p[i])   Prod[j] = Prod[j - 1] * j % rp[i];
            else    Prod[j] = Prod[j - 1];
        sol(n);
        int memo = yu, memoct = ct; yu = 1; ct = 0;
        for (int j = 1; j <= m; ++j)    sol(w[j]);
        yu = calc(yu, rp[nw]) * memo % rp[i];
        fakeans[i] = yu * quickpow(p[i], memoct - ct, rp[i]) % rp[i];
        yu = 1; ct = 0;
    }
    merge();
}