题目描述
一年一度的圣诞节快要来到了。每年的圣诞节小E都会收到许多礼物，当然他也会送出许多礼物。不同的人物在小E心目中的重要性不同，在小E心中分量越重的人，收到的礼物会越多。小E从商店中购买了n件礼物，打算送给m个人，其中送给第i个人礼物数量为wi。请你帮忙计算出送礼物的方案数（两个方案被认为是不同的，当且仅当存在某个人在这两种方案中收到的礼物不同）。由于方案数可能会很大，你只需要输出模P后的结果。
输入输出格式
输入格式：
输入的第一行包含一个正整数P，表示模；
第二行包含两个整整数n和m，分别表示小E从商店购买的礼物数和接受礼物的人数；
以下m行每行仅包含一个正整数wi，表示小E要送给第i个人的礼物数量。
输出格式：
若不存在可行方案，则输出“Impossible”，否则输出一个整数，表示模P后的方案数。
输入输出样例
100
4 2
1
2
12
100
2 2
1
2
Impossible
说明
【样例说明】
下面是对样例1的说明。
以“/”分割，“/”前后分别表示送给第一个人和第二个人的礼物编号。12种方案详情如下：
1/23 1/24 1/34
2/13 2/14 2/34
3/12 3/14 3/24
4/12 4/13 4/23
设P=p1^c1 * p2^c2 * p3^c3 * … *pt ^ ct，pi为质数。
对于15%的数据，n≤15，m≤5，pi^ci≤10^5；
在剩下的85%数据中，约有60%的数据满足t≤2，ci=1，pi≤10^5，约有30%的数据满足pi≤200。
对于100%的数据，1≤n≤10^9，1≤m≤5，1≤pi^ci≤10^5，1≤P≤10^9。
这个很好想啦，就是模数不是质数会很麻烦。
没错，这就是一道扩展卢卡斯+中国剩余定理的模板题。
大致的做法就是先把P质因数分解，然后得到在每个质因数^次数下的答案，最后用中国剩余定理合并。。
那么现在问题的关键变成了如何求n! mod p^k。
这个可以拆成三部分(准确的说是最多三部分)，
1:n!中是p倍数的可以直接拆出来，递归求(n/p)! mod p^k。
2.n!中可能有很多个完整的1~p^k-1,这部分可以先预处理一下(p^k-1)! (当然是不算p的倍数的阶乘，因为我们还要记录一下阶乘中p的次数最后一起处理)，然后快速幂一下。
3.n%(p^k)的部分可以直接乘预处理的了。
可能第一次写会比较费劲，，，，
(还有我为什么忘了中国剩余定理了23333,记住是∑a[i]*(P/mod[i])*inv(P/mod[i],mod[i]))
#include<bits/stdc++.h>
#define ll long long
#define maxn 100005
using namespace std;
int N,P,m,a[15],MOD;
int d[15],c[15],mo;
int D[15],jc[maxn];
int num,ans[15],tot; inline int add(int x,int y,const int ha){
x+=y;
if(x>=ha) return x-ha;
else return x;
} inline int ksm(int x,int y,const int ha){
int an=1;
for(;y;y>>=1,x=x*(ll)x%ha) if(y&1) an=an*(ll)x%ha;
return an;
} struct node{
int val,tmp;
node operator *(const node &U)const{
return (node){val*(ll)U.val%D[mo],tmp+U.tmp};
}
node operator /(const node &U)const{
return (node){val*(ll)ksm(U.val,D[mo]/d[mo]*(d[mo]-1)-1,D[mo])%D[mo],tmp-U.tmp};
}
}; inline void dvd(){
for(int i=2;i*(ll)i<=P;i++) if(!(P%i)){
d[++num]=i,D[num]=1;
while(!(P%i)) P/=i,D[num]*=i,c[num]++;
if(P==1) break;
}
if(P!=1) d[++num]=D[num]=P,c[num]=1;
} inline node getjc(int x){
node now=(node){1,0};
if(x>=d[mo]) now=now*getjc(x/d[mo]),now.tmp+=x/d[mo]; if(x>=D[mo]) now=now*(node){ksm(jc[D[mo]-1],x/D[mo],D[mo]),0}; now=now*(node){jc[x%D[mo]],0};
return now;
} inline node getC(int x,int y){
return getjc(x)/getjc(y)/getjc(x-y);
} inline void solve(int x){
jc[0]=1,mo=x;
for(int i=1;i<D[x];i++){
jc[i]=jc[i-1];
if(i%d[x]) jc[i]=jc[i]*(ll)i%D[x];
} ans[x]=1;
int lef=N;
node now;
for(int i=1;i<=m;i++){
now=getC(lef,a[i]);
// printf("%d %d %d %d\n",lef,a[i],now.val,now.tmp);
ans[x]=ans[x]*(ll)now.val%D[x]*(ll)ksm(d[x],now.tmp,D[x])%D[x];
lef-=a[i];
}
} inline int CRT(){
int an=0;
for(int i=1;i<=num;i++){
an=add(an,ksm(MOD/D[i],D[i]/d[i]*(d[i]-1)-1,D[i])*(ll)(MOD/D[i])%MOD*(ll)ans[i]%MOD,MOD);
}
return an;
} int main(){
scanf("%d%d%d",&P,&N,&m),MOD=P;
for(int i=1;i<=m;i++) scanf("%d",a+i),tot+=a[i]; if(tot>N){
puts("Impossible");
return 0;
} dvd(); for(int i=1;i<=num;i++) solve(i); printf("%d\n",CRT());
return 0;
}