Bsny最近公司運作不佳，本年度利潤才m元，但員工的獎金還是要發的，公司有n個員工，怎麽發獎金這個完全由老板Bsny自己決定。Bsny想要麽把這m元全發了，激勵一下員工，但具體怎麽分配方案有很多。比如m=1, n=2, 那麽可以員工1發1元，員工2發0元；也可以員工1發0元，員工2發1元，有兩種方案。

但其實，Bsny還是有點吝嗇的，他想這m元不一定全部作為獎金，可以部分留給自己，這樣的話，發獎金的方案數就更多了。還是以m=1, n=2為例子：

方案1：員工1發1元，員工2發0元

方案2：員工1發0元，員工2發1元

方案3：員工1發0元，員工2發0元

意味著老板Bsny發的獎金範圍為[0, m]。

好奇的Bsny想知道，給定n和m，他有多少種發獎金的方案？這個答案很大，所以再給定一個p，最終的答案取模p的余數.

輸入

第一行三個整數n, m, p。

輸出

僅一行，一個整數表示最終的答案取模p的余數。

樣例輸入

2 1 5

樣例輸出

3

對於p：設p=p1^c1 * p2^c2 * p3^c3 * … *pt ^ ct，pi為質數。

20%的數據：1 ≤ n, m≤ 15；

40%的數據：1≤n, m≤1000，p=10007；

60%的數據：保證t=1，ci=1，pi^ci≤10^5；

80%的數據：t≤2，ci=1，pi≤10^5；

100%的數據：1≤ n, m≤10^9，1≤pi^ci≤10^5，所有P不超過2^31-1。

這道題是多種數學方法的組合，代碼要求有點高。

首先答案是C(m,n+m),這個就用隔欄法吧，想像成m個實心球，拿到算1塊錢。這樣在中間放n-1個欄，分成n分，因為可以不拿錢，加上n個球，每個人得到的錢是分到的球數-1.。就是C(m+n-1,n-1)m from 0 to M;=C(n+m,m),畫個楊輝三角理解一下

然後就是求C(n+m,m)%P的值了

如果P是個質數的話，維護n!%p的值以及n!%p的值在%p意義下的逆元即可。

但P不是個質數，這時我們把P分解，

對於p：設p=p1^c1 * p2^c2 * p3^c3 * … *pt ^ ct，pi為質數。

如果我們知道x%(p1^c1)=a1,%(p2^c2)=a2...%(pt^ct)=at,就可以用中國剩余定理求x%p的值了

太巧妙了

參考LHQ的代碼

#include<iostream>
#include<cstdio>
#include<algorithm>
using namespace std;
#define ll long long
const ll MAXN=100005;
ll fac[MAXN],inv[MAXN],A[MAXN],M[MAXN];
ll p[MAXN],c[MAXN];
ll n,m,cntp,Mod;
ll power(ll a,ll b,ll p)
{
    //cout<<"power"<<endl;
    ll ret=1;
    while (b)
    {
        if (b&1)
        {
            ret*=a;
            if (p!=-1) ret%=p;
        }
        a*=a;
        if (p!=-1) a%=p;
        b>>=1;
    }
    return ret; 
}
ll exgcd(ll a,ll b,ll &x,ll &y)
{
    //cout<<"exgcd"<<endl;
    if (b==0) {x=1; y=0; return a;}
    ll d=exgcd(b,a%b,y,x);
    y-=a/b*x;
    return d;   
}
ll er(ll a,ll b)
{
    //cout<<"er"<<endl;
    ll x,y;
    exgcd(a,b,x,y);
    return (x%b+b)%b;
}
ll calcfac(ll n,ll &cnt,const ll &P,const ll &p)
{
    //cout<<"calfac"<<endl;
    if (n<p) return fac[n];
    ll seg=n/P,rem=n%P;
    ll ret=power(fac[P-1],seg,P);
    ret=ret*fac[rem]%P;
    cnt+=n/p;
    return ret*calcfac(n/p,cnt,P,p)%P;  
}
ll calcinv(ll n,ll &cnt,const ll &P,const ll &p)
{
    //cout<<"calcinv"<<endl;
    if (n<p) return inv[n];
    ll seg=n/P,rem=n%P;
    ll ret=power(inv[P-1],seg,P);
    ret=ret*inv[rem]%P;
    cnt-=n/p;
    return ret*calcinv(n/p,cnt,P,p)%P;  
}
void cal_c(ll n,ll m)
{
    //cout<<"cal_c"<<endl;
    for(ll i=1;i<=cntp;i++)
    {
        ll P=power(p[i],c[i],-1);
        fac[0]=1;
        for(ll j=1;j<P;j++)
        {
            fac[j]=fac[j-1];
            if(j%p[i]!=0)
            fac[j]=fac[j-1]*j%P;
        }
        inv[0]=1;
        for(ll j=1;j<P;j++)
        {
            inv[j]=inv[j-1];
            if(j%p[i]!=0)
            inv[j]=inv[j-1]*er(j,P)%P;
        }
        ll cnt=0;
        ll ret=calcfac(n,cnt,P,p[i]);
        ret=ret*calcinv(m,cnt,P,p[i])%P;
        ret=ret*calcinv(n-m,cnt,P,p[i])%P;
        ret=ret*power(p[i],cnt,P)%P;
        A[i]=ret; M[i]=P;     
    }
}
ll CRT(ll a[],ll m[],ll n)
{
    ll M=1,ans=0,x,y,Mi,d;
    for(ll i=1;i<=n;i++)
    M*=m[i];
    for(ll i=1;i<=n;i++)
    {
        Mi=M/m[i];
        d=exgcd(Mi,m[i],x,y);
        ans=(ans+Mi*x*a[i])%M;
    }
    if(ans<0) ans+=M;
    return ans;
}
int main()
{
    scanf("%lld %lld %lld",&n,&m,&Mod);
    for(ll i=2;i*i<=Mod;i++)
    if(Mod%i==0)
    {
        p[++cntp]=i;
        while(Mod%i==0){
            Mod/=i;c[cntp]++;
        }
    }
    if(Mod>1)p[++cntp]=Mod,c[cntp]=1;
    cal_c(n+m,n);
    cout<<CRT(A,M,cntp)<<endl;
}

還有中國剩余定理中兩兩互質的時候的求法。

對於n個方程x=Ai*y+Bi，令M=lcm(A1*A2*...*An)

ans=sigma(lcm(S|Ai!∈S)*Xi)其中Xi是第1個滿足=Ai*y+Bi的數。

理解一下就是算出每個條件對應的最小數並且不影響其他位子上的值

P1377發獎金