参考博客
exLucas：
求 \(C_n^m\bmod d\) ( \(d\) 不一定为质数)
1.
将 \(d\) 质因数分解为 \(d=p_1^{c_1}\times p_2^{c_2}\times\cdots\times p_k^{c_k}\)
\(\forall i,j\in[1,k]\) , \(p_i^{c_i}\) 与 \(p_j^{c_j}\) 互质，所以可以构造出如下同余方程：
\[\begin{cases}a_1\equiv C_n^m\ (\bmod \ p_1^{c_1})\\a_2\equiv C_n^m\ (\bmod \ p_2^{c_2})\\\cdots\\a^k\equiv C_n^m\ (\bmod \ p_k^{c_k})\end{cases} \]在求出所有的 \(a_i\) 后，就可以利用中国剩余定理求解出 \(C_n^m\) 。
所以问题转化为求出 \(C_n^m\bmod p^\alpha\) 的值。
2.
但是同余方程 \(ax\equiv 1(\bmod \ p)\) (乘法逆元) 有解的充要条件为 \(\gcd(a,p)=1\) 。
然而题目显然无法保证有解，所以无法直接求 \({\rm inv}_{m!}\) 和 \({\rm inv}_{(n-m)!}\)
所以可将原式转化为：
\[\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\times \frac{(n-m)!}{p^z}}\times p^{x-y-z}\bmod p^\alpha \]其中 \(x\) 为 \(n!\) 中包含 \(p\) 的个数，\(y,z\) 同理。
3.
若对于一个 \(n\) 可以求出 \(\frac{n!}{p^x}\) ，那么就可以求逆元了，所以现在要求 \(\frac{n!}{p^x}\bmod p^k\) 。
对 \(n!\) 进行变形可得：
\[\begin{array}{l}n!=1\times 2\times \cdots\times n\\ =(p\times 2p\times \cdots)\times\prod_{i=1,i\not\equiv 0(\bmod p)}^n i\\=p^{\lfloor\frac{n}{p}\rfloor}\times (\lfloor\frac{n}{p}\rfloor)!\times\prod_{i=1,i\not\equiv 0(\bmod p)}^n i\\=p^{\lfloor\frac{n}{p}\rfloor}\times (\lfloor\frac{n}{p}\rfloor)!\times(\prod_{i=1,i\not\equiv 0(\bmod p)}^{p^k} i)^{\lfloor\frac{n}{p^k}\rfloor}\times(\prod_{i=p^k\lfloor\frac{n}{p^k}\rfloor,i\not\equiv 0(\bmod p)}^ni)\end{array} \]设 \(f(n)=\frac{n!}{p^x}\ (f(0)=1)\)
\[\therefore f(n)=f(\lfloor\frac{n}{p}\rfloor)\times(\prod_{i=1,i\not\equiv 0(\bmod p)}^{p^k} i)^{\lfloor\frac{n}{p^k}\rfloor}\times(\prod_{i=p^k\lfloor\frac{n}{p^k}\rfloor,i\not\equiv 0(\bmod p)}^ni)\\ \]\[\therefore 原式=\frac{f(n)}{f(m)f(n-m)}p^{x-y-z}\bmod p^k \]这样就可以用 \(exgcd\) 求逆元了。
最后利用中国剩余定理得到答案。单次时间复杂度 \(O(p\log p)\)
代码
#include<bits/stdc++.h>
#define ll long long
using namespace std;
inline ll read(){
ll s=0,k=1;
char c=getchar();
while(c>'9'||c<'0'){
if(c=='-')k=-1;
c=getchar();
}
while(c>='0'&&c<='9'){
s=(s<<3)+(s<<1)+(c^48);
c=getchar();
}
return s*k;
}
ll n,m,p;
ll ksm(ll a,ll b,ll mod){
ll t=1;
for(;b;b>>=1,a=a*a%mod)
if(b&1) t=t*a%mod;
return t;
}
void exgcd(ll a,ll b,ll &x,ll &y){
if(!b){
x=1,y=0;
return ;
}
exgcd(b,a%b,x,y);
ll t=x;
x=y;
y=t-a/b*y;
}
ll inv(ll n,ll mod){
ll x,y;
exgcd(n,mod,x,y);
return (x%mod+mod)%mod;
}
ll CRT(int n,ll a[],ll P[]){
ll ans=0;
for(int i=1;i<=n;i++) ans=(ans+a[i]*(p/P[i])%p*inv(p/P[i],P[i])%p)%p;
return ans;
}
ll f(ll n,ll pi,ll pk){
if(!n) return 1;
ll t=1;
for(ll i=2;i<=pk;i++) if(i%pi) (t*=i)%=pk;
t=ksm(t,n/pk,pk);
for(ll i=n/pk*pk+1;i<=n;i++) if(i%pi) (t*=(i%pk))%=pk;
return t*f(n/pi,pi,pk)%pk;
}
ll C(ll n,ll m,ll pi,ll pk){
int k=0;
for(ll i=n;i;i/=pi) k+=i/pi;
for(ll i=m;i;i/=pi) k-=i/pi;
for(ll i=n-m;i;i/=pi) k-=i/pi;
return f(n,pi,pk)*inv(f(m,pi,pk),pk)%pk*inv(f(n-m,pi,pk),pk)%pk*ksm(pi,k,pk)%pk;
}
ll exLucas(ll n,ll m,ll p){
int cnt=0;
ll a[20],P[20],tmp=p;
for(int i=2;i*i<=p;i++)
if(tmp%i==0){
P[++cnt]=1;
while(tmp%i==0) P[cnt]*=i,tmp/=i;
a[cnt]=C(n,m,i,P[cnt]);
}
if(tmp>1) P[++cnt]=tmp,a[cnt]=C(n,m,tmp,tmp);
return CRT(cnt,a,P);
}
int main(){
n=read();m=read();p=read();
printf("%lld\n",exLucas(n,m,p));
return 0;
}