【JZOJ 3851】reward

Description

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$的余数.
对于$p$：设$p=p_1^{c_1} * p_2^{c_2} * p_3^{c_3} *\cdots *p_t^{c_t}$，$p_i$为质数。
对于100%的数据：$1≤n,m≤10^9，1≤pi^{ci}≤10^5$，所有p不超过$2^{31}-1$。

Analysis

想到了什么，SDOI2013方程！
抽象：$X_1+X_2+\cdots X_n=m,0\leq X_i\leq m$
形象：将$m$个球放入$n$个盒子，某个盒子可能为空。
为空？转化，把上式每个$X_i+1$，等式右边就加$n$。
这样又变成不能为空的情况了，对于一个$m$答案就是$C_{n+m-1}^{n-1}$
总的答案就是$\sum_{k=0}^{m}C_{n+k-1}^{n-1}=C_{n+m}^n$，至于为什么，可以自己画一个杨辉三角出来。
于是变成了标准的组合数取模。具体方法可以看上面的题目链接的题解。

Code

#include<cstdio>
#include<algorithm>
#define fo(i,a,b) for(int i=a;i<=b;i++)
using namespace std;
typedef long long ll;
const int N=100010;
int num,tot;
ll n,m,p,fac[N];
struct lyd
{
    ll x,y,z;
}c[N];
void prime(ll n)
{
    num=0;
    for(int i=2;i*i<=n;i++)
        if(n%i==0)
        {
            c[++num].x=i,c[num].y=0;
            for(;n%i==0;n/=i,c[num].y++);
        }
    if(n>1) c[++num].x=n,c[num].y=1;
}
ll qmi(ll x,ll n,ll mo)
{
    ll t=1;
    for(;n;n>>=1)
    {
        if(n&1) t=t*x%mo;
        x=x*x%mo;
    }
    return t;
}
void exgcd(ll a,ll b,ll &x,ll &y)
{
    if(!b)
    {
        x=1,y=0;
        return;
    }
    ll x1,y1;
    exgcd(b,a%b,x1,y1);
    x=y1,y=x1-a/b*y1;
}
ll ny(ll n,ll mo)
{
    ll x,y;
    exgcd(n,mo,x,y);
    return (x%mo+mo)%mo;
}
ll calcfac(ll n,int k,int lyd)
{
    if(n<c[k].x) return fac[n];
    ll mo=c[k].z;
    tot+=lyd*(n/c[k].x);
    return fac[n%mo]*qmi(fac[mo],n/mo,mo)%mo*calcfac(n/c[k].x,k,lyd)%mo;
}
ll CC(ll m,ll n,int k)
{
    ll mo=c[k].z;
    fac[0]=1;
    fo(i,1,mo)
        if(i%c[k].x==0) fac[i]=fac[i-1];
        else fac[i]=fac[i-1]*i%mo;
    tot=0;
    ll t=calcfac(m,k,1)*ny(calcfac(n,k,-1)*calcfac(m-n,k,-1)%mo,mo)%mo;
    return t*qmi(c[k].x,tot,mo)%mo;
}
ll C(ll m,ll n)
{
    ll ans=0;
    fo(i,1,num)
    {
        ll Mi=p/c[i].z,mo=c[i].z;
        (ans+=Mi*ny(Mi,mo)%p*CC(m,n,i)%p)%=p;
    }
    return ans;
}
int main()
{
    scanf("%lld %lld %lld",&n,&m,&p);
    prime(p);
    fo(i,1,num) c[i].z=qmi(c[i].x,c[i].y,p+1);
    printf("%lld\n",C(n+m,n));
    return 0;
}