bzoj2142 礼物（扩展Lucas+CRT）

给定n个物品，分给m个人，每个人拿到wi个礼物，问方案数mod P P不一定为质数

首先我们把剩下的礼物也分给一个人 答案明显不变 w[++m]=n-w1-w2-…-wm
那么显然答案就是$C_n^{w_1}C_{n-w_1}^{w_2}\cdots C_{n-w_1-...-w_{m-1}}^{w_m}$，也就是$\frac{n!}{w_1!w_2!\cdots w_m!}\%P$
P不是质数怎么办呢，我们设$P=\prod\limits_{i=1}^kp_i^{a_i}$,那么我们针对每一个$p_i^{a_i}$算一个余数，得到一个模线性方程组，然后再用中国剩余定理组合起来即可。（因为两两一定互质，所以普通CRT就可以啦）

现在我们的问题在于如何计算$\frac{n!}{w_1!w_2!\cdots w_m!}\%p_i^{a_i}$，因为阶乘可能和$p_i$不互质，就没有逆元，很gg。我们把分子分母都写成$xp_i^k$的形式，那么(x,pi)=1，可以算逆元了。剩下的若干pi的次幂可以直接减。

我们考虑如何计算互质的部分。去看这个典型的例子你就会啦：portal

我已经不会算复杂度了qaq

#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
#define ll long long
#define ld long double
#define inf 0x3f3f3f3f
#define N 110
inline char gc(){
    static char buf[1<<16],*S,*T;
    if(S==T){T=(S=buf)+fread(buf,1,1<<16,stdin);if(T==S) return EOF;}
    return *S++;
}
inline int read(){
    int x=0,f=1;char ch=gc();
    while(ch<'0'||ch>'9'){if(ch=='-')f=-1;ch=gc();}
    while(ch>='0'&&ch<='9') x=x*10+ch-'0',ch=gc();
    return x*f;
}
int P,nn,w[10],m,prime[N],mo[N],r[N],sum=0,tot=0;
inline void decomp(){
    int x=P;
    for(int i=2;i*i<=P;++i){
        if(x%i) continue;prime[++tot]=i;mo[tot]=1;
        while(x%i==0) x/=i,mo[tot]*=i;
    }if(x!=1) prime[++tot]=x,mo[tot]=x;
}
inline void exgcd(ll a,ll b,ll &x,ll &y){
    if(!b){x=1;y=0;return;}
    exgcd(b,a%b,x,y);ll t=x;x=y;y=t-a/b*y;
}
inline int ksm(ll x,int k,int mod){
    ll res=1;for(;k;k>>=1,x=x*x%mod) if(k&1) res=res*x%mod;return res;
}
inline int inv(int a,int mod){
    ll x,y;exgcd(a,mod,x,y);return (x%mod+mod)%mod;
}
inline int multi(int n,int p,int mod){//计算n! %mod 互质部分
    if(!n) return 1;ll res=1;
    for(int i=2;i<mod;++i) if(i%p) (res*=i)%=mod;
    res=ksm(res,n/mod,mod);
    for(int i=2;i<=n%mod;++i) if(i%p) (res*=i)%=mod;
    return res*multi(n/p,p,mod)%mod;
}
inline int calc(int n,int p,int mod){//n!/w1!w2!...wm! %mod
    ll res=multi(n,p,mod);
    for(int i=1;i<=m;++i){
        ll tmp=multi(w[i],p,mod);
        res=res*inv(tmp,mod)%mod;
    }int k=0;//计算非互质部分 
    for(ll i=p;i<=n;i*=p) k+=n/i;
    for(int i=1;i<=m;++i) for(ll j=p;j<=w[i];j*=p) k-=w[i]/j;
    return res*ksm(p,k,mod)%mod;
}
inline int CRT(int n,int mod){
    ll res=0;
    for(int i=1;i<=n;++i){
        ll a=P/mo[i],b=mo[i];
        res+=r[i]*a%mod*inv(a,b);res%=mod;
    }return res;
}
int main(){
//  freopen("a.in","r",stdin);
    P=read();nn=read();m=read();
    for(int i=1;i<=m;++i) w[i]=read(),sum+=w[i];
    if(sum>nn){puts("Impossible");return 0;}w[++m]=nn-sum;
    decomp();for(int i=1;i<=tot;++i) r[i]=calc(nn,prime[i],mo[i]);
    printf("%d\n",CRT(tot,P));
    return 0;
}