[bzoj 2142] 礼物

n 件有区别的物品给 m 个人, 第 i 个人 wi 件, 求方案数模 P, 或判断无解. (1≤n≤10^9, 1≤m≤5, 1≤pi^ci≤105)

题面没说P的范围. 可以认为是long long.

可重集的排列. 答案可以写成阶乘形式, 也可以用组合数表示. 组合数对素数取模有 Lucas 定理. 这里模不一定是素数, 所以分解成素数的幂之积, 分别求, 再用中国剩余定理合并答案.

然而 Lucas 定理不能处理模是素数的幂 (次数大于1) 的情形......有没有推广呢? 上网找了一下, 有所谓的扩展 Lucas 定理, 我不觉得它和 Lucas 定理有什么关联. 但方法还是不错的.

组合数对素数的幂取模

组合数取模麻烦的地方在于, 如果先对分子分母取模, 可能面临它们同时为 的情形. 但是, 如果模是素数的幂 , 我们知道只有 的倍数会贡献因子 . 把 提出来, 分子分母的次数作差. 如果差不小于 , 则答案为 ; 否则, 把约分后的 乘到分子上, 分母求逆元, 乘积即为所求.

考虑求 对 取模的结果, 并表示成 . 从 乘到 并取模是有循环节的: 它模等价于 . 预处理 表示 ( 除以不大于 的 的倍数) 模 . 把 从 中提出来, 则它们的乘积等于 . 是一个子问题, 递归下去求解. 的倍数提出来之后, 剩下的部分可以直接从 数组中取, 并配合快速幂.

设 fac(n, p, k) 返回二元组 (x, e), 伪代码为: fac(n, p, k) t = fac(n/p, p, k) return (t.x * f[p^k]^(n / p^k) * f[n mod p^k] mod p^k, t.e + n/p)

UPDATE Dec 24, 2018

高斯泛化的威尔逊定理

去年省选的时候这样写, 发现还是跑得太慢......去学习了 Tangjz 的 组合数取模, 发现:

Gauss's generalization of Wilson's theorem

因此上述代码中的 f[p^k] 实际上只会是 : 当 时, ; 当 或 时取 ; 当 且 时取 .

研究 的解的性质即可证明, 详细过程过几天写.

以下代码并未加入此优化.

中国剩余定理

令 , 两两互质, 则以下对应关系是一一映射:

说明: 可对 取模, 但 只能对 取模.

代码

可以用 Python 对拍, 但是 n 只能取到几百, 否则跑不出来.

#include <bits/stdc++.h>
#define Rep(i, a, b) for (int i = (a), i##_ = (b); i < i##_; ++i)
#define For(i, a, b) for (int i = (a), i##_ = (b); i <= i##_; ++i)
#define Down(i, a, b) for (int i = (a), i##_ = (b); i >= i##_; --i)

using namespace std;
typedef long long ll;
typedef pair<int, int> ii;

int cnt, mod[65];
ll w[5];
ii c[65];

void decompose(ll n)
{
	for (ll i = 2; i*i <= n; ++i)
		if (n % i == 0) {
			mod[cnt] = 1;
			c[cnt].first = i;
			do {
				mod[cnt] *= i;
				++c[cnt].second;
				n /= i;
			} while (n % i == 0);
			++cnt;
		}
	if (n > 1) {
		c[cnt].first = mod[cnt] = n;
		c[cnt].second = 1;
		++cnt;
	}
}

ll f[100001];

void init(ll p, ll m)
{
	f[0] = 1;
	For (i, 1, m)
		f[i] = f[i-1] * (i % p ? i : 1) % m;
}

ll fpm(ll x, ll n, ll m)
{
	ll y = 1;
	while (n) {
		if (n & 1)
			(y *= x) %= m;
		(x *= x) %= m;
		n >>= 1;
	}
	return y;
}

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

inline ll inv(ll x, ll m)
{
	ll y, t;
	exgcd(x, m, y, t); // yx + tm = 1
	return (y + m) % m;
}

ll fac(ll n, ll p, ll m, int& cnt)
{
	int t;
	return n ? (t = n/p, cnt += t, fpm(f[m], n/m, m) * f[n%m] % m * fac(t, p, m, cnt) % m) : 1;
}

int ans[65];

int main()
{
	ll p, n, m, s = 0;
	scanf("%lld%lld%lld", &p, &n, &m);
	Rep (i, 0, m) {
		scanf("%lld", &w[i]);
		s += w[i];
	}
	if (s > n) {
		puts("Impossible");
		return 0;
	}
	decompose(p);
	Rep (i, 0, cnt) {
		int a = 0, b = 0;
		ll t = c[i].first, M = mod[i];
		init(t, M);
		ll x = fac(n, t, M, a), y = fac(n-s, t, M, b);
		Rep (j, 0, m)
			(y *= fac(w[j], t, M, b)) %= M;
		ans[i] = a-b < c[i].second ? x * inv(y, M) % M * fpm(t, a-b, M) % M : 0;
	}
	ll r = 0;
	Rep (i, 0, cnt) if (ans[i]) {
		ll u = p/mod[i];
		r += u * inv(u, mod[i]) % p * ans[i] % p;
		r -= r >= p ? p : 0;
	}
	printf("%lld\n", r);
	return 0;
}