BZOJ 2142 - 礼物 | Sengxian's Blog

BZOJ 2142 - 礼物

#数学 #中国剩余定理 #组合数

Published on 2016-06-17

描述

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

分析

如果 $\sum w_i > n$ 无解输出 Impossible。
否则我们推下公式：

展开后裂项可得：

ans = \frac {n!}{w_1!w_2!w_3!\cdots w_n!(n-w_1-w_2-\cdots - w_n)!}

为了方便处理，如果 $\sum w_i < n$ 的话，我们补一个 $n - \sum w_i$ 即可当作 $\sum w_i = n$ 处理。

本题中 $P$ 并不是质数，处理这种非质数取模的问题，一般采用中国剩余定理（Chinese remainder theorem，CRT）。
将 $P$ 分解为 $P = p_1^{c_1}p_2^{c_2}\cdots p_n^{c_n}$ 的形式，可以知道 $p_i^{c_i}$ 两两互质，所以我们只要知道答案模 $p_i^{c_i}$ 的值，即可通过中国剩余定理合并为模 $P$ 的值。

所以我们只用考虑答案模 $p_i^{c_i}$ 的值，这个数也可能不是质数，而且阶乘高达 ${10}^9$ 级别，是根本无法直接算的，这需要我们进一步挖掘阶乘的性质。

我们先考虑 $c = 1$ 即 $p$ 就是质数的情况。注意到，即使 $p$ 是质数，这个式子仍然无法直接算！因为分母中可能存在 $p$ 这个因子，而 $p$ 对 $p$ 是没有逆元的！
我们考虑将分子分母变为为 $\frac {p^a \cdot x} {p^b \cdot y} \equiv p^{a - b} \cdot x \cdot y^{-1} \pmod {p^c}$ ，也就是将阶乘 $n!$ 分解为 $n! = p^a \cdot q$ 的形式。

定理：如果 $\sum m_i = n$ ，则 $n!$ 中素数 $p$ 的指数不会小于 $\prod m_i!$ 中素数 $p$ 的指数。

证明：对于 $n!$ 的分解式子中，素数 $p$ 的指数为：

\left\lfloor \frac n p\right\rfloor + \left\lfloor \frac n {p^2}\right\rfloor + \left\lfloor \frac n {p^3}\right\rfloor + ...

而 $\prod m_i!$ 的分解式子中，素数 $p$ 的指数为：

\sum \left\lfloor \frac {m_i} p\right\rfloor + \sum\left\lfloor \frac {m_i} {p^2}\right\rfloor + \sum\left\lfloor \frac {m_i} {p^3}\right\rfloor + ...

由于是向下取整，所以 $n!$ 中 $p$ 的个数不会小于 $\prod m_i!$ 中 $p$ 的个数，得证。这说明 $p^{a - b}$ 一定是有意义的。

现在，我们需要对 $n!$ 快速分解，得到 $n! = p^{a} \cdot q$ ，其中 $q$ 由于不含 $p$ 这个因子，是可以模 $p$ 的。

我们发现 $n! \bmod p$ 的计算过程是以 $p$ 为周期的，举例如下，当 $n = 10, p = 3$ 时：

\begin{aligned}n! &= 1 \times 2 \times 3 \times 4 \times 5 \times 6 \times 7 \times 8 \times 9 \times 10\\ &=(1 \times 2 \times 3) \times (4\times5 \times 6)\times(7 \times 8 \times 9) \times 10\\ &= (1 \times 2) \times (4 \times 5) \times (7 \times 8) \times 10 \times 3^3 \times (1 \times 2 \times 3) \end{aligned}

最后一步中的 $1 \times 2 \times 3$ 可递归处理，所以我们预处理出 $[1, p]$ 的阶乘，对于每个 $n!$ 即可 $O(\log_{p}n)$ 求解。

对于分母中的 $n!$ ，不需任何处理，只需要将最后分解出的 $n! = p^{a}\cdot q$ 中的 $q$ 用扩展欧几里得求出逆元即可。

