[数学系列 #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://pst.iorinn.moe/archives/exlucas.html
许可: https://pst.iorinn.moe/license.html若无特别说明,博客内的文章默认将采用 CC BY 4.0 许可协议 进行许可☆