[BZOJ2142]礼物-扩展lucas定理-中国剩余定理

礼物

Description

一年一度的圣诞节快要来到了。每年的圣诞节小E都会收到许多礼物，当然他也会送出许多礼物。不同的人物在小E心目中的重要性不同，在小E心中分量越重的人，收到的礼物会越多。小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

HINT

【样例说明】

下面是对样例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。

题不是很难，然而被Impossible坑了整整两个小时

(╯‵□′)╯︵┻━┻

思路:

首先可以得出求答案的式子:

ans=(nw[1])∗(n−w[1]w[2])∗(n−w[1]−w[2]w[3])...

拆式子化简？

确实可以拆，然而完全不需要。

考虑到模数不为质数，将模数拆成若干个质数的次幂分开计算，最后使用中国剩余定理合并答案。

考虑到模数依旧不是质数而是质数的某个次幂，采用扩展lucas定理计算组合数。

那么就做完了。无比的暴力。

唯一需要注意的一点是，记得判Impossible……

#include<bits/stdc++.h>

using namespace std;

typedef long long ll;
typedef pair<ll,ll> pr;
const ll Inf=1e18;

ll p,n,m;
ll w[19];
pr ans[10009],ori[10009];

inline ll read()
{
ll x=0;char ch=getchar();
while(ch<'0' || '9'<ch)ch=getchar();
while('0'<=ch && ch<='9')x=x*10+(ch^48),ch=getchar();
return x;
}

inline ll qpow(ll a,ll b,ll p=Inf)
{
ll ret=1;
while(b)
{
if(b&1)
ret=ret*a%p;
a=a*a%p;
b>>=1;
}
return ret;
}

inline ll phi(ll p,ll t)
{
return (p-1)*qpow(p,t-1);
}

inline ll inv(ll a,ll p,ll t,ll k)
{
ll phis=phi(p,t);
return qpow(a,phis-1,k);
}

inline ll fac(ll x,ll p,ll k)
{
if(!x)return 1;
ll ret=1;
for(ll i=1;i<=k;i++)
if(i%p)
ret=ret*i%k;
ret=qpow(ret,x/k,k);
for(ll i=x/k*k+1ll;i<=x;i++)
if(i%p)
ret=ret*i%k;
return ret*fac(x/p,p,k)%k;
}

inline ll c(ll n,ll m,ll p,ll t)
{
ll cnt=0,k=qpow(p,t);

for(ll i=n;i;i=i/p)cnt+=i/p;
for(ll i=m;i;i=i/p)cnt-=i/p;
for(ll i=n-m;i;i=i/p)cnt-=i/p;

ll ret=qpow(p,cnt,k)*fac(n,p,k)%k;
ret=ret*inv(fac(m,p,k),p,t,k)%k;
ret=ret*inv(fac(n-m,p,k),p,t,k)%k;
return ret;
}

inline ll calc(ll p,ll t)
{
ll ret=1,tot=n,k=qpow(p,t);
for(int pos=1;pos<=m;pos++)
{
ret=ret*c(tot,w[pos],p,t)%k;
tot-=w[pos];
}
return ret;
}

pr invv(ll x,ll y)
{
if(!y)return pr(1,0);
pr tmp=invv(y,x%y);
return pr(tmp.second,tmp.first-x/y*tmp.second);
}

inline ll crt(pr *aa,int n)
{
ll m=p,ret=0;
for(int i=1;i<=n;i++)
{
ll x=aa[i].first,a=aa[i].second;
ll mi=m/x,invs=(invv(mi,x).first%m+m)%m;
ret=(ret+a*mi*invs)%m;
}
return ret;
}

inline ll work(ll p)
{
int top=0;
for(ll i=2;i*i<=p;i++)
{
if(p%i==0)
{
int cnt=0;
while(p%i==0)
p/=i,cnt++;
ans[++top]=pr(qpow(i,cnt),calc(i,cnt));
}
}
if(p!=1)ans[++top]=pr(p,calc(p,1));
return crt(ans,top);
}

int main()
{
p=read();n=read();m=read();
ll sum=0;
for(int i=1;i<=m;i++)
w[i]=read(),sum+=w[i];
if(sum>n)
{
puts("Impossible");
return 0;
}
if(n>sum)w[++m]=n-sum;
printf("%lld\n",work(p));
return 0;
}