那么 $c_i \neq 1$ 怎么办？只需对原算法稍加修改即可。令 $\mathrm{fac}(i)$ 为除去 $p_i$ 的倍数外 $i$ 的阶乘。例如 $p_i = 3, c_i = 2$ ，那么 $\mathrm{fac}(10) = 1 \times 2 \times 4 \times 5 \times 7 \times 8 \times 10$ ，除去了 $3$ 的倍数 $3$ 、 $6$ 和 $9$ 。阶乘依然是以 $p_i^{c_i}$ 为周期的。在算法内部，统计因子 $p$ 的个数时，每次应该加 $\frac n {p_i}$ 而不是 $\frac n {p_i^{c_i}}$ ，递归处理时也应该是递归到 $\lfloor\frac n {p_i}\rfloor$ 。

详细参见代码。

代码

//  Created by Sengxian on 6/17/16.
//  Copyright (c) 2016年 Sengxian. All rights reserved.
//  BZOJ 2142 阶乘模任意数
#include <bits/stdc++.h>
using namespace std;

typedef long long ll;
inline int ReadInt() {
    static int n, ch;
    n = 0, ch = getchar();
    while (!isdigit(ch)) ch = getchar();
    while (isdigit(ch)) n = (n << 3) + (n << 1) + ch - '0', ch = getchar();
    return n;
}

void extgcd(ll a, ll b, ll &x, ll &y) {
    if (b == 0) x = 1, y = 0;
    else {
        extgcd(b, a % b, y, x);
        y -= (a / b) * x;
    }
}

ll china(int n, int *a, int *m) {
    ll M = 1, w, x, y, ans = 0;
    for (int i = 0; i < n; ++i) M *= m[i];
    for (int i = 0; i < n; ++i) {
        w = M / m[i];
        extgcd(w, m[i], x, y);
        (ans += w * x * a[i]) %= M;
    }
    return (ans + M) % M;
}

inline ll mod_pow(ll a, int b, int modu) {
    ll res = 1;
    while (b) {
        if (b & 1) (res *= a) %= modu;
        (a *= a) %= modu;
        b >>= 1;
    }
    return res;
}

const int maxn = 5 + 3, maxs = 6 + 3, maxp = 100000 + 3;
int n, p, m, w[maxn], a[maxs], mm[maxs];
int fact[maxp];

ll mod_fact(int n, const int &P, const int &p, ll &e) {
    if (n < p) return fact[n]; //注意不是 P 而是 p
    e += n / p;
    return mod_pow(fact[P - 1], n / P, P) * fact[n % P] % P * mod_fact(n / p, P, p, e) % P;
}

ll inv(ll a, ll p) {
    static ll x, y;
    extgcd(a, p, x, y);
    return (x % p + p) % p;
}

ll cal(int P, int p) {
    fact[0] = fact[1] = 1;
    for (int i = 2; i < P; ++i) fact[i] = (ll)fact[i - 1] * (i % p ? i : 1) % P;
    ll e1 = 0, e2 = 0;
    ll a = mod_fact(n, P, p, e1), b = 1;
    for (int i = 0; i < m; ++i) (b *= mod_fact(w[i], P, p, e2)) %= P;
    return a * inv(b, P) % P * mod_pow(p, e1 - e2, P) % P;
}


int main() {
#ifndef ONLINE_JUDGE
    freopen("test.in", "r", stdin);
#endif
    int sum = 0;
    p = ReadInt(), n = ReadInt(), m = ReadInt();
    for (int i = 0; i < m; ++i) sum += w[i] = ReadInt();
    if (sum > n) return puts("Impossible"), 0;
    else if (sum < n) w[m++] = n - sum;
    int cnt = 0;
    for (int i = 2; (ll)i * i <= p; ++i) if (p % i == 0) {
        int P = 1;
        while (p % i == 0) p /= i, P *= i;
        mm[cnt] = P, a[cnt++] = cal(P, i);
    }
    if (p > 1) mm[cnt] = p, a[cnt++] = cal(p, p);
    printf("%lld\n", china(cnt, a, mm));
    return 0;
}

BZOJ 1937 - [Shoi2004]Mst 最小生成树 BZOJ 2669 - [cqoi2012]局部极小值