【BZOJ】2142 礼物

【算法】中国剩余定理

【题意】给定n件物品分给m个人，每人分到wi件，求方案数%p。p不一定是素数。

【题解】

首先考虑n全排列然后按wi划分成m份，然后对于每份内都是全排列，除以wi!消除标号影响，注意剩余的(n-W)也视为一份。

所以ans=n!/(w₁!w₂!...w_m!(n-W)!)%p

也可以从排列组合公式方面考虑，即

ans=C(n,w1)*C(n-w1,w2)*C(n-w1-w2,w3)*...*C(n-w1-w2-...-w_(m-1),wm) mod P

      =n!/w1!/w2!/.../wm!/(n-W)! mod P （by POPOQQQ）

ans=n!/(w1!w2!...wn!(n-W)!) %p

到这里已经转化为经典问题：非素数模数下，中国剩余定理求阶乘及阶乘逆元。

根据唯一分解定理，P=p1^c1*...*pk^ck，对于每个模数p^c分别求答案，最后用中国剩余定理合并。

因为模数底数p不充分大，所以要先在分子和分母处找因子p，上下约去然后用快速幂额外加入答案，剩余部分正常运算。

现在问题是从n!和(n!)^(-1)中分离出p，着重考虑从n!中分离出p，逆元照做。

fac(x)表示x!除去1~x中p的倍数后的结果，则有：

n!=fac(n)*[p^(n/p)]*(n/p)!

fac(n)为分离后的部分，[p^(n/p)]为分离出来的p，(n/p)!为分离出p后剩余的倍数。下面一一解决。

<fac(n)>

对于n!，因为取模p^c且已经分离p，所以阶乘数列以p^c为周期（p^c~2*p^c-1取模后就是0至p^c-1）

那么预处理fac[0~p^c-1](不乘p的倍数)，fac(n)就是fac[n%pc]*{fac[p^c-1]^[n/(p^c)]}。

这部分答案为fac[n%pc]*power(fac[pc-1],n/pc) mod pc。

<[p^(n/p)]>累加进因子幂数，最后分子分母的因子幂数约去后用快速幂计算。

<(n/p)!>这部分又是阶乘，递归处理。

ll calc(ll x,ll p,ll pc)
{
    if(x<p)return fac[x];
    cnt+=x/p;
    return fac[x%pc]*calc(x/p,p,pc)%pc*power(fac[pc-1],x/pc,pc)%pc;
}

最后用中国剩余定理合并结果就可以了，注意ai是余数(即我们计算的结果)，Mi是除pc[i]外的模数之积，ti是Mi的逆元，最后对M取模。

#include<cstdio>
#include<algorithm>
#include<cstring>
#define ll long long
using namespace std;
const ll maxn=300010;
ll p[maxn],pc[maxn],w[maxn],cnt,tot,fac[maxn],n,m,pp,M[maxn];
void gcd(ll a,ll b,ll &x,ll &y)
{
    if(!b){x=1;y=0;}
    else {gcd(b,a%b,y,x);y-=x*(a/b);}
}
ll get_inv(ll x,ll mods)
{
    ll xx,yy;
    gcd(x,mods,xx,yy);
    return (((xx%mods)+mods)%mods);
}
ll power(ll x,ll k,ll mods)
{
    ll ans=1;
    while(k>0)
    {
        if(k&1)ans=(ans*x)%mods;
        x=(x*x)%mods;
        k>>=1;
    }
    return ans;
}
ll calc(ll x,ll p,ll pc)
{
    if(x<p)return fac[x];
    cnt+=x/p;
    return fac[x%pc]*calc(x/p,p,pc)%pc*power(fac[pc-1],x/pc,pc)%pc;//抄程序变量名错系列QAQ 
}
ll left(ll p,ll pc)
{
    fac[0]=1;
    for(int i=1;i<=pc-1;i++)fac[i]=(fac[i-1]*(i%p==0?1:i))%pc;
    cnt=0;
    ll up=calc(n,p,pc);
    ll upnum=cnt;
    for(int i=1;i<=pc-1;i++)fac[i]=(fac[i-1]*(i%p==0?1:get_inv(i,pc)))%pc;
    cnt=0;
    ll down=1;
    for(int i=1;i<=m;i++)down=(down*calc(w[i],p,pc))%pc;
    ll downnum=cnt;
    return up*down%pc*power(p,upnum-downnum,pc)%pc;
}
int main()
{
    scanf("%lld%lld%lld",&pp,&n,&m);
    ll sum=0;
    for(int i=1;i<=m;i++){scanf("%lld",&w[i]);sum+=w[i];}
    if(sum>n){printf("Impossible\n");return 0;} 
    tot=0;ll pop=pp;
    for(ll i=2;i*i<=pop;i++)
    {
        if(pop%i==0)p[++tot]=i,pc[tot]=1;
        while(pop%i==0)pc[tot]*=i,pop/=i;
    }
    if(pop>1)p[++tot]=pop,pc[tot]=pop;
    if(sum<n)w[++m]=n-sum;
    for(int i=1;i<=tot;i++)M[i]=pp/pc[i];
    sum=0;
    for(int i=1;i<=tot;i++)
    {
        sum=(sum+(left(p[i],pc[i])*get_inv(M[i],pc[i])%pp*M[i])%pp)%pp;//中国剩余定理合并时必须模M 
    }
    printf("%lld",((sum%pp)+pp)%pp);
    return 0;
}