浮点数与稳定计算

0.1 + 0.2 为什么不等于 0.3？同一个数学公式，为什么换个写法就能得到更准确的结果？要回答这些问题，需要同时看数怎样存储、每一步怎样舍入，以及问题本身对输入有多敏感。数值分析概览介绍了误差的基本分类；这里从实际计算展开。

符号、指数和有效数

Goldberg 的浮点运算论文（Oracle 授权重印）把浮点数写成符号、有效数（significand）和基数的指数幂。以 IEEE 754 binary64 的正规数为例：

$$x=(-1)^s(1+f)2^e,\qquad 0\le f<1.$$

64 位中，1 位存符号，11 位存指数，52 位存小数部分。正规数的二进制首位固定为 1，因此共有 53 位有效精度。若指数域存的是整数 $E$，则 $e=E-1023$；$E=0$ 和 $E=2047$ 留给零、次正规数、无穷与 NaN，不能套用这个正规数公式。

这不是小数点后固定保留多少位。在区间 $[2^e,2^{e+1})$ 内，正的正规数之间的间距是

$$\Delta=2^{e-52}.$$

数的量级每翻一倍，间距也翻一倍。例如，1 附近向上的间距是 $2^{-52}$，2 附近是 $2^{-51}$；到了 $10^{16}$ 附近，间距已经是 2，加上 1 就可能被舍掉。较小的浮点格式改变精度与范围；模型量化讨论的低位权重近似还涉及尺度和编码，不能只按总位数推断计算精度。

十进制显示与实际存值

Python 浮点数教程解释了 0.1 的二进制展开为何无限循环，以及 Python 为什么通常只显示能还原同一存值的最短十进制表示。显示成 0.1 不等于存下了精确的 $1/10$；增加显示位数也不会改变存值。

下面的例子使用 CPython 3.9 及以上版本的标准库（math.ulp 从 3.9 起提供），假定 float 是 binary64、运算采用最接近舍入，中点取偶数。第一个断言检查二进制基数和有效精度；以下各代码块都可独立运行。

import math
import sys
from fractions import Fraction

assert sys.float_info.radix == 2 and sys.float_info.mant_dig == 53
x = 0.1
print(repr(x))
print(format(x, ".17g"))
print(x.as_integer_ratio())
print(0.1 + 0.2, (0.1 + 0.2) == 0.3)
error = Fraction.from_float(x) - Fraction(1, 10)
print(error)
print(float(error), float(error / Fraction(1, 10)))
print(error / Fraction.from_float(math.ulp(x)))
print(math.ulp(1.0), math.ulp(2.0), math.ulp(1e16))
print(1e16 + 1.0 == 1e16)

0.1
0.10000000000000001
(3602879701896397, 36028797018963968)
0.30000000000000004 False
1/180143985094819840
5.551115123125783e-18 5.551115123125783e-17
2/5
2.220446049250313e-16 4.440892098500626e-16 2.0
True

as_integer_ratio() 给出实际存值的精确分数。把这个分数减去 $1/10$，得到 $1/180143985094819840$，约 $5.55\times10^{-18}$；每次算术运算还可能再引入舍入，因而 0.1 + 0.2 又比存下的 0.3 大了一步。

需要保留十进制输入语义时，可以从字符串构造 Decimal；需要精确有理数时，可以用 Fraction。先转成 float，再改用这两种表示，只能保留已经舍入的值，不能找回原输入。

绝对误差、相对误差与 ulp

设真实值为 $y$，计算值为 $\hat y$：

$$E_{\mathrm{abs}}=|\hat y-y|,\qquad E_{\mathrm{rel}}=\frac{|\hat y-y|}{|y|}\quad(y\ne0).$$

绝对误差与结果单位相同，相对误差则衡量损失的比例。真实值为零时，相对误差没有定义；接近零时，即使绝对误差很小，相对误差也可能很大。

ulp 是有效数末位所对应的数值单位，用来按浮点间距衡量误差。上面的 0.1 存值误差恰好是 math.ulp(0.1) 的 $2/5$。math.ulp 的定义对应正有限数向上相邻的间距（最大有限数另作处理）；在 2 的整数次幂处，向下间距与向上间距可能不同，使用 ulp 时要说清楚参照。

