计算C（n，m） % p，p不一定是质数
p=p1^k1 * p2^k2 * p3^k3 ………
我们可以求出C（n，m） ≡ ai mod pi^ki
对于方程组 x ≡ ai mod pi^ki
那么有C（n，m） ≡ x mod p
因为pi^ki 两两互质，所以如果已知ai，x可用中国剩余定理 计算
那么如何求出ai，即 C（n，m）% pi^ki
C(n,m)= n! / [ m! * (n-m)! ]
求出 n！，m！，（n-m）！，再求出 后两者的 逆元即可
如何求n！% pi^ki
假设 n=19 ，pi=3,ki=2
n=1*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19
=(1*2*4*5*7*8*10*11*13*14*16*17*19)*(3*6*8*12*15*18)
=(1*2*4*5*7*8)*(10*11*13*14*16*17)*19* 3^6 *(1*2*3*4*5*6)
将所有pi的倍数提出来
A、前两个括号 以pi^ki为一个周期，在模pi^ki意义下 余数相同，所以 只计算到pi^ki，再算几次幂即可
B、多出来的19，多出来的数不会超过pi^ki 个，枚举
C、3^6，pi 的倍数 先不管
D、最后一个括号，n/pi 下取整 的阶乘，递归下去即可
m！，（n-m）！ 可能在 模pi^ki 意义下没有逆元
所以把 其中pi的倍数提出来
就是说 我们算的是 (a*pi^x ) / (b*pi^y) % pi^ki
一般来说是 乘上 (b*pi^y) 的逆元，但不能保证互质，所以可能没有逆元
所以 把它变成 pi^(x-y) * (a/b) ，计算a/b 在模pi^ki 意义下的逆元
(a*pi^x ) / (b*pi^y) % pi^ki = pi^(x-y) * （a/b）* (a/b) ^(-1) % pi^ki
模板题：http://codeforces.com/problemset/gymProblem/100633/J
#include<cstdio> using namespace std; typedef long long LL; int mod; int Pow(int a,LL b,int m) { int res=1; for(;b;a=1LL*a*a%m,b>>=1) if(b&1) res=1LL*res*a%m; return res; } int get_fac(LL n,int pk,int p) { if(!n) return 1; int ans=1; if(n/pk) { for(int i=1;i<=pk;++i) if(i%p) ans=1LL*ans*i%pk; ans=Pow(ans,n/pk,pk); } int m=n%pk; for(int i=1;i<=m;++i) if(i%p) ans=1LL*ans*i%pk; return 1LL*ans*get_fac(n/p,pk,p)%pk; } void exgcd(int a,int b,int &x,int &y) { if(!b) { x=1; y=0; return; } exgcd(b,a%b,y,x); y-=a/b*x; } int get_inv(int a,int b) { int x,y; exgcd(a,b,x,y); x=(x%b+b)%b; return x; } int exLucas(LL n,LL m,int pk,int p) { int fn=get_fac(n,pk,p); int fm=get_fac(m,pk,p); int fnm=get_fac(n-m,pk,p); LL k=0; for(LL i=n;i;i/=p) k+=i/p; for(LL i=m;i;i/=p) k-=i/p; for(LL i=n-m;i;i/=p) k-=i/p; int ans=1LL*fn*get_inv(fm,pk)*get_inv(fnm,pk)%pk*Pow(p,k,pk)%pk; return ans=1LL*ans*(mod/pk)*get_inv(mod/pk,pk)%mod; } int main() { LL n,m; int ans=0; int pk; scanf("%I64d%I64d%d",&n,&m,&mod); for(int t=mod,i=2;i<=mod;++i) if(!(t%i)) { pk=1; while(!(t%i)) pk*=i,t/=i; ans+=exLucas(n,m,pk,i); ans-=ans >=mod ? mod : 0; } printf("%d",ans); }
