【BZOJ2142】礼物 组合数+CRT

【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

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。

题解：答案很简单，$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;
}