首页 > 其他好文 > 详细

P1377发奖金

时间：2017-08-27 17:14:01 阅读：361 评论：0 收藏：0 [点我收藏+]

标签：数学输入 rem while algorithm lld 中国剩余定理分配很多

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

样例输出

对于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;
}