參考：http://blog.csdn.net/wzq_qwq/article/details/46709471
首先推組合數，設sum為每個人禮物數的和，那麽答案為
\[ ( C_{n}^{sum}C_{sum}^{w[1]}c_{sum-w[1]}^{w[2]}... \]
設w[0]=n-sum，然後化簡成階乘的形式：
\[ \frac{n!}{w[0]!w[1]!...w[n]!} \]
註意到這裏p不是質數，所以把p拆成質數的方相乘的形式，最後用中國剩余定理合並即可
然後現在的問題是怎麽快速求出階乘
假設當前的質數的方為p=3那麽1x2x3x4x5x6x7x8x9x10x11=1x2x4x5x7x8x10x11x 3x(1x2x3)，註意到後面又是一個階乘，但是範圍更小，所以可以遞歸來做，然後前面乘的3被模消去了

#include<iostream>
#include<cstdio>
using namespace std;
const int N=100005;
long long P,n,m,w[10],p[N],cnt[N],mod[N],tot,sum,a[N];
struct qwe
{
    int a,b;
};
void exgcd(long long a,long long b,long long &x,long long &y,long long &d)
{
    if(!b)
    {
        x=1;
        y=0;
        d=a;
        return
 
;
    }
    exgcd(b,a%b,y,x,d);
    y=y-a/b*x;
}
long long china()
{
    long long d,x=0,y;
    for(int i=1;i<=tot;i++)
    {
        long long r=P/mod[i];
        exgcd(mod[i],r,d,y,d);
        x=(x+r*y*a[i])%P;
    }
    return (x+P)%P;
}
long long ksm(long long a,long long b,long long mod)
{
    long
 
 long r=1ll;
    while(b)
    {
        if(b&1)
            r=r*a%mod;
        a=a*a%mod;
        b>>=1;
    }
    return r;
}
long long inv(long long a,long long b)
{
    long long x,y,d;
    exgcd(a,b,x,y,d);
    return (x%b+b)%b;
}
qwe fac(long long k,long long n)
{
    qwe r;
    if(!n)
    {
        r.a=0,r.b=1;
        return r;
    }
    long long x=n/p[k],y=n/mod[k],ans=1ll;
    if(y)
    {
        for(int i=2;i<mod[k];i++)
            if(i%p[k]!=0)
                ans=ans*i%mod[k];
        ans=ksm(ans,y,mod[k]);
    }
    for(int i=y*mod[k]+1;i<=n;i++)
        if(i%p[k]!=0)
            ans=ans*i%mod[k];
    qwe tmp=fac(k,x);
    r.a=x+tmp.a,r.b=ans*tmp.b%P;
    return r;
}
long long clc(int k,long long n,long long m)
{
    if(n<m)
        return 0;
    qwe a=fac(k,n),b=fac(k,m),c=fac(k,n-m);
    return ksm(p[k],a.a-b.a-c.a,mod[k])*a.b%mod[k]*inv(b.b,mod[k])%mod[k]*inv(c.b,mod[k])%mod[k];
}
long long wk(long long n,long long m)
{
    for(int i=1;i<=tot;i++)
        a[i]=clc(i,n,m);
    return china();
}
int main()
{
    scanf("%lld%lld%lld",&P,&n,&m);
    for(int i=1;i<=m;i++)
        scanf("%lld",&w[i]),sum+=w[i];
    int x=P;
    for(int i=2;i*i<=x;i++)
        if(x%i==0)
        {
            p[++tot]=i;
            mod[tot]=1;
            while(x%i==0)
            {
                x/=i;
                cnt[tot]++;
                mod[tot]*=i;
            }
        }
    if(x>1)
    {
        p[++tot]=x;
        mod[tot]=x;
        cnt[tot]=1;
    }
    if(sum>n)
    {
        puts("Impossible");
        return 0;
    }
    long long ans=wk(n,sum)%P;
    for(int i=1;i<=m;i++)
    {
        ans=ans*wk(sum,w[i])%P;
        sum-=w[i];
    }
    printf("%lld\n",ans);
    return 0;
}

bzoj 2142: 禮物【中國剩余定理+組合數學】