【bzoj2142】礼物 组合数学+中国剩余定理

手抖看了一眼boss题，就入坑了
C(n,m)%P怎么求？
P=p1^c1 * p2^c2 * p3^c3 * … *pt ^ ct，pi为质数
Pi=pi^ci
因为P太大了，而且不是素数，所以无从下手
但是我们发现可以用中国剩余定理来解决，先算出模Pi的值，最后用中国剩余定理合并。
因为Pi只包含一个质因子，所以问题简化了许多。
C(n,m)=n!/m!(n-m)!
有一种思路，算出n!%Pi和m!和(n-m)!在%Pi意义下的逆元，然后直接乘起来
很可惜，这样是不行的
因为m!和(n-m)!有可能与Pi不是互质的，也就是有可能没有逆元
我们发现，其实分子上也有许多pi因子，那不妨把分字分母的pi因子拿出来一起处理了
问题转化成两个问题
1、计算n!中pi因子的个数
2、计算n!去掉所有的pi因子后%Pi的值


问题1：
n!中有至少一个pi因子的数的个数n/pi
n!中有至少两个pi因子的数的个数n/pi/pi
……
以此类推，可以递归的计算。


问题2：
以19为例
19!%9=(1*2*3*4*5*6*7*8*9*10*11*12*13*14*15*16*17*18*19) %9
     =(1*2*4*5*7*8*10*11*13*14*16*17*19)*3^6*(1*2*3*4*5*6) %9
先把有pi因子的数提到右边，剩下的数放到左边
容易发现，左边的数的乘积，其实有一个长度为Pi的循环节，于是我们单独计算出一个循环节的乘积后快速幂即可。

右边的部分，把pi因子提出来后，发现又是阶乘的形式，所以我们可以递归来解决了。

#include<cstdio>
#include<cstring>
#include<cstdlib>
#include<cmath>
#include<algorithm>
#include<iostream>

using namespace std;

int seq[50],tot,mod[50],inv[50];
long long sum[50];
int n,m;
long long Mod;
int fac[50][100010];
int w[10];
long long ans=1;

long long power(long long x,int y,int mod)
{
	long long ans=1;
	while (y)
	{
		if (y&1) ans=ans*x%mod;
		x=x*x%mod;
		y>>=1;
	}
	return ans;
}

void decompose()
{
	long long n=Mod;
	for (int i=2;(long long)i*i<=n;i++)
	  if (n%i==0)
	  {
	  	seq[++tot]=i;mod[tot]=1;
	  	while (n%i==0) mod[tot]*=i,n/=i;
	  }
	if (n!=1) seq[++tot]=n,mod[tot]=n; 
}

int cal(int x,int i)//除掉seq[i]因子的x!%mod[i]的值 
{
	int ans=1;
	if (x<seq[i]) ans=fac[i][x];
	else
	{
		ans=power(fac[i][mod[i]-1],x/mod[i],mod[i]);
		ans=(long long)ans*fac[i][x%mod[i]]%mod[i];
		ans=(long long)ans*cal(x/seq[i],i)%mod[i];
	}
	return ans;
}

int calc(int x,int i)
{
	int ans=0;
	if (x<seq[i]) return 0;
	else
	{
		ans+=x/seq[i];
		ans+=calc(x/seq[i],i);
	}
	return ans;
}

int C(int x,int y,int i)
{
	int num=calc(x,i)-calc(x-y,i)-calc(y,i);
	int ans=(long long)cal(x,i)*power(cal(x-y,i),mod[i]/seq[i]*(seq[i]-1)-1,mod[i])*power(cal(y,i),mod[i]/seq[i]*(seq[i]-1)-1,mod[i])*power(seq[i],num,mod[i])%mod[i];
	return ans;
}

long long C(int x,int y)
{
	for (int i=1;i<=tot;i++) 
	  inv[i]=C(x,y,i);
	long long ans=0;
	for (int i=1;i<=tot;i++) ans+=(long long)inv[i]*sum[i]*power(sum[i],mod[i]/seq[i]*(seq[i]-1)-1,mod[i])%Mod;
	return ans;
}

int main()
{
	//freopen("gift.in","r",stdin);
	//freopen("gift.out","w",stdout);
	scanf("%lld%d%d",&Mod,&n,&m);
	long long cnt=0;
	for (int i=1;i<=m;i++) {scanf("%d",&w[i]);cnt+=w[i];}
	if (cnt>n)
	{
		printf("Impossible\n");
		return 0;
	}
	decompose();
	for (int i=1;i<=tot;i++)
	{
		sum[i]=Mod/mod[i];
		fac[i][0]=1;
		for (int j=1;j<=100000;j++)
		  if (j%seq[i]!=0) fac[i][j]=(long long)fac[i][j-1]*j%mod[i];
		  else fac[i][j]=fac[i][j-1];
	}
	for (int i=1;i<=m;i++) ans=ans*C(n,w[i])%Mod,n-=w[i];
	printf("%lld\n",ans);
	return 0;
}