思路与算法

可行时方案数为n!/(r!∏wi!)，r=n−Σwi。把P分解为互质素数幂q=p^c。每个q下，分别计算阶乘中p的次数和去除p因子的单位部分；单位部分用长度q的前缀乘积表与递归n/p得到，再乘逆元和p的指数差。最后CRT合并各q余数。

正确性说明

依次选择每人的礼物是不同元素的不相交子集，剩余r件不送，因此多项式系数精确计数。对每个p，勒让德和给出p指数；去p阶乘把非p倍数按q整块计积，而p倍数去掉一个p后递归为(n/p)!的同类部分。分母单位与q互质，可求逆，乘回p指数差恰恢复该整数系数模q；CRT唯一恢复模P余数。Σwi>n无方案，P1可行时余数0。

复杂度

试除P为O(sqrt P)，每个素幂表初始化O(q)，每个阶乘去p递归O(log_p n)层、每层模幂O(log n)，总O(sqrt P+Σq+(m+1)Σlog_p n·log n)，空间O(max q+m)。

易错点

礼物不同，不能只数人数分配；允许剩余。不能在合数模数下直接逆n!，需分开p因子。模数总体上界P1e9为公开修订，原每素幂≤1e5仍保留；Impossible区分可行但余数0。

参考代码（C++11）

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
ll powmod(ll a, ll b, ll m) {
    ll z = 1 % m;
    for (; b; b >>= 1, a = a * a % m)
        if (b & 1)
            z = z * a % m;
    return z;
}
ll eg(ll a, ll b, ll &x, ll &y) {
    if (!b) {
        x = 1;
        y = 0;
        return a;
    }
    ll u, v, g = eg(b, a % b, u, v);
    x = v;
    y = u - a / b * v;
    return g;
}
ll inv(ll a, ll m) {
    ll x, y;
    eg(a, m, x, y);
    return (x % m + m) % m;
}
int main() {
    ll P, n;
    int m;
    cin >> P >> n >> m;
    vector<ll> w(m);
    ll sum = 0;
    for (auto &x : w) {
        cin >> x;
        sum += x;
    }
    if (sum > n) {
        cout << "Impossible\n";
        return 0;
    }
    if (P == 1) {
        cout << 0 << '\n';
        return 0;
    }
    w.push_back(n - sum);
    ll tmp = P, ans = 0;
    vector<pair<ll, ll>> factors;
    for (ll p = 2; p * p <= tmp; p++)
        if (tmp % p == 0) {
            ll q = 1;
            while (tmp % p == 0) {
                tmp /= p;
                q *= p;
            }
            factors.push_back({p, q});
        }
    if (tmp > 1)
        factors.push_back({tmp, tmp});
    for (auto item : factors) {
        ll p = item.first, q = item.second;
        vector<ll> f(q + 1, 1);
        for (ll i = 1; i <= q; i++)
            f[i] = i % p ? f[i - 1] * i % q : f[i - 1];
        function<ll(ll)> unit = [&](ll a) -> ll {
            if (!a)
                return 1;
            return powmod(f[q], a / q, q) * f[a % q] % q * unit(a / p) % q;
        };
        auto val = [&](ll a) {
            ll z = 0;
            while (a) {
                a /= p;
                z += a;
            }
            return z;
        };
        ll e = val(n), z = unit(n);
        for (ll a : w) {
            e -= val(a);
            z = z * inv(unit(a), q) % q;
        }
        z = z * powmod(p, e, q) % q;
        ll part = P / q;
        ans = (ans + z * part % P * inv(part % q, q)) % P;
    }
    cout << ans << '\n';
}