组合数取模

求：

$$C^n_m \pmod {mo}$$
$n,m,mo<10^9$
$P_i<10^9$（不一定为质数）

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

根据中国剩余定理，可得：
设$mo=\prod p_i^{k_i}，P_i=p_i^{k_i}$
当$X \equiv a_i \pmod {p_i}$有唯一解时，有性质
设$M_i=mo/P_i，t_i=M_i^{-1} \pmod {P_i}$

$$X=\sum a_i*M_i*t_i \pmod {mo}$$
很显然，只要算出每一个的 $C^n_m \mod {P_i}$，再用中国剩余定理组合起来即可，
$$C^n_m=\dfrac{m!}{n!(m-n)!}$$
对于每一个的（X!），我们可以把它写成 $a*(p_i)^b$,这样即可以处理那些对 $P_i$没有逆元的数，也可以在约分的时候愉快的把许多的 $p_i$约掉，
在 $1*2*3.....*X$中,时 $p_i$的倍数的数有 $X/p_i$个，我们把 $X/p_i$，向下递归的同时，也可以获得 $X/p_i$个b（乘积中多了 $X/p_i$个 $p_i$）,
剩下的不是 $p_i$的倍数的数，我们把它们分成 $（(1*2*3.....*(P_i-1))*((P_i+1)*...*(2*P_i-1))*.....*(剩下的)）$(注意这些都是 $P_i$)
显然，这些数每一组的乘积( $\mod P_i$)都是一定的，于是搞出一组后就可以用快速幂愉快的全部搞出来，
然后我们又可以愉快的发现，一堆东西可以预处理( $P_i$不会太大)，
于是。。。。。
愉快的搞定

---

标在这里：

#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 IE
ll 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;
}