扩展卢卡斯

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

根据唯一分解定理

[p = p_1^{a_1}p_2^{a_2}dots p_k^{a_k} ]

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

[left{egin{array}{ll}C_n^m mod p_1^{a_1}\C_n^m mod p_2^{a_2}\dots\end{array} ight. ]

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

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

此时问题等价于

[frac{n!}{m!(n-m)!}mod p^k ]

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

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

转换式子

[frac{frac{n!}{p^x}}{frac{m!}{p^y}frac{(n-m)!}{p^z}}p^{x - y -z}mod p^k ]

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

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

此时问题等价于

[frac{n!}{p^x}mod p^k ]

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

那么原式(= p^{lfloorfrac{n}{p} floor}(1⋅2⋅3dots)(1⋅2⋅dots) \ = p^{lfloorfrac{n}{p} floor}({lfloorfrac{n}{p} floor})!prod_{i = 1,i otequiv0(mod p)}^ni)

而(mod p)是由循环节的

[= p^{lfloorfrac{n}{p} floor}({lfloorfrac{n}{p} floor})!(prod_{i = 1,i otequiv0(mod p)}^{p^k}i)^{frac{n}{p^k}}(prod_{i = p^klfloorfrac{n}{p^k} floor,i otequiv0(mod p)}^n i) ]

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

((prod_{i = p^klfloorfrac{n}{p^k} floor,i otequiv0(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) = f(lfloorfrac{n}{p} floor)(prod_{i = 1,i otequiv0(mod p)}^{p^k}i)^{lfloorfrac{n}{p^k} floor}(prod_{i = p^klfloorfrac{n}{p^k} floor,i otequiv0(mod p)}^n i)\f(0) = 1 ]

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

那么原式转换为

[frac{f(n)}{f(m)f(n-m)}p^{x-y-z} mod p^k ]

(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))满足递推式

[g(n) = lfloorfrac{n}{p} floor + g(lfloorfrac{n}{p} floor)\g(n) = 0(n < p) ]

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

答案就是

[frac{f(n)}{f(m)f(n-m)}p^{g(x) - g(y) - g(z)} mod p^k ]

模板

时间复杂度(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
", exLucas(n, m, p));
    return 0;
}

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