组合数取模

求：

Cnm(modmo)

n,m,mo<109
Pi<109（不一定为质数）

p不是质数，(〃＞皿＜)，好恶心

根据中国剩余定理，可得：
设mo=∏pkii，Pi=pkii
当X≡ai(modpi)有唯一解时，有性质
设Mi=mo/Pi，ti=M−1i(modPi)

X=∑ai∗Mi∗ti(modmo)

很显然，只要算出每一个的CnmmodPi，再用中国剩余定理组合起来即可，

Cnm=m!n!(m−n)!

对于每一个的（X!），我们可以把它写成a∗(pi)b,这样即可以处理那些对Pi没有逆元的数，也可以在约分的时候愉快的把许多的pi约掉，
在 1∗2∗3.....∗X中,时pi的倍数的数有X/pi个，我们把X/pi，向下递归的同时，也可以获得X/pi个b（乘积中多了X/pi个pi）,
剩下的不是pi的倍数的数，我们把它们分成（(1∗2∗3.....∗(Pi−1))∗((Pi+1)∗...∗(2∗Pi−1))∗.....∗(剩下的)）(注意这些都是Pi)
显然，这些数每一组的乘积(modPi)都是一定的，于是搞出一组后就可以用快速幂愉快的全部搞出来，
然后我们又可以愉快的发现，一堆东西可以预处理(Pi不会太大)，
于是。。。。。
愉快的搞定

---

标在这里：

#include<cstdio>#include<cstdlib>#include<cstring>#include<algorithm>#define fo(i,a,b) for(int i=a;i<=b;i++)using namespace std;typedef long long ll;typedef double db; ll read(ll &n){    char ch;int q=1;n=0;    for(ch=getchar();(ch<'0' || ch>'9')&&(ch!='-');ch=getchar());    if(ch=='-')q=-1,ch=getchar();    for(;ch>='0' && ch<='9';ch=getchar())n=n*10+ch-48;    n*=q;return n;}ll P,n,n1,n2,m,ans,mo;ll p[100][4],p1[100][15000],p0;//0 is p,1 is k,2 is P^k,3 is IEll a[10];struct qww{ll a,k;};ll exgcd(ll a,ll b,ll &x,ll &y){    if(!b){x=1,y=0;return a;}    ll zhi=exgcd(b,a%b,x,y);    ll t=x;x=y,y=t-a/b*y;    return zhi;}ll IE(ll a,ll mo){    ll x=0,y=0;    exgcd(a,mo,x,y);    return (x%mo+mo)%mo;}ll ksm(ll a,ll w,ll mo){    ll ans=1;    while(w)    {if(w%2)(ans*=a)%=mo;(a*=a)%=mo,w>>=1;}    return ans;}qww Ch(ll q,int mn){    qww ans;ans.a=1;ans.k=0;ll mo=p[mn][2];    if(q<2)return ans;    ans=Ch(q/p[mn][0],mn);    ans.k+=q/p[mn][0];    (ans.a*=p1[mn][q%mo]*ksm(p1[mn][mo-1],q/mo,mo)%mo)%=mo;    return ans;}ll C(ll m,ll n){    ll ans=0;    fo(i,1,p0)    {        qww x,y,z;        x=Ch(m,i);        y=Ch(n,i);        z=Ch(m-n,i);        x.k-=y.k+z.k;x.a=(x.a*IE(y.a*z.a,p[i][2]))%mo;        (ans+=x.a*p[i][3]%mo*ksm(p[i][0],x.k,mo)%mo*(mo/p[i][2]))%=mo;    }    return ans;}int main(){    ll q,w,_;    read(n),read(m),mo=q=read(P);    for(int i=2;i*i<=q;i++)        if(!(q%i))        {            p[++p0][0]=i,p[p0][2]=1;            while(!(q%i))q/=i,p[p0][1]++,p[p0][2]*=i;            p[p0][3]=IE(mo/p[p0][2],p[p0][2]);            p1[p0][0]=1;            fo(j,1,p[p0][2]-1)p1[p0][j]=p1[p0][j-1]*((j%p[p0][0])?j:1)%p[p0][2];        }    if(q>1)    {        p[++p0][0]=p[p0][2]=q,p[p0][1]=1,p[p0][3]=IE(mo/p[p0][2],p[p0][2]);        p1[p0][0]=1;        fo(j,1,p[p0][2]-1)p1[p0][j]=p1[p0][j-1]*((j%p[p0][0])?j:1)%p[p0][2];    }    printf("%lld\n",C(m,n));    return 0;}