扩展卢卡斯定理

计算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);
}