前阵子有人拿一段统计脚本问我,为什么一个只算几十万项的循环跑得那么慢。我打开一看,里面把1 + 1/2 + 1/4 + ...这种等比级数用 for 循环一项一项往上垒,还一直加到了 10 的负十几次方才肯停。改成闭式公式之后,同样的结果,时间从几秒掉到微秒级。这件事让我意识到,级数求和这件事看着是中学课本里的老内容,但真正在工程和算法场景里把它用对、用快、用准的人并不多。大部分人记住的是"等差用首项加末项乘项数除以二""等比用一减公比分之……",可真拿到一个形如∑ k·x^k或者∑ 1/(k(k+2))的式子,就卡在那儿了。
这篇内容就是把我这些年整理的级数求和八个常用公式和配套的级数求和法完整摊开讲一遍。不只是列公式,更重要的是讲清楚每个公式为什么成立、在什么条件下才能用、用错会出什么问题。适合三类人看:正在准备考试需要把公式串起来的学生、写代码时遇到累加性能问题的开发者、做数据分析或数值计算需要估算截断误差的从业者。哪怕你只记得住两三个公式,把判断顺序和边界条件搞明白,收益也比死记八个公式大得多。
1. 八个公式各管哪一摊,先给一张全景图
很多人学级数求和的方式是"遇到一个记一个",结果是记了十几个公式,用的时候还是不知道该掏哪一个。真正有效的组织方式是按照式子的结构特征分类,因为求和法的本质就是"识别结构 → 匹配变形手段"。结构认对了,公式自己就浮出来了。
1.1 八个公式的分工,按结构而不是按难度排
我习惯把这八个公式分成三组来看。第一组是线性结构组,包括等差数列求和、等比数列求和、无穷等比级数求和,这三条处理的是项与项之间成固定比例或固定差的关系,是使用频率最高的三条。第二组是幂和组,包括自然数平方和、自然数立方和,它们处理的是∑ k^p这种幂次累加,往上还能递推到四次、五次。第三组是变形技巧组,包括裂项相消、差比数列(等差乘等比)求和、p 级数判敛与巴塞尔和,这一组不是"套公式",而是"套手法",需要你先对通项做一次改造。
把八个公式列成清单是这样的:
| 编号 | 名称 | 公式 | 关键前提 |
|---|---|---|---|
| 1 | 等差数列求和 | S = n(a₁ + aₙ)/2 | 相邻项差恒定 |
| 2 | 等比数列求和 | S = a₁(1 − qⁿ)/(1 − q) | q ≠ 1,q = 1 时 S = n·a₁ |
| 3 | 无穷等比级数 | S = a₁/(1 − q) | 必须严格满足 |q| < 1 |
| 4 | 自然数平方和 | ∑k² = n(n+1)(2n+1)/6 | n 为正整数 |
| 5 | 自然数立方和 | ∑k³ = [n(n+1)/2]² | 结果恰为前 n 项和的平方 |
| 6 | 裂项相消 | ∑1/[k(k+1)] = 1 − 1/(n+1) | 通项能拆成两个相邻结构的差 |
| 7 | 差比数列求和 | ∑k·x^k = x[1 − (n+1)xⁿ + nxⁿ⁺¹]/(1−x)² | x ≠ 1,通常 |x| < 1 |
| 8 | p 级数与巴塞尔和 | ∑1/k² = π²/6;∑1/k^p 在 p > 1 时收敛 | 判敛看 p 与 1 的大小关系 |
这张表看起来简单,但每一行的"关键前提"那一列才是真正值钱的地方。我见过太多次把第 3 条用在q = 2上的情况,也见过有人拿第 4 条去算∑(2k+1)²却忘了先把括号展开。
1.2 三条最容易踩的边界,提前说清楚
第一条边界是收敛性。有限项求和只要项数确定,怎么算都有结果;但一旦上限变成无穷,结果是否存在就完全取决于通项衰减的速度。第 3 条、第 8 条都属于无穷级数,用之前必须先确认收敛。一个简单的判断经验:通项趋于零是必要条件,但不充分——调和级数∑1/k的通项趋近于零,可它的和是发散的,部分和大约以ln n + 0.5772的速度增长。这个反例我建议每个人都记牢,它是区分"会套公式"和"真懂级数"的分水岭。
第二条边界是公比为 1 的分支。第 2 条公式分母是1 − q,q = 1时直接除零。很多人在写代码时把这条公式封成一个函数,忘了加这个分支,结果程序在某些输入下直接抛异常或者返回 inf。正确做法是在函数入口就判断,q == 1时返回n * a1。
第三条边界是整数运算的中间溢出。第 4 条n(n+1)(2n+1)/6,如果 n 取到 10⁶ 量级,n(n+1)(2n+1)大约是 2×10¹⁸,已经超出 32 位有符号整数的范围。在 C、Java 这类语言里,必须用 64 位整型,或者调整运算顺序先做除法,这一点后面第 5 章会展开。
记住一件事:公式的适用条件不是注释里的客套话,而是能不能用的开关。先验条件,再代数字。
2. 等差与等比:两条地基公式的推导脉络和失效点
这两条公式几乎所有人在中学就背过,但真正能讲清楚"为什么是这个形式"的人比例并不高。而理解推导过程的价值在于,当你遇到变体(比如每隔两项取一项、或者带系数的等差)时,能自己推出来而不是到处搜。
2.1 等差数列:倒序相加消掉的是变量而不是数字
设首项 a₁,末项 aₙ,项数 n,公差 d。把和写成两遍,一遍正序一遍倒序:
S = a₁ + (a₁+d) + ... + (a₁+(n-1)d) S = aₙ + (aₙ-d) + ... + (aₙ-(n-1)d)两式逐项相加,每一对的和都是a₁ + aₙ,一共有 n 对,于是2S = n(a₁ + aₙ)。这个手法的精妙之处在于:它并没有去求每一项的具体值,而是利用"对称位置的两项之和恒定"这个性质,把 n 个未知量压缩成了一个常量。这也是整个级数求和领域的核心思想——寻找结构中的不变量,用它替换掉逐项累加。
把这个思路稍微推广一下,就能处理"等差但隔项"的情况。比如求1 + 3 + 5 + ... + (2n−1),首项 1,末项 2n−1,项数 n,代入得n(1 + 2n − 1)/2 = n²。这个结果非常有用:前 n 个奇数之和恰好是 n 的平方。我在做复杂度分析时经常用到它,比如某个循环的内层执行次数是奇数序列,一眼就能写出总量是 n²。
需要注意的是项数的计算。很多人在这里出错,比如求3 + 7 + 11 + ... + 99,公差是 4,项数不是(99−3)/4而是(99−3)/4 + 1 = 25。少加这个 1 是最高频的错误之一,我个人的习惯是写完项数之后回头验算一次末项:a₁ + (n−1)d是不是等于给定的末项,如果不等,说明项数算错了。
2.2 等比数列:错位相减与那个不能忘的 q = 1 分支
等比数列求和的标准推导是错位相减法。设S = a₁ + a₁q + a₁q² + ... + a₁qⁿ⁻¹,两边同乘 q 得qS = a₁q + a₁q² + ... + a₁qⁿ,两式相减,中间所有项两两抵消,剩下(1−q)S = a₁ − a₁qⁿ,于是得到公式。
这个推导里有一个隐含的假设:相减的时候我们没有除以零,也就是说q ≠ 1。当q = 1时,每一行都相同,减出来是零等于零,推不出任何信息。此时必须单独处理,S = n·a₁。这不是数学上的"特例补丁",而是两种完全不同的结构——公比为 1 的等比数列其实就是常数数列。
我在实际写代码时,会把这两个分支写在同一个函数里:
def geo_sum(a1, q, n): if n <= 0: return 0.0 if abs(q - 1.0) < 1e-12: # 浮点比较不能直接用 == return a1 * n return a1 * (1 - q ** n) / (1 - q)这里abs(q - 1.0) < 1e-12而不是q == 1.0,是因为浮点数在多次运算后很难精确等于 1。这一点在纯数学题里无所谓,但在工程代码里会真实地导致除零或者巨大误差。我曾经处理过一个信号衰减的计算,参数是从配置文件读的浮点数,理论值应该是 0.5,但因为前面做了一次归一化,实际值变成了 0.49999999999999994,q ** n在 n 很大时下溢成 0,结果反倒正常;但如果那个参数被误设成 1.0000000001,程序就会得到一个荒谬的巨大值而不报错。这类"静默错误"比直接崩溃更危险。
另外提一句,第 2 条公式在|q| < 1且 n 很大时,qⁿ会趋近于零,此时可以直接用第 3 条的无穷版本a₁/(1−q)。这个近似什么时候能用?经验法则是看qⁿ相对于结果是否能忽略。比如 q = 0.5,n = 50 时qⁿ ≈ 8.9×10⁻¹⁶,已经低于双精度浮点的有效位数,可以直接丢弃。但 q = 0.99,n = 50 时qⁿ ≈ 0.605,跟 1 是同量级的,绝对不能近似。判断标准不能拍脑袋,得实际算一下。
3. 幂和公式:平方和、立方和长成那样的原因
∑k² = n(n+1)(2n+1)/6这条公式,绝大多数人是背下来的,很少有人知道它是怎么来的。但如果你能理解它的推导路径,就能自己推出∑k³、∑k⁴,甚至写代码自动生成任意次幂的求和公式。这才是"级数求和法"里"法"字的含义——掌握方法比掌握结论重要得多。
3.1 平方和:差分法把三次方拆成二次方
推导的核心恒等式是(k+1)³ − k³ = 3k² + 3k + 1。对这个式子从 k = 1 加到 k = n,左边是一个裂项相消的结构,所有中间项抵消,只剩(n+1)³ − 1;右边拆成三项:3∑k² + 3∑k + n。于是:
(n+1)³ − 1 = 3∑k² + 3·n(n+1)/2 + n解出∑k²就得到n(n+1)(2n+1)/6。整个过程只用了两个已知结论:裂项相消和等差数列求和。这就是差分法的通用套路——用高一次的幂的差分,反解出低一次幂的累加。
把 k 从 1 换成 0 或者从 2 开始,结果会变,但变的只是边界项。我在实际计算中遇到过不少"从 0 开始"的场景,比如数组下标从 0 计数的累加,此时∑_{k=0}^{n} k²和∑_{k=1}^{n} k²是一样的,因为 0² = 0,不影响。但如果是从 2 开始,就得减去 1。这个细节听起来很笨,但我在 code review 时见过至少三次因为下标起点搞错导致统计结果偏小的情况。
顺带说一个容易被忽略的用途:平方和公式可以用来估算循环的总执行次数。比如一个双层循环,外层 i 从 1 到 n,内层 j 从 1 到 i,那么内层总执行次数就是∑i = n(n+1)/2,是 O(n²) 级别。如果是三层嵌套且每层上限递增,总次数就是∑i² = n(n+1)(2n+1)/6,仍是 O(n³) 但系数是 1/3。做性能预估的时候,这个精确系数比大 O 记号有用得多。
3.2 立方和与前 n 项和的平方,这个巧合值得琢磨
∑k³ = [n(n+1)/2]²这个结果第一次看到会觉得很神奇——立方和竟然等于和的平方。它并不是巧合,而是有一个很直观的几何解释:把 1³ 看成边长为 1 的立方体,2³ 看成 8 个边长为 1 的小立方体,整个构造可以拼成一个边长为n(n+1)/2的大正方形。这个"无字证明"我第一次看到时就记住了,比任何代数推导都牢固。
顺着差分法的思路,∑k³可以从(k+1)⁴ − k⁴ = 4k³ + 6k² + 4k + 1推出来。推的过程中会再次用到∑k²和∑k的结果,这形成了一个递推链条:知道∑k就能推∑k²,知道∑k²就能推∑k³,以此类推。这就引出了下面这个话题。
3.3 更高次幂和:递推公式与什么时候该上程序
一般地,∑k^p是一个关于 n 的 p+1 次多项式,这个结论可以用数学归纳法证明,也是伯努利数的来源。手推 p = 4 还行,∑k⁴ = n(n+1)(2n+1)(3n²+3n−1)/30,再往上系数就变得很丑陋了。这时候更实际的做法是写一个小脚本用符号计算自动生成,而不是继续背。
# 概念示意:用差分法递推求解 ∑k^p 的多项式 # 实际工程中通常直接调用符号计算库 import sympy as sp n = sp.symbols('n', positive=True) k = sp.symbols('k', integer=True, positive=True) for p in range(1, 6): expr = sp.summation(k**p, (k, 1, n)) print(f"p={p}:", sp.factor(expr))跑出来的结果可以直接抄进文档里。我在做数值积分误差估计的时候就这么干过——需要知道某个离散化方案的高阶项系数,手推太慢且容易错,直接生成然后检查几个 n 值验证一下就完事。
需要强调的是,∑k^p的多项式性质只对正整数幂成立。一旦 p 变成负数,比如∑1/k或∑1/k²,性质完全变了:前者发散,后者收敛到 π²/6,都不是多项式。这个分界线必须划清楚。
4. 裂项相消与无穷级数:把求和转化为相减的手法
前三条公式处理的是"结构已知"的数列,而现实中大量遇到的式子结构并不规整,比如1/(1×2) + 1/(2×3) + ...。这类问题的解法不是套公式,而是改造通项。裂项相消就是最经典的一种改造,它把加法问题转换成减法问题,然后再用抵消来降维。
4.1 裂项拆分的三类常见套路
第一类是相邻因子差。1/[k(k+1)] = 1/k − 1/(k+1),从 1 加到 n 之后,中间所有项正负抵消,只剩下1 − 1/(n+1)。这个结果有个很好的性质:当 n 趋于无穷时,和趋近于 1。也就是说无穷级数∑1/[k(k+1)]收敛到 1,而且收敛速度是 1/n 级别,不算快。
第二类是间隔因子。1/[k(k+d)] = (1/d)·(1/k − 1/(k+d))。这里的关键是那个1/d的系数,很多人会漏掉。比如∑1/[k(k+2)],正确答案是(1/2)(1 + 1/2 − 1/(n+1) − 1/(n+2)),而不是1 − 1/(n+2)。我见过太多人在这里出错,因为 d = 1 的情况太典型了,形成了思维惯性。
第三类是根式有理化。1/(√k + √(k+1)) = √(k+1) − √k,从 1 加到 n 之后中间抵消,得到√(n+1) − 1。这类题在竞赛里出现的频率很高,识别特征是分母里有两个相邻的根号相加。
把这三类整理成一张对照表,用的时候照着找结构:
| 通项形态 | 拆分结果 | 求和后剩余项 |
|---|---|---|
| 1/[k(k+1)] | 1/k − 1/(k+1) | 首项减末项 |
| 1/[k(k+d)] | (1/d)[1/k − 1/(k+d)] | 前 d 项减后 d 项 |
| 1/(√k + √(k+1)) | √(k+1) − √k | √(n+1) − √1 |
| 1/[(2k−1)(2k+1)] | (1/2)[1/(2k−1) − 1/(2k+1)] | 首项减末项,步长为 2 |
最后一行那种步长为 2 的情况要特别注意:抵消之后剩下的不是一项,而是首尾各一项的一半,容易算错。
4.2 无穷等比级数:收敛条件是硬门槛不是建议
第 3 条公式S = a₁/(1−q)是我认为最容易被滥用的公式。它的前提是|q| < 1,这个条件不是"最好满足",而是"不满足就完全错"。举个反例:1 + 2 + 4 + 8 + ...,公比是 2,如果硬套公式会得到1/(1−2) = −1。一个全为正数的级数求和得到负数,这显然是荒谬的,但因为公式算起来太顺,很多人不会去回头检查合理性。
|q| = 1时同样不能用。q = 1时是常数数列,和发散;q = −1时是1 − 1 + 1 − 1 + ...,部分和在 0 和 1 之间振荡,按柯西意义也不收敛(虽然在某些求和意义下可以赋予 1/2 的值,但那属于完全不同的理论框架,工程里不要这么用)。
那么怎么快速判断|q| ≥ 1的等比级数该给什么答案?答案是"发散",没有有限和。在程序里遇到这种情况,应当抛出异常或者返回一个明确的标记,而不是返回一个数字。这一点在写数值库时尤其重要——错误必须大声喊出来,不能悄悄返回一个错误答案。
顺带提一个高频变体:∑ k·x^k型(等差乘等比),这就是第 7 条公式。它的推导同样用错位相减法,只不过乘完 x 相减之后,剩下的不是常数项而是另一个等比数列。结果形式是x[1 − (n+1)xⁿ + nxⁿ⁺¹]/(1−x)²。当|x| < 1且 n 趋于无穷时,后两项都趋于零,结果简化为x/(1−x)²。这个无穷版本在概率论里会用到,计算几何分布的期望时就是它。
4.3 p 级数与巴塞尔问题:先判敛再谈计算
∑1/n^p这类级数的收敛性完全由 p 决定:p > 1 收敛,p ≤ 1 发散。这个结论非常干脆,也极其重要,因为它是很多判敛法的参照标准。判断一个复杂级数的收敛性时,常用的思路就是找一个 p 级数做比较,看谁衰减得快。
最有名的特例是 p = 2,也就是巴塞尔问题,∑1/n² = π²/6 ≈ 1.644934。这个结果为什么会出现 π,直观上很难想象,因为式子全是整数倒数。欧拉的证明思路是把正弦函数展开成无穷乘积,再把泰勒展开的系数对上,这个过程需要用一点技巧,我这里不展开。但值得记住的是结论本身:收敛速度是 n 的倒数级别,要精确到 6 位有效数字,大概需要加到一百万项。这个数字在写代码估算截断误差时非常有用。
至于 p = 1 的调和级数,很多人误以为它收敛,因为通项确实趋近于零。但它的部分和Hₙ ≈ ln n + γ,其中 γ 是欧拉-马歇罗尼常数约 0.5772。n = 10⁶ 时Hₙ ≈ 14.39,还在缓慢增长。也就是说,这个级数的和可以是任意大的正数,只是增长得极其缓慢。在工程中如果需要判断"加到第多少项才超过某个阈值",这个近似式可以直接给出答案,不用真的循环。
5. 落到代码:误差、判敛、溢出三座山
把公式写在纸上是一回事,让它在程序里稳定工作完全是另一回事。我见过太多代码,数学推导完全正确,但实现出来结果差了百分之几甚至直接崩掉。这一章就讲实打实的问题。
5.1 整数型求和的溢出与类型选择
前面提过,n(n+1)(2n+1)的增长是立方级的。n = 10⁴ 时大约是 2×10¹²,已经超过 32 位有符号整型的上限约 2.1×10⁹。n = 10⁶ 时约 2×10¹⁸,还在 64 位上限约 9.2×10¹⁸ 之内,但已经不远了。n = 10⁷ 时就会溢出 64 位。
防御的办法有两种。一是用大整数类型,比如 Python 的 int 天然支持任意精度,Java 用 BigInteger,C++ 用__int128或者自己实现。二是调整运算顺序提前约分。后一种更快也更省资源,具体做法是先用 6 去除 n、n+1、2n+1 中能被整除的那个。
def sum_of_squares(n: int) -> int: # 提前约分,减小中间值,Python 大整数下也更快 a, b, c = n, n + 1, 2 * n + 1 for factor in (2, 3): for i in range(3): if factor == 2 and [a, b, c][i] % 2 == 0: [a, b, c][i] //= 2 break if factor == 3 and [a, b, c][i] % 3 == 0: [a, b, c][i] //= 3 break return a * b * c这段代码写得有点绕,实际用的时候可以写得更直白一些。核心思想就是:6 = 2 × 3,三个连续整数的乘积一定能被 6 整除,所以一定存在一个因子可以约掉,先约再乘能显著降低溢出的风险。
关于浮点数和整数的选择,我的经验是:只要结果在整数范围内有意义,优先用整数算,最后再转浮点。因为整数运算是精确的,浮点运算每一步都可能引入误差。
5.2 无穷级数截断:停止条件怎么定才靠谱
计算无穷级数只能截断,问题是加到第几项才算"够"。这个判断不能凭感觉,得有误差上界。
对交错级数,有莱布尼茨定理:截断误差小于第一个被舍弃的项的绝对值。比如算∑(−1)ⁿ/(2n+1) = π/4,加到第 n 项时误差不超过1/(2n+3)。要 10⁻⁶ 精度,解1/(2n+3) < 10⁻⁶,n 大约 50 万。这个判据非常干净,可以直接写进循环条件。
对正项级数,误差估计需要用到积分判别法。如果通项是单调递减的,那么尾巴∑_{k=n+1}^{∞} a_k大约在∫_{n}^{∞} f(x)dx和∫_{n+1}^{∞} f(x)dx之间。比如∑1/k²,尾巴大约是 1/n,要 10⁻⁶ 精度,n 需要 10⁶。这跟前面说的结果一致。
实践中还有一个更简单粗暴的做法:设定一个最小增量阈值,当新加的项小于阈值时停止。这个做法在工作上没问题,但风险是如果收敛极慢,可能加了很久都到不了阈值;如果项有正有负,还可能误判。所以我一般会同时设置一个最大迭代次数上限做保护,避免死循环。
def sum_infinite(term_func, tol=1e-12, max_iter=10_000_000): total = 0.0 for k in range(1, max_iter + 1): t = term_func(k) total += t if abs(t) < tol: return total, k raise RuntimeError(f"未在 {max_iter} 次内收敛,当前和 {total}")这个函数的两个返回值我特意都返回了,因为迭代次数本身是重要信息——如果某个输入下迭代次数从几十跳到了几百万,说明那里可能有数值问题。
5.3 数值稳定性:当加法本身成为误差来源
有个反直觉的事实:大量浮点数相加,误差可能会累积得相当可观,即使每一项本身都很小。原因是浮点数的有效位数有限,当一个大的数和一个小很多的数相加时,小的那部分可能会被直接"吃掉"。比如1e16 + 1.0在双精度下结果仍然是1e16,因为 1.0 相对太小了。
在求和几百万项时,这个效应会累积。解决办法是Kahan 求和算法,用一次额外的减法把被舍入掉的部分暂存下来,下次加回去。代价是每项多两次浮点运算,收益通常是有效位数能多出一两位。
def kahan_sum(items): total = 0.0 comp = 0.0 # 补偿量:记录被舍入掉的部分 for x in items: y = x - comp t = total + y comp = (t - total) - y total = t return total还有一个更实用的技巧是改变求和顺序。对于收敛缓慢的正项级数,可以先加大的项再加小的项,或者用分块求和(把序列分成若干块先各自求和,再把块和加起来)。我在处理一个信号累加的场景时用过这个方法,同样一百万项,直接顺序加得到的结果和分块加差在第六位有效数字上,对于需要长期累积的场景,这个差距不能忽略。
一个经验判断:如果你的结果对第六位以后的数字敏感,就应该考虑用 Kahan 或者更高精度的算法;如果只关心量级,普通加法足够。
6. 实战选型:拿到一个求和式之后怎么走
前面五章讲了公式和实现,但真正落到具体问题上,最难的一步是"我该从哪儿下手"。这一章把我自己的判断流程写出来,还包括几个我真实踩过的坑。
6.1 我的判断顺序:判敛、定结构、选工具
第一步永远是看上限。上限是确定的数 n,还是无穷?如果是无穷,立刻进入判敛流程。判敛的快速路径是:先看通项是否趋于零,不趋于零直接发散;趋于零的话,找有没有 p 级数的结构或者能用比较判别法。
第二步是识别结构。把通项写出来,对照前面那张八个公式的表,看它落在哪一类。这里有个技巧:先看分子分母的次数差。如果通项能写成多项式的形式,大概率是幂和类;如果有明显的乘积因子在分母上,考虑裂项;如果通项里同时有 n 和指数,考虑差比数列。
第三步是判断能不能改造。有些式子乍看不属于任何一类,但稍微变形就归位了。比如∑(2k+1)²,展开成4k² + 4k + 1之后就变成了三个已知级数的组合。再比如∑k(k+1),展开成k² + k,同样是组合。展开再组合是我用得最多的手法,因为它把陌生的式子变成了熟悉式子的线性组合。
剩下的情况,也就是通项既不规则又不收敛快或没有闭式解,才是真的需要写循环逐项累加。这种时候反而简单了,重点转移到误差控制上。
6.2 几个我真实踩过的坑
第一个坑是把有限和公式用在无穷级数上而不检查收敛。我曾经写过一个计算几何分布期望的函数,直接用x/(1−x)²这个公式,结果调用方传了一个 x = 1.2 的参数进来,函数返回了一个负数。更糟的是这负数被塞进后续的统计里,直到最后的结果明显不合理才被发现。教训是:函数入口必须做参数校验,不满足前提就直接报错,而不是返回一个数字让错误传播下去。
第二个坑是浮点比较用了精确等于。前面提过 q = 1 的判断,我在早期的代码里写过if q == 1,在上百次测试里都没出问题,因为测试数据都是整数或者简单的分数。上线之后用户传进来一个通过计算得到的参数,理论上应该等于 1,但浮点运算后变成了 0.9999999999999999,程序走进了错误的分支,结果差了十万倍。这种 bug 极难复现,因为它在不同的机器和编译器下表现可能不同。
第三个坑是忽略了项的符号交替。有一次算一个交错级数,为了保证精度我设了很小的阈值 1e-15,结果迭代次数爆炸。后来才想明白,交错级数的项会逐渐变小但可能变得非常慢,而莱布尼茨定理告诉我们误差只跟下一项的大小有关,所以完全可以用更宽松的阈值。调整阈值之后速度提升了三个数量级,精度反而没有下降。理解误差的来源比无脑调小阈值有效得多。
第四个坑是把有限项的公式当成恒等式去化简。比如∑_{k=1}^{n} 1/[k(k+1)] = 1 − 1/(n+1),有人会误以为∑1/k加上∑1/(k+1)也能这样合并,实际上这两个都是发散的,不能这么操作。发散的级数不支持任意的加减重组,这是级数理论里一个必须守住的底线。
6.3 一份可以直接抄的检查清单
我在每次做涉及级数求和的实现时,会过一遍下面这几条。你也可以把它贴在自己的代码注释里:
- 上限是有限还是无穷?无穷的话收敛性确认了吗?
- 公式里的所有分母,有没有可能为零的分支?
- 中间结果的最大值是多少,会不会超出所用类型的范围?
- 浮点比较是否都用容差而不是精确相等?
- 如果用了截断,误差上界怎么算,有没有可验证的判据?
- 如果结果可以手工验证,我验算了前两三项吗?
最后一条尤其重要。任何级数求和的结果,取 n = 1、2、3 代进去都对得上,几乎可以排除大半的实现错误。这个习惯我是从一次惨痛的经历里学来的——当时花了一下午排查一个复杂的求和函数,最后发现问题在通项本身写错了,如果一开始就代 n = 1 验一下,五分钟就能发现。
我个人在实际使用中发现,这八个公式里真正高频的其实只有四五个,剩下的是在特定场景下才会用到的。所以如果你时间有限,优先级应该是:裂项相消 > 等比数列求和 > 平方和 > 无穷等比级数 > 其他。裂项排第一是因为它的应用面最广,很多看起来不像裂项的式子,经过一步构造之后都能用;无穷等比级数排得靠后,是因为它的失效后果最严重,需要更多的前置检查。
再分享一个小技巧:遇到一个陌生的求和式,先把它前三四项手算出来,看看数值有没有规律。如果前三项是 1, 3, 6, 10 这种三角形数,那通项大概率是 k(k+1)/2;如果是 1, 4, 9, 16,就是 k²。先看数字规律,再决定用什么公式,这条路比死盯通项更快,而且不容易选错方向。