binary64 · ulp

浮点数：实数在代码里的残影

前四页造出的 $\mathbb{R}$ 是一个不可数集，每个元素要么是无穷不循环的小数展开，要么是一个 Cauchy 列的等价类。这两种东西都装不进有限的存储。代码里那个叫 number 的类型是另一种对象，本页量清它到底是什么，以及由此产生的几条必须知道的行为。

1 · 可表示数的有限性

IEEE 754 的双精度格式（binary64）用 $64$ 个二进制位存一个数： $1$ 位符号、 $11$ 位阶码、 $52$ 位尾数。位数有限，取值就有限。

定义 1.1（可表示数） binary64 的每个有限取值都可写成 $\pm m \cdot 2^{e}$ ，其中 $m$ 是不超过 $2^{53}$ 的非负整数、 $e$ 是整数。这样的数都是有理数，且分母是 $2$ 的幂。

由此立刻得到两条结论。第一，任何无理数都不是可表示数： $\sqrt 2$ 、 $\pi$ 、 $e$ 在代码里存的一律是某个邻近的二进制有理数。第二，可表示数的个数是可以数清的——实测 FINITE_DOUBLE_COUNT 为 $18437736874454810623$ ，约 $1.84 \times 10^{19}$ 个，含 $\pm$ 无穷与 NaN 之外的全部取值。

这些数在数轴上的分布极不均匀。实跑 countDoublesIn：区间 $[1, 2)$ 内有 $2^{52} = 4503599627370496$ 个可表示数，而 $[0, 1)$ 内有 $4607182418800017408$ 个，是前者的 $1023$ 倍，占全部有限可表示数的 $24.988\%$ 。四分之一的双精度值挤在 $0$ 与 $1$ 之间，另有四分之一挤在 $-1$ 与 $0$ 之间；剩下一半铺满整条数轴的其余部分。

图 1-1 ·

\log_2

相邻间距随

\log_2

取值的阶梯，以及选中值与它左右两个邻居的精确十进制写法。可点选不同量级，读数栏给出该处的间距与本量级内的可表示数个数。

2 · 十分之一在二进制下的循环节

$0.1$ 这个十进制小数不是可表示数。理由与有理数与小数展开 §3 完全同源，只是换了底： $1/10$ 的分母含因子 $5$ ，而 $2$ 的幂里没有 $5$ ，所以 $1/10$ 在二进制下是循环小数。实跑 binaryExpansion(1, 10) 得到 $0.0\overline{0011}$ ，前置段一位、循环节四位——那个 $4$ 就是 $2$ 在模 $5$ 下的乘法阶。

存进 double 时循环节被截断，得到的是最靠近 $1/10$ 的那个可表示数。它的精确值写得出来，因为分母是 $2$ 的幂：

\texttt{0.1} = 7205759403792794 \times 2^{-56} = 0.1000000000000000055511151231257827021181583404541015625

$55$ 位小数，比 $1/10$ 略大。同理 0.2 存进去的值也略大，两者相加再舍入，落到的可表示数不是存 0.3 时落到的那一个：

例 2.1 0.1 + 0.2 的精确值是 $0.3000000000000000444089209850062616169452667236328125$ ，而 0.3 的精确值是 $0.299999999999999988897769753748434595763683319091796875$ 。两个数的位模式相差 $1$ ，也就是相邻的两个可表示数，所以 0.1 + 0.2 !== 0.3。打印时看不出来是因为 Number.prototype.toString 输出的是能唯一还原该值的最短十进制串。

图 2-1 · 几个算式的打印值与精确十进制值的对照，以及

1/q

的二进制展开（横线标出循环节）。可切换算式与分母，读数栏给出差值折合多少个间距。

同一个原因让加法失去结合律。实测 (0.1 + 0.2) + 0.3 得 $0.6000000000000001$ ，而 0.1 + (0.2 + 0.3) 得 $0.6$ ：每一步各自舍入，先后次序不同则舍入的位置不同。把 0.1 累加十次得到的是 $0.9999999999999999$ ，累加一百次得到 $9.99999999999998$ ，与 $10$ 相差 $1.95 \times 10^{-14}$ 。误差随累加次数增长，这是数值求和要先排序或走 Kahan 补偿的原因。

