P1377发奖金

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的数。

理解一下就是算出每个条件对应的最小数并且不影响其他位子上的值