「一本通 6.6 练习 8」礼物

Description

原题来自：BZOJ 2142

一年一度的圣诞节快要来到了。每年的圣诞节小 E 都会收到许多礼物，当然他也会送出许多礼物。不同的人物在小 E 心目中的重要性不同，在小 E 心中分量越重的人，收到的礼物会越多。

小 E 从商店中购买了 n 件礼物，打算送给 m 个人，其中送给第 iii 个人礼物数量为 w_i​i​​。请你帮忙计算出送礼物的方案数（两个方案被认为是不同的，当且仅当存在某个人在这两种方案中收到的礼物不同）。由于方案数可能会很大，你只需要输出模 PPP 后的结果。

Input

输入的第一行包含一个正整数 P，表示模数；

第二行包含两个正整数 n 和 m，分别表示小 E 从商店购买的礼物数和接受礼物的人数；

以下 mmm 行每行仅包含一个正整数 w_ii​​，表示小 E 要送给第 iii 个人的礼物数量。

Output

若不存在可行方案，则输出 Impossible，否则输出一个整数，表示模 P 后的方案数。

Sample Input

100
4 2
1
2

Sample Output

12

Sample Explaination

12 种方案详情如下： {1}{2,3},{1}{2,4},{1}{3,4},{2}{1,3},{2}{1,4},{2}{3,4},{3}{1,2},{3}{1,4},{3}{2,4},{4}{1,2},{4}{1,3},{4}{2,3}

 Hint

设 P=p₁^c₁×p₂^c₂×p₃^c₃×⋯×p_t^c_i p_i为质数1​c​1​​​​×p​2​c​2​​​​×p​3​c​3​​​​×⋯×p​t​c​t​​​​，pip_ip​i​​ 为质数。

对于 100%的数据，1≤n≤10⁹,1≤m≤5,1≤p_i^c_i≤10⁵⁹​​,1≤m≤5,1≤p​i​c​i​​​​≤10​5​​。

最近一直在做CRT，乍一看还以为又是一个水题，实际上并没有那么简单。

是时候来总结一下中国剩余定理CRT，并且开启奇妙的扩展卢卡斯定理ex_lucas了。

审题，题意很简单，要求 Π_(i=1->m)C(n-∑_(j=1->i-1)w_j,w_j)%P

组合数取模，很容易想到卢卡斯定理。

然而根据题意，P并不一定是质数，普通卢卡斯定理只适用于质数。

怎么办呢？那么你觉得我所说的ex_lucas是干什么用的呢？

如果我们把P质因数分解，不就得到了质数吗？

我们只需要求出ans%p_i^k_i 然后再CRT解关于它们的余数方程，就好了呀。

　　简单介绍CRT：

　　　　对于余数方程xΞa_i(mod p_i) 其中p_i彼此互质。

　　　　设P=∏p_i，则最小自然数解x=(∑P/p_i*a_i*s_i)%P;其中s_i表示方程(P/p_i)*x+p_i*y=1的解。

　　　　证明略。而对于pi不互质的情况需要ex_CRT，十分恶心不再这里介绍了。

//tim[i]即为pi，aans[i]即为a[i]，prr即为P
int CRT(){
    int ans=0,x,y;
    for(int i=1;i<=nom;++i){
        ex_gcd(prr/tim[i],tim[i],x,y);
        x=(x%tim[i]+tim[i])%tim[i];
        ans=(ans+prr/tim[i]*aans[i]%prr*x%prr)%prr;
    }
    return ans;
}

接下来我们继续考虑求ans%p_i^k_i，下一个问题还是组合数如何对非质数取模。

然而现在这个非质数不再是普通的合数了，而是一个质数的幂，即要求C(n,m)%p^k。

那么过程中求阶乘时不能很轻易地求出它的逆元了，因为它本身就可能含有p这个因子。

那还不简单？把阶乘的因子p都提出来不就好了嘛？

则C(n,m)%p^k=(n! / p^k1)*inv(m! / p^k2)*inv((n-m)! / p^k3)*p^k1+k2+k3

现在inv内部的东西与p互质了，可以求逆元了。

接下来再来考虑如何求(n! / p^k1)%p^k，除个p^k1不是什么难事，先讨论n！%p^k，如22！%3²

