多彩编程 多彩编程MZPH · CODE BLOG
ARTICLE DETAIL

文章详情

深耕前端与后端开发技术的一线实战笔记与踩坑复盘。

数论算法入门:从欧拉筛到Lucas定理的完整逻辑链

数论算法入门:从欧拉筛到Lucas定理的完整逻辑链 很多人在入门数论算法的时候被卡住的地方往往不是某一个定理本身而是这些知识点之间的“衔接感”为什么要先学欧拉筛快速幂和它有什么关系裴蜀定理与拓展欧几里得是不是同一件事乘法逆元为什么要单独拉出来讲Lucas定理又是解决什么场景的。其实这一串内容是一条完整的逻辑链先用欧拉筛批量拿到素数表用快速幂处理指数运算用扩展欧几里得解同余方程再用乘法逆元处理模意义下的除法最后用Lucas定理处理大组合数取模。把它们串起来理解比孤立地背模板有效得多。如果你是刚接触算法竞赛、或者正在刷数论专题但总感觉“知识点都见过、做题还是不会”这篇内容就是为你准备的。我会从每个知识点的“为什么存在”讲起给出可以直接抄进代码库的模板再补上那些文档里不会写、但实操时一定会踩的细节坑。1. 欧拉筛让每个合数被它的最小质因子筛掉1.1 为什么不用埃氏筛埃氏筛的核心思想很简单从小到大枚举质数然后把它的倍数全部标记为合数。比如筛到质数2就把4、6、8、10……全部划掉筛到3就把9、12、15……划掉。问题在于一个合数会被重复标记很多次12既会被2标记又会被3标记。虽然可以通过优化减少重复理论上界仍是 O(n log log n)在n为10^7量级时性能已经吃紧。欧拉筛线性筛的思路是只让每个合数被它的最小质因子标记一次。这是它名字叫“线性筛”的由来——每个数在遍历过程中最多被操作一次总复杂度严格 O(n)。1.2 欧拉筛核心逻辑先说代码模板语言用C因为竞赛里最常用const int N 1e7 5; vectorint primes; bool notPrime[N]; void eulerSieve(int n) { for (int i 2; i n; i) { if (!notPrime[i]) { primes.push_back(i); } for (int p : primes) { if (i * p n) break; notPrime[i * p] true; if (i % p 0) break; // 核心保证每个合数只被最小的质因子筛掉 } } }这个模板看起来简单但里面有两个细节决定了它为什么是“线性”的第一内层循环的退出条件是i * p n而不是i n / p。前者在i和p都很大的时候可能溢出int建议写成if (i n / p) break;避免整数溢出。这是我踩过的真实坑n取到1e8时i * p会超过int范围变成负数筛出来一堆错误结果。第二if (i % p 0) break;是灵魂。当i能被p整除说明p是i的最小质因子此时如果继续用更大的质因子去标记就会造成重复。直接break终止循环保证每个合数在唯一的时机被标记它的最小质因子乘以某个数。1.3 欧拉筛的常见坑欧拉筛虽然快但有几个容易翻车的点数组初始化notPrime[0]和notPrime[1]要手动置为true否则0和1会被误判为质数。这个错误非常隐蔽因为很多题的数据范围从2开始小数据测不出来。质数表是否需要包含n题目要求筛 [2, n] 还是 [2, n) 写之前先确认边界否则多筛一个或少筛一个都是错。空间占用bool数组虽然每个元素只用1字节但n1e8时仍需约100MB内存这在部分内存限制256MB的题目里很危险。替代方案是用bitset或者只在需要小范围质数时用区间筛。2. 快速幂把指数拆成二进制来加速2.1 从朴素循环到二进制拆分计算a^b mod m朴素写法是一个循环乘b次复杂度 O(b)。当b达到10^9甚至10^18时这个循环是跑不完的。快速幂的思路是利用指数的二进制表示比如计算a^11因为11的二进制是1011也就是11 2^3 2^1 2^0所以a^11 a^(2^3) * a^(2^1) * a^(2^0)把a不断平方a^1, a^2, a^4, a^8……然后根据二进制位决定乘不乘进结果。这样只需要循环log2(b)次约等于30次循环就能处理10^9级别的指数。2.2 两种实现写法递归版本最容易理解long long fastPow(long long a, long long b, long long m) { if (b 0) return 1 % m; long long t fastPow(a, b / 2, m); t t * t % m; if (b 1) t t * a % m; return t; }迭代版本更推荐避免递归栈开销也更好背long long fastPow(long long a, long long b, long long m) { long long res 1 % m; while (b) { if (b 1) res res * a % m; a a * a % m; b 1; } return res; }为什么用1 % m而不是直接1因为如果m1任何数取模都是0直接返回1就错了。这个边界在竞赛里经常出现属于典型的“样例测不出来但提交就WA”的坑。2.3 指数为负的情况普通快速幂处理的是非负指数。如果指数是负数表示求它的逆元这依赖第4节讲的内容。在代码实现时我会做一个封装long long modPow(long long a, long long b, long long m) { if (b 0) { // 前提a与m互质 return fastPow(modInverse(a, m), -b, m); } return fastPow(a, b, m); }这个封装在以后写组合数、概率题时会很方便。3. 裴蜀定理与拓展欧几里得算法3.1 裴蜀定理先给结论裴蜀定理描述了一个非常基础的事实对于任意整数 a、b方程ax by gcd(a, b)一定存在整数解。再扩展一步方程ax by c有整数解当且仅当gcd(a, b)能整除 c。这是判断不定方程是否有解的最快方法。你可以把gcd(a, b)理解为“a和b能组合出来的最小正增量”所有能被组合出来的数都是它的倍数。3.2 拓展欧几里得的推导过程扩展欧几里得算法一般缩写为exgcd做的事情是不仅求出 gcd(a, b)还顺带求出满足ax by gcd(a, b)的一组整数解。它的推导用到了欧几里得算法的递归过程。假设我们递归调用了exgcd(b, a % b)得到了b * x1 (a % b) * y1 gcd(b, a % b)因为gcd(a, b) gcd(b, a % b)所以我们想从这个已知解推出a * x b * y gcd(a, b)的解。关键在于把a % b展开a % b a - floor(a / b) * b代入上面的等式b * x1 (a - floor(a / b) * b) * y1 gcd整理后a * y1 b * (x1 - floor(a / b) * y1) gcd对比a * x b * y gcd得到x y1 y x1 - floor(a / b) * y1这就是拓展欧几里得的核心递推公式。模板代码long long exgcd(long long a, long long b, long long x, long long y) { if (b 0) { x 1; y 0; return a; } long long g exgcd(b, a % b, y, x); y - a / b * x; return g; }这里用了一个C技巧递归调用时把y和x的位置交换这样返回后y x1、x y1自然成立再补上y - a / b * x就完成了递推。很多教材里写的是临时变量版本而这个交换写法更简洁实战中推荐。3.3 从一组解到通解exgcd只给我们一组特解但题目往往需要更多信息要么是最小非负解要么是某个范围内的解。假设已经得到特解(x0, y0)那么通解是x x0 k * (b / g) y y0 - k * (a / g)其中g gcd(a, b)k为任意整数。为什么周期是b / g而不是b因为要保持a * x b * y不变x每增加一个b / ga * x增加的量是a * (b / g) (a * b) / g lcm(a, b)这是a和b的公共倍数因此y需要减少(lcm / b) a / g来抵消。求x的最小非负解常用的代码是long long g exgcd(a, b, x, y); b_g b / g; x (x % b_g b_g) % b_g;这个(x % mod mod) % mod的写法在数论里到处都是用来确保结果落在 [0, mod) 区间。4. 同余方程与乘法逆元4.1 用extgcd解同余方程线性同余方程长这样a * x ≡ b (mod m)它等价于存在整数k使得a * x - b k * m移项得a * x - m * k b这就转化成了不定方程a * x m * t b其中 t -k。根据裴蜀定理有解的充要条件是gcd(a, m) | b。具体解法分两步用exgcd求出a * x m * t gcd(a, m)的一组解(x0, t0)两边同时乘以b / gcd(a, m)得到原方程的一组特解4.2 乘法逆元是什么在模意义下除法没有直接定义。a / b在模m下是什么意思可以理解为“乘以b的逆元”。b的逆元指的是一个整数b^{-1}使得b * b^{-1} ≡ 1 (mod m)。一旦有了逆元a / b就可以写成a * b^{-1} mod m。这解决了组合数计算中模除法的问题。求逆元有三种常用方法各有适用场景方法一费马小定理适用于m为质数当m是质数p且a不是p的倍数时a^(p-1) ≡ 1 (mod p)所以a * a^(p-2) ≡ 1 (mod p)即逆元就是a^(p-2) mod p。long long inv fastPow(a, p - 2, p);这是三行代码里最快的解法前提是p必须是质数。方法二扩展欧几里得适用于gcd(a, m) 1解方程a * x m * y 1得到的x就是逆元long long x, y; exgcd(a, m, x, y); long long inv (x % m m) % m;这个方法不要求m是质数只需要a与m互质实用性更广。方法三线性递推适用于预处理1到n的逆元如果题目需要连续计算[1, n]所有数的逆元单个求是O(n log m)用递推可以做到 O(n)inv[1] 1; inv[i] m - (m / i) * inv[m % i] % m;这个公式的正确性推导用的是替换法把m拆成i * (m / i) m % i在模m下等于0然后变形出递推式。注意m必须是质数才能保证每个数都可逆且递推公式对非质数m不成立。4.3 逆元的应用场景逆元最典型的应用是组合数取模。计算C(n, k) n! / (k! * (n - k)!)时分母在模意义下不能直接除先预处理阶乘数组f和逆元数组inv然后C(n, k) f[n] * inv[k] % m * inv[n - k] % m;这样组合数计算变成三次乘法和两次取模速度极快是组合数学题的标准打法。5. Lucas定理大组合数取模救场5.1 Lucas定理解决什么问题直接算C(n, m) mod p如果n和m非常大比如10^18而p相对较小比如10^6以下预处理阶乘数组根本开不下。Lucas定理说的是把n和m分别写成p进制表示n n_k * p^k n_{k-1} * p^{k-1} ... n_0 m m_k * p^k m_{k-1} * p^{k-1} ... m_0那么C(n, m) ≡ C(n_k, m_k) * C(n_{k-1}, m_{k-1}) * ... * C(n_0, m_0) (mod p)也就是说大组合数取模可以被拆成若干个“位数上的小组合数”的乘积而每个小组合数的n、m都不超过p-1在预处理的阶乘表范围内。5.2 模板实现long long lucas(long long n, long long m, int p) { if (m 0) return 1; return C(n % p, m % p, p) * lucas(n / p, m / p, p) % p; }其中C(n, m, p)用预处理阶乘逆元实现long long C(int n, int m, int p) { if (m n) return 0; return fac[n] * inv[m] % p * inv[n - m] % p; }需要注意fac和inv数组的长度只开到p因为Lucas递归中每次调用C的参数都小于p所以预处理的规模不会比p大多少。实操中还有一个小优化当n和m的量级差不多时直接判断m 0提前返回能少递归好几层。6. 常见问题排查与经验总结6.1 问题速查表这几块内容在实战中遇到WA或TLE最常出问题的点如下问题现象可能原因排查办法筛完的质数表不完整数组下标越界或n边界不对打印前20个质数核对快速幂结果总为0取了模数为1的极端值检查边界1 % mexgcd返回的g为负数输入存在负数gcd处理不一致先取绝对值符号最后处理逆元结果不正确前提条件gcd(a, m) ! 1先判断互质再求逆Lucas定理超时未剪枝m0或预处理数组过大加提前返回缩小数组6.2 我的几条实践心得第一模板不要背要会推导。exgcd的结论如果忘了公式现场推一遍也就五分钟但背错一个系数就是半小时的debug时间。尤其注意求y - a / b * x那一步很多教材写法不同必须理解为什么这么写。第二先判边界再写代码。数论题的WA往往出在“0、1、n1”这些极小边界上。我的习惯是先写几个极端case的单测再跑大样例能省下很多提交罚时。第三乘号两边都要考虑溢出。a * a % m在a接近10^18时会爆long long。必要时用__int128过渡或者手写快速乘把乘法拆成多次加法取模。第四能预处理就预处理。如果题目是多组询问逆元和阶乘宁可一次性算好存下来也不要每次重新求。O(n)预处理 O(1)查询的组合数方案比单次计算快一个数量级。我个人最深的体会是这七个知识点单独看都不难但把它们组合起来以后你才能看懂组合数、概率DP、容斥原理这些更上层的内容。数论是一环扣一环的体系把这条逻辑链打通了后面学什么都会顺畅很多。
返回列表