在正规数范围内，正确的最接近舍入带来至多半个局部 ulp 的误差。对非零、没有上溢或下溢的正规运算结果，常用模型是

$$\operatorname{fl}(a\circ b)=(a\circ b)(1+\delta), \qquad |\delta|\le u=2^{-53},$$

其中 $\circ$ 是基本算术运算。这里的单位舍入误差 $u$ 是 sys.float_info.epsilon（1 向上相邻的间距 $2^{-52}$）的一半。这个单步界不保证整个程序的相对误差也在 $u$ 内。

舍入、上溢、下溢与特殊值

最接近舍入选择距离精确结果最近的存值；恰好在中点时，取有效数末位为偶数的那一个。1e16 + 1 就处在中点，前面的运行结果显示它舍回了 1e16。

情况	binary64 中的含义	对计算的影响
上溢	结果超出有限范围；最大有限值约为 $1.798\times10^{308}$	在最接近舍入下，足够大的结果变成带符号的无穷；语言接口也可能抛异常
次正规数	满足 $0<\lvert x\rvert<2^{-1022}$（幅值小于最小正正规数），用隐含首位 0 表示	渐进下溢保留小值，间距固定为 $2^{-1074}$，相对精度逐渐降低
舍入到零	小到无法保留为非零次正规数	最小正次正规数的一半在取偶数规则下舍入为零
无穷	inf 与 -inf 是特殊值	可表示范围溢出或某些无界结果；inf - inf 产生 NaN
NaN	“非数”，可来自无效运算	与自身也不相等；用 math.isnan 检查，用 math.isfinite 检查结果是否有限

次正规数本身不一定有误差。零也有正负两种编码。Goldberg 的历史论文详细解释了渐进下溢与特殊值；实际语言接口的行为需另看定义。Python math 的异常约定说明，某些无效运算和上溢会抛异常，而不是直接返回 IEEE 特殊值：

import math
import sys

print(sys.float_info.min)
tiny = math.ulp(0.0)
print(tiny, tiny / 2)
print(1e308 * 1e308)
nan = float("inf") - float("inf")
print(nan, nan == nan, math.isnan(nan))
for operation in (
    lambda: 1.0 / 0.0,
    lambda: math.sqrt(-1.0),
    lambda: math.exp(1000.0),
):
    try:
        operation()
    except (ZeroDivisionError, ValueError, OverflowError) as exc:
        print(type(exc).__name__)

2.2250738585072014e-308
5e-324 0.0
inf
nan False True
ZeroDivisionError
ValueError
OverflowError

5e-324 是最小正次正规数的短十进制显示，不是精确的十进制值。把中间量缩放到合适范围，通常比等到它变成零或无穷后再补救有效；仅检查最终结果有限，也不能证明它准确。

消去误差：换一个等价公式

两个接近的数相减，可能把已有误差放大到结果的主要部分。减法本身甚至可以精确；损失可能已经发生在两个被减数的计算中。Goldberg 在消去误差一节区分了这种情况与精确操作数相减。

考虑 $x$ 接近零时的 $g(x)=\sqrt{1+x}-1$。对实数 $x\ge-1$，乘以共轭项可得

$$g(x)=\frac{(\sqrt{1+x}-1)(\sqrt{1+x}+1)}{\sqrt{1+x}+1} =\frac{x}{\sqrt{1+x}+1}.$$

取 $x=10^{-16}$，直接式先把 1 + x 舍入成 1，再减 1 就得到零。改写式把小量 $x$ 留在分子里，分母接近 2，即使分母舍入为 2，仍可保留约 $5\times10^{-17}$ 的结果。

import math
from decimal import Decimal, localcontext

