BZOJ2142 礼物 扩展lucas 快速幂 数论

时间:2023-03-09 20:05:46
BZOJ2142 礼物 扩展lucas 快速幂 数论

原文链接http://www.cnblogs.com/zhouzhendong/p/8110015.html


题目传送门 - BZOJ2142


题意概括

  小E购买了n件礼物,送给m个人,送给第i个人礼物数量为wi。计算出送礼物的方案数模P后的结果。
  设P=p1^c1 * p2^c2 * p3^c3 * … *pt ^ ct,pi为质数。
  对于100%的数据,1≤n≤10^9,1≤m≤5,1≤pi^ci≤10^5。

题解

  首先,我们可以列出答案:
  ans=∑1<=i<=n C(n,n-∑1<=j<i w[j])
  然后算那个C(a,b)显然用n!预处理不行。
  然后我们有一个叫ex_lucas的算法可以搞定他。具体自行百度。

代码

#include <cstring>
#include <cstdio>
#include <algorithm>
#include <cstdlib>
#include <cmath>
using namespace std;
typedef long long LL;
LL P,n,m,w[5];
LL px[30],py[30],cnt;
bool check(){
int tot=0;
for (int i=1;i<=m;i++)
tot+=w[i];
return tot<=n;
}
void divide_prime(LL x){
cnt=0;
for (LL i=2;x>1&&i*i<=x;i++)
if (x%i==0){
px[++cnt]=i;
while (x%i==0)
x/=i,py[cnt]++;
}
if (x>1)
px[++cnt]=x,py[cnt]=1;
}
LL Pow(LL x,LL y,LL mod){
if (y==0)
return 1LL;
LL xx=Pow(x,y/2,mod);
xx=xx*xx%mod;
if (y&1LL)
xx=xx*x%mod;
return xx;
}
void ex_gcd(LL a,LL b,LL &x,LL &y){
if (!b)
x=1,y=0;
else
ex_gcd(b,a%b,y,x),y-=a/b*x;
}
LL Inv(LL X,LL mod){
if (!X)
return 0;
LL a=X,b=mod,x,y;
ex_gcd(a,b,x,y);
x=(x%b+b)%b;
return x;
}
LL ex_lucas(LL n,LL pi,LL pk){
if (!n)
return 1LL;
LL ans=1;
for (LL i=2;i<=pk;i++)
if (i%pi)
ans=ans*i%pk;
ans=Pow(ans,n/pk,pk);
for (LL i=2;i<=n%pk;i++)
if (i%pi)
ans=ans*i%pk;
return ans*ex_lucas(n/pi,pi,pk)%pk;
}
LL C(LL n,LL m,LL pi,LL pk){
if (m>n)
return 0;
LL a=ex_lucas(n,pi,pk),b=ex_lucas(m,pi,pk),c=ex_lucas(n-m,pi,pk);
LL k=0,ans;
for (LL i=n;i;i/=pi,k+=i);
for (LL i=m;i;i/=pi,k-=i);
for (LL i=n-m;i;i/=pi,k-=i);
ans=a*Inv(b,pk)%pk*Inv(c,pk)%pk*Pow(pi,k,pk)%pk;
return ans*(P/pk)%P*Inv(P/pk,pk)%P;
}
LL C(LL n,LL m){
LL ans=0;
for (int i=1;i<=cnt;i++)
ans=(ans+C(n,m,px[i],Pow(px[i],py[i],P+1)))%P;
return ans;
}
int main(){
scanf("%lld%lld%lld",&P,&n,&m);
for (int i=1;i<=m;i++)
scanf("%lld",&w[i]);
if (!check()){
puts("Impossible");
return 0;
}
divide_prime(P);
LL ans=1;
for (int i=1;i<=m;i++)
ans=ans*C(n,w[i])%P,n-=w[i];
printf("%lld",ans);
return 0;
}