[数学系列 #3] 扩展Lucas

问题:求$\binom{n}{m} \bmod p$,$p\leq 100000$,$p$不一定为质数。

考虑到唯一分解:$p=\prod {p_i}^{k_i}$

则问题可以分解成同余方程组:$x\equiv \binom{n}{m} \pmod{{p_i}^{k_i}}$ ,所求即为 $x$ 。

现仅需考虑如何求得$\frac{n!}{m!(n-m)!} \bmod {p_i}^{k_i}$

进一步考虑如何求得$n! \bmod {p_i}^{k_i}$

我们考虑将$n!$中的$p_i$因子剔除考虑。

则:$22! = (1\times 2\times 4\times 5\times 7\times 8)\times (10\times 11\times 13\times 14\times 16\times 17)\times (19\times 20\times 22)\times 3^6\times (1\times 2\times 3\times 4\times 5\times 6\times 7)$

发现前两组的取模结果相同,第三组冗余可以暴力算,最后一组即为$\left \lfloor \frac{22}{3} \right \rfloor!$

其中所有的涉及逆元计算皆为互质所以一定有解,exgcd即可。

最后再考虑$n!$中的$p_i$因子个数$f(n)$,则有递推式:$$f(n)=f(\left \lfloor \frac{n}{p_i} \right \rfloor)+\left \lfloor \frac{n}{p_i} \right \rfloor$$

同时需要注意建议全程使用long long计算。虽然尝试非常小心地使用了int但还是不知道哪里溢出了。

Code

// Code by ajcxsu
// Problem: exlucas

#include<bits/stdc++.h>
using namespace std;

typedef int _int;
#define int long long

void exgcd(int a, int b, int &x, int &y) {
    if(!b) x=1, y=0;
    else exgcd(b, a%b, y, x), y-=a/b*x;
}
int inv(int a, int mo) {
    int x, y; exgcd(a, mo, x, y);
    return (x%mo+mo)%mo;
}
int qpow(int x, int y, int mo) {
    int ret=1;
    while(y) {
        if(y&1) ret=1ll*x*ret%mo;
        x=1ll*x*x%mo, y>>=1;
    }
    return ret;
}

int f(int x, int p, int k) {
    if(x<=1) return 1;
    int ret=1;
    for(int i=2;i<=p;i++) if(i%k) ret=1ll*ret*i%p;
    ret=qpow(ret, x/p, p);
    for(int i=2;i<=x%p;i++) if(i%k) ret=1ll*ret*i%p; // <= 注意
    return 1ll*ret*f(x/k, p, k)%p;
}

int x, y, p, c, rp, ans;
int g(int x, int p) { int ret=0; for(int i=x/p;i;i/=p) ret+=i; return ret; }
void work(int ki) {
    int na, tc=0;
    int m=1, M, t;
    while(p%ki==0) m*=ki, p/=ki;
    M=rp/m, t=inv(M, m);
    na=1ll*f(x, m, ki)*inv(f(y, m, ki), m)%m*inv(f(x-y, m, ki), m)%m;
    tc=g(x, ki)-g(y, ki)-g(x-y, ki);
    na=1ll*na*qpow(ki, tc, m)%m, ans=(ans+1ll*na*M%rp*t%rp)%rp;
}

_int main() {
    ios::sync_with_stdio(false), cin.tie(0);
    cin>>x>>y>>p, rp=p;
    for(int i=2;i*i<=p;i++) if(p%i==0) work(i);
    if(p>1) work(p);
    cout<<ans<<endl;
    return 0;
}

本文链接:https://acxblog.site/archives/exlucas.html
文章采用知识共享署名-非商业性使用 4.0 国际许可协议进行许可。