上周有个朋友发来一道题让我帮忙看给一个 n最大 10^7求 1 到 n 所有数的约数个数之和。他说自己写了一个 O(n√n) 的版本n 到 10^5 就卡死了问我是不是该上个什么高级数据结构。我看了眼他的代码第一反应不是“优化”而是“换数学视角”。这题本质是统计有多少对 (a,b) 满足 a×b ≤ n把枚举约数改成枚举倍数一个调和级数就解决了。这种事情我碰到过太多次不是题难而是大家把数学公式当成了课本上的装饰品。这篇文章我想把几个在算法里“低调但极其好用”的数学公式系统地聊一遍。我会讲清楚三件事每个公式是什么、为什么在算法里有用、用的时候有哪些坑。内容覆盖调和级数、斯特林公式、费马小定理、贝叶斯公式和拉格朗日插值。适合正在刷题的人、准备面试的人、写工程代码的人以及所有对“算法背后的数学直觉”感兴趣的人。1. 从一道“约数统计”题说起公式到底被谁忽视了1.1 一个常见到让人忽略的暴力写法题目是这样的给定 n求 1 到 n 所有数的约数个数之和。大部分人的第一版代码长这样对每个 i用从 1 到 √i 的循环数约数。单看一个数√i 并不大但所有数加起来就是 Σ√i约等于 (2/3)n^1.5。n 到 10^7 时这个量级是十亿次级别本地跑完得几十秒放到在线评测系统里基本就是超时。我当时问那个朋友你有没有想过约数这个东西是“成对出现”的一个正整数 i 如果有一个约数 d那么一定存在一个配对约数 i/d。所以“1 到 n 的所有约数”本质上就是所有满足 a×b ≤ n 的有序数对 (a,b)。这句话看起来像废话但把它写下来问题的复杂度就变了。1.2 为什么第一反应不是数学而是“加个缓存”很多人的第一反应是用筛法预处理每个数的约数表然后累加或者用记忆化缓存已经算过的约数个数。这些做法不是不行而是没有抓住问题的本质。筛法预处理约数个数确实可以做到 O(n log n)但前提是你已经知道“枚举倍数”这个套路如果不知道你还是会陷入“对每个数单独分解”的思维定式。我的观点是数据结构是战术数学是战略。遇到计数类问题先别急着写代码先在草稿纸上把问题翻译成一个数学表达式。很多时候表达式一变形复杂度就自动降下来了。这也是为什么我特别想写这篇文章——数学公式不是拿来背的是拿来对问题做“变形”的。1.3 这篇文章要聊的五个公式一句话版为了让你有个整体印象我先把后面要展开的几个公式列在这里调和级数H_n 1 1/2 ... 1/n ≈ ln n γ用来估计算法里“枚举倍数”这类循环的真实复杂度。斯特林公式n! ≈ √(2πn)(n/e)^n用来估算阶乘量级、组合数大小和概率。费马小定理a^(p-1) ≡ 1 (mod p)用来在模素数意义下求逆元本质是“取模世界里的除法”。贝叶斯公式P(A|B) P(B|A)P(A)/P(B)用来做概率推理也是很多机器学习算法的底层框架。拉格朗日插值给定 k1 个点构造一个 k 次多项式穿过所有点用来在稀疏样本之间做精确估计。这五个公式的共同特点是教科书里都有但算法课不怎么强调刷题时也容易被跳过。但它们在关键时刻能省掉一整层循环或者把一个看似无法下手的题变成套公式。2. 调和级数复杂度分析里那个假装是 log 的隐藏常数2.1 从计数约数到枚举倍数双层循环的复杂度真相回到刚才的约数统计问题。既然约数是成对出现的我们反过来考虑对于每个正整数 a它能作为多少个数的约数答案是 floor(n/a)也就是 a 的倍数在 1 到 n 中的个数。于是总数就是Σ_{a1}^{n} floor(n/a) ≈ n × H_n ≈ n ln n γn这个式子看着抽象我换成代码你就懂了cnt 0 for a in range(1, n 1): for b in range(a, n 1, a): # 步长是a也就是枚举a的倍数 cnt 1内层循环的次数是 n/a所以总次数是 n/1 n/2 ... n/n。这是典型的调和级数求和。很多人看到双重循环就默认是 O(n²)但这个内层循环的步长不是 1而是 i总次数只有 n ln n。举一个具体的数字n 10000 时Σ floor(n/i) 大约是 93648。你自己心里想一下一个双层循环只执行不到十万次这和你印象里的“O(n²) 一定爆炸”完全不同。所以以后看到内层步长是外层变量的循环先算一算 Σ n/i而不是直接写个 O(n²) 的帽子扣上去。2.2 埃氏筛的 O(n log log n) 到底从哪来埃氏筛是另一个调和级数思想的经典案例。你要筛出 1 到 n 的所有素数做法是从 2 开始把每个素数的倍数都标记为合数。表面上的循环是is_prime [True] * (n 1) is_prime[0] is_prime[1] False for i in range(2, int(n ** 0.5) 1): if is_prime[i]: for j in range(i * i, n 1, i): is_prime[j] False内层循环的次数是 (n - i²) / i对所有素数 i 求和。这里面用到一个数论结论不超过 n 的素数倒数和大约是 ln ln n 加一个常数。把 n 提出来总标记次数就是 n × ln ln n 的量级所以复杂度写作 O(n log log n)。这个 log log n 非常反直觉。ln ln 10^12 大约是 3.3也就是说 n 到了一万亿这个复杂度的“隐常数”也才三倍多。这也是为什么埃氏筛在实际跑起来快到离谱。如果你只看“两层循环”就认定它是 O(n√n)那你根本没有理解筛法的本质。理解调和级数之后你会发现所谓的“线性筛”和“埃氏筛”之间的差距并没有你想象的那么大。2.3 一段能直接跑的调和级数代码回到约数统计问题最直接的写法就是用倍数枚举统计每个数的约数个数def divisor_count_sum(n): cnt [0] * (n 1) for d in range(1, n 1): for multiple in range(d, n 1, d): cnt[multiple] 1 return sum(cnt[1:])这个代码的时间复杂度是 O(n log n)空间是 O(n)。n 取 10^7 时内层总执行次数约等于 1.6×10^8在 C 里可以在一秒多跑完Python 会慢一些。如果你想进一步优化到 O(n)可以用线性筛维护每个数的最小素因子然后递推约数个数但那是另一种思路不是这篇文章的重点。这里有一个工程层面的坑上面的代码虽然复杂度可以接受但内存访问是跳跃的。尤其是 n 很大的时候每次 cnt[multiple] 1 都会触碰到不同的 cache line实际耗时可能比理论复杂度高不少。如果你在性能敏感场景里写这种代码可以考虑分块处理或者用数组填充配合连续访问的写法。不过在算法题里n ≤ 10^7 这个量级通常不需要为了 cache 优化折腾知道有这回事就行。2.4 记忆点看到 n/i 的求和第一反应是 H_n我自己用调和级数最多的场景不是约数统计而是分析各种“均摊复杂度”。比如分块处理、跳跃表更新、并查集按秩合并的势能分析背后都有调和级数的影子。一个特别有用的直觉是ln n 的值远比你想象的小。ln 1000 ≈ 6.9ln 10^6 ≈ 13.8ln 10^9 ≈ 20.7。所以 n log n 在 n10^6 时只是 1400 万次操作这在绝大多数评测环境下都能接受。以后看到形如“for i in 1..n: for j in step i:”的代码不要急着喊 O(n²)。先在草稿纸上写一个 Σ n/i再想想它等于 nH_n。光是这一个习惯就能帮你少写很多“过度优化”的代码。3. 斯特林公式阶乘量级的“目测尺”3.1 为什么需要一把尺子阶乘是组合数学和概率论里最常见的量但它的增长速度非常反直觉。你背过 20! 2432902008176640000大约 2.4×10^18已经超过 64 位整数的上限100! 更是大到 9.33×10^157。在算法题里阶乘经常不是作为最终结果出现而是作为中间量参与概率估算、组合数比较、期望计算。如果每一步都精确计算很快就会溢出或者超时。这时候你需要一把“目测尺”。斯特林公式就是这把尺子n! ≈ √(2πn) × (n / e)^n它可以精确到什么程度相对误差大约是 1/(12n)。n10 时误差约 0.8%n100 时误差约 0.08%n1000 时误差约 0.008%。如果你只是要量级这个精度完全够用。3.2 公式长什么样误差有多大斯特林公式最常见的写法是n! ≈ √(2πn) (n/e)^n更精确的版本还包含一个修正项n! ≈ √(2πn) (n/e)^n × (1 1/(12n) 1/(288n²) - ...)在算法题里通常只需要取前三项的对数形式就够了。为什么因为很多题目要的不是阶乘本身而是“多少位”或者“概率是否大于某个阈值”这时候对数域是唯一不会溢出的选择。判断一个数的十进制位数本质是求 log10然后取整数加一。结合斯特林公式log10(n!) n log10 n - n log10 e 0.5 log10(2πn)我来演示一个实际计算100! 的位数。log10 100 2n log10 n 200log10 e ≈ 0.4343n log10 e ≈ 43.430.5 log10(2π×100) ≈ 0.5 × 2.798 ≈ 1.399合计 ≈ 200 - 43.43 1.399 ≈ 157.97位数 158而 100! 真实值就是 158 位数。你看不用算完整乘积两三行心算就能确定位数。这种能力在估算复杂度、估算搜索空间时特别管用。3.3 用对数形式估算位数和概率来看一个稍微复杂一点的例子估算从 52 张牌里随机抽 13 张的组合数 C(52,13) 的量级。先用对数形式展开log10 C(52,13) log10(52!) - log10(13!) - log10(39!)用斯特林公式分别算得到的结果约等于 11.6 亿不对实际上 C(52,13) ≈ 6.35×10^11log10 约 11.80。用斯特林公式估算误差在 0.1 以内完全够用。另一个典型场景是“生日悖论”。要估算 50 个人里面存在两个人生日相同的概率可以算“都不相同”的概率P (365/365)(364/365)...(316/365)。这个乘积直接用浮点算会非常小容易下溢但取 log 之后就是一个很漂亮的求和。用 log 方法算出的概率大约 0.97和精确值几乎一致。3.4 代码实现与精度陷阱在工程代码里最好不要手写斯特林公式的浮点部分直接用标准库的对数伽马函数最稳。Python 里是math.lgammaC 里是std::lgamma。import math def log10_fact(n): # math.lgamma(n1) 返回 ln(n!) return math.lgamma(n 1) / math.log(10) def approx_choose_log10(n, k): return log10_fact(n) - log10_fact(k) - log10_fact(n - k)这里有一个很容易踩的坑lgamma返回的是自然对数不是常用对数。如果你想直接得到十进制位数记得除以ln(10) ≈ 2.302585。另外lgamma对非常小的 n 也有定义但 n0 或 n1 时直接返回 0不需要走近似公式。还有一个很多人忽略的点斯特林公式是“近似”不是“精确”。如果你拿它去求一个精度要求达到 1e-15 的浮点结果那你会被误差坑到。它的定位是量级估算和对数域计算不是替代精确计算。在算法题里需要精确阶乘取模时应该预处理阶乘表或使用质因数分解而不是用斯特林公式硬凑。4. 费马小定理模运算里的“除法替身”4.1 取模世界里的除法为什么不能直接做在整数世界里a/b 就是 a 乘以 b 的倒数。但在取模世界里分数没有直接意义。比如在模 7 意义下2/3 到底是多少你没法直接说 0.666...因为模算术只接受整数。正确的做法是找一个整数 x使得 3x ≡ 2 (mod 7)。如果模数是素数那么基本上每个非零元素都有乘法逆元除法就能变成“乘以逆元”。问题变成了怎么高效求一个数在模 p 意义下的逆元最常用的方法之一就是费马小定理。4.2 费马小定理怎么来的费马小定理说的是如果 p 是素数a 是任意不被 p 整除的整数那么a^(p-1) ≡ 1 (mod p)这个定理的证明非常漂亮。考虑集合 {a, 2a, 3a, ..., (p-1)a}在模 p 意义下这些数恰好是 1 到 p-1 的一个排列。把它们全部乘起来得到a^(p-1) × (p-1)! ≡ (p-1)! (mod p)因为 (p-1)! 在模 p 下不等于 0可以约掉就得到了 a^(p-1) ≡ 1。利用这个结论a × a^(p-2) ≡ a^(p-1) ≡ 1所以 a 的逆元就是 a^(p-2) mod p。这是求逆元最简单的方法因为它不需要任何前置知识只要会快速幂。4.3 快速幂求逆元代码与边界快速幂的实现几乎是模板级的def mod_pow(a, e, p): a % p r 1 while e: if e 1: r r * a % p a a * a % p e 1 return r def mod_inv(a, p): # 要求 p 是素数且 a 和 p 互质 return mod_pow(a, p - 2, p)这个代码在 p 是素数且 a % p ! 0 的时候是正确且高效的复杂度是 O(log p)。但在 C 里写的时候要注意乘法溢出p 如果是 1e97两个数相乘大约是 1e18刚好卡在 64 位整数边界内但如果直接用 int 相乘一定溢出。我的习惯是直接用long long保守一点可以用__int128做中间乘法。有一个隐藏很深的条件费马小定理要求 a 不被 p 整除。如果 a % p 0那么 a 在模 p 下没有逆元因为这个数乘以任何整数模 p 都等于 0不可能等于 1。很多题会故意在这种边界上设坑尤其是当输入数据量大、你没注意取模时。4.4 实际场景组合数取模和哈希设计费马小定理最常见的工程场景是组合数取模。要计算 C(n,k) mod p其中 p 是 1e97 这种大素数流程是预处理阶乘表和阶乘逆元表然后用公式 fac[n] * inv_fac[k] % p * inv_fac[n-k] % p。求 inv_fac 的时候不需要对每个阶乘都做一次快速幂只需要对最大的那个做一次快速幂然后逆推fac [1] * (n 1) for i in range(1, n 1): fac[i] fac[i - 1] * i % p inv_fac [1] * (n 1) inv_fac[n] mod_inv(fac[n], p) for i in range(n, 0, -1): inv_fac[i - 1] inv_fac[i] * i % p原理是 fac[i-1] fac[i] / i所以逆元也满足 inv_fac[i-1] inv_fac[i] * i。除了组合数取模费马小定理还常用于哈希表设计中的“除法哈希”和概率算法的随机数生成。它和扩展欧几里得是两种并列的求逆元方案模数是素数时选费马模数是合数时用扩展欧几里得。不要混着来。5. 贝叶斯公式被“最大似然”挤出视野的先验5.1 一个反直觉的概率题先来做一道选择题某种罕见病在人群中的发病率是 0.1%。有一种检测手段如果一个人真的患病检测阳性的概率是 99%如果一个人没患病检测阳性的概率是 1%假阳性率。现在有一个人的检测结果是阳性请问他真正患病的概率是多少一半以上的人会脱口而出“99%”因为检测准确率听起来很高。但实际上答案是大约 9%。计算过程就是贝叶斯公式P(患病 | 阳性) P(阳性 | 患病) × P(患病) / P(阳性)其中 P(阳性) 要展开成两部分患病且阳性 未患病且假阳性。也就是0.99 × 0.001 / (0.99 × 0.001 0.01 × 0.999) ≈ 0.090看到没一个 99% 准确的检测阳性后真实患病概率只有 9%。原因是发病率太低了假阳性人数远远超过真实患者人数。这道题在算法面试里出现过不止一次它考的根本不是概率计算能力而是你有没有“先验”这个概念。5.2 贝叶斯公式拆解先验、似然、后验贝叶斯公式的标准形式是P(A|B) P(B|A) × P(A) / P(B)三个关键角色P(A) 是“先验概率”指在看到数据 B 之前你对事件 A 的判断。P(B|A) 是“似然”指在 A 为真的前提下观察到 B 的概率。P(A|B) 是“后验概率”指观察到 B 之后你对事件 A 的更新判断。很多机器学习初学者只盯着似然 P(B|A)认为“哪个假设能让数据出现的概率最大就选哪个”。这就是最大似然估计。但最大似然有一个盲区它忽略了先验。如果先验非常低即使似然很高后验也可能很低。正则化在数学上就是给参数加了一个先验分布比如 L2 正则等价于假设参数服从均值为 0 的高斯分布。很多博客讲正则化只知道“惩罚复杂度”不知道它其实是贝叶斯框架的自然推论。这就是被忽视的地方。5.3 朴素贝叶斯分类的最小实现贝叶斯公式在工程里最常见的落地是朴素贝叶斯分类器。它的核心假设是在给定类别 c 的条件下各个特征之间相互独立。于是P(c | x) ∝ P(c) × Π P(x_i | c)以文本分类为例特征是“单词是否出现”类别是“垃圾邮件”和“正常邮件”。训练时统计每个词在两个类别里的出现概率以及两类的先验 P(垃圾) 和 P(正常)。预测时对一篇文章计算每个词在两类下的概率乘积再乘上先验。代码实现上有一个关键细节要用 log 加法和而不是概率连乘。因为一篇文章有成百上千个单词每个词概率都在 0.01 到 0.99 之间连乘之后很容易下溢成 0。取 log 之后乘法变成加法数值稳定性大幅提升。import math def predict(words, log_p_class, log_p_word, vocab): scores {} for c in log_p_class: s log_p_class[c] for w in words: s log_p_word[c].get(w, log_p_word[c][__default__]) scores[c] s return max(scores, keyscores.get)5.4 被忽视的三个工程陷阱第一个陷阱是零概率问题。如果某个词在训练集的某个类别里一次都没出现概率算出来是 0连乘后整个后验变成 0。解决方案是拉普拉斯平滑分子加 1分母加词汇表大小。第二个陷阱是特征独立性假设几乎永远不成立。真实文本里“免费”和“点击”高度相关但朴素贝叶斯假装它们没关系。即便如此它在很多场景下依然表现良好这种“错误假设但可用”的现象本身就是工程里值得习惯的。第三个陷阱是类不平衡。如果垃圾邮件只占 1%先验就会把后验压得很低。这时单纯的贝叶斯分类阈值可能需要校准而不是机械地选择后验最大的类别。如果你在算法题里看到概率题不要急着列树状图。先写一个贝叶斯公式把分母展开把先验写进去很多题的答案会自动浮现出来。6. 拉格朗日插值给定点太少时如何让曲线“合理地”穿过它们6.1 插值不是拟合数学给了一个精确解拉格朗日插值解决的是这样一个问题给定 k1 个点 (x_0, y_0), (x_1, y_1), ..., (x_k, y_k)并且所有 x_i 互不相同我们要找一个不超过 k 次的多项式让它在每个给定点上都取到对应的 y 值。这和“拟合”不一样拟合是找一条近似贴近数据的曲线插值是要求严格穿过所有已知点。为什么算法里需要这东西因为有些题目会故意只给你很少的信息比如给出某个函数在 0, 1, 2, ..., k 这些点上的取值让你求这个函数在 n 处的值并且 n 可以非常大大到没法直接递推。如果一个函数本质上是多项式那你只需要 k1 个点就能完全确定它然后直接算 n 处的值。这就是拉格朗日插值的威力。6.2 公式的直觉与 O(k²) 实现拉格朗日插值的公式长这样L(x) Σ_{i0}^{k} y_i × Π_{j≠i} (x - x_j) / (x_i - x_j)看着吓人但直觉非常朴素。对每一个 i构造一个“只有在 x_i 处等于 1在其他所有已知点处等于 0”的基函数然后把它乘上 y_i再加起来。这样在任意一个已知点 x_i 处只有第 i 项起作用其他项都是 0所以多项式必然精确穿过所有已知点。一个最直观的实现是两重循环def interpolate(xs, ys, x): total 0 k len(xs) for i in range(k): term ys[i] for j in range(k): if i ! j: term * (x - xs[j]) / (xs[i] - xs[j]) total term return total复杂度是 O(k²)。对大多数题目 k 在几千以内这个复杂度完全能接受。6.3 优化到 O(k)前缀后缀积技巧如果你要在大量 x 上反复求值O(k²) 可能不够。优化方法也不复杂预先算好分母 denom[i] Π_{j≠i} (x_i - x_j)以及分子中 x 相关的部分。关键是分子可以拆成“前缀积 × 后缀积”pref[j] Π_{t0}^{j-1} (x - x_t)suf[j] Π_{tj1}^{k-1} (x - x_t)那么第 i 项分子就是 pref[i] × suf[i]。有了前缀后缀积计算一次 L(x) 的复杂度从 O(k²) 降到 O(k)。这里还有一个很大的坑当 x 恰好等于某个已知点 x_i 时分子分母都变成 0直接套公式会得到 0/0。正确的做法是先判断 x 是否等于某个 xs[i]如果是直接返回 ys[i]。很多人在这个 corner case 上翻车。6.4 高次插值的坑多项式很美但别乱用拉格朗日插值在数学上是精确解但在浮点世界里要小心。高次多项式本身很平滑但当你用它在区间外做外推时结果可能剧烈震荡。经典的反例是龙格现象对一个均匀取点的函数做高次插值次数越高在端点附近的误差反而越大。原因不是插值公式错了而是高次多项式“太灵活”为了穿过所有点它可以剧烈弯曲。工程上如果你处理的是实验数据不要用高次插值去猜测两点之间的曲线。更好的做法是分段低次插值比如三次样条。算法题里之所以敢用高次插值是因为很多时候取模运算让一切变成精确的整数算术不存在浮点误差问题多项式次数再高也不怕。这一点要分清楚数学公式是好工具但用错了场景再优美的公式也会变成灾难。自然数幂和是一个很好的正面例子。S_k(n) 1^k 2^k ... n^k 是关于 n 的 k1 次多项式所以只需要算出前 k2 个点然后用拉格朗日插值求 S_k(n) 的模值。即使 n 到了 10^18只要 k 不大插值都能在 O(k) 时间内算完。这种题用朴素循环肯定超时用贝努利数又要背公式拉格朗日插值是最容易理解也最容易写对的路线。最后分享一个我自己的习惯。学算法这几年我学到最值钱的一件事就是遇到题先问自己有没有数学结构可以依赖再想数据结构怎么配合。遇到复杂度看不懂的循环写个 n1000 的样例数一下执行次数再和 H_n、ln n 对一对遇到取模除法第一反应是问模数是不是素数而不是闷头实现一个除法函数遇到概率题先写 P(A|B) ...把分母展开把先验写进去遇到给点求值先想拉格朗日插值能不能派上用场。这几个动作几乎成了我的条件反射。顺便说一句手边可以存几个常用常数ln 2 ≈ 0.693欧拉常数 γ ≈ 0.577log10 e ≈ 0.4343ln 10 ≈ 2.303。别小看它们很多复杂度估算和手算题最后都落到这几个数上。公式本身不难难的是你愿不愿意在写循环之前先停下来做一分钟数学。