bzoj 2142: 礼物(组合数取模终极版) 组合数学+中国剩余定理+exgcd

题意

一年一度的圣诞节快要来到了。每年的圣诞节小E都会收到许多礼物，当然他也会送出许多礼物。不同的人物在小E心目中的重要性不同，在小E心中分量越重的人，收到的礼物会越多。小E从商店中购买了n件礼物，打算送给m个人，其中送给第i个人礼物数量为wi。请你帮忙计算出送礼物的方案数（两个方案被认为是不同的，当且仅当存在某个人在这两种方案中收到的礼物不同）。由于方案数可能会很大，你只需要输出模P后的结果。
m<=5,n,p<=10^9，p不一定为素数。

分析

显然答案就等于若干个组合数搞起来，这里就不推了。
因为p不是素数，所以显然不能用lucas定理来搞，那么我们就考虑把p分解质因数，分开求组合数然后用天朝剩余定理合并即可。
那么对于一个组合数Cmnmodpa=n!m!(n−m)!
考虑如何求n!

复制一波别人的题解：
c(n,m)我们把他看做三个阶乘，分别为n!,m!,(n-m)!
但是如果直接用逆元的话并不可以，因为有可能在这些阶乘中出现取模的数，这样的话就直接变为0了。所以我们要把这些数提取出来。
引用一名大牛的方法
n!=t1∗piu1
m!=t2∗piu2
(n−m)!=t3∗piu3
所以t1,t2,t3我们怎么求呢？
假设19!,pi=3（popoqqq大爷举的栗子）
19!=1*2…*19
19!=1∗2∗4∗5∗7∗8∗10∗11∗13∗14∗16∗17∗19∗36∗(1∗2∗3∗4∗5∗6)
所以后面的123456又是个子问题可以递归求解，而这个前面的乱七八糟一堆东西我们只需要把它乘到一起就好了！并且我们提取出来了6个pi，现在继续递归下去的话，乱七八糟的东西再次乘到一起，pi继续提取，然后继续递归就搞定n！了！
现在我们只需要t1∗niyuan(t2)∗niyuan(t3)∗piu1−u2−u3
这样我们就可以搞定所有的c(n,m)最终乘到一起就好了！

代码

#include<iostream>#include<cstdio>#include<cstdlib>#include<cstring>#include<algorithm>#include<cmath>using namespace std;typedef long long LL;LL n,m,p,a[20];int tot;LL ksm(LL x,LL y,LL mo){    LL ans=1;    while (y)    {        if (y&1) ans=ans*x%mo;        x=x*x%mo;y>>=1;    }    return ans;}LL work(LL p,LL pa,LL n,LL op){    if (n<=1) return 1;    tot+=op*(n/p);    LL x=1;    for (int i=1;i<=pa;i++)        if (i%p>0) x=x*i%pa;    x=ksm(x,n/pa,pa);    for (int i=1;i<=n%pa;i++)        if (i%p>0) x=x*i%pa;    return x*work(p,pa,n/p,op)%pa;}LL gcd(LL x,LL y){    if (!y) return x;    else return gcd(y,x%y);}void exgcd(LL a,LL b,LL &x,LL &y){    if (!b)    {        x=1;y=0;        return;    }    exgcd(b,a%b,x,y);    LL tmp=x;    x=y;    y=tmp-(a/b)*y;}void merge(LL &a1,LL &m1,LL a2,LL m2){    LL a=m1,b=m2,c=a2-a1,x,y;    exgcd(a,b,x,y);    LL d=gcd(a,b);c/=d;    x*=c;y*=c;    a1=(x*m1%(m1*m2)+a1+m1*m2)%(m1*m2);m1*=m2;}LL get_ny(LL x,LL p,LL mo){    return ksm(x,mo/p*(p-1)-1,mo);}LL solve(LL n,LL m){    LL p1=sqrt(p),pp=p,a=-1,mo=-1;    for (int i=2;i<=p1;i++)        if (pp%i==0)        {            LL w=1;            while (pp%i==0) pp/=i,w*=i;            tot=0;            LL x=work(i,w,n,1),y=work(i,w,m,-1),z=work(i,w,n-m,-1);            x=x*get_ny(y,i,w)%w*get_ny(z,i,w)%w*ksm(i,tot,w)%w;            if (a==-1) a=x,mo=w;            else merge(a,mo,x,w);        }    if (pp>1)    {        LL w=pp;        tot=0;        LL x=work(pp,w,n,1),y=work(pp,w,m,-1),z=work(pp,w,n-m,-1);        x=x*get_ny(y,pp,w)%w*get_ny(z,pp,w)%w*ksm(pp,tot,w)%w;        if (a==-1) a=x,mo=w;        else merge(a,mo,x,w);    }    return a;}int main(){    scanf("%lld",&p);    scanf("%lld%lld",&n,&m);    LL sum=0;    for (int i=1;i<=m;i++)    {        scanf("%lld",&a[i]);        sum+=a[i];    }    if (sum>n)    {        puts("Impossible");        return 0;    }    LL ans=1;    for (int i=1;i<=m;i++)    {        ans=ans*solve(n,a[i])%p;        n-=a[i];    }    printf("%lld",ans);    return 0;}