x = 1e-16
naive = math.sqrt(1.0 + x) - 1.0
stable = x / (math.sqrt(1.0 + x) + 1.0)
with localcontext() as ctx:
    ctx.prec = 80
    dx = Decimal.from_float(x)
    reference = dx / ((Decimal(1) + dx).sqrt() + Decimal(1))
    print(f"reference = {reference:.18e}")
    for name, value in [("naive", naive), ("stable", stable)]:
        absolute = abs(Decimal.from_float(value) - reference)
        if reference != 0:
            relative = absolute / abs(reference)
            print(f"{name} = {value:.18e}; relative error = {relative:.3e}")
        else:
            print(f"{name} = {value:.18e}; absolute error = {absolute:.3e}")

reference = 4.999999999999999770e-17
naive = 0.000000000000000000e+00; relative error = 1.000e+0
stable = 4.999999999999999895e-17; relative error = 2.500e-17

参考值用 80 位十进制精度计算，是高精度近似。Decimal.from_float(x) 保留同一个二进制输入，因此这里比较的是算法误差，没有混入十进制输入与二进制输入之间的差别。直接式相对误差为 1，改写式约为 $2.5\times10^{-17}$。独立地看，展开式 $g(x)=x/2-x^2/8+O(x^3)$ 也给出约 $5\times10^{-17}$，第二项约为 $-1.25\times10^{-33}$。

把 x 改成 0 时，真实结果也为零，代码改报绝对误差。改写式避免了零附近的消去，但不能避免结果本身的下溢：最小正次正规输入的结果约为其一半，在 binary64 中仍舍入为零。

同样的思路适用于专门设计的库函数：log1p(x)在零附近准确计算 $\log(1+x)$；expm1(x)避免直接计算 $e^x-1$ 的消去。使用时仍须满足各函数的实数定义域。这些改写也会影响自动微分所经过的中间量；自动微分本身不会消除浮点误差。

求和顺序与补偿求和

实数加法满足结合律，浮点加法通常不满足。取 $10^{16}$、十个 1 和 $-10^{16}$，精确和为 10。逐项从左往右加，每个 1 都可能在大部分和中被舍掉；先抵消两个大数则能保住它们。

Goldberg 的求和误差分析解释了 Kahan 补偿求和：除了部分和，还保留一次加法丢失的低位信息，供下一次加法修正。下面显式写出逐项求和，避免依赖不同 Python 版本中内置 sum 的实现。

import math

def sequential(values):
    total = 0.0
    for value in values:
        total += value
    return total

def kahan(values):
    total = correction = 0.0
    for value in values:
        adjusted = value - correction
        new_total = total + adjusted
        correction = (new_total - total) - adjusted
        total = new_total
    return total

values = [1e16] + [1.0] * 10 + [-1e16]
print(sequential(values))
print(sequential([1e16, -1e16] + [1.0] * 10))
print(kahan(values), math.fsum(values))
hard = [1e16, 1.0, -1e16]
print(kahan(hard), math.fsum(hard))

0.0
10.0
10.0 10.0
0.0 1.0

correction = (new_total - total) - adjusted 估计这次加法相对于实际加入量的偏差；下一步先减去这项补偿。在十个 1 的例子里，它能回收丢失的贡献。但最后一行说明，普通 Kahan 对 [1e16, 1.0, -1e16] 仍会丢掉 1，不能把补偿求和当成精确算术。

Python 的 math.fsum保留多个中间部分和，两个例子都得到精确和。它的精度依赖 IEEE 算术保证与通常的中点取偶数舍入；文档还指出，某些非 Windows 构建的扩展精度加法可能造成二次舍入，影响末位。

对 $n\ge1$ 个有限浮点输入，将其存值视为精确输入，令 $S=\sum_i x_i$。以下逐项求和的界只计加法舍入，不包含输入转换误差；假设没有上溢或下溢，非零结果满足前述单步模型，精确为零的加法取 $\delta=0$。每个输入沿计算链最多经历 $n-1$ 次舍入，把这些乘法因子展开，可得到绝对误差界：

$$|\hat S-S|\le\gamma_{n-1}\sum_{i=1}^n|x_i|, \qquad \gamma_k=\frac{ku}{1-ku}\quad(ku<1).$$

