exlucas

exlucas 用来求如下式子
$$C_n^m\mathrm{mod}p$$
其中 $p$ 不一定为质数。

---

lucas 定理用于求解 $p$ 为质数的情况。对于其他情况，容易想到对 $p$ 进行质因数分解
$$p=\prod _i{p}_i^{k_i}$$
设 $C_n^m\mathrm{mod}p=ans$，由中国剩余定理可得模线性方程组：
$$\{\begin{array}{l}ans\equiv a_1(\mathrm{mod}{p}_1^{k_1})\\ ans\equiv a_2(\mathrm{mod}{p}_2^{k_2})\\ \dots \end{array}$$
若能快速求出 $a$，就能解出 $ans$，得到答案。

现在问题转化为求
$$C_n^m\mathrm{mod}{p}^k$$
其中 $p$ 为质数。

---

运用组合数公式可得上式为
$$\frac{n!}{m!(n-m)!}\mathrm{mod}{p}^k$$
考虑到分子分母都可能含有 $p$ 这个质因子，所以不能直接用逆元求。

所以我们同时提出 $n!,m!,(n-m)!$ 所有的质因子 $p$，可得
$$\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\frac{(n-m)!}{p^z}}\cdot p^{x-y-z}\mathrm{mod}{p}^k$$

现在分子分母都与 $p$ 互质，可以用逆元了。但关键是要求 $\frac{n!}{p^x}\mathrm{mod}{p}^k$。

先考虑求 $n!\mathrm{mod}{p}^k$ 吧，好求一点。

这里我们提出 $1$ 至 $n$ 中 $p$ 的倍数，容易得到它们的乘积为 $p^{\lfloor \frac{n}{p}\rfloor }\cdot \lfloor \frac{n}{p}\rfloor !$，因为它们都能提出一个 $p$，然后剩下 $1$ 到 $\frac{n}{p}$，乘在一起就是 $\lfloor \frac{n}{p}\rfloor !$。

即
$$n!\equiv p^{\lfloor \frac{n}{p}\rfloor }\cdot \lfloor \frac{n}{p}\rfloor !\cdot \left(\prod _{i=1,i\mathrm{mod}p\ne 0}i\right)(\mathrm{mod}{p}^k)$$

众所周知，$x\equiv x+p^k(\mathrm{mod}{p}^k)$，所以上式后半部分的连乘一定存在若干长度为 $p^k$ 的循环，循环中的数的乘积相等。

而 $n$ 里面有几个循环呢？为 $\lfloor \frac{n}{p^k}\rfloor$。

还有 $n$ 的末尾不一定是一整个循环，所以要暴力乘上它们。

令函数 $g(n)=n!$，则有
$$g(n)\equiv p^{\lfloor \frac{n}{p}\rfloor }\cdot g\left(\lfloor \frac{n}{p}\rfloor \right)\cdot {\left(\prod _{i=1,i\mathrm{mod}p\ne 0}^{p^k-1}i\right)}^{\lfloor \frac{n}{p^k}\rfloor }\cdot \left(\prod _{i=1,i\mathrm{mod}p\ne 0}^{n\mathrm{mod}{p}^k}i\right)(\mathrm{mod}{p}^k)$$
在代码实现可使用递归的方式求得 $g(n)$，注意边界 $g(0)=1$。

返回来看 $\frac{n!}{p^x}$，发现上面的式子中只有 $p^{\lfloor \frac{n}{p}\rfloor }$ 才是 $p$ 的倍数，而 ${\left(\prod _{i=1,i\mathrm{mod}p\ne 0}^{p^k-1}i\right)}^{\lfloor \frac{n}{p^k}\rfloor }\cdot \left(\prod _{i=1,i\mathrm{mod}p\ne 0}^{n\mathrm{mod}{p}^k}i\right)$ 都不是 $p$ 的倍数，$f\left(\lfloor \frac{n}{p}\rfloor \right)$ 递归下去也是类似的情况，所以 $p^x$ 全被约掉了。这样，令 $f(n)=\frac{n!}{p^x}$，我们就可以得到一个优美的式子：
$$f(n)\equiv f\left(\lfloor \frac{n}{p}\rfloor \right)\cdot {\left(\prod _{i=1,i\mathrm{mod}p\ne 0}^{p^k-1}i\right)}^{\lfloor \frac{n}{p^k}\rfloor }\cdot \left(\prod _{i=1,i\mathrm{mod}p\ne 0}^{n\mathrm{mod}{p}^k}i\right)(\mathrm{mod}{p}^k)$$
求 $f(n)$ 的代码实现如下：

