首页
登录 | 注册

[BZOJ2142]礼物-扩展lucas定理-中国剩余定理

礼物

Description

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

Input

输入的第一行包含一个正整数P,表示模;
第二行包含两个整整数n和m,分别表示小E从商店购买的礼物数和接受礼物的人数;
以下m行每行仅包含一个正整数wi,表示小E要送给第i个人的礼物数量。

Output

若不存在可行方案,则输出“Impossible”,否则输出一个整数,表示模P后的方案数。

Sample Input

100
4 2
1
2

Sample Output

12

HINT

【样例说明】

下面是对样例1的说明。

以“/”分割,“/”前后分别表示送给第一个人和第二个人的礼物编号。12种方案详情如下:

1/23 1/24 1/34

2/13 2/14 2/34

3/12 3/14 3/24

4/12 4/13 4/23

【数据规模和约定】

设P=p1^c1 * p2^c2 * p3^c3 * … *pt ^ ct,pi为质数。

对于100%的数据,1≤n≤109,1≤m≤5,1≤pi^ci≤10^5。


题不是很难,然而被Impossible坑了整整两个小时

(╯‵□′)╯︵┻━┻


思路:
首先可以得出求答案的式子:
ans=(nw[1])(nw[1]w[2])(nw[1]w[2]w[3])...

拆式子化简?
确实可以拆,然而完全不需要。

考虑到模数不为质数,将模数拆成若干个质数的次幂分开计算,最后使用中国剩余定理合并答案。
考虑到模数依旧不是质数而是质数的某个次幂,采用扩展lucas定理计算组合数。
那么就做完了。无比的暴力。

唯一需要注意的一点是,记得判Impossible……

#include<bits/stdc++.h>

using namespace std;

typedef long long ll;
typedef pair<ll,ll> pr;
const ll Inf=1e18;

ll p,n,m;
ll w[19];
pr ans[10009],ori[10009];

inline ll read()
{
    ll x=0;char ch=getchar();
    while(ch<'0' || '9'<ch)ch=getchar();
    while('0'<=ch && ch<='9')x=x*10+(ch^48),ch=getchar();
    return x;
}

inline ll qpow(ll a,ll b,ll p=Inf)
{
    ll ret=1;
    while(b)
    {
        if(b&1)
            ret=ret*a%p;
        a=a*a%p;
        b>>=1;
    }
    return ret;
}

inline ll phi(ll p,ll t)
{
    return (p-1)*qpow(p,t-1);
}

inline ll inv(ll a,ll p,ll t,ll k)
{
    ll phis=phi(p,t);
    return qpow(a,phis-1,k);
}

inline ll fac(ll x,ll p,ll k)
{
    if(!x)return 1;
    ll ret=1;
    for(ll i=1;i<=k;i++)
        if(i%p)
            ret=ret*i%k;
    ret=qpow(ret,x/k,k);
    for(ll i=x/k*k+1ll;i<=x;i++)
        if(i%p)
            ret=ret*i%k;
    return ret*fac(x/p,p,k)%k;
}

inline ll c(ll n,ll m,ll p,ll t)
{
    ll cnt=0,k=qpow(p,t);

    for(ll i=n;i;i=i/p)cnt+=i/p;
    for(ll i=m;i;i=i/p)cnt-=i/p;
    for(ll i=n-m;i;i=i/p)cnt-=i/p;

    ll ret=qpow(p,cnt,k)*fac(n,p,k)%k;
    ret=ret*inv(fac(m,p,k),p,t,k)%k;
    ret=ret*inv(fac(n-m,p,k),p,t,k)%k;
    return ret;
}

inline ll calc(ll p,ll t)
{
    ll ret=1,tot=n,k=qpow(p,t);
    for(int pos=1;pos<=m;pos++)
    {
        ret=ret*c(tot,w[pos],p,t)%k;
        tot-=w[pos];
    }
    return ret;
}

pr invv(ll x,ll y)
{
    if(!y)return pr(1,0);
    pr tmp=invv(y,x%y);
    return pr(tmp.second,tmp.first-x/y*tmp.second);
}

inline ll crt(pr *aa,int n)
{
    ll m=p,ret=0;
    for(int i=1;i<=n;i++)
    {
        ll x=aa[i].first,a=aa[i].second;
        ll mi=m/x,invs=(invv(mi,x).first%m+m)%m;
        ret=(ret+a*mi*invs)%m;
    }
    return ret;
}

inline ll work(ll p)
{
    int top=0;
    for(ll i=2;i*i<=p;i++)
    {
        if(p%i==0)
        {
            int cnt=0;
            while(p%i==0)
                p/=i,cnt++;
            ans[++top]=pr(qpow(i,cnt),calc(i,cnt));
        }
    }
    if(p!=1)ans[++top]=pr(p,calc(p,1));
    return crt(ans,top);
}

int main()
{
    p=read();n=read();m=read();
    ll sum=0;
    for(int i=1;i<=m;i++)
        w[i]=read(),sum+=w[i];
    if(sum>n)
    {
        puts("Impossible");
        return 0;
    }
    if(n>sum)w[++m]=n-sum;
    printf("%lld\n",work(p));
    return 0;
}


2020 jeepxie.net webmaster#jeepxie.net
10 q. 0.009 s.
京ICP备10005923号