恶心东西爬、、、
我们要求解一个 \(\binom{n}{m}\mod M\)，\(M\) 是不太大的正整数，\(n,m\) 是可能比较大的正整数。
首先我们分解 \(M=\prod_{i=1}^kp_i^{x_i}\)，我们对于每一个 \(i\in[1,k]\) 求出 \(\binom{n}{m}\mod p_i^{x_i}\)，然后就会组成一个方程组，\(Ans\equiv\binom{n}{m}\pmod p_i^{x_i}\)，因为模数显然互质，所以此时的 \(Ans\) 可以通过 CRT 求解。
问题是怎么求这个 \(\binom{n}{m}\mod p^k\) 呢，我们知道有 \(a\) 在 \(\mod p\) 意义下有逆元的充要条件是 \(a\bot p\)，这个地方就显得不一定有逆元了，不能直接求解，我们考虑变形：
\[\binom{n}{m}\equiv\frac{n!}{m!(n-m)!}\equiv\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\frac{(n-m)!}{p^z}}p^{x-y-z}\pmod {p^k} \]然后求解 \(\frac{n!}{p^x}\pmod {p^k}\)，这个 \(x\) 尽量往大了取，先挑出 \(p\) 的倍数直接做。
\[\frac{n!}{p^x}\equiv p^{⌊\frac{n}{p}⌋}(⌊\frac{n}{p}⌋!)\prod_{i 不是 p 的倍数，i=1}^{n}i\pmod {p^k} \]后面那一段其实是有一个长度为 \(p^k\) 的循环节的，手玩可得，好像是一个威尔逊定理的推论，所以我们可以拆开写 >w<
\[\frac{n!}{p^x}\equiv p^{⌊\frac{n}{p}⌋}(⌊\frac{n}{p}⌋!)(\prod_{i 不是 p 的倍数，i=1}^{p^k} i)^{⌊\frac{n}{p^k}⌋}(\prod_{i 不是 p 的倍数,i=p^k⌊\frac{n}{p^k}⌋}^{n}i)\pmod {p^k} \]写成递归的函数形式，那么 \(f(n)\equiv\frac{n!}{p^x}\pmod p\)，有下柿：
\[f(n)=f(⌊\frac{n}{p}⌋)(\prod_{i 不是 p 的倍数，i=1}^{p^k} i)^{⌊\frac{n}{p^k}⌋}(\prod_{i 不是 p 的倍数,i=p^k⌊\frac{n}{p^k}⌋}^{n}i)\pmod {p^k} \]直接递归暴力做就行了谢谢喵，放回去原式，发现 \(\frac{\frac{n!}{p^x}}{\frac{m!}{p^y}\frac{(n-m)!}{p^z}}p^{x-y-z}\) 前面的分数部分可以做了，那么还差一个 \(p^{x-y-z}\) 的指数，这个东西显然是 \(f(n)\) 中省去的 \(p^{⌊\frac{n}{p}⌋}\) 的指数部分，可以直接累加这一部分。
然后终于他妈的做完了 /lh
#include<bits/stdc++.h>
#define int long long
#define up(i,l,r) for(int i=l; i<=r; ++i)
#define dn(i,r,l) for(int i=r; i>=l; --i)
using namespace std;
int exgcd(int a,int b,int &x,int &y) {
if(!b) return y=0, x=1, a;
int ans=exgcd(b,a%b,y,x);
return y-=a/b*x, ans;
}
int inv(int a,int p) {
int x, y; exgcd(a,p,x,y);
return (x%p+p)%p;
}
int ksm(int a,int b,int p) {
int res=1%p;
for( ; b; b>>=1) {
if(b&1) res=res*a%p;
a=a*a%p;
}
return res;
}
int fac(int n,int pi,int pk) {
if(!n) return 1;
int res=1;
up(i,2,pk) if(i%pi) res=res*i%pk;
res=ksm(res,n/pk,pk);
up(i,2,n%pk) if(i%pi) res=res*i%pk;
return res*fac(n/pi,pi,pk)%pk;
}
int C(int n,int m,int pi,int pk) {
int above=fac(n,pi,pk), l=fac(m,pi,pk), r=fac(n-m,pi,pk), k=0;
for(int i=n; i; i/=pi) k+=i/pi;
for(int i=m; i; i/=pi) k-=i/pi;
for(int i=n-m; i; i/=pi) k-=i/pi;
return above*inv(l,pk)%pk*inv(r,pk)%pk*ksm(pi,k,pk)%pk;
}
int exlucas(int n,int m,int p) {
int ans=0, x=p, t=sqrt(x);
up(i,2,t) if(x%i==0) {
int mul=1;
while(x%i==0) x/=i, mul*=i;
ans=(ans+C(n,m,i,mul)*(p/mul)%p*inv(p/mul,mul)%p)%p;
}
if(x>1) ans=(ans+C(n,m,x,x)*(p/x)%p*inv(p/x,x)%p)%p;
return ans;
}
signed main() {
ios::sync_with_stdio(0);
cin.tie(0);
int n, m, p;
cin >> n >> m >> p;
cout << exlucas(n,m,p);
return 0;
}