扩展卢卡斯

$ C_n^m \mod p\ ,p不一定为质数$

根据唯一分解定理

得到\(k\)个互质的\(p_i^{a_i}\)，满足

最后用中国剩余定理合并求得最小\(x\)解即可

那么问题转换为求\(C_n^m \mod p^k,[p∈prime]\)

此时问题等价于

那么只需要求\(m!(n-m)!\)对\(p^k\)的逆元就可以了

但是可能不存在逆元，因为\(p^k\)不一定为质数，而对于非质数的模采用扩展欧几里得求解\(\frac{1}{a}\mod p\)，但需要满足\(gcd(a,p) = 1\)

转换式子

其中\(x\)为\(n!\)中\(p\)因子的个数，\(y\)为\(m!\)中\(p\)因子的个数，\(z\)为\((n-m)!\)中\(p\)因子的个数

那么如果对于一个n，可以求出\(\frac{n!}{p^x}\)，就可以求出逆元了

此时问题等价于

对\(n! = 1⋅2⋅3...⋅n = (p⋅2p⋅3p\dots)(1⋅2\dots)\)

那么原式\(= p^{\lfloor\frac{n}{p}\rfloor}(1⋅2⋅3\dots)(1⋅2⋅\dots) \\ = p^{\lfloor\frac{n}{p}\rfloor}({\lfloor\frac{n}{p}\rfloor})!\prod_{i = 1,i\not\equiv0(\mod p)}^ni\)

而\(mod\ p\)是由循环节的

其中\((\prod_{i = 1,i\not\equiv0(\mod p)}^{p^k}i)\)是循环节\(1∼p\)中所有无p因子数的乘积。

\((\prod_{i = p^k\lfloor\frac{n}{p^k}\rfloor,i\not\equiv0(\mod p)}^n i)\)是余项的乘积

比如$22! \mod 3^2 $

\(22! \equiv (3⋅6⋅9⋅12⋅15⋅18⋅21)(1⋅2⋅4⋅\dots)(\mod 3^2) \\ \equiv 3^7(1⋅2⋅3⋅4⋅5⋅6⋅7)(1⋅2⋅4⋅5⋅7⋅8)(10⋅11⋅13⋅14⋅16⋅17)(19⋅20⋅22)(\mod 3^2)\\ \equiv 3^7(1⋅2⋅3⋅4⋅5⋅6⋅7)(1⋅2⋅4⋅5⋅7⋅8)^2(19⋅20⋅22)(\mod 3^2)\)

设\(f(n) = \frac{n!}{p^x}\)

那么求\(f(n)\)的时间复杂度是\(O(log_pn)\)

那么原式转换为

\(f(m)，f(n-m)\)一定与\(p^k\)互质，就可以用扩展欧几里得求逆元了

而对于\(p^{x - y - z}\)

比如求\(f(n) = \frac{n!}{p^x}\)里的\(x\)

设\(g(n)\)表示\(n!\)中有多少个p的因子，也就是\(x\)

而\(g(n)\)满足递推式

就可以用时间复杂度为\(O(log_pn )\)求出

答案就是

模板

时间复杂度\(O(plogp)\)

所以\(p ≤1e6\),n和m可以很大

#include <iostream>
#include <cstdio>
#define ll long long
using namespace std;
void ex_gcd(ll a, ll b, ll &d, ll &x, ll &y){
    if(!b){
        d = a, x = 1, y = 0;
        return;
    }
    ex_gcd(b, a % b, d, y, x);
    y -= x * (a / b);
}
ll inv(ll a, ll p){//求a在p模下的逆元
    ll d, x, y;
    ex_gcd(a, p, d, x, y);
    return (x % p + p) % p;//保证有解了
}
ll pow(ll a, ll b, ll p){
    ll ans = 1; a %= p;
    while(b){
        if(b & 1)ans = ans * a % p;
        b >>= 1;
        a = a * a % p;
    }
    return ans;
}
ll crt(ll a, ll M, ll m){//x = a(mod m)
    return inv(M / m, m) * (M / m) * a % M;
}
ll fac(ll n, ll p, ll pk){//f(n)
    if(!n)return 1;
    ll ans = 1;
    for(ll i = 1; i <= pk; i++)//循环节
        if(i % p) ans = ans * i % pk;
    ans = pow(ans, n / pk, pk);
    for(ll i = pk * (n / pk); i <= n; i++)//余项
        if(i % p) ans = ans * (i % pk) % pk;
    return ans * fac(n / p, p, pk) % pk;
}
ll C(ll n, ll m, ll p, ll pk){
    ll N = fac(n, p, pk), M = fac(m, p, pk), Z = fac(n - m, p, pk);
    ll sum = 0;
    for(ll i = n; i; i = i / p) sum += i / p;//g(n)
    for(ll i = m; i; i = i / p) sum -= i / p;//g(m)
    for(ll i = n - m; i; i = i / p) sum -= i / p;//g(n-m)
    return N * pow(p, sum, pk) % pk * inv(M, pk) % pk * inv(Z, pk) % pk;
}
ll exLucas(ll n, ll m, ll p){//C(n,m) % p
    ll ans = 0;
    ll t = p;
    for(ll i = 2; i * i <= p; i++){
        ll k = 1;
        while(t % i == 0) k *= i,t /= i;
        ans = (ans + crt(C(n, m, i, k), p, k)) % p;
    }
    if(t > 1)
        ans = (ans + crt(C(n, m, t, t), p, t)) % p;
    return ans;
}
int main(){
    ll n, m, p;
    scanf("%lld%lld%lld", &n, &m, &p);
    printf("%lld\n", exLucas(n, m, p));
    return 0;
}

对于\(p\)很大的数，可以先将p分解成几个较小的质数乘积，再利用中国剩余定理求解