当 $ku$ 很小时，$\gamma_k$ 约为 $ku$。若 $S\ne0$，相对误差界还要乘上 $\sum_i|x_i|/|S|$；正负项严重抵消时，这个因子会很大。

对于大批量数据，平衡的两两求和把加法组织成树，将最长链条从 $n-1$ 缩短为 $\lceil\log_2 n\rceil$；在相同舍入条件下，界中的 $\gamma_{n-1}$ 可以换成对应树深的 $\gamma_k$。先加小量也常能减少它被大部分和吞掉的机会。不过，混合正负数时，没有一种简单排序能保证最佳结果。改变并行归约顺序也可能改变末位。选择顺序、补偿或更高精度时，要结合数的量级和所需误差；NumPy 笔记给出数组归约与比较的使用方法。

条件数、后向稳定性与容差

条件数衡量数学问题对输入扰动的敏感程度。对 $y=a-b\ne0$，若每个输入的相对扰动至多为 $\varepsilon$，由三角不等式可得

$$|\Delta y|\le\varepsilon(|a|+|b|),\qquad \frac{|\Delta y|}{|y|}\le\kappa\varepsilon, \qquad \kappa=\frac{|a|+|b|}{|a-b|}.$$

把 $a=1.000001$、$b=1.000000$ 视为精确的实数输入，差为 $10^{-6}$，$\kappa=2000001$。输入相对不确定性为 $10^{-8}$ 时，输出绝对不确定性的界是 $2.000001\times10^{-8}$，相对界约为 2%。即便计算完全精确，也不能消除这种输入不确定性。

后向稳定性考察算法的结果能否看作附近输入问题的精确解：若计算 $f(x)$ 得到 $\hat y=f(x+\Delta x)$，且输入扰动按指定范数或逐分量尺度很小，后向误差就小。所谓“后向稳定”要求在适用的输入范围内得到与舍入精度相称的扰动界，不能只给单次结果找到任意一个附近输入。前向误差则直接比较 $\hat y$ 与 $f(x)$；当扰动很小、$f$ 可微且 $f(x)\ne0$ 时，一阶关系约为“相对前向误差 $\lesssim$ 条件数 × 相对后向误差”。病态问题即使用后向稳定算法，也可能有明显的前向误差。

前面的平方根例子则有不同来源：对 $x>0$，其相对条件数

$$\left|\frac{xg'(x)}{g(x)}\right| =\frac{\sqrt{1+x}+1}{2\sqrt{1+x}} \longrightarrow1\quad(x\to0^+).$$

问题本身在这里条件良好，直接式却丢掉了结果。这个区分能帮助判断该改善输入，还是改写算法。

容差应由误差预算决定：输入的不确定性、离散化或近似误差，以及算法的舍入误差。绝对容差要有结果的单位，相对容差要对应有意义的非零尺度；机器精度只描述表示，不能代替这个预算。

math.isclose对有限值使用

$$|a-b|\le\max\bigl(\text{rel\_tol}\max(|a|,|b|),\text{abs\_tol}\bigr).$$

import math

print(math.isclose(0.1 + 0.2, 0.3, rel_tol=0.0, abs_tol=math.ulp(0.3)))
print(math.isclose(1e-12, 0.0, rel_tol=1e-9))
print(math.isclose(1e-12, 0.0, rel_tol=1e-9, abs_tol=1e-11))

True
False
True

第一行允许一个以存下的 0.3 为参照的 ulp 差异，适合这个已分析的短运算。后两行演示零附近必须提供有意义的绝对容差；1e-11 只是演示阈值，不是通用默认值。NaN 与任何值都不接近，无穷只与同号无穷接近。数组比较应按 NumPy 的接口定义选择容差，不能直接把 math.isclose 的公式套到另一个 API。

实际检查时，先确定输入与输出尺度，再排查非有限中间量和消去，选择可靠的求和或专用函数，最后用精确解、高精度参考或可证明的界检验误差。两个算法彼此接近，只说明它们一致，还需要独立依据才能判断准确。

我的阅读 · 打开过的笔记会记在当前浏览器里。