算法学习笔记：Lucas 定理/exLucas
\(\text{Update 2025/4/3}\)：修正了一些笔误。感谢 @Lonelyhacker 指出错误。
一些记号
对于质数 \(p\)，我们用 \(v_p(n)\) 表示 \(n\) 的质因数分解中 \(p\) 的次数，用 \((n)_p\) 表示将 \(n\) 中的所有 \(p\) 因子除去后的数。可以得到
Lucas 定理
Lucas 定理常用于求大组合数取模一个较小的质数。
我们经常通过 \(O(n)\) 递推预处理阶乘和阶乘的逆元，来实现 \(O(1)\) 回答组合数，这种方法的局限性在于 \(n\) 较大时会炸，或者当模数较小时无法保证分母的阶乘部分存在逆元。Lucas 定理就可以很好地解决这两个问题。
Lucas 定理：对于质数 \(p\)，有
特别地，我们规定当 \(a<b\) 时 \(\binom{a}{b}=0\)。
为了证明 Lucas 定理，我们先引出一个引理。
引理：对于质数 \(p\)，\((x+y)^p\equiv x^p+y^p\pmod{p}\)。
证明：用二项式定理对其展开，考察各项系数 \(\binom{p}{i}=\frac{p!}{i!(p-i)!}\)，当 \(1\leq i<p\) 时，分母中显然不含 \(p\) 因子，而分子中含有 \(p\) 因子，因此此时 \(\binom{p}{i}\equiv 0\pmod{p}\)，而当 \(i=0\) 或 \(i=p\) 时，\(\binom{p}{i}\equiv 1\pmod{p}\)。综上，\((x+y)^p\) 在 \(\bmod{p}\) 意义下只有 \(x^p\) 和 \(y^p\) 处系数不为 \(0\)。引理得证。
Lucas 定理的证明：根据引理，\((1+x)^p\equiv 1+x^p\pmod{p}\)。令 \(n=qp+r(0\leq r<p)\)，则
\((1+x^p)^q\) 展开后的各项次数都是 \(p\) 的倍数，而 \((1+x)^r\) 展开后的各项次数都 \(<p\)。根据多项式乘法，\([x^m](1+x)^n\) 的系数由两个多项式中次数和为 \(m\) 的项贡献而得，这种表示方法是唯一的，即 \(m=\lfloor\frac{m}{p}\rfloor p+(m\bmod{p})\)，因此对应系数相乘就得到了 Lucas 定理的形式，得证。
Lucas 定理的一种等价理解就是把 \(n,m\) 都表示成 \(p\) 进制数，那么 \(\binom{n}{m}\) 就是其对应数位的组合数之积。
对于 \(\binom{n\bmod{p}}{m\bmod{p}}\bmod{p}\)，上下部分都 \(<p\)，我们可以正常递推 \(O(p)\) 预处理阶乘和阶乘的逆元 \(O(1)\) 回答，而对于 \(\binom{\lfloor\frac{n}{p}\rfloor}{\lfloor\frac{m}{p}\rfloor}\bmod{p}\)，这是一个子问题，递归求解即可，边界为 \(m=0\) 时返回 \(1\)。每次递归规模变为原来的 \(\frac{1}{p}\)，因此 Lucas 定理可以以 \(O(\log_p(n))\) 的复杂度回答一个组合数。
代码：
int comb(int n, int m) {
if (n < m) return 0;
return 1ll * fac[n] * ivfac[m] % p * ivfac[n - m] % p;
}
int lucas(int n, int m) {
if (n < m) return 0;
if (!m) return 1;
return 1ll * comb(n % p, m % p) * lucas(n / p, m / p) % p;
}
Wilson 定理
在讲解 exLucas 前，我们需要引出 Wilson 定理。
Wilson 定理：\((p-1)!\equiv -1\pmod{p}\) 当且仅当 \(p\) 为质数。
证明：逆元关系是相互的，我们考虑 \([1,p-1]\) 中的数能否两两配对成逆元。考察逆元为本身的情况，即 \(x^2\equiv 1\pmod{p}\)，由 Lagrange 定理，该方程至多有 \(2\) 个根，我们容易直接给出这 \(2\) 个根：\(x\equiv\pm 1\pmod{p}\)。因此对于 \([2,p-2]\) 中的数，逆元成对出现，因此它们的乘积 \(\bmod{p}\) 的结果为 \(1\)，再乘上 \(1\) 和 \(-1\)，就得到了乘积 \(\bmod{p}\) 的结果为 \(-1\)。
再来讨论 \(p\notin \mathbb{P}\) 时的情况。这时考虑反证，假设此时有 \((p-1)!\equiv -1\pmod{p}\)，即存在整数 \(t\) 使得 \((p-1)!=tp-1\)，而因为 \(p\) 是合数，因此必定存在质数 \(q\mid p\)，即 \(p=tq(t\in \mathbb{Z})\)，因此 \((p-1)!=tkq-1\equiv -1\pmod{q}\)，但是因为 \(q<p\)，所以 \(q\) 必然包含在 \([1,p-1]\) 中，则 \((p-1)!\equiv 0\pmod{q}\)，矛盾。因此当 \(p\) 为合数时一定有 \((p-1)!\not\equiv -1\pmod{p}\)。得证。
不过，对于我们更有用的是，Wilson 定理在模数为质数幂 \(p^k\) 时的推广。
Wilson 定理的推广：对于质数 \(p\) 和正整数 \(k\)，
证明：我们还是考虑逆元是否能成对匹配，实际上我们就是在考虑 \(x^2\equiv 1\pmod{p^k}\)，此时 Lagrange 定理不再适用，我们无法保证除了 \(x\equiv \pm 1\pmod{p}\) 外没有其它的解。
先讨论 \(p>2\) 的情况。对上面的同余方程进行因式分解得到 \(x^2-1\equiv (x-1)(x+1)\equiv 0\pmod{p^k}\)。此时 \(\gcd(x-1,x+1)\mid 2\)，因此 \(x-1\) 和 \(x+1\) 不可能同时为 \(p\) 的倍数，必然有 \(x\equiv 1\pmod{p^k}\) 或者 \(x\equiv -1\pmod{p^k}\)。于是这种情况依然有 \((p^k!)_p\equiv -1 \pmod{p}\)。
再来讨论 \(p=2\) 的情况。此时 \(x-1\) 和 \(x+1\) 可能同时为 \(2\) 的倍数，考虑这种情况，设 \(x-1=2y\)，\(x+1=2(y+1)\)，将同余方程两侧同时除掉一个 \(4\)，得到 \(y(y+1)\equiv 0\pmod{2^{k-2}}\)，此时必然有 \(y\equiv \pm 1\pmod{2^{k-2}}\)，即 \(x\equiv 2^{k-1}\pm 1\pmod{p^k}\)。因此，当 \(p=2\) 时，\(x\) 有 \(4\) 种可能的解：\(x\equiv \pm 1\pmod{2^k}\) 或 \(x\equiv 2^{k-1}\pm 1\pmod{p^k}\)。而当 \(k=1\) 时，\(4\) 个解重合；当 \(k=2\) 时，两对解重合。即当 \(p=2\land k\leq 2\) 时依然有 \((p^k!)_p\equiv -1 \pmod{p}\)，而当 \(p=2\land k\geq 3\) 时，\((p^k!)_p\equiv 1 \pmod{p}\)。得证。
Wilson 定理即为该推广在 \(k=1\) 时满足 \((p!)_p=(p-1)!\) 时的特殊情况。
exLucas
Luogu P4720 【模板】扩展卢卡斯定理/exLucas
题意：给出 \(n,m,t\)，求 \(\binom{n}{m}\bmod{t}\)，\(1\leq m\leq n\leq 10^{18}\)，\(2\leq t\leq 10^6\)，不保证 \(t\) 是质数。
我们尝试快速求解组合数对任意模数取模的结果。
对于特定模数推广到任意模数的问题，一种很套路的做法就是考虑模数的质因数分解 \(t=\prod_{k=1}^kp_i^{c_i}\)，倘若我们能快速算出模数是质数幂 \(p_i^{c_i}\) 的结果 \(a_i\)，我们就能列出一个同余方程组：
模数两两互质，我们就可以 CRT 求出答案了。
因此对于组合数对任意模数取模的问题，我们将其特化成对质数幂 \(p^k\) 取模的问题，即求解
直接用定义展开组合数，得到
一个很直接的想法就是求出分母的逆元，我们无法保证其与 \(p^k\) 互质，即无法保证逆元存在，因此考虑除去各个阶乘的所有 \(p\) 因子：
此时分母必定和 \(p^k\) 互质，我们就可以用 exgcd 求出它的逆元了。
先考虑如何快速计算 \(v_p(n!)\)。这个就是阶乘的质数因子计数，是经典的。具体来说，先考察 \([1,n]\) 中所有 \(p\) 的倍数，它们至少各贡献 \(1\) 个 \(p\) 因子，总计贡献 \(\lfloor\frac{n}{p}\rfloor\) 个 \(p\) 因子；再考察 \([1,n]\) 中所有 \(p^2\) 的倍数，它们至少各贡献 \(2\) 个 \(p\) 因子，而其中的 \(1\) 个已经被前面计算了，所以总计贡献 \(\lfloor\frac{n}{p^2}\rfloor\) 个 \(p\) 因子……以此类推。可以得出
于是我们就可以 \(O(\log_p(n))\) 地计算出 \(v_p(n!)\)。这部分代码如下：
int cntp(ll x, int p) {
int res = 0;
for (; x; x /= p) res += x / p;
return res;
}
再考虑如何快速计算 \((n!)_p\)。先来看 \([1,n]\) 中所有与 \(p\) 互质的数的乘积，令
由于是在 \(\bmod{p^k}\) 意义下，我们有 \(f(1,p^k)=f(p^k+1,2p^k)=\cdots\)，因此我们将 \([1,n]\) 中的数每 \(p^k\) 个为一组，每组与 \(p\) 互质的数的乘积相等，即
注意到 \(f(1,p^k)\) 就是 \((p^k!)_p\)，可以使用 Wilson 定理的推广得出。
上面的式子只考虑了所有与 \(p\) 互质的数的乘积，再来看所有 \(p\) 的倍数把 \(p\) 因子除尽后的乘积，令
显然有 \((kp)_p=(k)_p\)，因此我们可以把乘积内的 \(p\) 去掉，得到
这是一个子问题，我们可以递归求解！
把 \(f(1,n)\) 和 \(g(n)\) 两部分乘起来，我们就得到了
我们可以 \(O(p^k)\) 预处理出所有的 \(f(1,i)(1\leq i< p^k)\)，做到单次询问 \(O(\log_p(n))\)。
这部分的代码：
f[0] = 1;
for (int i = 1; i < pp; ++i) f[i] = i % p ? 1ll * f[i - 1] * i % pp : f[i - 1];
int calcp(ll x, int p, int c) {
int pp = qpow(p, c), res = 1;
bool neg = p > 2 || c <= 2;
for (; x; x /= p) {
if ((x / pp) & neg) res = pp - res;
res = 1ll * res * f[x % pp] % pp;
}
return res;
}
至此，我们已经能快速算出 \(v_p(n!)\) 和 \((n!)_p\) 了，于是计算 \(\binom{n}{m}\bmod{p^k}\) 的问题也就迎刃而解了。
综上，我们得到了预处理 \(O(\sqrt{t}+\sum p_i)\)，单次询问 \(O(\log_p{n})\) 的 exLucas 算法。
注意实现中，在 \(v_p(n!)-v_p(m!)-v_p((n-m)!)\geq k\) 时可以直接返回 \(0\)。
exLucas 的代码：
ll n, m;
int p;
int sz, pr[C], cnt[C], f[P];
int qpow(int a, int b) {
int res = 1;
for (; b; b >>= 1) {
if (b & 1) res = 1ll * res * a;
a = 1ll * a * a;
}
return res;
}
ll exgcd(ll a, ll b, ll &x, ll &y) {
if (!b) { x = 1, y = 0; return a; }
ll res = exgcd(b, a % b, x, y), tmp = x;
x = y, y = tmp - a / b * y;
return res;
}
ll calc_inv(ll x, int p) {
ll t1, t2; exgcd(x, p, t1, t2);
t1 %= p; if (t1 < 0) t1 += p;
return t1;
}
void expand(int x) {
for (int i = 2; i * i <= x; ++i) {
if (x % i) continue;
pr[++sz] = i;
while (!(x % i)) ++cnt[sz], x /= i;
}
if (x != 1) pr[++sz] = x, cnt[sz] = 1;
}
int cntp(ll x, int p) {
int res = 0;
for (; x; x /= p) res += x / p;
return res;
}
int calcp(ll x, int p, int c, int pp) {
int res = 1;
bool neg = p > 2 || c <= 2;
for (; x; x /= p) {
if ((x / pp) & neg) res = pp - res;
res = 1ll * res * f[x % pp] % pp;
}
return res;
}
int lucas_pk(ll n, ll m, int p, int c, int pp) {
int x = cntp(n, p) - cntp(m, p) - cntp(n - m, p);
if (x >= c) return 0;
f[0] = 1;
for (int i = 1; i < pp; ++i) f[i] = i % p ? 1ll * f[i - 1] * i % pp : f[i - 1];
int res = qpow(p, x);
int fn = calcp(n, p, c, pp), fm = calcp(m, p, c, pp), fnm = calcp(n - m, p, c, pp);
res = 1ll * res * fn % pp * calc_inv(fm, pp) % pp * calc_inv(fnm, pp) % pp;
return res;
}
int ex_lucas(ll n, ll m, int p) {
int res = 0;
for (int i = 1; i <= sz; ++i) {
int pp = qpow(pr[i], cnt[i]);
int mp = p / pp;
ll inv = calc_inv(mp, pp);
int v = lucas_pk(n, m, pr[i], cnt[i], pp);
res = (res + 1ll * inv * mp % p * v % p) % p;
}
return res;
}
浙公网安备 33010602011771号