【BZOJ2142】礼物（扩展Lucas）

Description

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

---

Solution

参考博客

容易看出答案是：
$$Ans=\left(\genfrac{}{}{0ex}{}{n}{w_1}\right)\left(\genfrac{}{}{0ex}{}{n-w_1}{w_2}\right)\left(\genfrac{}{}{0ex}{}{n-w_1-w_2}{w_3}\right)\dots$$

$p$没有保证是质数，所以考虑对于每一个组合数用扩展$Lucas$定理求出。

将$p$表示成${\prod }_i{p}_i^{k_i}$，我们可以对于所有不同的质因子$p_i$求出$\left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}{p}_i^{k_i}$，然后用CRT合并，所以现在主要是要求$\left(\genfrac{}{}{0ex}{}{n}{m}\right)\mathrm{m}\mathrm{o}\mathrm{d}{p}_i^{k_i}$。

我们首先只考虑对于$n!$的计算，$m!$、$(n-m)!$同理。
像网上大多数博客一样，举这个例子：
假设$n=22$，$p_i=3$，$k_i=2$
那么$n!=1\times 2\times \cdots \times 22$
然后将其中是3的倍数的数提出来：
$(1\times 2\times 4\times 5\times 7\times 8\times 10\times 11\times 13\times 14\times 16\times 17\times 19\times 20\times 22)\times 3^6\times (1\times 2\times 3\times 4\times 5\times 6\times 7)$
然后发现这个式子可以分成三部分：
1、$p_i^{k_i}$,这个可以直接快速幂
2、对于阶乘，我们可以递归求解
3、对于前面的缺了$3$的倍数的阶乘，我们发现以$p_i^{k_i}$为周期，它们是同余的！
即$(1\times 2\times 4\times 5\times 7\times 8)\equiv (10\times 11\times 13\times 14\times 16\times 17)(\mathrm{m}\mathrm{o}\mathrm{d}{3}^2)$
注意$p_i^{k_i}$要放在最后算，因为$m!$、$(n-m)!$如果包含$p_i^{k_i}$，就求不了逆元了。
注：计算$n!$中质因子$p_i$的个数公式为：$x={\sum }_{j=1}^{\infty }\lfloor \frac{n}{p_i^j}\rfloor$

---

Code

/************************************************
 * Au: Hany01
 * Date: May 27th, 2018
 * Prob: [BZOJ2142] 礼物
 * Email: hany01@foxmail.com
************************************************/

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;
typedef pair<int, int> PII;
#define File(a) freopen(a".in", "r", stdin), freopen(a".out", "w", stdout)
#define rep(i, j) for (register int i = 0, i##_end_ = (j); i < i##_end_; ++ i)
#define For(i, j, k) for (register int i = (j), i##_end_ = (k); i <= i##_end_; ++ i)
#define Fordown(i, j, k) for (register int i = (j), i##_end_ = (k); i >= i##_end_; -- i)
#define Set(a, b) memset(a, b, sizeof(a))
#define Cpy(a, b) memcpy(a, b, sizeof(a))
#define x first
#define y second
#define pb(a) push_back(a)
#define mp(a, b) make_pair(a, b)
#define ALL(a) (a).begin(), (a).end()
#define SZ(a) ((int)(a).size())
#define INF (0x3f3f3f3f)
#define INF1 (2139062143)
#define debug(...) fprintf(stderr, __VA_ARGS__)
#define y1 wozenmezhemecaia

template <typename T> inline bool chkmax(T &a, T b) { return a < b ? a = b, 1 : 0; }
template <typename T> inline bool chkmin(T &a, T b) { return b < a ? a = b, 1 : 0; }

int p;

inline int read()
{
	register int _, __; register char c_;
    for (_ = 0, __ = 1, c_ = getchar(); c_ < '0' || c_ > '9'; c_ = getchar()) if (c_ == '-') __ = -1;
    for ( ; c_ >= '0' && c_ <= '9'; c_ = getchar()) _ = (_ << 1) + (_ << 3) + (c_ ^ 48);
    return _ * __;
}

inline LL Pow(LL a, LL b, LL Mod)
{
	static int Ans;
	for (Ans = 1; b; b >>= 1, (a *= a) %= Mod)
		if (b & 1) (Ans *= a) %= Mod;
	return Ans;
}

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

LL fac(LL n, LL p, LL pk)
{
	if (!n) return 1;
	LL res = 1;
	for (register LL i = 2; i <= pk; ++ i)
		if (i % p) (res *= i) %= pk;
	res = Pow(res, n / pk, pk);
	for (register LL i = 2; i <= n % pk; ++ i)
		if (i % p) (res *= i) %= pk;
	return res * fac(n / p, p, pk) % pk;
}

inline LL inv(LL n, LL Mod)
{
	static LL x, y, t;
	gcd(n, Mod, x, y);
	t = ((x % Mod) + Mod) % Mod;
	return t;
}

inline LL C(LL n, LL m, LL p, LL k, LL pk)
{
	if (n < m) return 0;
	register LL t1 = fac(n, p, pk), t2 = fac(m, p, pk), t3 = fac(n - m, p, pk), cnt = 0;
	for (register LL i = n; i; i /= p) cnt += i / p;
	for (register LL i = m; i; i /= p) cnt -= i / p;
	for (register LL i = n - m; i; i /= p) cnt -= i / p;
	return t1 * inv(t2, pk) % pk * inv(t3, pk) % pk * Pow(p, cnt, pk) % pk;
}

inline LL CRT(LL c, LL m) { return c * inv(p / m, m) % p * (p / m) % p; }

inline LL exLucas(LL n, LL m)
{
	LL Ans = 0, tmp = p;
	for (register int i = 2; i * i <= tmp; ++ i)
		if (!(tmp % i)) {
			LL cnt = 0, prod = 1;
			while (!(tmp % i)) tmp /= i, prod *= i, ++ cnt;
			(Ans += CRT(C(n, m, i, cnt, prod), prod)) %= p;
		}
	if (tmp > 1) (Ans += CRT(C(n, m, tmp, 1, tmp), tmp)) %= p;
	return Ans;
}

int main()
{
#ifdef hany01
	File("bzoj2142");
#endif

	static int n, m, w[7], sumw = 0;
	static LL Ans = 1;

	p = read(), n = read(), m = read();
	For(i, 1, m) sumw += w[i] = read();

	if (sumw > n) { puts("Impossible"); return 0; }

	For(i, 1, m) (Ans *= exLucas(n, w[i])) %= p, n -= w[i];
	printf("%lld\n", Ans);

    return 0;
}
//重见金英人未见。相思一夜天涯远。
//    -- 晏几道《蝶恋花·黄菊开时伤聚散》