一尺之槌，日取其半，1075日而竭

今天的故事源自数月前我在 r/learnpython 上回答的一个问题。以下是对该问题的转述：

将 64 位浮点数 1.0 不断除以 2，多少次后恰好变为 0？

x, i = 1.0, 0
while x:
    x /= 2.0
    i += 1
print(i)  # 1075

fn main() {
    let mut x = 1.0f64;
    let mut i = 0;
    while x != 0. {
        x /= 2.0;
        i += 1;
    }
    println!("{}", i);  // 1075
}

let x = 1,
  i = 0
while (x) {
  x /= 2
  i++
}
console.log(i) // 1075

古语云：“一尺之槌，日取其半，万世不竭”。但很显然，此话不适用于这个问题。浮点数的精度是有限的，当 x / 2 足够小时，浮点数终无法表示，便会得到结果 0。如以上代码所示，循环将终结在第 1075 次。

1075 次。这表明 64 位浮点能表示的最小数在 $2^{-1075}$ 这个量级，再小的数超出了表示范围，也就归于零了。但 1075 这个数很奇怪，它离 1024 很近，但又多出来一些，有零有整的，这是为什么呢？

这一切还要从浮点数的分类说起。

浮点数的分类

精通 Javascript 的读者会记得，JS 的数字体系中有一些奇怪的值，如 ±Infinity、NaN、±0。这些值行为怪异，有着和“正常值”与众不同的运算法则。

这不是因为 JS 标新立异，恰恰相反，这是因为 JS 遵循了 IEEE 754 标准。JS 采用 IEEE 754 标准中的 64 位浮点数（以下简称 FP64）表示数字，而 ±Infinity、NaN 等是标准中定义好的一系列特殊值，用来表示浮点运算中的诸多异常结果，如除以零、上下溢出等等。

这不免让人好奇：FP64 中一共有哪几类值呢？除了无穷大、零和 NaN 外，其他的都是正常值吗？

Rust 提供了一个非常有用的函数 f64::classify()，顾名思义即为浮点数作分类。该函数返回的结果是 core::num::FpCategory，其定义如下：

pub enum FpCategory {
    Nan,
    Infinite,
    Zero,
    Subnormal,
    Normal,
}

可以看到浮点数被分成五类，其中有三类 Infinite、Zero 和 Nan 为我们熟知，分别代表 无穷大、零和 NaN 三种情况。但刨去这些特殊值，还有两个类别 Subnormal 和 Normal。它们分别代表什么呢？

要了解二者的区别，我们要先知道 FP64 的二进制表示是如何构成的。

FP64 的二进制表示

64 位浮点数，顾名思义，即采用 64 个比特位表示一个实数 x。这 64 位由三部分组成：符号 s、指数 e 和底数 m，分别占据 1、11 和 52 位长度。下图详细展示了 -0.9375 的二进制表示，用不同颜色标出了这三个部分：

我们常用科学计数法 $\pm a \times p^n$ 表示大数。FP64 也是类似的，只不过是在二进制意义下。FP64 使用如下公式，将 s, e, m 组合成实数 x：

$$ \begin{align} x = (-1)^s \times 2^{e-1023} \times \left(1+m / 2 ^{52}\right) \label{eq:normal} \end{align} $$

不难发现其与科学计数法的相似性。FP64 将底数的范围固定为 $[1,2)$，最后的 52 位 m 正是底数的小数部分；中间的 11 位 e 取值范围为 $[1, 2046]$，用于合成指数；而最高位的符号位 s 顾名思义，即指示 x 的正负性。在后文中，为了叙述方便，我们也会将 x 的 FP64 表示中的这三部分记为 $\text{fp}_s(x)$、$\text{fp}_e(x)$ 和 $\text{fp}_m(x)$。

细心的读者可能又会奇怪了，为什么 e 的范围是 $[1, 2046]$ 呢？e 有 11 位，理论上取值范围应该是 $[0, 2047]$，还有两个值 $e=0$ 和 $e=2047$ 去哪了？

