扩展卢卡斯定理学习笔记

扩展卢卡斯定理学习笔记

【模板】扩展卢卡斯

前置知识

1. \(exgcd\)
2. 中国剩余定理
3. 一定数学能力

卢卡斯定理是用来求\(C_n^m \bmod{p}\)的工具,若\(p\)为质数,则\(C_n^m\equiv C_{\lfloor{\frac{n}{p}}\rfloor}^{\lfloor{\frac{m}{p}}\rfloor}*C_{n\%p}^{m\%p}(\bmod \ p)\)

但是当\(p\)不是质数的时候,就需要用到扩展卢卡斯定理了.下面简单说一下过程:

1

将模数\(P\)分解成\(\prod p_i^{k_i}\)的形式,也就是将\(P\)表示成若干个质数的指数的乘积形式.

2

设答案为\(ans\),则有\[ans\equiv C^{m}_{n}(\mod \ p_i ^{k_i})\],也就是说,我们只需要求解\(C_n^m\)在模某个质数的几次幂下的结果,然后将这些结果用中国剩余定理合并就可以求出答案.

3

下面用\(p\)表示\(p_i\),\(k\)表示\(k_i\).
\[C_n^m=\frac{n!}{m!*(n-m)!}\]
将\(n!,m!,(n-m)!\)中含有\(p_i\)的因子提出来,将上面的式子进行变形,则有:
\[C_n^m\equiv \frac{\frac{n!}{p^{a1}}}{\frac{m!}{p^{a2}}*\frac{(n-m)!}{p^{a3}}}*p^{a1-a2-a3}(\mod p^k)\]
那么此时式子的左半边都是与\(p^k\)互质的,可以直接用\(exgcd\)求逆元,右边可以快速幂解决.

4

现在我们需要解决\[\frac{n!}{p^{a1}}\],根据别人得出的式子,我们可以知道:
\[n!\equiv p^{\lfloor \frac{n}{p} \rfloor}*\lfloor \frac{n}{p} \rfloor!*(\prod_{i=1}^{p^k}num_i)^{\lfloor \frac{n}{p^k} \rfloor}*(\prod_{i=1}^{n\%p^k}num_i)(\mod \ p^k)\]
其中\(num_i\)表示不含\(p\)因子的数字.
其中\((\prod_{i=1}^{p^k}num_i)^{\lfloor \frac{n}{p^k} \rfloor}\)可以直接枚举\(1\)到\(p^k\)的数字乘起来,然后次方做一遍快速幂.
\((\prod_{i=1}^{n\%p^k}num_i)\)直接枚举乘起来,\(\lfloor \frac{n}{p} \rfloor!\)递归求解,边界是当\(n=0\)时,返回\(1\).
因为我们要求\(\frac{n!}{p^{a1}}\),所以在乘的时候可以直接把式子第一部分的\(p^{\lfloor \frac{n}{p} \rfloor}\)直接忽略掉.

int fac(int n, int pi, int pk){
    if(!n) return 1;
    int res = 1;
    for(int i = 2; i < pk; i++)
        if(i % pi) (res *= i) %= pk;
    res = qpow(res, n/pk, pk);
    for(int i = 2; i <= n%pk; i++)
        if(i % pi) (res *= i) %= pk;
    return fac(n/pi, pi, pk)*res%pk;
}

5

根据\[C_n^m\equiv \frac{\frac{n!}{p^{a1}}}{\frac{m!}{p^{a2}}*\frac{(n-m)!}{p^{a3}}}*p^{a1-a2-a3}(\mod p^k)\],计算出了左边式子的上下部分,就可以直接用\(exgcd\)求出下面两个结果在模\(p^k\)意义下的逆元,然后乘起来.

inline int C(int n, int m, int pi, int pk){
    if(n < m) return 0;
    int r1 = fac(n, pi, pk), r2 = fac(m, pi, pk);
    int r3 = fac(n-m, pi, pk), cnt = 0;
    for(int i = n; i; i /= pi) cnt += i/pi;
    for(int i = m; i; i /= pi) cnt -= i/pi;
    for(int i = n-m; i; i /= pi) cnt -= i/pi;
    return r1*inv(r2, pk)%pk*inv(r3, pk)%pk*qpow(pi, cnt, pk)%pk;
}

