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

2024-10-10 | Tech

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

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

Python

Rust

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 𝑎 \times 𝑝 𝑛$ 表示大数。FP64 也是类似的，只不过是在二进制意义下。FP64 使用如下公式，将 s, e, m 组合成实数 x：

$𝑥 = (- 1) 𝑠 \times 2 𝑒 - 1023 \times (1 + 𝑚 / 252) (1)$

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

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

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

$𝑒 \in [1, 2046]$ 时， $𝑥$ 为 Normal；
$𝑒 = 2047, 𝑚 = 0$ 时， $𝑥$ 为 ±Infinity；
$𝑒 = 2047, 𝑚 \neq 0$ 时， $𝑥$ 为 NaN；
$𝑒 = 0, 𝑚 = 0$ 时， $𝑥$ 为 ±0；
$𝑒 = 0, 𝑚 \neq 0$ 时， $𝑥$ 为 Subnormal。

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

#Subnormal FP

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

$𝑥 = (- 1) 𝑠 \times 2 - 1022 \times (𝑚 / 252) (2)$

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

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

由于 $f p 𝑒 (𝑥)$ 固定为了 $0$ ，表达式 $(2)$ 中不再出现字母 $𝑒$ ，而是直接用 $2 - 1022$ ；
底数的范围不再是 $[1, 2)$ ，而是 $[0, 1)$ ，因此表达式中 $(2)$ 的 $1 + 𝑚 / 252$ 变成了 $𝑚 / 252$ 。

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

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

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

#解密 1075

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

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

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

$f p 𝑠 (𝑥 0) = 0, f p 𝑒 (𝑥 0) = 1023, f p 𝑚 (𝑥 0) = 0$

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

$𝑥 1 = 2 - 1 ⇝ f p 𝑠 (𝑥 1) = 0, f p 𝑒 (𝑥 1) = 1022, f p 𝑚 (𝑥 1) = 0 𝑥 2 = 2 - 2 ⇝ f p 𝑠 (𝑥 2) = 0, f p 𝑒 (𝑥 2) = 1021, f p 𝑚 (𝑥 2) = 0 \dots 𝑥 1022 = 2 - 1022 ⇝ f p 𝑠 (𝑥 1022) = 0, f p 𝑒 (𝑥 1022) = 1, f p 𝑚 (𝑥 1022) = 0$

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

现在我们需要另一条规则：如果 $𝑥$ 和 $𝑥 / 2$ 都是 Subnormal FP，则做除法相当于将 m 右移 1 位，即 $f p 𝑚 (𝑥 / 2) = f p 𝑚 (𝑥) 𝑠 ℎ 𝑟 1$ 。这一规则同样可以从 Subnormal FP 的定义 $(2)$ 中直接得到。运用这条规则，我们有：

$𝑥 1023 = 2 - 1023 ⇝ f p 𝑠 (𝑥 1023) = 0, f p 𝑒 (𝑥 1023) = 0, f p 𝑚 (𝑥 1023) = 251 𝑥 1024 = 2 - 1024 ⇝ f p 𝑠 (𝑥 1024) = 0, f p 𝑒 (𝑥 1024) = 0, f p 𝑚 (𝑥 1024) = 250 \dots 𝑥 1074 = 2 - 1074 ⇝ f p 𝑠 (𝑥 1074) = 0, f p 𝑒 (𝑥 1074) = 0, f p 𝑚 (𝑥 1074) = 20$

再往后，即使是 Subnormal FP 也无法表示如此小的数，于是 $𝑥 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
链接：http://i.hsfzxjy.site/repeatly-divide-one-until-zero/
许可：CC BY-NC-ND 4.0.
著作权归作者所有。本文不允许被用作商业用途，非商业转载请注明出处。

Rust Floating-Point Subnormal FP IEEE 754