BZOJ 2142 礼物 组合数学 CRT 中国剩余定理

2142: 礼物

Time Limit: 10 Sec  Memory Limit: 259 MB
Submit: 1450  Solved: 593
[Submit][Status][Discuss]

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


【样例说明】
下面是对样例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。

HINT

Source

Solution

根据题目大意可以列出式子C(n,sum)*C(sum, w[1])*C(sum-w[1], w[2])··········

如果p是质数的话，显然就是Lucas的裸题。

但现在p不是质数了，要怎么办呢？把p分解成一个个pi^ci，这样的话就保证了两两互质，CRT求解。

那怎么算C(x,y)%(pi^ci)呢？把C(x,y)转化成x!、y!和(x-y)!。

由于需要求逆元，而n!不一定与pi^ci互质，我们就需要把n!分解为k*(pi^x)，其中，k与pi互质。

对于当前的n!，我们发现可以把它整理为一些模(pi^ci)下的循环串乘上(pi^(n/pi))再乘上(n/pi)!

举个栗子，17!%(3^2) = (1*2*4*5*7*8*10*11*13*14*16*17)*(3^5)*(1*2*3*4*5)

很显然，栗子最后的那部分就是5!%3这个显然可以递归下去。

对于最前的那部分，我们可以发现，它是由一些循环节（模3^2下）组成的，且均不含有质因数3。

仔细分析，就是[1,3^2]中与质因数3互质的数。求出循环节之后就可以直接用快速幂完成。

对于中间的部分，我们收集起来，算组合数的时候统一再做一次快速幂。

本题主要的问题是处理好逆元的问题，把n!分解成k*(pi^x)，含pi的部分单独算，k与pi^ci互质，可以直接求解逆元。

Code

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <string>
#include <algorithm>

using namespace std;

#define REP(i, a, b) for (int i = (a), i##_end_ = (b); i <= i##_end_; ++i)
#define fi first
#define se second
#define mp make_pair
typedef long long LL;
typedef pair<LL, LL> pa;
const int maxn = 1e5;
LL p, a[105];
LL prime[maxn+10], p_cnt;
bool isNotPrime[maxn+10];
LL p_num[105], p_mod[105], p_sum[105], cnt;

void prepare()
{
    REP(i, 2, maxn)
    {
        if (!isNotPrime[i]) prime[++p_cnt] = i;
        REP(j, 1, p_cnt)
        {
            if (i*prime[j] > maxn) break ;
            isNotPrime[i*prime[j]] = 1;
            if (i%prime[j] == 0) break ; 
        }
    }
    LL temp = p;
    REP(i, 1, p_cnt)
    {
        if (prime[i] > temp || temp == 1) break ;
        if (temp%prime[i] == 0)
        {
            p_num[++cnt] = prime[i], p_mod[cnt] = 1, p_sum[cnt] = 0;
            while (temp%prime[i] == 0)
                p_mod[cnt] *= prime[i], p_sum[cnt] ++, temp /= prime[i];
        }
    }
}

LL power(LL x, LL y, LL MOD)
{
    LL ret = 1;
    while (y > 0)
    {
        if (y&1) ret = (ret*x)%MOD;
        x = (x*x)%MOD;
        y >>= 1;
    }
    return ret;
}

pa calc(int k, LL num)
{
    if (num == 0) return mp(0, 1);
    LL x = num/p_num[k], y = num/p_mod[k];
    LL ans = 1;
    if (y)
    {
        LL t_num = 1;
        REP(i, 1, p_mod[k]-1)
            if (i%p_num[k] != 0) t_num = (t_num*i)%p_mod[k];
        ans = power(t_num, y, p_mod[k]);
    }
    if (num%p_mod[k] != 0)
    {
        REP(i, 1, num%p_mod[k])
            if (i%p_num[k] != 0) ans = (ans*i)%p_mod[k];
    }
    pa temp = calc(k, x);
    return mp(temp.fi+x, temp.se*ans%p_mod[k]);
}

LL ex_gcd(LL aa, LL bb, LL &x, LL &y)
{
    if (bb == 0)
    {
        x = 1, y = 0;
        return aa;
    }
    LL ret = ex_gcd(bb, aa%bb, x, y);
    LL temp = x;
    x = y, y = temp-(aa/bb)*y;
    return ret;
}

LL inv(LL k, LL MOD)
{
    LL x, y;
    ex_gcd(k, MOD, x, y);
    x = (x%MOD+MOD)%MOD;
    return x;
}

LL CRT()
{
    LL ret = 0;
    REP(i, 1, cnt)
    {
        LL Mi = p/p_mod[i];
        ret = (ret+((Mi*inv(Mi, p_mod[i]))%p*a[i])%p)%p;
    }
    ret = (ret%p+p)%p;
    return ret;
}

LL work(LL x, LL y)
{
    REP(i, 1, cnt)
    {
        pa t1 = calc(i, x), t2 = calc(i, y), t3 = calc(i, x-y);
        a[i] = power(p_num[i], t1.fi-t2.fi-t3.fi, p_mod[i]);
        a[i] = (((a[i]*t1.se)%p_mod[i]*inv(t2.se, p_mod[i]))%p_mod[i]*inv(t3.se, p_mod[i]))%p_mod[i];
    }
    return CRT();
}

int main()
{
//    freopen("a.in", "r", stdin);
//    freopen("a.out", "w", stdout);
    LL n, m, w[15], sum = 0;
    scanf("%lld %lld %lld", &p, &n, &m);
    prepare();
    REP(i, 1, m) scanf("%lld", &w[i]), sum += w[i];
    if (sum > n) { puts("Impossible"); return 0; }
    LL ans = work(n, sum);
    REP(i, 1, m)
        ans = (ans*work(sum, w[i]))%p, sum -= w[i];//, printf("%lld
", ans);
    ans = (ans+p)%p;
    printf("%lld
", ans);
    return 0;
}