添加浮点小数：趣味与盈利

# 为了乐趣和利润添加浮点小数数 来源：https://blog.vero.site/post/float 2026-08-30（2527字）分类：数学 (https://blog.vero.site/category/math)，计算机科学 (https://blog.vero.site/category/cs) 许多人知道不应该使用大多数编程语言中的浮点数进行小数计算，例如涉及美元和美分的计算。这是因为小数无法用这种浮点数精确表示，因此会遇到舍入误差。 著名的是，使用IEEE双精度浮点数 (https://en.wikipedia.org/wiki/Double-precision_floating-point_format)，`0.1 + 0.2`不等于`0.3`，而是`0.30000000000000004`。 另一方面，我经常使用Python REPL来加总小数金额以制作收据。当然，我的风险较低；我是在汇总1到100美元范围内的金额，并知道在将总和复制到其他地方之前手动舍弃微小的误差。但实际上，这些误差通常根本不会出现。 如果我们对所有不超过1.00的0.01倍数进行两两求和，并查看打印的结果是偏大（蓝色）、偏小（红色）还是正确，就会得到一个很酷的模式： 图1：浮点数如何影响两个不超过1.00的0.01倍数的求和 这个模式从何而来？ ### 浮点数简介 一个双精度浮点数 \(x\) 由1个符号位、11个指数位和52个小数位（按此顺序排列）组成，总共64位。 符号位提供符号，\(+\)或\(-\)。指数位表示一个介于-1022和1023（含）之间的整数 \(E\)。小数位表示一个小于 \(2^{52}\) 的非负整数 \(F\)，对应于有效数字 \(1 + F/2^{52} \in [1, 2)\)；这个始终存在的 \(1\) 称为“隐藏位”。浮点数的值是 \(x = \pm 2^E(1 + F/2^{52})\)。 因此，最后一个小数位的有效值是 \(2^{E-52}\)，这个量被称为 \(x\) 的**ulp**（“最后一位的单位”）。除了当 \(F = 0\) 且 \(x\) 恰好是2的幂时一侧的一个边界情况外，最接近 \(x\) 的两个浮点数与它相差恰好一个ulp。示例： ``` "隐藏位"52位 float(0.1) = 0.00011001100110011001100110011001100110011001100110011010（二进制） ulp(float(0.1)) = 0.00000000000000000000000000000000000000000000000000000001（二进制） ``` 这个描述对我们的目的来说已经足够，但忽略了其他许多情况：0、次正规数、无穷大和NaN；它们使用了我没有描述的两种指数位值。为简单起见，我也不考虑其他精度的浮点数。 ### 加法剖析 举个例子（参考例如 qntm (https://qntm.org/notpointthree)），让我们逐步分析当你在Python REPL中输入`0.1 + 0.2`时会发生什么。1 (https://blog.vero.site/post/float#fn1) 首先，Python表达式`0.1`求值为某个浮点数：具体来说，根据IEEE标准的要求，是最接近0.1的浮点数。我们将这个*精确*数字写作 \(\text{float}(0.1)\)。2 (https://blog.vero.site/post/float#fn2) 类似地，Python表达式`0.2`求值为另一个浮点数，\(\text{float}(0.2)\)。 现在，我们可以想象`\+`分两步求值。首先，Python计算*精确*的和 \(\text{float}(0.1) + \text{float}(0.2)\)。然后，将其舍入到最接近的浮点数。得到的值是 \(\text{float}(\text{float}(0.1) + \text{float}(0.2))\)。（这不是它字面上的工作方式——第一步中的精确和从未在任何地方具体化——但在数学上是准确的。） 实际上，这里有一个微妙之处，直到我详细推导时才注意到：\(\text{float}(0.1) + \text{float}(0.2)\) 与它的两个最近浮点数*等距*！发生这种情况时，Python舍入到具有偶数有效数字的浮点数，在本例中是向上舍入。3 (https://blog.vero.site/post/float#fn3) 最后，要打印这个值，Python必须将其转换为十进制数。这种转换非常微妙，并且IEEE标准没有精确规定！尽管 \(\text{float}(\text{float}(0.1) + \text{float}(0.2))\) 不完全等于0.3，但它非常接近，因此将其显示为“0.3”并非不合理。另一种选择是精确显示，如“0.3000000000000000444089209850062616169452667236328125”。也可以想象在各种舍入量后显示：0.300000000000000044，或0.30000000000000004441，等等。甚至可以考虑，比如0.30000000000000005，因为仍然成立 \(\text{float}(\text{float}(0.1) + \text{float}(0.2)) = \text{float}(0.30000000000000005)\)；也就是说，Python表达式`0.1 + 0.2 == 0.30000000000000005`为真。这种转换的微妙之处体现在，直到Python 3.1中的一个特定补丁 (https://bugs.python.org/issue1580) 之前，如果你在REPL中输入`1.1`，Python会将其作为`1.1000000000000001`回显给你。 关于此十进制转换应如何工作的标准描述由Steele和White在1990年 (https://dl.acm.org/doi/10.1145/93548.93559) 形式化4 (https://blog.vero.site/post/float#fn4)，他们提出了三个标准5 (https://blog.vero.site/post/float#fn5)： 1. 十进制表示应可往返转换：如果你再次输入它，应该得到相同的浮点数。这排除了输出“0.3”。 2. 在满足标准1的前提下，十进制表示应尽可能短。这排除了像“0.30000000000000004441”这样的输出。 3. 在满足标准1和2的前提下，十进制表示应尽可能接近浮点数。这排除了像“0.30000000000000005”这样的输出。 遵循这个算法，我们可以理解为什么`0.1 + 0.2 == 0.30000000000000004`。 ### 推广 让我们考虑一般情况：假设你要尝试将精确的正十进制量 \(a\) 和 \(b\) 相加，其精确和为 \(c\)。嗯，不是完全一般。我们假设这些值是“合理的美元金额”，为正数且小于7万亿美元；我认为这应该足以涵盖我需要提交的收据。略高于此（\(2^{46}= 70,368,744,177,664\)），浮点数变得比美分的倍数更稀疏，这不太好。 问题是，\(\text{float}(\text{float}(a) + \text{float}(b))\) 与 \(\text{float}(c)\) 相比如何？ 令它们的差为 \[\Delta := \text{float}(\text{float}(a) + \text{float}(b)) - \text{float}(c)。\] 我们可以通过引入误差函数 \(\text{error}(x) := \text{float}(x) - x\) 来推理它。然后，我们可以将 \(\Delta\) 重写为 \[\begin{aligned}\Delta ={} &\text{error}(a) + \text{error}(b) \\&+ \text{error}(\text{float}(a) + \text{float}(b)) - \text{error}(c)。\qquad(*)\end{aligned}\] 我们可以对每个误差项进行限定。回想一下，浮点数的*ulp*（“最后一位的单位”）是最后一位的值。为了轻微滥用符号，我们将允许自己写 \(\text{ulp}(x)\)，即使 \(x\) 无法精确表示为浮点数，并理解这意味着 \(\text{ulp}(\text{float}(x))\)。所以 \(\text{float}(x) \pm \text{ulp}(x)\) 也是浮点数6 (https://blog.vero.site/post/float#fn6)，它们离 \(x\) 的距离必须不小于 \(\text{float}(x)\) 本身（否则，\(\text{float}(x)\) 会求值为更接近的值）；这意味着对于所有（合理的）\(x\)，我们有 \[\|\text{error}(x)\| \leq \frac{\text{ulp}(\text{float}(x))}{2}。\] 此外，等号仅在 \(x\) 恰好位于两个最接近浮点数的中间时成立，如果 \(x\) 是合理的金额，这是不可能成立的。7 (https://blog.vero.site/post/float#fn7) 我们可以将这个界限逐项应用于 \((*)\)，得出 \(\|\Delta\| < 2\text{ulp}(c)\)。此外，因为 \(\Delta\) 是 \(c\) 附近两个浮点数的差，所以它是 \(\text{ulp}(c)\) 的倍数。8 (https://blog.vero.site/post/float#fn8) 由此我们得出结论 \(\Delta \in \{-\text{ulp}(c), 0, +\text{ulp}(c)\}\)——也就是说，结果最多偏离答案1个ulp。 然而，这里有一个推导可以产生 \(\|\Delta\|\) 更紧的中间界限：不失一般性地假设 \(a \leq b\)。那么，\(\text{float}(c) + \text{ulp}(c) - \text{float}(b)\) 是一个可表示的浮点数，因为结果的ulp ≤ \(b\) 和 \(c\) 的ulp。因此，至少这是一个可用的 \(a\) 的近似值。并且它是一个高估： \[\begin{aligned}a &= c - b \\ &\leq \text{float}(c) + \frac{\text{ulp}(c)}{2} - \text{float}(b) + \frac{\text{ulp}(c)}{2} \\ &= \text{float}(c) - \text{float}(b) + \text{ulp}(c)。\end{aligned}\] 因此，\[\begin{aligned}\text{float}(a) &\leq \text{float}(c) - \text{float}(b) + \text{ulp}(c) \\ \text{float}(a) + \text{float}(b) &\leq \text{float}(c) + \text{ulp}(c)。\end{aligned}\] 从这个不等式中减去 \(a + b = c\)，我们得到 \[\text{error}(a) + \text{error}(b) \leq \text{error}(c) + \text{ulp}(c)。\] 从另一侧也适用相同的界限。结果是，如果我们令 \[\delta := \text{error}(a) + \text{error}(b) - \text{error}(c),\] 我们有 \[-1 \leq \frac{\delta}{\text{ulp}(c)} \leq 1。\] 如前所述，我们知道 \(\|\delta - \Delta\| \leq \text{ulp}(c)/2\)，并且同样因为 \(\Delta\) 是 \(\text{ulp}(c)\) 的倍数，我们看到 \(\Delta \in \{-\text{ulp}(c), 0, +\text{ulp}(c)\}\)。 如果我们绘制 \(\delta/\text{ulp}(c)\) 的热力图，我们看到的可以描述为图1的一个更连续的版本： 图2：对两个不超过1.00的0.01倍数求和产生的组合浮点误差 我们现在可以理解图1是图2的“舍入”版本，棋盘模式出现在由于舍入到具有偶数有效数字的浮点数而导致 \(\delta\) 恰好为 \(\pm\text{ulp}(c)/2\) 的区域： 图3：具有奇数浮点有效数字的0.01倍数 并且，我们可以将图2解释为三个 \(\text{error}\) 函数副本“干涉”的结果：一个水平方向，一个垂直方向，一个对角线方向（尽管分母在变化）。 图4：将图2分解为各项 唯一剩下的问题是，为什么 \(\text{error}(x)\) 看起来像那样？ ### 一维误差 图5：不超过2.00的0.01倍数处的误差函数 首先我们观察到，每当 \(x\) 是2的精确幂时，\(\text{error}(x) = 0\)。在这样两个幂之间，让我们比较 \(\text{error}(x)\) 和 \(\text{error}(x+0.01)\)。我们有 \(\text{ulp}(x) = \text{ulp}(x+0.01)\)，所以 \(\text{float}(x + 0.01) \equiv 0 \equiv \text{float}(x) \bmod \text{ulp}(x)\)，所以 \[\text{error}(x + 0.01) \equiv \text{error}(x) - 0.01 \bmod \text{ulp}(x);\] 也就是说，\(\text{error}(x)\) 是一个模 \(\text{ulp}(x)\) 下“公差为-0.01的等差数列”。因此，该函数的环绕行为导致了我们之前图形中的周期性模式。 让我们关注图2的右下象限，\([0.5, 1] \times [0.5, 1]\)。在这个区域，我们可以计算出 \(\text{ulp}(0.5) = 2^{-53}\) 和 \(\text{ulp}(1) = 2^{-52}\)，然后得到 \[\begin{aligned}\frac{0.01 \bmod \text{ulp}(0.5)}{\text{ulp}(0.5)} &\equiv 0.92\equiv -0.08 \bmod 1 \\ \frac{0.01 \bmod \text{ulp}(1)}{\text{ulp}(1)} &\equiv 0.96\equiv -0.04 \bmod 1,\end{aligned}\] 两者都“接近0”。因为 \(0.08 \approx 1/12\)，当 \(x \in [0.5, 1]\) 时，\(\text{error}(x)\) 在环绕之前有12个“步长”；因为 \(0.04 = 1/25\)，当 \(x \in [1, 2]\) 时，\(\text{error}(x)\) 在环绕之前有25个“步长”。 为了更好地理解这个模式，我们可以计算出 \[\frac{0.01 \bmod 2^{-n}}{2^{-n}} = 0.01 \times 2^n \bmod 1 = \frac{2^n \bmod 100}{100}。\] 实际上，双精度浮点数中小数位的数量52是一个很好的巧合，使得 \(2^{52}\) 在模100下“接近0”；这就是为什么误差函数没有那么多次环绕，所以我们有平滑的区域。如果我们将图表扩展到 \(a, b \in [0, 2]\)，使得 \(c\) 可以达到 \([2, 4]\)，我们会看到更混乱的棋盘格和对角线，因为 \[\frac{0.01 \bmod \text{ulp}(2)}{\text{ulp}(2)} \equiv 0.48\bmod 1\] 并且 \(\text{error}(x)\) 在 \([2, 4]\) 中大约每隔一步环绕一次，这会与 \(c\) 的有效数字的奇偶性以更复杂的方式相互干涉。 图6：浮点数如何影响两个不超过2.00的0.01倍数的求和 图7：对两个不超过2.00的0.01倍数求和产生的组合浮点误差 图8：不超过4.00的具有奇数浮点有效数字的0.01倍数 图9：不超过4.00的0.01倍数处的误差函数 ### 附录：浮点数中的无误差变换 （这可能比正文更实用） 你实际上如何在计算机上计算像 \(\text{error}(x) = \text{float}(x) - x\) 这样的函数，例如，生成本文中的图表？显然，你不能在你正在研究的相同浮点格式中直接计算它。在那个格式中，\(\text{float}\) 是恒等函数；当你试图表达 \(x\) 时，误差已经产生了。 概念上最简单的方法是使用某种精确的有理数算术，比如Python的`fractions`。对于我的初步探索，我使用了我自己的 Noulith (https://github.com/betaveros/noulith)（在随意地为其`rational`类型添加了一堆特性和错误修复之后...）。 然而，事实证明，有一系列间接方法可以在不离开浮点格式的情况下处理这类误差。我相信这些技术被称为“无误差变换”。 **2Sum** (https://en.wikipedia.org/wiki/2Sum)（Møller，1965年）：从 \(a\) 和 \(b\)，计算 \(s\) 和 \(t\)，使得 \(a +_\text{float} b = s\) 且 \(a + b = s + t\) 精确成立。 ``` def two_sum(a: float, b: float) -> tuple[float, float]: s = a + b bb = s - a return s, (a - (s - bb)) + (b - bb) ``` **Veltkamp分割**9 (https://blog.vero.site/post/float#fn9)：从 \(a\)，计算 \(h\) 和 \(\ell\)，使得 \(a = h + \ell\) 精确成立，并且 \(h\) 和 \(\ell\) 最多有26个有效位（隐藏位之后）。这很有用，因为以浮点数相乘两个这样的浮点数是精确的。（你可以通过更改魔术常量来重新分配 \(h\) 和 \(\ell\) 之间的有效位数。） ``` def veltkamp_split(a: float) -> tuple[float, float]: c = ((1 << 27) | 1) * a hi = c - (c - a) return hi, a - hi ``` **Dekker乘积**10 (https://blog.vero.site/post/float#fn10)：从 \(a\) 和 \(b\)，计算 \(p\) 和 \(r\)，使得 \(a \times_\text{