bzoj2142 礼物
原题地址:http://www.lydsy.com/JudgeOnline/problem.php?id=2142
题意:
小E从商店中购买了n件礼物,打算送给m个人,其中送给第i个人礼物数量为wi。请你帮忙计算出送礼物的方案数(两个方案被认为是不同的,当且仅当存在某个人在这两种方案中收到的礼物不同)。
由于方案数可能会很大,你只需要输出模P后的结果。
数据范围
设P=p1^c1 * p2^c2 * p3^c3 * … *pt ^ ct,pi为质数。
对于100%的数据,1≤n≤109,1≤m≤5,1≤pi^ci≤10^5。
题解:
答案肯定是
C(n,sum)∗C(sum,w1)∗C(sum−w1,w2)∗...C(sum−...wn−1,wn)modP
考虑求C(n,m)%P,
这个 P=p1^c1 * p2^c2 * p3^c3 * … *pt ^ ct,pi为质数。
想到分开求C(n,m)% pi^ci 再用CRT合并。
C(n,m)=n!/(m!(n-m)!)
我们不能直接求分子的逆元,因为分子与pi^ci并不互质。
于是考虑把n!,m!,(n-m)!的pi都提出来。
就转化成:
n!=p^a * t1
m!=p^b * t2
(n-m)!=p^c *t3
那么C(n,m) mod p^k 转化为:
t1∗inv(t2)∗inv(t3)∗pa−b−c
a,b,c可以递归求解 f(x)=f(⌊x/p⌋) +⌊x/p⌋
或者就那样一直除下去。
那么t1,t2,t3 呢 ?
对于1~n ,由于那些 模
p
不为0的数在模
例如:
(1∗2∗4∗5∗7)≡(10∗11∗13∗14∗16)(modpiki)
于是先求 pk 长度内的乘积,快速幂一下,再求剩下的 n%(p^k)个数中模 p 不为0的数的积,乘起来得到t。
至此,所有需要的求出,按照
再用CRT合并即可。
代码:
#include<cstdio>
#include<iostream>
#include<cstring>
#include<algorithm>
#define LL long long
using namespace std;
int P,N,M,p[100005],ai[100005],mi[100005],w[10],cnt=0;
LL sum=0;
void init(int x)
{
for(int i=2;i<=x;i++)
{
if(x%i==0)
{
p[++cnt]=i;
mi[cnt]=i;
x=x/i;
while(x%i==0)
{
mi[cnt]=mi[cnt]*i;
x=x/i;
}
}
if(x==1) break;
}
}
int modpow(int A,int B,int mod)
{
int ans=1; int base=A;
for(;B;B>>=1)
{
if(B&1) ans=(1LL*ans*base)%mod;
base=(1LL*base*base)%mod;
}
return ans;
}
int cal(int x,int k)
{
if(x==0) return 1;
int ret=1;
for(int i=1;i<mi[k];i++) //不含p的因子的一个完整循环节
if(i%p[k]!=0) ret=(1LL*ret*i)%mi[k];
ret=modpow(ret,x/mi[k],mi[k]);
for(int i=1;i<=x%mi[k];i++) //一个不完整循环节
if(i%p[k]!=0) ret=(1LL*ret*i)%mi[k];
return (1LL*ret*cal(x/p[k],k))%mi[k];
}
void exgcd(int a,int b,int &x,int &y) // a*x+b*y=1
{
if(b==0)
{
x=1; y=0; return;
}
int x0,y0; //x*a+y*b=x0*b+y0*(a-(a/b)*b)
exgcd(b,a%b,x0,y0);
x=y0;
y=x0-(a/b)*y0;
}
int inv(int a,int b)
{
int x,y;
exgcd(a,b,x,y);
return (x+b)%b;
}
void exlucas(int n,int m,int k)
{
if(n<m) {ai[k]=0; return;}
int t1=cal(n,k); int t2=cal(m,k); int t3=cal(n-m,k);
int c=0;
for(int i=n;i;i/=p[k]) c+=(i/p[k]);
for(int i=m;i;i/=p[k]) c-=(i/p[k]);
for(int i=n-m;i;i/=p[k]) c-=(i/p[k]);
int po=modpow(p[k],c,mi[k]);
int iv=(1LL*inv(t2,mi[k])*inv(t3,mi[k]))%mi[k];
int ans=(1LL*po*iv)%mi[k];
ans=(1LL*ans*t1)%mi[k];
ai[k]=ans;
}
int CRT()
{
int ans=0;
for(int i=1;i<=cnt;i++)
{
int Mi=P/mi[i];
int Ri=inv(Mi,mi[i]);
ans=(ans+(1LL*((1LL*Mi*Ri)%P)*ai[i])%P)%P;
}
return ans;
}
LL getans(int n,int m) //C(n,m)
{
for(int i=1;i<=cnt;i++) // 所有 p^k
exlucas(n,m,i);
return CRT();
}
int main()
{
scanf("%d%d%d",&P,&N,&M);
for(int i=1;i<=M;i++)
{
scanf("%d",&w[i]); sum+=w[i];
}
if(sum>N)
{
printf("Impossible\n");
return 0;
}
init(P);
LL ans=1;
ans=(1LL*ans*getans(N,sum))%P;
for(int i=1;i<=M;i++)
{
ans=(1LL*ans*getans(sum,w[i]))%P;
sum-=w[i];
}
printf("%I64d\n",ans);
return 0;
}