威尔逊/扩展卢卡斯简明教程
Dog_E
·
·
算法·理论
卢卡斯定理/威尔逊定理
威尔逊定理
前置知识：
1) 唯一分解定理
2) 线性同余方程
3) 乘法逆元
4) （扩展）中国剩余定理
定义
若 p 为质数,则有：
(p-1)!\equiv -1\pmod p.
证明
当 p=2 时，命题显然成立。
当 p\geq 3，此时要证 \mathbf{Z}_p 中所有非零元素的积为 \overline{-1}.(\mathbf{Z}_p 是 \mod p 下的剩余系)
但是要注意 $a$ 和 $a^{-1}$ 可能相等。若 $a=a^{-1}$，则 $a^2\equiv 1\pmod p$，即
$$
0\equiv a^2-1\equiv (a+1)(a-1)\pmod p
$$
从而 $a\equiv 1\pmod p$ 或 $a\equiv -1\pmod p$.
形象的说，仅有 $\mathbf{Z}_p$ 中的两个元素 $\{\overline{1},\overline{p-1}\}$ 和它们的逆元相等，进而 $\mathbf{Z}_p$ 中所有非零元素的积为 $\overline{-1}$.
### 推广
对于整数 $n$，令 $(n!)_p$ 表示所有小于等于 $n$ 但不能被 $p$ 整除的正整数的乘积，即：
$$(n!)_p=\frac{n!}{ \lfloor \frac{n}{p} \rfloor !p^{\lfloor \frac{n}{p} \rfloor}}$$
Wilson 定理指出 $(p!)_p=(p-1)!\equiv -1\pmod p$ 且可被推广至模素数 $p$ 的**幂次**。
该推广定理为接下来的`扩展卢卡斯定理`提供了理论基础。
### 应用
在某些情况下，出现以某个素数 $p$ 为模的公式，并且分**子分母中包含阶乘**，比如经常出现在**排列组合数计算**中。
只有当阶乘同时出现在分数的分子和分母中时，这个问题才有意义。否则，后续项 $p!$ 将减少为**零**。但在分数中，因子 $p$ 可以抵消，结果将是非零。
因此，要计算 $n! \bmod p$，而不考虑阶乘中出现因数 $p$。写下素因子分解，去掉所有因子 $p$，并计算乘积模。
用 $(n!)_p$ 表示这个修改的因子。例如：
$$
(7!)_p \equiv 1 \cdot 2 \cdot \underbrace{1}_{3} \cdot 4 \cdot 5 \underbrace{2}_{6} \cdot 7 \equiv 2 \bmod 3
$$
这种修正的阶乘，可用于快速计算各种带取模和组合数的公式的值。
#### 计算余数的算法
计算上述去掉因子 $p$ 的阶乘模 $p$ 的余数。
$$
\begin{aligned}
(n!)_p &= 1 \cdot 2 \cdot 3 \cdot \ldots \cdot (p-2) \cdot (p-1) \cdot \underbrace{1}_{p} \cdot (p+1) \cdot (p+2) \cdot \ldots \cdot (2p-1) \cdot \underbrace{2}_{2p} \\
&\quad \cdot (2p+1) \cdot \ldots \cdot (p^2-1) \cdot \underbrace{1}_{p^2} \cdot (p^2 +1) \cdot \ldots \cdot n \pmod{p} \\\\
&= 1 \cdot 2 \cdot 3 \cdot \ldots \cdot (p-2) \cdot (p-1) \cdot \underbrace{1}_{p} \cdot 1 \cdot 2 \cdot \ldots \cdot (p-1) \cdot \underbrace{2}_{2p} \cdot 1 \cdot 2 \\
&\quad \cdot \ldots \cdot (p-1) \cdot \underbrace{1}_{p^2} \cdot 1 \cdot 2 \cdot \ldots \cdot (n \bmod p) \pmod{p}
\end{aligned}
$$
可以清楚地看到，除了最后一个块外，阶乘被划分为几个长度相同的块。
$$
\begin{aligned}
(n!)_p&= \underbrace{1 \cdot 2 \cdot 3 \cdot \ldots \cdot (p-2) \cdot (p-1) \cdot 1}_{1\text{st}} \cdot \underbrace{1 \cdot 2 \cdot 3 \cdot \ldots \cdot (p-2) \cdot (p-1) \cdot 2}_{2\text{nd}} \cdot \ldots \\\\
& \cdot \underbrace{1 \cdot 2 \cdot 3 \cdot \ldots \cdot (p-2) \cdot (p-1) \cdot 1}_{p\text{th}} \cdot \ldots \cdot \quad \underbrace{1 \cdot 2 \cdot \cdot \ldots \cdot (n \bmod p)}_{\text{tail}} \pmod{p}.
\end{aligned}
$$
块的主要部分 $(p-1)!\ \mathrm{mod}\ p$ 很容易计算，可以应用 Wilson 定理：
$$
(p-1)!\equiv -1\pmod p
$$
总共有 $\lfloor \frac{n}{p} \rfloor$ 个块，因此需要将 $\lfloor \frac{n}{p} \rfloor$ 写到 $-1$ 的指数上。可以注意到结果在 $-1$ 和 $1$ 之间切换，因此只需要查看指数的奇偶性，如果是奇数，则乘以 $-1$。除了使用乘法，也可以从结果中减去 $p$.
最后一个部分块的值可以在 $O(p)$ 的时间复杂度单独计算。
这只留下每个块的最后一个元素。如果隐藏已处理的元素，可以看到以下模式：
$$
(n!)_p = \underbrace{ \ldots \cdot 1 } \cdot \underbrace{ \ldots \cdot 2} \cdot \ldots \cdot \underbrace{ \ldots \cdot (p-1)} \cdot \underbrace{ \ldots \cdot 1 } \cdot \underbrace{ \ldots \cdot 1} \cdot \underbrace{ \ldots \cdot 2} \cdots
$$
这也是一个修正的阶乘，只是维数小得多。它是：
$$
\left(\left\lfloor \frac{n}{p} \right\rfloor !\right)_p
$$
因此，在计算修改的阶乘 $(n!)_p$ 中，执行了 $O(p)$ 个操作，剩下的是计算 $(\lfloor \frac{n}{p} \rfloor !)_p$，于是有了递归，递归深度为 $O(\log_p n)$，因此算法的总时间复杂度为 $O(p \log_p n)$.
如果预先计算阶乘 $0!, 1!, 2!, \dots, (p-1)!$ 模 $p$，那么时间复杂度为 $O(\log_p n)$.
### 评价
适用于阶乘 $n!$ 比较大，而模数 $p$ 比较小的情况，并且 $p \in {Prime}$ 限制条件较多。
如果 模数 $p$ **比较大**，并且 $p \in {Prime}$ 可以使用 `NTT 快速阶乘`复杂度 $\Theta(\sqrt n \log n)$.
在竞赛中用处不是很大，但是是初等数论四大定理之一，是后续定理的前置理论基础。
## 扩展卢卡斯定理
### 定义
卢卡斯定理内容如下：对于质数 $p$，有
$$
C_{n}^{m}\bmod p = C_{\left\lfloor n/p \right\rfloor}^{\left\lfloor m/p\right\rfloor}\cdot C_{n\bmod p}^{m\bmod p}\bmod p
$$
### 扩展
卢卡斯定理中对于模数 $p$ 要求必须为**素数**，那么对于 $p$ 不是素数的情况，就需要用到 `EXLucas` 定理。
### 应用
考虑令 $p=p_1^{k_1} p_2^{k_2} p_3^{k_3} \dots p_{cnt}^{k_{cnt}}
原问题转化为求： \begin{cases} C_n^m \bmod p_1^{k_1} \\ C_n^m \bmod p_2^{k_2} \\ C_n^m \bmod p_3^{k_3} \\ \dotsm \end{cases}
的一个解，使用中国剩余定理合并即可。
针对其中一个式子 C_n^m \bmod p^{k} 即 \frac {n!} {m!(n-m)!} \bmod p^k ，不能直接求得逆元，因为 p^{k} 不保证为质数， m!(n-m)! 可能不存在逆元。
设
f(n)=n! \bmod p^k
类似地，下列算式也提取了公因子p,使得剩余部分可以**使用逆元**。
$$
\begin{aligned}
n! &= 1 \times 2 \times 3 \cdots n \\
&= (p \cdot 2p \cdot 3p \cdots) \cdot [1 \times 2 \times 3 \cdots (p-1) \cdot (p + 1) \cdot (p + 2) \cdots] \\
&= p^{\lfloor \frac{n}{p} \rfloor} (1 \times 2 \times 3 \cdots) \cdot [1 \times 2 \times 3 \cdots (p-1)] \cdot [(p + 1) \cdot (p + 2) \cdots) \\
&= \large p^{\lfloor \frac{n}{p} \rfloor} {\lfloor\frac{n}{p}\rfloor}! \prod_{p \nmid i}^n i \\
&= \large p^{\lfloor \frac{n}{p} \rfloor} {\lfloor \frac{n}{p}\rfloor}! (\prod_{p \nmid i}^{p^k} i)^{\lfloor \frac{n}{p^k} \rfloor} \textit{rest} \\
&= \large p^{\lfloor \frac{n}{p} \rfloor} f({\lfloor \frac{n}{p} \rfloor}) (\prod_{p \nmid i}^{p^k} i)^{\lfloor \frac{n}{p^k} \rfloor} \textit{rest} \\
\end{aligned}
$$
其中 $\text{rest}$ 是指最后一段 $p^k$ 的余项， $\displaystyle \large (\prod_{p \nmid i}^{p^k} i)$ 是在 $\bmod p^k$ 下的一个循环节，这样的循环节有 ${\lfloor \frac{n}{p^k} \rfloor}$ 个。
每进行一次上述操作，我们都可以从阶乘提取出 $p^{\lfloor \frac n p \rfloor}$ 这个因子，并把原问题变为一个规模更小的问题 $f({\lfloor \frac{n}{p} \rfloor})$ ，并在递归中累加因子的次数 $\lfloor \frac n p \rfloor$ ，记为 $h(n)$ ，最后总计 $h(n)$ 就是 $n!$ 中 $p$ 这个因子的数量。
至此，$n!$ 里的 $p$ 的因子已经被全部分离出来。
所以 $C_n^m \bmod p^{k}$即为：
$$ \large \frac {\frac {n!}{p^{h(n)}}} {\frac {m!}{p^{h(m)}}\frac {(n-m)!}{p^{h(n-m)}}} p^{h(n)-h(m)-h(n-m)} \bmod p^k
$$
即
$$
\frac{f(n)}{f(m)f(n-m)} p^{h(n)-h(m)-h(n-m)} \bmod p^k
$$
由于 $f(n)$ 已经是被约掉 $p$ 因子的，一定有逆元，并
使用 `CRT/ EXCRT` 合并所有的算数得到答案。
####代码
`````` C++ {.line-numbers}
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef __int128_t i128;
inline ll qpow(ll a, ll b, const ll& mod) {
if (!b)
return 1;
ll res = qpow(a, b / 2, mod);
if (b % 2)
return res * res % mod * a % mod;
return res * res % mod;
}
void EXgcd(i128 a, i128 b, i128& d, i128& x, i128& y) {
if (!b) {
d = a, x = 1, y = 0;
} else {
EXgcd(b, a % b, d, y, x);
y -= x * (a / b);
}
}
namespace EXCRT {
ll a[100010], m[100010];
ll EXCRT(const int n, const ll* a = a, const ll* m = m) {
i128 A = a[1], M = m[1];
A = (A % M + M) % M; // 避免出现负解
for (int i = 2; i <= n; i++) {
i128 p, q, d;
EXgcd(M, m[i], d, p, q);
if ((a[i] - A) % d)
return -1;
i128 bi = m[i] / d;
p *= (a[i] - A) / d;
p = (p % bi + bi) % bi; // 最小正整数
A = M * p + A;
M = M / d * m[i];
}
return A;
}
} // namespace EXCRT
/**
* @brief 瓶颈在模数MOD分解的因子上,复杂度的上界由单个因子p^k决定
* @param n 阶乘n
* @param Prime 模MOD分解出来的质因数,保证Prime是质数
* @param mod 质因数的幂次Prime^k
* @param h 返回n!里包含Prime的次方数
*/
i128 Wilson(ll n, ll Prime, ll mod, i128& h) {
if (n == 0)
return 1;
i128 res = 1, tmp = 1, hh = 0;
res *= Wilson(n / Prime, Prime, mod, hh);
res = (res % mod + mod) % mod;
h += n / Prime + hh; // 存Prime的次方
for (ll i = 1; i < mod; ++i) { // 计算循环节
if (i % Prime)
tmp = tmp * i % mod;
}
res *= qpow(tmp, n / mod, mod);
res = (res % mod + mod) % mod;
for (ll i = n / mod * mod; i <= n; i++) { // 循环节尾剩余的部分
if (i % Prime) {
res = res * i % mod;
}
}
return res;
}
/**
* @brief 瓶颈在模数MOD分解的因子上,复杂度的上界由单个因子p^k决定
*/
ll EXLucas(ll n, ll m, ll MOD) {
ll Prime[100] = {}, k[100] = {};
ll p = MOD, cnt = 0;
for (ll i = 2; i * i <= p; i++) { // 分解质因数
if (p % i == 0) {
Prime[++cnt] = i;
while (p % i == 0) {
p /= i, k[cnt]++;
}
}
}
if (p != 1)
Prime[++cnt] = p, k[cnt]++;
for (int i = 1; i <= cnt; i++) {
i128 exp = 0, mod = qpow(Prime[i], k[i], MOD + 1); // 不要取模
i128 res = 1, inv, y, h, d;
res *= Wilson(n, Prime[i], mod, h = 0), exp += h;
// 把Prime的因数全部提取出去,保证逆元有解
EXgcd(Wilson(m, Prime[i], mod, h = 0), mod, d, inv, y);
inv = (inv % mod + mod) % mod;
res *= inv;
res = (res % mod + mod) % mod;
exp -= h;
EXgcd(Wilson(n - m, Prime[i], mod, h = 0), mod, d, inv, y);
inv = (inv % mod + mod) % mod;
res *= inv;
res = (res % mod + mod) % mod;
exp -= h;
EXCRT::a[i] = res * qpow(Prime[i], exp, mod) % mod;
EXCRT::m[i] = mod;
}
return EXCRT ::EXCRT(cnt);
}
int main() {
ios::sync_with_stdio(0);
ll n, m, p;
cin >> n >> m >> p;
cout << EXLucas(n, m, p) << "\n";
return 0;
}
``````
### 评价
金牌题，一般不用到，用到不一般。
若不考虑 `EXCRT` 的复杂度，通过预处理 $\frac{n!}{n以内的p的所有倍数的乘积}\bmod{p}$，可以使时间复杂度优化至单次 $\Theta(p + \log p)$。而如果 p 是固定的，我们在一开始就可以对 p 进行分解，并进行预处理，可以达到总复杂度 $\Theta(p + T\log p)$.
该算法的瓶颈在模数 $p$ 分解的因子上,复杂度的上界由因子 $Prime^k$ 决定。
**综合得：该算法比较适合阶乘数较大，但是模 $p$ 分解出来的质数较小的情况。**
**另外，`EXLucas` 可以进行改造，实现 `EXWilson` 的功能(计算带 $(n!)_p$ 的方程)**
**特殊情况：** 一般在 **模数固定**的情况下，建议预处理模的因数。多次计算 `EXWilson` 循环节的代价较大，也建议预处理阶乘。
**实例：**
[[SDOI2010] 古代猪文](https://www.luogu.com.cn/problem/P2480)
>
> 给出 $G,n$（$1 \leq G,n \leq 10^9$），求：
>
>$$G^{\sum_{k\mid n}C_{n}^{k}} \bmod 999~911~659$$
>
> **代码**
>`````` C++ {.line-numbers}
>#include <bits/stdc++.h>
>using namespace std;
>typedef long long ll;
>#define MMOD (ll)999911659
>
>const ll Prime[5] = {0, 2, 3, 4679, 35617};
>ll f[5][40000];
>
>void init() {
> for (int i = 1; i <= 4; i++) {
> f[i][0] = 1;
> for (int j = 1; j <= Prime[i]; j++) {
> f[i][j] = f[i][j - 1] * j % Prime[i];
> }
> }
>}
>
>ll qpow(ll a, ll b, const ll& mod) {
> if (!b)
> return 1;
> ll res = qpow(a, b / 2, mod);
> if (b % 2)
> return res * res % mod * a % mod;
> return res * res % mod;
>}
>
>void EXgcd(ll a, ll b, ll& d, ll& x, ll& y) {
> if (!b) {
> d = a, x = 1, y = 0;
> } else {
> EXgcd(b, a % b, d, y, x);
> y -= x * (a / b);
> }
>}
>
>namespace EXCRT {
>ll a[10], m[10];
>
>ll EXCRT(const ll n, const ll* a = a, const ll* m = m) {
> ll A = a[1], M = m[1];
> A = (A % M + M) % M; // 避免出现负解
>
> for (int i = 2; i <= n; i++) {
> ll p, q, d;
> EXgcd(M, m[i], d, p, q);
> if ((a[i] - A) % d)
> return -1;
> ll bi = m[i] / d;
>
> p *= (a[i] - A) / d;
> p = (p % bi + bi) % bi; // 最小正整数
> A = M * p + A;
> M = M / d * m[i];
> }
> return A;
>}
>} // namespace EXCRT
>
>ll Wilson(ll n, ll Prime, ll Case, ll& h) {
> if (n == 0)
> return 1;
> ll res = 1, tmp = 1, hh = 0;
>
> res *= Wilson(n / Prime, Prime, Case, hh);
> res = (res % Prime + Prime) % Prime;
> h += n / Prime + hh; // 存Prime的次方
>
> tmp = f[Case][Prime - 1]; // 循环节
>
> res *= qpow(tmp, n / Prime, Prime);
> res = (res % Prime + Prime) % Prime;
>
> res = res * f[Case][n % Prime] % Prime; // 循环余项
>
> return res;
>}
>
>ll EXLucas(ll n, ll m, ll MOD) {
> ll cnt = 4;
>
> for (ll i = 1; i <= cnt; i++) {
> ll exp = 0, mod = Prime[i];
> ll res = 1, inv, y, h, d;
>
> res *= Wilson(n, Prime[i], i, h = 0), exp += h;
>
> // 把Prime的因数全部提取出去,保证逆元有解
> EXgcd(Wilson(m, Prime[i], i, h = 0), mod, d, inv, y);
> inv = (inv % mod + mod) % mod;
> res *= inv;
> res = (res % mod + mod) % mod;
> exp -= h;
>
> EXgcd(Wilson(n - m, Prime[i], i, h = 0), mod, d, inv, y);
> inv = (inv % mod + mod) % mod;
> res *= inv;
> res = (res % mod + mod) % mod;
> exp -= h;
>
> EXCRT::a[i] = res * qpow(Prime[i], exp, mod) % mod;
> EXCRT::m[i] = mod;
> }
> return EXCRT ::EXCRT(cnt);
>}
>
>int main() {
> ios::sync_with_stdio(0);
> init();
>
> ll n, g, sum = 0;
> cin >> n >> g;
> if (g % MMOD == 0) {
> cout << "0\n";
> return 0;
> }
> ll Eu_MOD = MMOD - 1;
> for (ll i = 1; i * i <= n; i++) {
> if (n % i == 0) {
> sum = (sum + EXLucas(n, i, Eu_MOD) % Eu_MOD) % Eu_MOD;
> if (i * i != n)
> sum = (sum + EXLucas(n, n / i, Eu_MOD) % Eu_MOD) % Eu_MOD;
> }
> }
> cout << qpow(g, sum, MMOD) << "\n";
> return 0;
>}
>``````
>
## 习题
- [P3807【模板】卢卡斯定理](https://www.luogu.com.cn/problem/P3807)
- [P2480 [SDOI2010] 古代猪文 卢卡斯定理](https://loj.ac/problem/10229)
- [P4720【模板】扩展卢卡斯](https://www.luogu.com.cn/problem/P4720)
- [P4345 [SHOI2015] 超能粒子炮·改](https://www.luogu.com.cn/problem/P4345)
- [Ceizenpok’s formula](http://codeforces.com/gym/100633/problem/J)