事实上，这两个值被 IEEE 754 特意空了出来，用于表示一些特殊的值。比如当 $s = e = m = 0$ 时，x 的值为 +0。更具体些，根据 e 和 m 的不同取值，FP64 的值可以分为以下五类：

• $e \in [1, 2046]$ 时，$x$ 为 Normal；
• $e=2047, m = 0$ 时，$x$ 为 ±Infinity；
• $e=2047, m \not = 0$ 时，$x$ 为 NaN；
• $e = 0, m = 0$ 时，$x$ 为 ±0；
• $e = 0, m \not = 0$ 时，$x$ 为 Subnormal。

本文只关注最后一类 Subnormal 浮点数，以下简称 Subnormal FP。

Subnormal FP

Subnormal FP，顾名思义，即比 Normal FP 还要小的一类数。事实上，它们的计算也有别于 Normal FP：

$$ \begin{equation} x = (-1)^s \times 2^{-1022} \times (m / 2^{52}) \label{eq:subnormal} \end{equation} $$

下图给出了一个 Subnormal 浮点数 $x=-3\times 2^{-1074}$ 的例子

Subnormal FP 与 Normal FP 的区别主要体现在以下几点：

• 由于 $\text{fp}_e(x)$ 固定为了 $0$，表达式 $\eqref{eq:subnormal}$ 中不再出现字母 $e$，而是直接用 $2^{-1022}$；
• 底数的范围不再是 $[1,2)$，而是 $[0,1)$，因此表达式中 $\eqref{eq:subnormal}$ 的 $1+m/2^{52}$ 变成了 $m/2^{52}$。

容易看出，由于 $2^{-1022}$ 这一项的存在，Subnormal FP 所表示的数绝对值都很小。这实际上弥补了 Normal FP 和 0 之间的空隙，使得 FP64 能够表示更多的数。我们还能看到，这种设计是融洽且符合直觉的：

• Normal FP 的最小值与 Subnormal FP 的最大值紧密相连，分别是 $2^{-1022}$ 和 $2^{-1022} \times (1 - 2^{-52})$；
• 零可以看作是 Subnormal FP 的特例。事实上，把 $s = m = 0$ 带入 $\eqref{eq:subnormal}$，正好可以得到 $0$。

尽管都能表示实数，但由于 Normal FP 和 Subnormal FP 的计算方式不同，IEEE 754 还是正式地将它们区分开来，并置于两个独立的类别中。

解密 1075

有了 Normal FP 和 Subnormal FP 的前置知识，我们终于可以解释为什么开头的代码会恰好在 1075 次后结束。本节先尝试用一个思想实验解释这一现象，然后用 Rust 代码验证这一点。

在思想实验前，我们先给出一条显而易见的规则：如果 $x$ 和 $x/2$ 都是 Normal FP，则除法发生前后，s 和 m 都不变，只有 e 减一，即 $\text{fp}_e(x/2)=\text{fp}_e(x)-1$。读者可以结合 Normal FP 的定义 $\eqref{eq:normal}$ 自行验证这一点。

现在考察程序的起点，也就是 $x_0 = 1$ 的二进制表示：

$$\text{fp}_s(x_0) = 0, \text{fp}_e(x_0) = 1023, \text{fp}_m(x_0) = 0$$

根据前述规则，当 x 还在 Normal FP 的表示范围内时，每次除以二相当于将指数减一，直至指数达到 0。此过程可持续 1022 次：

$$ \begin{aligned} x_1 &= 2^{-1} &\leadsto &\quad \text{fp}_s(x_1) = 0, \text{fp}_e(x_1) = 1022, \text{fp}_m(x_1) = 0 \\ x_2 &= 2^{-2} &\leadsto &\quad \text{fp}_s(x_2) = 0, \text{fp}_e(x_2) = 1021, \text{fp}_m(x_2) = 0 \\ &&\ldots &\\ x_{1022} &= 2^{-1022} &\leadsto &\quad \text{fp}_s(x_{1022}) = 0, \text{fp}_e(x_{1022}) = 1, \text{fp}_m(x_{1022}) = 0 \\ \end{aligned} $$

