Lucas&&Exlucas
Lucas和Exlucas可以求模p意义下大数的组合数。
先考虑p为质数的情况，那么直接上Lucas定理即可。
Lucas 定理基本内容：
[C_n^m=C_{n mod p}^{m mod p}*C_{n/p}^{m/p} (mod p) p是质数
]
对于Lucas的实现直接递归处理即可，注意特判n＜m的情况。
如果p不为质数，那么用Exlucas求解。
实际上，Exlucas和Lucas定理没有关系……
由于p不为质数，所以可以考虑把p质因数分解：
[p=p_1^{c_1}*p_2^{c_2}...p_t^{c_t}
]
那么可以问题可以转化成求解下列每一个式子的值，然后用CRT合并答案即可：
[C_n^m mod p_1^{c_1}\
C_n^m mod p_2^{c_2}\
...\
C_n^m mod p_t^{c_t}\
]
继续转化，把组合数写成阶乘形式,相当于要求：
[frac{n!}{m!·(n-m)!} mod p^k
]
发现直接做不好做，逆元问题都不好解决，考虑把左边每一项中的因子p给提出来:
[frac{frac{n!}{p^{a_1}}}{frac{m!}{p^{a_2}}*frac{(n-m)!}{p^{a3}}}*p^{a^1-a^2-a^3} mod p^k
]
这样逆元就可以直接用Exgcd求了。
那么现在重点解决的问题又变成了：
[n! mod p^k 其中n还要除去所有的因子p
]
举个例子：n=22，p=3，k=2
把其中所有p（也就是3）的倍数提取出来，得到：
[22!=3^7×(1×2×3×4×5×6×7)×(1×2×4×5×7×8)×(10×11×13×14×16×17)×(19×20×22) ]
可以发现：
[(1×2×4×5×7×8)equiv(10×11×13×14×16×17) (mod 3)
]
所以对于这种一段一段的可以直接处理，最后求一下它（n/p^k）次方即可。
对于3的次方不需要考虑，因为在外面会乘上来。
那么为什么不把3,6彻底分解呢？
这是为了递归的方便，可以发现存在一项7！，也就是（n/p）！，可以递归处理。
所以在一层层递归中，每次每个含有因数p的数提且仅提出一个p，那么可以实现递归，且最后可以把所有p都提出!
然后就可以开心的求出组合数了O(∩_∩)O！
IL int qpow(int x,int p,int mod) {
RG int ans=1;
for(;p;p>>=1,x=x*x%mod)
if(p&1) ans=ans*x%mod;
return ans;
}
IL int exgcd(int &x,int &y,int a,int b) {
if(!b) return x=1,y=0,a;
RG int gcd=exgcd(x,y,b,a%b),z=x;
x=y,y=z-y*(a/b);
return gcd;
}
IL int Fac(int a,int p,int pk) {
if(a==0) return 1;
RG int i,ans=1;
for(i=1;i<pk;++i)
if(i%p) ans=ans*i%pk;
ans=qpow(ans,a/pk,pk);
for(i=1;i<=a%pk;++i)
if(i%p) ans=ans*i%pk;
return ans*Fac(a/p,p,pk)%pk;
}
IL int C(int a,int b,int p,int pk) {
if(a<b) return 0;
RG int i,x,y,k=0,f1=Fac(a,p,pk),f2=Fac(b,p,pk),f3=Fac(a-b,p,pk);
exgcd(x,y,f2,pk),x=(x+pk)%pk,f2=x;
exgcd(x,y,f3,pk),x=(x+pk)%pk,f3=x;
for(i=a/p;i;i/=p) k+=i;
for(i=b/p;i;i/=p) k-=i;
for(i=(a-b)/p;i;i/=p) k-=i;
return f1*f2%pk*f3%pk*qpow(p,k,pk)%pk;
}
IL int CRT() {
RG int i,x,y,ans=0;
for(i=1;i<=tot;++i) {
exgcd(x,y,p/Pk[i],Pk[i]),x=(x+Pk[i])%Pk[i];
ans=(ans+a[i]%p*(p/Pk[i])%p*x%p)%p;
}
return ans;
}
signed main()
{
RG int i,j,x;
n=gi(),m=gi(),p=gi();
for(i=2,x=p;i*i<=x;++i) {
if(x%i==0) {
pr[++tot]=i;
while(x%i==0) ++cnt[tot],x/=i;
}
}
if(x>1) pr[++tot]=x,cnt[tot]=1;
for(i=1;i<=tot;++i) {
for(j=1,x=1;j<=cnt[i];++j) x=x*pr[i];
a[i]=C(n,m,pr[i],Pk[i]=x);
}
printf("%lld
",CRT());
return 0;
}