【BZOJ2142】礼物
Description
小E从商店中购买了n件礼物，打算送给m个人，其中送给第i个人礼物数量为wi。请你帮忙计算出送礼物的方案数（两个方案被认为是不同的，当且仅当存在某个人在这两种方案中收到的礼物不同）。由于方案数可能会很大，你只需要输出模P后的结果。
Input
输入的第一行包含一个正整数P，表示模；
第二行包含两个整整数n和m，分别表示小E从商店购买的礼物数和接受礼物的人数；
以下m行每行仅包含一个正整数wi，表示小E要送给第i个人的礼物数量。
Output
若不存在可行方案，则输出“Impossible”，否则输出一个整数，表示模P后的方案数。
Sample Input
100
4 2
1
2
4 2
1
2
Sample Output
12
【样例说明】
下面是对样例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为质数。
对于100%的数据，1≤n≤109，1≤m≤5，1≤pi^ci≤10^5。
【样例说明】
下面是对样例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为质数。
对于100%的数据，1≤n≤109，1≤m≤5，1≤pi^ci≤10^5。
题解：答案很简单，$ans=C_{tot}^{w1}C_{tot-w1}^{w2}C_{tot-w1-w2}^{w3}...$。
但是问题来了，首先模数既不是质数，也不是若干独立质数的乘积，而是质数的幂次之积。我们还是先用中国剩余定理，将模数变成质数的幂次，但此时普通的lucas定理无法使用，我们怎么求组合数呢？
我们考虑将组合数变成阶乘相除的形式，再将分子和分母中的阶乘都化成$a imes p^b$的形式（a对$p^c$取模），最后答案的a等于分子的a乘上分母的a的逆元，答案的b等于分子的b-分母的b。那么我们考虑如何将一个数的阶乘化成$a imes p^b$的形式。
以$p=5，p^c=25$为例，然后将n的阶乘中5的倍数都提出来考虑：
$n!=(1 imes2 imes3 imes4 imes6 imes7 imes...) imes5^{lfloorfrac n 5 floor} imes(lfloorfrac n 5 floor)!$
那么如何处理前面的那坨东西呢？我们可以先预处理出$1...24$的阶乘（不计算p的倍数），即为jc，然后前面的那堆东西就变成了$jc[24]^{lfloor frac n {25} floor} imes jc[n\%25]$。
对于后面的那个阶乘，我们可以递归处理下去，每次n的大小会除以p，所以复杂度是可以接受的~
最后用CRT合并，求逆元时用exgcd即可，注意判Impossible的情况。
#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
typedef long long ll;
const int maxn=100010;
ll n,m,num,P,PP,PC,Pri,ans;
ll w[10],jc[maxn];
inline ll pm(ll x,ll y)
{
ll ret=1;
while(y)
{
if(y&1) ret=ret*x%PP;
x=x*x%PP,y>>=1;
}
return ret;
}
void exgcd(ll a,ll b,ll &x,ll &y)
{
if(!b)
{
x=1,y=0;
return ;
}
exgcd(b,a%b,x,y);
ll t=x;
x=y,y=t-a/b*x;
}
inline ll ine(ll a)
{
ll x,y;
exgcd(a,PP,x,y);
return x;
}
struct node
{
ll x,y;
node() {x=1,y=0;}
node(ll a,ll b) {x=a,y=b;}
node operator + (const node &a) const {return node(x*a.x%PP,y+a.y);}
node operator * (const ll &a) const {return node(pm(x,a),y*a);}
node operator - (const node &a) const {return node(x*(ine(a.x)+PP)%PP,y-a.y);}
};
inline node getjc(ll x)
{
if(x<P) return node(jc[x],0);
return node(pm(jc[PP-1],x/PP)*jc[x%PP]%PP,x/P)+getjc(x/P);
}
inline ll solve()
{
ll i,sum=n;
for(jc[0]=1,i=1;i<PP;i++)
{
jc[i]=jc[i-1];
if(i%P) jc[i]=jc[i]*i%PP;
}
node a,b;
for(i=1;i<=m;i++) b=b+getjc(w[i]),sum-=w[i];
a=getjc(n),b=b+getjc(sum);
a=a-b;
return a.x*pm(P,a.y);
}
inline ll CRT()
{
ll i,t=Pri;
for(i=2;i*i<=t;i++)
{
if(t%i==0)
{
P=i,PP=1,PC=0;
while(t%i==0) t/=i,PP*=i,PC++;
ans=(ans+ine(Pri/PP)*solve()*(Pri/PP))%Pri;
}
}
if(t!=1)
{
P=PP=t,PC=1;
ans=(ans+ine(Pri/PP)*solve()*(Pri/PP))%Pri;
}
return (ans+Pri)%Pri;
}
int main()
{
scanf("%lld%lld%lld",&Pri,&n,&m);
ll tmp=0;
for(ll i=1;i<=m;i++) scanf("%lld",&w[i]),tmp+=w[i];
if(tmp>n)
{
printf("Impossible");
return 0;
}
printf("%lld",CRT());
return 0;
}