此时若再进一步，我们将得到 $x_{1023} = 2^{-1023}$，但这个数已经离开了 Normal FP 的范围，进入到 Subnormal FP 的领域。

现在我们需要另一条规则：如果 $x$ 和 $x/2$ 都是 Subnormal FP，则做除法相当于将 m 右移 1 位，即 $\text{fp}_m(x/2)=\text{fp}_m(x)\ shr\ 1$。这一规则同样可以从 Subnormal FP 的定义 $\eqref{eq:subnormal}$ 中直接得到。运用这条规则，我们有：

$$ \begin{aligned} x_{1023} &= 2^{-1023} &\leadsto &\quad \text{fp}_s(x_{1023}) = 0, \text{fp}_e(x_{1023}) = 0, \text{fp}_m(x_{1023}) = 2^{51} \\ x_{1024} &= 2^{-1024} &\leadsto &\quad \text{fp}_s(x_{1024}) = 0, \text{fp}_e(x_{1024}) = 0, \text{fp}_m(x_{1024}) = 2^{50} \\ &&\ldots &\\ x_{1074} &= 2^{-1074} &\leadsto &\quad \text{fp}_s(x_{1074}) = 0, \text{fp}_e(x_{1074}) = 0, \text{fp}_m(x_{1074}) = 2^0 \\ \end{aligned} $$

再往后，即使是 Subnormal FP 也无法表示如此小的数，于是 $x_{1075}$ 终于归零。在此期间，我们经历了 1022 个 Normal FP，52 个 Subnormal FP，以及最后 1 个 0，共 1075 个值。我们可以用如下 Rust 代码验证这一点：

fn main() {
    let mut counter = std::collections::HashMap::new();
    let mut x = 1f64;
    while x != 0. {
        x /= 2.0f64;
        *counter.entry(format!("{:?}", x.classify())).or_insert(0) += 1;
    }
    dbg!(counter);
}
// Output:
// [src/main.rs:8:5] counter = {
//     "Zero": 1,
//     "Normal": 1022,
//     "Subnormal": 52,
// }

Subnormal FP 的支持性问题

最后，本文就 Subnormal FP 再作进一步讨论。

尽管 IEEE 754 定义了 Subnormal FP，这并不意味着它在所有场景下都被支持。在实际中，Subnormal FP 的支持可能面临如下限制：

• 硬件本身不支持。这可能为了简化硬件设计，或者为了提高性能而舍弃了 Subnormal FP 的支持；
• 该功能被编程性地关闭了。某些硬件即便支持 Subnormal FP，其运算性能可能远不如 Normal FP。因此在一些场景下，为了提高浮点运算性能，程序会选择关闭 Subnormal FP 特性。

例如，在 x86 架构中，可通过设置 CSR 寄存器的 FTZ 标志位来控制 Subnormal FP 的支持。

fn main() {
    unsafe {
        use core::arch::x86_64::{_mm_getcsr, _mm_setcsr, _MM_FLUSH_ZERO_ON};
        _mm_setcsr(_mm_getcsr() | _MM_FLUSH_ZERO_ON);
    }
    let mut x = 1.0f64;
    let mut i = 0;
    while x != 0. {
        x /= 2.0;
        i += 1;
    }
    println!("{}", i); // 1023
}

在上面的代码中，我们通过 _mm_setcsr 函数设置了 FTZ 标志位。该标志位的出现指示 CPU：当运算结果得到 Subnormal FP 时，直接将其截断为 0。因此，上述代码会在 1023 次后结束，而不是 1075 次。

以上的例子表明，Subnormal FP 不是一个可靠的特性。在实际应用中，可能因为硬件的限制或是某些库的设置，浮点运算无法得到 Subnormal FP 值。通常来说，我们不提倡依赖 Subnormal FP 进行计算以获得更高的精度。如果确实需要 Subnormal FP，我们应当在程序中显式地开启或检查这一特性，以确保计算的正确性。

作者：hsfzxjy
链接：
许可：CC BY-NC-ND 4.0.
著作权归作者所有。本文不允许被用作商业用途，非商业转载请注明出处。