「ExLucas」学习笔记

2022-10-28 20:07:49

「ExLucas」学习笔记

前置芝士

  • 中国剩余定理 \(CRT\)

  • \(Lucas\) 定理

  • \(ExGCD\)

  • 亿点点数学知识

给龙蝶打波广告

Lucas 定理

\(C^m_n = C^{m\% mod}_{n\% mod} \times C^{\frac{m}{mod}}_{\frac{n}{mod}}\)

适用条件

  • 给出的数据范围较大（无法用线性求出）

  • 模数很烂的时候（会使阶乘中出现 \(0\)）

  • \(mod\) 必须为质数

证明

证明很恶心，略。

模板

某谷P4720

#include <cstdio>

#include <cstring>

#include <iostream>

#include <algorithm>

#define int long long

#define DEBUG puts ("emmmm");

using namespace std;

const int maxn = 1e5 + 50, INF = 0x3f3f3f3f;

inline int read () {

	register int x = 0, w = 1;

	char ch = getchar();

	for (; ch < '0' || ch > '9'; ch = getchar()) if (ch == '-') w = -1;

	for (; ch >= '0' && ch <= '9'; ch = getchar()) x = x * 10 + ch - '0';

	return x * w;

}

int T, n, m, mod;

int jc[maxn];

inline void Init () {

	jc[1] = 1;

	for (register int i = 2; i <= n + m; i ++) {

		jc[i] = jc[i - 1] * i % mod;

	}

}

inline int qpow (register int a, register int b) {

	register int ans = 1;

	while (b) {

		if (b & 1) ans = ans * a % mod;

		a = a * a % mod;

		b >>= 1;

	}

	return ans;

}

inline int C (register int a, register int b) {

	if (a == 0 || b == 0 || a == b) return 1;

	if (a < b) return 0;

	return jc[a] * qpow (jc[a - b], mod - 2) % mod * qpow (jc[b], mod - 2) % mod;

}

inline int Lucas (register int a, register int b) {

	if (a == 0 || b == 0) return 1;

	return C (a % mod, b % mod) * Lucas (a / mod, b / mod) % mod;

}

signed main () {

	T = read();

	while (T --) {

		n = read(), m = read(), mod = read();

		Init ();

		printf ("%lld\n", Lucas (n + m, n));

	}

	return 0;

}

扩展 Lucas

若题目中给出的 \(mod\) 不能保证是质数，当我们在求的时候，还是会出现 \(0\) 的情况，\(ExLuacs\) 就是来解决这种问题的。

STEP1

对于一个非质数 \(p\)，我们可以将其进行质因数分解，化成 \(\prod_ip_i^{k_i}\) 的形式。

我们就可以将原式子 \(C^m_n(mod \; p)\) 化成若干个同余方程：

