bzoj2142 礼物

稍有数学常识的人都知道，答案是

$$\frac{n!}{\prod_{i=1}^ma_i!(n-\sum_{i=1}^ma_i)!}$$
$p$不是质数这一点也比较好解决，分解质因数分别求解以后用中国剩余定理合并就可以了。
现在关键是求$x! \% p^q$。
把$1..x$中$p$的倍数单独拎出来，可以得到$p^ww!$，$w!$可以递归计算。剩下的部分都和$p$互质，可以暴力计算。每$p^q$个数看成一组，显然每一组答案都相同。$O(p^q)$暴力计算出这个数，快速幂之后再乘上单独的一部分就可以了。
最后我们得到的答案是$x*p^y$，其中$x$和$p$互质。做除法的时候，互质部分直接求逆元，的幂只需要指数相减。

#include<cstdio>
#include<cmath>
#include<algorithm>
using namespace std;
#define LL long long
int p[50],mod[50],a[10],res[50],tot;
int pow(int b,int k,int p)
{
    int ret=1;
    for (;k;k>>=1,b=(LL)b*b%p)
        if (k&1) ret=(LL)ret*b%p;
    return ret;
}
void cal(int x,int k,int &t,int &q)
{
    if (x<p[k])
    {
        q=0;
        t=1;
        for (int i=2;i<=x;i++) t=(LL)t*i%mod[k];
        return;
    }
    int y=1,z,t1,q1;
    for (int i=2;i<mod[k];i++)
        if (i%p[k]) y=(LL)y*i%mod[k];
    t=pow(y,x/mod[k],mod[k]);
    for (int i=mod[k]*(x/mod[k])+1;i<=x;i++)
        if (i%p[k]) t=(LL)t*i%mod[k];
    q=x/p[k];
    cal(q,k,t1,q1);
    t=(LL)t*t1%mod[k];
    q+=q1;
}
void exgcd(int a,int b,int &x,int &y)
{
    if (!b)
    {
        x=1;
        y=0;
    }
    else
    {
        exgcd(b,a%b,y,x);
        y-=a/b*x;
    }
}
int main()
{
    int tp,P,N,M,lim,s=0,x,y,w,z,q,ans=0;
    scanf("%d%d%d",&P,&N,&M);
    for (int i=1;i<=M;i++) scanf("%d",&a[i]);
    lim=sqrt(P+0.5);
    tp=P;
    for (int i=2;i<=lim;i++)
        if (P%i==0)
        {
            tot++;
            p[tot]=i;
            mod[tot]=1;
            while (P%i==0)
            {
                P/=i;
                mod[tot]*=i;
            }
        }
    if (P>1)
    {
        tot++;
        mod[tot]=p[tot]=P;
    }
    P=tp;
    for (int i=1;i<=M;i++) s+=a[i];
    if (s>N)
    {
        printf("Impossible\n");
        return 0;
    }
    if (s<N) a[++M]=N-s;
    for (int i=1;i<=tot;i++)
    {
        cal(N,i,x,q);
        for (int j=1;j<=M;j++)
        {
            cal(a[j],i,y,w);
            q-=w;
            exgcd(y,mod[i],w,z);
            w=(w%mod[i]+mod[i])%mod[i];
            x=(LL)x*w%mod[i];
        }
        x=(LL)x*pow(p[i],q,mod[i])%mod[i];
        exgcd(P/mod[i],mod[i],y,w);
        y=(y%P+P)%P;
        ans=(ans+(LL)P/mod[i]*y%P*x)%P;
    }
    printf("%d\n",ans);
}