6

再根据\[ans\equiv C^{m}_{n}(\mod \ p ^{k})\],我们可以用中国剩余定理合并这些结果(这里我用的是扩展中国剩余定理)

inline int exlucas(int n, int m, int p){
    int pi, pk, res = 0, x, y, gcd, c, M = 1;
    for(int i = 2; i*i <= p; i++){
        if(p % i == 0){
            pi = i, pk = 1;
            while(p % i == 0) p /= i, pk *= i;
            gcd = exgcd(M, pk, x, y), c = Mod((C(n, m, pi, pk)-res)%pk+pk, pk);
            x = Mod(x*(c/gcd)%pk+pk, pk);
            res += M*x, M *= pk/gcd, res %= M;
        }
    }
    if(p > 1){ // 最后分解质因数后可能存在一个大于sqrt(p)的大质数
        pi = pk = p; gcd = exgcd(M, pk, x, y);
        c = Mod((C(n, m, pi, pk)-res)%pk+pk, pk);
        x = Mod(x*(c/gcd)%pk+pk, pk);
        res += M*x, M *= pk/gcd, res %= M;
    }
    return res;
}

下面放一下完整代码

#include<bits/stdc++.h>
using namespace std;
typedef int _int;
#define int long long

int n, m, p;

int exgcd(int a, int b, int &x, int &y){
    if(!b){ x = 1, y = 0; return a; }
    int gcd = exgcd(b, a%b, x, y), tmp;
    tmp = x, x = y, y = tmp-a/b*y;
    return gcd;
}

inline int Mod(int x, int mod){ return x > mod ? x-mod : x; }

inline int inv(int a, int mod){
    int x, y; exgcd(a, mod, x, y);
    return (x%mod+mod)%mod;
}

inline int qpow(int x, int n, int mod){
    int res = 1;
    for(; n; x = x*x%mod, n >>= 1)
        if(n & 1) (res *= x) %= mod;
    return res;
}

int fac(int n, int pi, int pk){
    if(!n) return 1;
    int res = 1;
    for(int i = 2; i < pk; i++)
        if(i % pi) (res *= i) %= pk;
    res = qpow(res, n/pk, pk);
    for(int i = 2; i <= n%pk; i++)
        if(i % pi) (res *= i) %= pk;
    return fac(n/pi, pi, pk)*res%pk;
}

inline int C(int n, int m, int pi, int pk){
    if(n < m) return 0;
    int r1 = fac(n, pi, pk), r2 = fac(m, pi, pk);
    int r3 = fac(n-m, pi, pk), cnt = 0;
    for(int i = n; i; i /= pi) cnt += i/pi;
    for(int i = m; i; i /= pi) cnt -= i/pi;
    for(int i = n-m; i; i /= pi) cnt -= i/pi;
    return r1*inv(r2, pk)%pk*inv(r3, pk)%pk*qpow(pi, cnt, pk)%pk;
}

inline int exlucas(int n, int m, int p){
    int pi, pk, res = 0, x, y, gcd, c, M = 1;
    for(int i = 2; i*i <= p; i++){
        if(p % i == 0){
            pi = i, pk = 1;
            while(p % i == 0) p /= i, pk *= i;
            gcd = exgcd(M, pk, x, y), c = Mod((C(n, m, pi, pk)-res)%pk+pk, pk);
            x = Mod(x*(c/gcd)%pk+pk, pk);
            res += M*x, M *= pk/gcd, res %= M;
        }
    }
    if(p > 1){
        pi = pk = p; gcd = exgcd(M, pk, x, y);
        c = Mod((C(n, m, pi, pk)-res)%pk+pk, pk);
        x = Mod(x*(c/gcd)%pk+pk, pk);
        res += M*x, M *= pk/gcd, res %= M;
    }
    return res;
}

_int main(){
    cin >> n >> m >> p;
    cout << exlucas(n, m, p) << endl;
    return 0;
}

转载于:https://www.cnblogs.com/BCOI/p/10368553.html