【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;
}