reward

题目大意

球的个数在0..M中，将球分到N个盒子里，每个盒子可以没有，问方案数，答案模P。
N,M,P<=10^9，将P分成互不互质的数的乘积，且这些数均可表示为质数的幂，最大的数<=10^5。

隔板问题

我们可以看作将M个球放入N+1个盒子里，由于每个盒子可以没有，所以方案数为$C_{N+M}^N$
考虑该如何进行组合数取模呢？

中国剩余定理

将P分解成互不互质的数的乘积，且这些数均可表示为质数的幂。那么我们可以列出模方程，于是可以利用中国剩余定理解出模P剩余系下的唯一解，即为组合数模P的值。于是，我们现在需要做的，是如何快速算组合数模一个质数的幂。

快速算阶乘

设p为要模的数，pp为p唯一的质因子。
我们可以统计最终结果有多少个pp，计为cnt。然后计算分子除pp外的积为A，分母除pp外的积为B，答案即为$A*B^{-1}*pp^{cnt}$模p的值。
现在我们考虑如何快速统计阶乘中有多少个pp以及阶乘除pp外的积。
举个例子，p=9，pp=3,n=20
$1*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19*20$
我们发现，每9个为一个周期，即
$1*2*4*5*7*8 Mod 9=10*11*13*14*16*17 Mod 9=19*20*22*23*25*26 Mod 9等$
我们预处理fac[i]表示$i！$中排除pp的倍数的积，可知一个周期值为fac[p-1]。
对于上面那个例子等价于
$(1*2*4*5*7*8)^2*3*6*9*12*15*18*19*20$
余下部分为fac[n%p]。
然后只剩下
$3*6*9*12*15*18$
因此cnt+=n/pp。剩余部分可以变为
$1*2*3*4*5*6$
继续递归处理即可。附代码理解：

ll calcfac(ll n,ll p,ll pp){
    if (n<pp) return fac[n];
    ll t=quicksortmi(fac[p-1],n/p,p);
    t=t*fac[n%p]%p;
    cnt+=n/pp;
    t=t*calcfac(n/pp,p,pp)%p;
    return t;
}

因此本题解决。

参考程序

#include<cstdio>
#include<algorithm>
#define fo(i,a,b) for(i=a;i<=b;i++)
using namespace std;
typedef long long ll;
ll a[20],b[20],c[20],d[20],e[20],pri[32000+10],fac[100000+10],num[100000+10];
bool bz[32000+10];
ll i,j,k,l,t,n,m,p,pp,top,xx,yy,ans,cnt;
ll quicksortmi(ll x,ll y,ll p){
    if (!y) return 1;
    if (y==1) return x%p;
    ll t=quicksortmi(x,y/2,p);
    t=t*t%p;
    if (y%2) t=t*(x%p)%p;
    return t;
}
void gcd(ll a,ll b){
    if (!b){
        xx=1;
        yy=0;
    }
    else{
        gcd(b,a%b);
        swap(xx,yy);
        yy-=xx*(a/b);
    }
}
ll getny(ll x,ll y){
    gcd(x,y);
    xx=(xx%y+y)%y;
    return xx;
}
ll calcfac(ll n,ll p,ll pp){
    if (n<pp) return fac[n];
    ll t=quicksortmi(fac[p-1],n/p,p);
    t=t*fac[n%p]%p;
    cnt+=n/pp;
    t=t*calcfac(n/pp,p,pp)%p;
    return t;
}
ll calc(ll x,ll y,ll p,ll pp){
    ll i,j;
    fac[0]=1;num[0]=0;
    fo(i,1,p-1){
        num[i]=0;
        j=i;
        while (j%pp==0){
            num[i]++;
            j/=pp;
        }
        num[i]+=num[i-1];
        if (i%pp==0) fac[i]=fac[i-1];
        else fac[i]=fac[i-1]*i%p;
    }
    cnt=0;
    ll A=calcfac(y,p,pp);
    ll tot=cnt;
    cnt=0;
    ll B=calcfac(x,p,pp);
    B=B*calcfac(y-x,p,pp)%p;
    B=getny(B,p);
    return A*B%p*quicksortmi(pp,tot-cnt,p)%p;
}
int main(){
    //freopen("E:/wzd/11.28/yj.txt","w",stdout);
    scanf("%lld%lld%lld",&n,&m,&p);
    fo(i,2,32000){
        if (!bz[i]) pri[++k]=i;
        fo(j,1,k){
            if (pri[j]*i>32000) break;
            bz[i*pri[j]]=1;
            if (i%pri[j]==0) break;
        }
    }
    pp=p;
    fo(i,1,k){
        if (pp%pri[i]==0){
            d[++top]=1;e[top]=pri[i];
            while (pp%pri[i]==0){
                d[top]*=pri[i];
                pp/=pri[i];
            }
        }
    }
    if (pp>1){
        d[++top]=pp;
        e[top]=d[top];
    }
    fo(i,1,top) c[i]=p/d[i];
    fo(i,1,top) a[i]=calc(n,n+m,d[i],e[i]);
    fo(i,1,top) b[i]=getny(c[i],d[i]);
    fo(i,1,top) ans=(ans+a[i]*b[i]%p*c[i]%p)%p;
    printf("%lld\n",ans);
    return 0;
}