3 · 间距与判等容差

相邻两个可表示数之间的距离叫 ulp（unit in the last place）。它不是常数：尾数位数固定，绝对间距随量级翻倍。

几个量级处的相邻间距，由 nextUp 减自身实测得出。
取值	相邻间距	说明
$0.1$	$1.3878 \times 10^{-17}$	比 $1$ 处的间距还小一个量级
$1$	$2.220446049250313 \times 10^{-16}$	等于 `Number.EPSILON`，即 $2^{-52}$
$10^9$	$1.1920928955078125 \times 10^{-7}$	绝对容差 $10^{-9}$ 在此处已经不起作用
$2^{53}$	$2$	此后整数不再全部可表示
$10^{300}$	$1.4870 \times 10^{284}$	间距本身已是天文数字

Number.EPSILON 是 $1$ 处的间距，不是「浮点的误差上限」。这个误读会直接写出错的判等：绝对容差 $10^{-9}$ 在 $1$ 附近够用，到 $10^9$ 处已经小于该处的间距 $1.19 \times 10^{-7}$ ，判等退化成 ===。反方向的坑同样存在：相对容差在 $0$ 附近失效，因为两个数的量级都趋于零时相对差可以任意大。

建议 · 判等取绝对与相对两条容差的较宽者：|a − b| <= absTol || |a − b| <= relTol * max(|a|, |b|)。绝对容差管 $0$ 附近，相对容差管大量级。closeTo 就是这么实现的，单测里 $10^{16}$ 与 $10^{16} + 2$ 这一对（相邻的两个可表示数）在绝对容差 $10^{-9}$ 下判为不等、在相对容差 $10^{-12}$ 下判为相等，正好把两条各自的适用范围划出来。

$2^{53}$ 这条界值得单独记。超过它之后间距是 $2$ ，2 ** 53 + 1 === 2 ** 53 成立，奇数一律存不下。这是 JavaScript 里 Number.MAX_SAFE_INTEGER 取 $2^{53} - 1$ 的来历，也是数据库的 64 位整数 ID 传到前端要用字符串的来历。

4 · 精确算术的代价

要精确，就得放弃「一个数一个格子」这件事。用 BigInt 的分子与分母表示有理数，加减乘除全部精确：本系列的 core/completeness.ts 就是这么做的，二分 $64$ 步之后端点仍是精确分数，分母 $20$ 位十进制数，一位不差。

代价有两处。一是分母会长：那 $64$ 步里分母按 $2$ 的幂增长，位数与步数成正比，运算时间随之上升。二是这条路只覆盖有理数——开方与超越函数出不来。 $\sqrt 2$ 本身写不成分数（见无理数：反证、递降与超越），所以精确有理算术对 $\sqrt 2$ 能给出的只是一列越来越窄的区间，给不出那个数。

于是「实数」在代码里只有两种落法：要么放弃精确，取一个邻近的二进制有理数，接受每步舍入；要么放弃「得到一个数」，只保留一个可以任意收窄的区间与一套区间上的运算。前者是 double，后者是符号计算与区间算术。两者都不是 $\mathbb{R}$ ，而 $\mathbb{R}$ 的完备性正来自那些装不下的对象。

5 · 参考文献

Double-precision floating-point format. Wikipedia. binary64 的位布局、取值范围与特殊值。https://en.wikipedia.org/wiki/Double-precision_floating-point_format
Unit in the last place. Wikipedia. ulp 的定义与几种常见约定。https://en.wikipedia.org/wiki/Unit_in_the_last_place
Machine epsilon. Wikipedia. 机器 epsilon 的两种定义，以及它同 ulp 的关系。https://en.wikipedia.org/wiki/Machine_epsilon
Goldberg, D. (1991). What every computer scientist should know about floating-point arithmetic. ACM Computing Surveys, 23(1), 5–48.
Number.MAX_SAFE_INTEGER. MDN. $2^{53} - 1$ 这条界的来历与它的实际后果。https://developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Global_Objects/Number/MAX_SAFE_INTEGER