卢卡斯定理(Lucas定理)

搞懂Lucas定理！

Lucas 定理内容如下：对于质数p，有C(n, m) ≡ C(n/p,m/p) * C(n%p, m%p) (mod p)其中n/p, m/p下取整即下面的图片
![[Pasted image 20240429091355.png]]
卢卡斯定理 - OI Wiki (oi-wiki.org)
非常的简单啊
因为是用于组合数的，具体应用的分析和时间复杂度在这里[[排列组合#利用Lucas定理求解]]

代码

如果不计算求组合数的时间复杂度的话，时间复杂度为O(logn)

int lucas(int a, int b)
{
	if (a < p && b < p) return C(a, b); // a>=b 所以主要是判断a < p时终止
	return (LL)C(a % p, b % p) * lucas(a/p, b/p) % p;
}

在上面主要是a < p时终止，那么如果b == 0怎么办?首先无论a是否为0，当b == 0时它都一定为1，因为不选物品也是一种方案，不可能出现0，而1对于之前求得没有影响，也多算不了几次，所以不要紧。

证明

这里推荐一篇文章卢卡斯定理证明 - onlyblues - 博客园 (cnblogs.com)里面的证明非常严谨，这里我采用感性的思想证明一下
这里写明文章中第一种证明方式

因为我不会LaTeX语法，所以符号看起来会很别扭，可以用笔写下来看看是什么样
要证明C(a, b) ≡ C(a/p,b/p) * C(a%p, b%p) (mod p)其中`a/p, b/p`下取整

使用的公式
(x + y)^p ≡ x^p + y^p (mod p)
(x + y) ^ p^a ≡ x ^ p^a + y ^ p^a (mod p), p^a即p的a次方

证明
(x + y)^p 可以用二项式定理分解，其中每一项为C(p, i) * x^i * y^(p-i)
因为p是质数，所有比他小的正整数(除了1)都无法通过除的方式消去它
那么当i在[1, p - 1]时可知C(p, i)一定可以化为p * k的形式，k为任意正整数
那么此时C(p, i) % p = 0, 那么C(p, i) * x^i * y^(p-i) % p = 0
所以这些项一定为0
由此可知(x + y)^p ≡ C(p, 0) * y^p + C(p, p) * x^p (mod p)
即(x + y)^p ≡ y^p + x^p (mod p)

而(x + y)^p^a ≡ x^p^a + y^p^a (mod p)可以通过同样的方式证明这里省去

接下来开始真正的证明

首先我们可以把C(a, b)中的a，b分解为p进制的形式，即
a = ak*p^k + ak-1 * p^(k - 1) + ... + a0*p^0, 其中ai为非负整数
b = bk*p^k + bk-1 * p^(k - 1) + ... + b0*p^0, 其中bi为非负整数



由(x + y)^p ≡ x^p + y^p (mod p)可得
(1 + x)^p ≡ 1 + x^p (mod p)
由(x + y)^p^a ≡ x^p^a + y^p^a (mod p)可得
(1 + x)^p^a ≡ 1 + x^p^a (mod p)

因为(1 + x)^a = (1 + x)^(ak*p^k + ak-1 * p^(k - 1) + ... + a0*p^0)
= (1+x)^(ak*p^k) * (1+x)^(ak-1 * p^k-1) * ... * (1+x)^a0*p^0
≡ (1 + x^p^k)^ak * (1 + x^p^(k-1))^ak-1 * ... * (1 + x^p^0)^a0 (mod p)

即(1 + x)^a ≡ (1 + x^p^k)^ak * (1 + x^p^(k-1))^ak-1 * ... * (1 + x^p^0)^a0 (mod p) 设为(1)

我们要求C(a, b)，而C(a, b)恰好等于(1 + x)^a中含x^b项的系数
而x^b = x^(bk*p^k + bk-1 * p^k-1 + ... + b0*p^0)
      = x^(bk*p^k) * x^(bk-1 * p^k-1) * ... * x^b0*p^0 设为(2)
      
这时候看看(2)中的x^(bk*p^k)和(1)中的(1 + x^p^k)^ak, 发现了什么?
可以发现前者在后者的系数为C(ak, bk)项中
那么就可以把这两个联系起来了，类比一下
由(1)和(2)可得C(a, b) ≡ C(ak, bk) * C(ak-1, bk-1) * ... * C(a0, b0) (mod p) ,设为(3)
这时候证明就差不多了

因为
a/p = ak*p^k + ak-1 * p^(k - 1) + ... + a1*p^1
b/p = bk*p^k + bk-1 * p^(k - 1) + ... + b1*p^1,
a0 = a % p
b0 = b % p
所以把C(a, b)换成C(a/p, b/p)可发现(3)中的C(a0, b0)就没了
于是
C(a/p, b/p) ≡ (ak, bk) * C(ak-1, bk-1) * ... * C(a1, b1) (mod p)
C(a%p, b%p) = C(a0, b0) 
那么C(a, b) ≡ C(a/p, b/p) * C(a%p, b%p) (mod p)

证毕
上面有个小问题，就是a >= b是知道的
但是怎么保证ai >= bi, 从而满足C(ai, bi)成立呢?
事实上是存在ai < bi这种情况的, 但是，我们可以想想
a >= b是一定的，如果ai < bi, 那么a就一定比b多含有p的因子, 
因为一定有一个ak > bk, 如果这样在最后的阶乘上a!要比b!至少多出p^k的因子, k > i >=0 
k >= 1, 且因为ai < bi, ak > bk, (a - b)!的结果不可能含有足够多的p因子来抵消掉a!上的p^k, 那么此时C(a, b) % p == 0
而在代码中因为遇到这种情况我们在求组合数的时候会返回一个0, 最后结果也就是0
因此，这里虽然是理想化了，但不影响正确性，也不影响证明。

可以自己试试
欢迎指出错误