「ExLucas」學習筆記
阿新 • • 發佈：2020-10-02
# 「ExLucas」學習筆記
## 前置芝士
* 中國剩餘定理 $CRT$
* $Lucas$ 定理
* $ExGCD$
* ~~億點點~~數學知識
[給龍蝶打波廣告](https://www.cnblogs.com/DZN2004/p/13663114.html)
## Lucas 定理
$C^m_n = C^{m\% mod}_{n\% mod} \times C^{\frac{m}{mod}}_{\frac{n}{mod}}$
### 適用條件
* 給出的資料範圍較大（無法用線性求出）
* 模數很爛的時候（會使階乘中出現 $0$）
* $mod$ 必須為質數
### 證明
[被某人強迫打的廣告](https://www.cnblogs.com/rui-4825/p/13340278.html)
### 模板
[某谷P3807](https://www.luogu.com.cn/problem/P3807)
```c++
#include
#include
#include
#include
#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$ 的，直接不計算即可。
```c++
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](https://www.luogu.com.cn/problem/P4720)
```c++
#include
#include
#include
#include
#include
#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](https://www.luogu.com.cn/problem/P2183)
思路很簡單，就是沒取一個 $w[i]$，總數就得減小，依次用 $ExLucas$ 求組合數即可。
#### 程式碼
```c++
#include
#include
#include
#include
#include
#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;
}
```