ll f(ll n,ll p,ll mod)
{
    if(!n) return 1;
    ll ans=1;
    for(int i=1;i<mod;i++) if(i%p) ans=ans*i%mod;
    ans=ksm(ans,n/mod,mod);
    for(int i=1;i<=n%mod;i++) if(i%p) ans=ans*i%mod;
    return ans*f(n/p,p,mod)%mod;
}

---

求得 $f(n)$ 后，$C_n^m\equiv \frac{f(n)}{f(m)f(n-m)}\cdot p^{x-y-z}(\mathrm{mod}{p}^k)$，如何求 $p^{x-y-z}$ 呢？

其实和上面求 $f(n)$ 的思路类似，求 $f(n)$ 时在 $g(n)$ 的基础上把 $p^{\lfloor \frac{n}{p}\rfloor }$ 约掉了，而 $p^{\lfloor \frac{n}{p}\rfloor }$ 正是想要求的 $p^x$ 的一部分。在递归的过程中，$p^x$ 就等于 $p^{\lfloor \frac{n}{p}\rfloor }+p^{\lfloor \frac{\lfloor \frac{n}{p}\rfloor }{p}\rfloor }+\dots$。

下方为求 $C_n^m\mathrm{mod}{p}^k$ 的代码，结合代码方便理解。

ll c(ll n,ll m,ll p,ll mod)
{
    if(n<m) return 0;
    ll f1=f(n,p,mod),f2=f(m,p,mod),f3=f(n-m,p,mod),x=0;
    for(ll i=n;i;i/=p) x+=i/p;
    for(ll i=m;i;i/=p) x-=i/p;
    for(ll i=n-m;i;i/=p) x-=i/p;
    return f1*inv(f2,mod)%mod*inv(f3,mod)%mod*ksm(p,x,mod)%mod;
}

---

至此，我们已经可以解决原问题了。

题目：P4720 【模板】扩展卢卡斯定理/exLucas

完整代码如下

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
ll n,m,mod;
void exgcd(ll a,ll b,ll &d,ll &x,ll &y)
{
    if(!b) d=a,x=1,y=0;
    else exgcd(b,a%b,d,y,x),y-=a/b*x;
}
ll inv(ll n,ll mod)
{
    ll d,x,y;
    exgcd(n,mod,d,x,y);
    return (x%mod+mod)%mod;
}
ll ksm(ll a,ll b,ll mod)
{
    ll ans=1;
    while(b){
        if(b&1) ans=ans*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return ans;
}
ll f(ll n,ll p,ll mod)
{
    if(!n) return 1;
    ll ans=1;
    for(int i=1;i<mod;i++) if(i%p) ans=ans*i%mod;
    ans=ksm(ans,n/mod,mod);
    for(int i=1;i<=n%mod;i++) if(i%p) ans=ans*i%mod;
    return ans*f(n/p,p,mod)%mod;
}
ll c(ll n,ll m,ll p,ll mod)
{
    if(n<m) return 0;
    ll f1=f(n,p,mod),f2=f(m,p,mod),f3=f(n-m,p,mod),x=0;
    for(ll i=n;i;i/=p) x+=i/p;
    for(ll i=m;i;i/=p) x-=i/p;
    for(ll i=n-m;i;i/=p) x-=i/p;
    return f1*inv(f2,mod)%mod*inv(f3,mod)%mod*ksm(p,x,mod)%mod;
}
ll C(ll n,ll m,ll mod)
{
    ll czn=mod,ans=0;
    for(ll i=2;i*i<=czn;i++){
        if(czn%i==0){
            ll x=1;
            while(czn%i==0) czn/=i,x*=i;
            ans=(ans+c(n,m,i,x)*(mod/x)%mod*inv(mod/x,x))%mod;
        }
    }
    if(czn>1) ans=(ans+c(n,m,czn,czn)*(mod/czn)%mod*inv(mod/czn,czn))%mod;
    return ans;
}
int main()
{
    scanf("%lld%lld%lld",&n,&m,&mod);
    printf("%lld",C(n,m,mod));
}