题目地址：洛谷：【P2183】[国家集训队]礼物 – 洛谷、BZOJ：Problem 2142. — 礼物

题目描述

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

输入输出格式

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

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

输入输出样例

输入样例#1：

100
4 2
1
2

输出样例#1：

12

输入样例#2：

100
2 2
1
2

输出样例#2：

Impossible

说明

【样例说明】
下面是对样例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为质数。
对于15%的数据，n≤15，m≤5，pi^ci≤10^5；
在剩下的85%数据中，约有60%的数据满足t≤2，ci=1，pi≤10^5，约有30%的数据满足pi≤200。
对于100%的数据，1≤n≤10^9，1≤m≤5，1≤pi^ci≤10^5，1≤P≤10^9。

题解

本题需要的数学姿势有：数学笔记：数论（Lucas、CRT） | KSkun’s Blog
参考资料：题解 P2183 【[国家集训队]礼物】 – 没名字可被用的博客 – 洛谷博客
无解仅当\sum_{i=1}^m w_i < n。
有解的情况下，要求的就是\mathrm{C}_n^{w_1} \times \mathrm{C}_{n-w_1}^{w_2} \times \cdots \bmod P，对于每个组合数，应用扩展Lucas定理求解即可。

代码

// Code by KSkun, 2018/4
#include <cstdio>

typedef long long LL;

struct Num {
    LL num, pcnt; // pcnt: p因子数量，num: 提取p因子后乘起来的积
    Num(LL num = 0, LL pcnt = 0): num(num), pcnt(pcnt) {}
};

inline char fgc() {
    static char buf[100000], *p1 = buf, *p2 = buf;
    return p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 100000, stdin), p1 == p2) ? EOF 
        : *p1++;
}

inline LL readint() {
    register LL res = 0, neg = 1;
    char c = fgc();
    while(c < '0' || c > '9') {
        if(c == '-') neg = -1;
        c = fgc();
    }
    while(c >= '0' && c <= '9') {
        res = res * 10 + c - '0';
        c = fgc();
    }
    return res * neg;
}

const int MAXN = 100005, MAXM = 10;
const LL INF = 1e15;

inline LL fpow(LL n, LL k, LL p) {
    LL res = 1; n %= p;
    while(k) {
        if(k & 1) res = res * n % p;
        n = n * n % p;
        k >>= 1;
    }
    return res;
}

LL P, n, m, w[MAXM];

// calculate prime divisors in x
LL divi[MAXN], dcnt[MAXN], dpow[MAXN], dtot;

inline void calDivisor(LL x) {
    dtot = 0;
    for(LL i = 2; i * i <= P; i++) {
        if(x % i == 0) {
            divi[++dtot] = i;
            while(x % i == 0) {
                x /= i; dcnt[dtot]++;
            }
            dpow[dtot] = fpow(i, dcnt[dtot], INF);
        }
    }
    if(x != 1) {
        dtot++; 
        divi[dtot] = dpow[dtot] = x; dcnt[dtot] = 1;
    }
}

// calculate x! mod p^k
inline Num calFactorial(LL x, LL id) {
    if(x == 1 || x == 0) return Num(1, 0); // 1! = 0! = 1
    LL pk = dpow[id], res = 1;
    LL pnum = 0;
    for(LL i = x; i; i /= divi[id]) pnum += i / divi[id]; // 计算p因子的数量
    Num nxt = calFactorial(x / divi[id], id); // 递归计算提取了一个p因子的p倍数组成的阶乘
    res = res * nxt.num % pk; // 合并答案
    if(x >= pk) { // 如果有多段超过p^k的，预处理应用快速幂
        LL fac = 1;
        for(LL i = 1; i < pk; i++) {
            if(i % divi[id] == 0) continue;
            fac = fac * i % pk;
        }
        res = res * fpow(fac, x / pk, pk) % pk;
    }
    for(LL i = x; i >= 1; i--) { // 非整段p^k，暴力求解
        if(i % pk == 0) break;
        if(i % divi[id] == 0) continue;
        res = res * i % pk;
    }
    return Num(res, pnum);
}

inline LL exgcd(LL a, LL b, LL &x, LL &y) {
    if(!b) {
        x = 1; y = 0; return a;
    }
    LL res = exgcd(b, a % b, x, y);
    LL t = x; x = y; y = t - a / b * y;
    return res;
}

inline LL calInverse(LL n, LL p) {
    LL x, y;
    exgcd(n, p, x, y);
    return (x % p + p) % p;
}

LL crtm[MAXN];

// CRT求解组合数
inline LL crt() {
    LL res = 0;
    for(int i = 1; i <= dtot; i++) {
        res = (res + crtm[i] * P / dpow[i] % P * calInverse(P / dpow[i], dpow[i]) % P) % P;
    }
    return res;
}

int main() {
    P = readint(); n = readint(); m = readint();
    LL wsum = 0;
    for(int i = 1; i <= m; i++) {
        w[i] = readint();
        wsum += w[i];
    }
    if(wsum > n) {
        puts("Impossible"); return 0;
    }
    calDivisor(P);
    LL ans = 1;
    for(int i = 1; i <= m; i++) {
        for(int j = 1; j <= dtot; j++) {
            Num N = calFactorial(n, j), 
                M = calFactorial(w[i], j), 
                NM = calFactorial(n - w[i], j);
            // 分别代表n!、m!、(n-m)!
            N.pcnt -= M.pcnt + NM.pcnt;
            if(N.pcnt >= dcnt[j]) { // 如果p因子更多，说明能整除，模意义下为0
                crtm[j] = 0;
            } else {
                crtm[j] = N.num * calInverse(M.num, dpow[j]) % dpow[j] 
                    * calInverse(NM.num, dpow[j]) % dpow[j] 
                    * fpow(divi[j], N.pcnt, dpow[j]) % dpow[j]; // 把p因子乘进去
            }
        }
        ans = ans * crt() % P;
        n -= w[i];
    }
    printf("%lld", ans);
    return 0;
}