\(\left\{\begin{matrix}
C^m_n \equiv b_1 (mod \; p_1^{k_1})\\
C^m_n \equiv b_2 (mod \; p_2^{k_2})\\
C^m_n \equiv b_3 (mod \; p_3^{k_3})\\
......\\
C^m_n \equiv b_i (mod \; p_i^{k_i})
\end{matrix}\right.\)

这样最后用 \(CRT\) 求出 \(C^m_n\) 即可。

STEP2

  • 现在问题变成了如何求每个 \(b_i\) 。

\(b_i = C^m_n (mod \; p_i ^ {k_i}) = \frac{n!}{m! \times (n - m)!} (mod \; p_i ^ {k_i})\)

但是我们会发现 \(p_i ^ {k_i}\) 仍不是质数， \(m!\) 和 \((n - m)!\) 的逆元仍求不出来。

所以我们将 \(n!\) 和 \(m!\) 和 \((n - m)!\) 中的所有质因子 \(p_i\) 都提出来，化成：

\(\frac{\frac{n!}{p_i^{k_1}}}{\frac{m!}{p_i^{k_2}} \times \frac{(n - m)!}{p_i^{k_3}}} \times p_i^{k_1-k_2-k_3}\)

这样分母上的就可以求出逆元了。

STEP3

  • 现在问题变成了如何求每个 \(\frac{n!}{p_i^{k_1}}\)

举个栗子！！

例如 \(n=22,p=3,k=2\)

\(n!=22\times 21\times 20\times 19\times 18\times 17\times 16\times 15\times 14\times 13\times 12\times 11\times 10\times 9\times 8\times 7\times 6\times 5\times 4\times 3\times 2\times 1\)

\(=3^7\times(1\times 2\times 3\times 4\times 5\times 6\times 7) \times (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)\)

我们会发现这个式子由三部分组成：

  • \(3^7\) 为 \(p^{\frac{n!}{p}}\)

  • \(7!\) 可以继续递归下去求解

  • 可以看出是在 \((mod \; 9)\) 意义下是一个循环节，长度为 \(\frac{n}{p_i^{k_i}}\)，类似 \(19\times 20\times 22\) 这样剩下的直接暴力求即可。

但是我们会发现第一部分会被原式子的分母消掉，所以不用计算，对于剩下的包含质因子 \(p_i\) 的，直接不计算即可。

inline int Calc (register int n, register int p, register int pk) {

	if (n == 0) return 1;

	register int ans = 1;

	for (register int i = 1; i <= pk; i ++) { // 每个循环节

		if (i % p) ans = ans * i % pk;

	}

	ans = qpow (ans, n / pk, pk); // 计算所有的循环节

	for (register int i = 1; i <= n % pk; i ++) { // 乘下剩下的

		if (i % p) ans = ans * i % pk;

	}

	return ans * Calc (n / p, p, pk) % pk;

}

最后

现在我们已经将所有要用的东西都求出来了，最后直接倒着退回去即可。

代码

某谷P4720

#include <cstdio>

#include <cstring>

#include <iostream>

#include <algorithm>

#include <cmath>

#define int long long

#define DEBUG puts ("emmmm")

const int maxn = 1e5 + 50, INF = 0x3f3f3f3f;

using namespace std;

inline int read () {

	register int x = 0, w = 1;

	char ch = getchar ();

	for (; ch < '0' || ch > '9'; ch = getchar ()) if (ch == '-') w = -1;

	for (; ch >= '0' && ch <= '9'; ch = getchar ()) x = x * 10 + ch - '0';

	return x * w;

}

int n, m, p, tot;

int b[maxn], c[maxn], d[maxn];

inline int qpow (register int a, register int b, register int mod) {

	register int ans = 1;

	while (b) {

		if (b & 1) ans = ans * a % mod;

		a = a * a % mod;

		b >>= 1;

	}

	return ans;

}

inline int ExGCD (register int a, register int b, register int &x, register int &y) {

	if (b == 0) {

		x = 1, y = 0;

		return a;

	}

	register int d = ExGCD (b, a % b, x, y);

	register int tmp = x;

	x = y;

	y = tmp - (a / b) * y;

	return d;

}

inline int Inv (register int a, register int mod) { // 利用扩展欧几里德求逆元

	register int x = 0, y = 0;

	ExGCD (a, mod, x, y);

	return (x % mod + mod) % mod;

}

inline int Calc (register int n, register int p, register int pk) {

	if (n == 0) return 1;

	register int ans = 1;

	for (register int i = 1; i <= pk; i ++) { // 每个循环节

		if (i % p) ans = ans * i % pk;

	}

	ans = qpow (ans, n / pk, pk); // 计算所有的循环节

	for (register int i = 1; i <= n % pk; i ++) { // 乘下剩下的

		if (i % p) ans = ans * i % pk;

	}

	return ans * Calc (n / p, p, pk) % pk;

}

inline int C (register int n, register int m, register int p, register int pk) {

	if (n == 0 || m == 0 || n == m) return 1;

	if (n < m) return 0;

	register int nn = Calc (n, p, pk), mm = Calc (m, p, pk), nm = Calc (n - m, p, pk), cnt = 0, k = n - m;

	while (n) n /= p, cnt += n;

	while (m) m /= p, cnt -= m;

	while (k) k /= p, cnt -= k;

	return nn * Inv (mm, pk) % pk * Inv (nm, pk) % pk * qpow (p, cnt, pk) % pk;

}

inline int CRT () { // 中国剩余定理

	register int M = 1, ans = 0;

	for (register int i = 1; i <= tot; i ++) {

		M *= c[i];

	}

	for (register int i = 1; i <= tot; i ++) {

		d[i] = Inv (M / c[i], c[i]);

	}

	for (register int i = 1; i <= tot; i ++) {

		ans += b[i] * (M / c[i])  * d[i];

	}

	return (ans % M + M) % M;

}

inline void ExLucas (register int n, register int m, register int p) {

	register int tmp = sqrt (p);

	for (register int i = 2; i <= tmp && p >= 1; i ++) { // 将p拆分质因数

		register int pk = 1;

		while (p % i == 0) p /= i, pk *= i;

		if (pk > 1) {

			b[++ tot] = C (n, m, i, pk), c[tot] = pk;

		}

	}

	if (p > 1) b[++ tot] = C (n, m, p, p), c[tot] = p;

	printf ("%lld\n", CRT ());

}

signed main () {

	n = read(), m = read(), p = read();

	ExLucas (n, m, p);

	return 0;

}

例题

[国家集训队]礼物

某谷P2183

思路很简单，就是没取一个 \(w[i]\)，总数就得减小，依次用 \(ExLucas\) 求组合数即可。

代码

#include <cstdio>

#include <cstring>

#include <iostream>

#include <algorithm>

#include <cmath>

#define int long long

#define DEBUG puts ("emmmm")

const int maxn = 1e5 + 50, INF = 0x3f3f3f3f;

using namespace std;

inline int read () {

	register int x = 0, w = 1;

	char ch = getchar ();

	for (; ch < '0' || ch > '9'; ch = getchar ()) if (ch == '-') w = -1;

	for (; ch >= '0' && ch <= '9'; ch = getchar ()) x = x * 10 + ch - '0';

	return x * w;

}

int n, m, p, totw, tot, ans = 1;

int w[maxn];

int b[maxn], c[maxn], d[maxn];

inline void Init () {

	memset (b, 0, sizeof b);

	memset (c, 0, sizeof c);

	memset (d, 0, sizeof d);

	tot = 0;

}

inline int qpow (register int a, register int b, register int mod) {

	register int ans = 1;

	while (b) {

		if (b & 1) ans = ans * a % mod;

		a = a * a % mod;

		b >>= 1;

	}

	return ans;

}

inline int ExGCD (register int a, register int b, register int &x, register int &y) {

	if (b == 0) {

		x = 1, y = 0;

		return a;

	}

	register int d = ExGCD (b, a % b, x, y);

	register int tmp = x;

	x = y;

	y = tmp - (a / b) * y;

	return d;

}

inline int Inv (register int a, register int mod) {

	register int x = 0, y = 0;

	ExGCD (a, mod, x, y);

	return (x % mod + mod) % mod;

}

inline int Calc (register int n, register int p, register int pk) {

	if (n == 0) return 1;

	register int ans = 1;

	for (register int i = 1; i <= pk; i ++) {

		if (i % p) ans = ans * i % pk;

	}

	ans = qpow (ans, n / pk, pk);

	for (register int i = 1; i <= n % pk; i ++) {

		if (i % p) ans = ans * i % pk;

	}

	return ans * Calc (n / p, p, pk) % pk;

}

inline int C (register int n, register int m, register int p, register int pk) {

	if (n == 0 || m == 0 || n == m) return 1;

	if (n < m) return 0;

	register int nn = Calc (n, p, pk), mm = Calc (m, p, pk), nm = Calc (n - m, p, pk), cnt = 0, k = n - m;

	while (n) n /= p, cnt += n;

	while (m) m /= p, cnt -= m;

	while (k) k /= p, cnt -= k;

	return nn * Inv (mm, pk) % pk * Inv (nm, pk) % pk * qpow (p, cnt, pk) % pk;

}

inline int CRT () {

	register int M = 1, ans = 0;

	for (register int i = 1; i <= tot; i ++) {

		M *= c[i];

	}

	for (register int i = 1; i <= tot; i ++) {

		d[i] = Inv (M / c[i], c[i]);

	}

	for (register int i = 1; i <= tot; i ++) {

		ans += b[i] * (M / c[i]) * d[i];

	}

	return (ans % M + M) % M;

}

inline int ExLucas (register int n, register int m, register int p) {

	Init ();

	register int tmp = sqrt (p);

	for (register int i = 2; i <= tmp && p > 1; i ++) {

		register int pk = 1;

		while (p % i == 0) p /= i, pk *= i;

		b[++ tot] = C (n, m, i, pk);

		c[tot] = pk;

	}

	if (p > 1) {

		b[++ tot] = C (n, m, p, p);

		c[tot] = p;

	}

	return CRT ();

}

signed main () {

	p = read(), n = read(), m = read();

	for (register int i = 1; i <= m; i ++) {

		w[i] = read();

		totw += w[i];

	}

	if (totw > n) {

		puts ("Impossible");

	} else {

		register int sum = n;

		for (register int i = 1; i <= m; i ++) {

			ans = (ans * ExLucas (sum, w[i], p)) % p;

			sum -= w[i];

		}

		printf ("%lld\n", ans);

	}

	return 0;

}
上一篇：Java for LeetCode 116 Populating Next Right Pointers in Each Node


下一篇：[LeetCode] 117. Populating Next Right Pointers in Each Node II 每个节点的右向指针 II