bzoj2142 礼物 ( 扩展Lucas )

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,w_1)* C(sum-w_1,w_2)*...C(sum-...w_{n-1}，w_n) \mod P$

考虑求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) * p^{a-b-c}$
a，b，c可以递归求解 f(x)=f(⌊x/p⌋) +⌊x/p⌋
或者就那样一直除下去。

那么t1，t2，t3 呢 ?
对于1~n ,由于那些 模$p$不为0的数在模$p^k$意义下以$p^k$为周期。
例如：
$(1*2*4*5*7)\equiv(10*11*13*14*16)\pmod {{p_i}^{k_i}}$

于是先求$p^k$长度内的乘积，快速幂一下，再求剩下的 n%(p^k)个数中模$p$不为0的数的积，乘起来得到t。

至此，所有需要的求出，按照
$C(n,m) mod \ p^k =t1 * inv (t2) * inv (t3) * p^{a-b-c}$
再用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;
}