=1×2×3×4×5×6×7×8×9×10×11×12×13×14×15×16×17×18×19×20×21×22

提出含p=3的项

=3⁷×(1×2×3×4×5×6×7)×(1×2×4×5×7×8×10×11×13×14×16×17×19×20×22)

蓝色部分：是3⁷，即p^n/p，快速幂。

绿色部分：是7！，即(n/p)！，递归解决。

红色部分：是与p互质的数的乘积，其中以p^k为分段可以发现前两各完整的段里有

　　　　(1×2×4×5×7×8)Ξ(10×11×13×14×16×17) (mod p^k)

　　　　这样就有了一个循环节，快速幂解决。多余部分暴力乘出来就可以。

int fac(int n,int p,int pt,int ans=1){
    if(!n)return 1;
    for(int i=1;i<pt;++i)if(i%p)ans=ans*i%pt;
    ans=pow(ans,n/pt,pt);
    for(int i=1;i<=n%pt;++i)if(i%p)ans=ans*i%pt;
    return ans*fac(n/p,p,pt)%pt;
}

逐步回代，解决啦。

#include<cstdio>
#include<cmath>
#define int long long
using namespace std;
int m,pr,sqr,modd[100005],tim[100005],aans[100005],p[100005],nop,num[6],nom,lef[6],prr;bool np[100005];
int pow(int b,int t,int p,int ans=1){
    while(t){
        if(t&1)ans=ans*b%p;
        b=b*b%p;
        t/=2;
    }
    return ans;
}
int fac(int n,int p,int pt,int ans=1){
    if(!n)return 1;
    for(int i=1;i<pt;++i)if(i%p)ans=ans*i%pt;
    ans=pow(ans,n/pt,pt);
    for(int i=1;i<=n%pt;++i)if(i%p)ans=ans*i%pt;
    return ans*fac(n/p,p,pt)%pt;
}
void ex_gcd(int a,int b,int &x,int &y){
    if(!b){x=1;y=0;return;}
    ex_gcd(b,a%b,x,y);
    int res=x;x=y;y=res-a/b*y;
}
int inv(int n,int p){
    int x,y;
    ex_gcd(n,p,x,y);
    return (x%p+p)%p;
}
int c(int b,int t,int p,int pt){
    if(b<t)return 0;int cnt=0;
    for(int i=b;i;i/=p)cnt+=i/p;
    for(int i=t;i;i/=p)cnt-=i/p;
    for(int i=b-t;i;i/=p)cnt-=i/p;
    return fac(b,p,pt)*inv(fac(t,p,pt),pt)%pt*inv(fac(b-t,p,pt),pt)%pt*pow(p,cnt,pt)%pt;
}
int CRT(){
    int ans=0,x,y;
    for(int i=1;i<=nom;++i){
        ex_gcd(prr/tim[i],tim[i],x,y);
        x=(x%tim[i]+tim[i])%tim[i];
        ans=(ans+prr/tim[i]*aans[i]%prr*x%prr)%prr;
    }
    return ans;
}
int ex_lucas(){
    for(int i=1;i<=nop;++i){
        if(pr%p[i]==0)modd[++nom]=p[i],tim[nom]=1;
        while(pr%p[i]==0)pr/=p[i],tim[nom]*=p[i];
    }
    for(int ii=1;ii<=nom;++ii){
        aans[ii]=1;
        for(int i=1;i<=m;++i)aans[ii]=(aans[ii]*c(lef[i],num[i],modd[ii],tim[ii]))%tim[ii];
    }
    return CRT();
}
signed main(){
    scanf("%lld%lld%lld",&pr,&lef[1],&m);prr=pr;
    for(int i=1;i<=m;++i)scanf("%lld",&num[i]),lef[i+1]=lef[i]-num[i];
    if(lef[m+1]<0){puts("Impossible");return 0;}
    for(int i=2;i<=100000;++i){
        if(!np[i])p[++nop]=i;
        for(int j=1;j<=nop&&i*p[j]<=100000;++j)
            if(i%p[j])np[i*p[j]]=1;
            else{np[i*p[j]]=1;break;}
    }
    printf("%lld\n",ex_lucas());
}