同一份代码，两台机器，两个答案

double s = 0.0;
for (int i = 0; i < 10000000; i++) s += 0.1;
printf("%.17g\n", s);

本机实测

下面这些都是在写这本书的这台机器上真跑出来的（clang 14，x86-64 macOS）：

$ cc -O0 fma.c -o t && ./t
a*b+c   = -4.9303806576313238e-32
(x+y)+z = 1
x+(y+z) = 0
sum 0.1 x1e7 = 999999.99983897537

$ cc -O2 -ffast-math fma.c -o t && ./t
a*b+c   = -4.9303806576313238e-32
(x+y)+z = 1
x+(y+z) = 1                        ★ 变了
sum 0.1 x1e7 = 1000000.0000223046  ★ 变了

$ cc -O2 -ffp-contract=off fma.c -o t && ./t
a*b+c   = 0                        ★ 又变了

答案是 C，而且 D 也对——-ffast-math 那次给出的 1000000.0000223046 比 -O0 那次的 999999.99983897537 更接近正确答案。

这一点很值得强调：优化改变结果，方向是随机的。它有时更准，有时更差，而你无法预测是哪一种。

顺带一个交叉核对：-O0 那个 999999.99983897537，和第 4 章里 Node.js 跑出来的、以及 Python 跑出来的，逐位相同。三套完全不同的运行时，只要严格按 IEEE 754 顺序执行，就一定给出同一个位模式。这就是标准的价值，也正是 -ffast-math 放弃掉的东西。

为什么会变：四个来源

① 加法不满足结合律

volatile double x = 1e16, y = -1e16, z = 1.0;

(x + y) + z  =  0 + 1  =  1
x + (y + z)  =  1e16 + (-1e16 + 1)  =  1e16 + (-1e16)  =  0
                                        ↑ 那个 1 在 -1e16 面前被抹掉了

这是这一章的根。实数加法满足结合律，浮点加法不满足。而编译器为了向量化和指令级并行，必须重排加法顺序——重排就换答案。

-ffast-math 的作用就是「允许编译器假装浮点加法满足结合律」。上面那个例子里，加了这个开关之后 x+(y+z) 被编译器直接重排成了 (x+y)+z，输出从 0 变成 1。

② 并行归约的宽度

不用编译器，光是「用几路累加器」就够了。同一份 65536 个数的数组，同一个求和逻辑：

六种宽度，六个不同的答案。注意偏移是带符号的，而且分布在精确值两侧——1 路差 −13 格，16 路差 +13 格，同样远、方向相反。也就是说，你没法靠「哪个宽度更准」来选：它们各自的误差是不同顺序的舍入攒出来的，方向纯属偶然。而「宽度」这个东西，在真实系统里由 SIMD 位宽、线程数、GPU 的 block 大小、cuDNN 挑的 kernel 决定——这些都不是你写在代码里的。换一张显卡、换一个 batch size、换一个 cuDNN 版本，宽度就可能变。

③ FMA：一次舍入还是两次

fma(a, b, c) 算 a*b + c，中间的乘积不舍入，整个式子只舍一次。这比分开算更准，但结果不同：

a = 1 + 2⁻⁵²,  b = 1 − 2⁻⁵²,  c = −1

a × b 的精确值 = 1 − 2⁻¹⁰⁴
double 存不下 ⇒ 抹成 1
(a*b) + c = 1 − 1 = 0                      ← 分开算，两次舍入

fma(a, b, c) = −2⁻¹⁰⁴ = −4.9303806576313238e−32   ← 一次舍入，精确

★ 本机默认（-ffp-contract=on）：编译器把 a*b+c 自动变成 fma ⇒ 得到 −4.93e−32
★ 加 -ffp-contract=off                              ⇒ 得到 0

注意这里有多刁钻：同一行源码 a*b + c，编译器默认就会把它换成 FMA（C 标准允许这种「收缩」）。所以「有没有 FMA 硬件」「编译器版本」「优化级别」三者中的任何一个变化，都会让这一行给出不同的数。

④ 数学库函数不保证正确舍入

IEEE 754 要求加减乘除和开方必须正确舍入（第 8 章），但没有要求 sin、exp、pow、log。于是 glibc、macOS 的 libm、MSVC、Intel MKL、CUDA 的实现各自可以差最后一两个 ulp。

这就是为什么「同一个 Python 脚本，Linux 和 Mac 上结果差一点」——大部分时候差的不是 Python，是底下的 libm。

-ffast-math 到底放弃了什么

这个开关不是「稍微不那么严格」，它是一整包假设。GCC 和 clang 下它等价于同时打开：

第二行最危险，值得单独说：-ffinite-math-only 会让 isnan() 失效。因为编译器被告知「不可能有 NaN」，于是它可以把这个判断优化成常量假。你的 NaN 保护代码被编译器删掉了，而它什么也不会说。

更麻烦的是：这个开关会传染。如果一个共享库是用 -ffast-math 编译的，在某些平台上它会在加载时修改进程级的 FPU 控制字（打开 flush-to-zero），于是你的程序里所有的次正规数都被冲成 0——包括不是你编译的那些代码。

#include <stdio.h>
#include <math.h>
int main(void) {
  volatile double a = 1.0 + 0x1p-52, b = 1.0 - 0x1p-52, c = -1.0;
  printf("a*b+c   = %.17g\n", a * b + c);
  printf("fma     = %.17g\n", fma(a, b, c));
  volatile double x = 1e16, y = -1e16, z = 1.0;
  printf("(x+y)+z = %.17g\n", (x + y) + z);
  printf("x+(y+z) = %.17g\n", x + (y + z));
  double s = 0.0;
  for (int i = 0; i < 10000000; i++) s += 0.1;
  printf("sum 0.1 x1e7 = %.17g\n", s);
  return 0;
}

$ cc -O0 fp.c -o t -lm && ./t
a*b+c   = -4.9303806576313238e-32
fma     = -4.9303806576313238e-32
(x+y)+z = 1
x+(y+z) = 0
sum 0.1 x1e7 = 999999.99983897537

$ cc -O2 -ffast-math fp.c -o t -lm && ./t
x+(y+z) = 1                          ★ 结合律被「修好」了
sum 0.1 x1e7 = 1000000.0000223046    ★ 循环被向量化，累加顺序变了

$ cc -O2 -ffp-contract=off fp.c -o t -lm && ./t
a*b+c   = 0                          ★ 不许收缩成 FMA
fma     = -4.9303806576313238e-32    （显式调 fma 不受影响）

这一章的一句话

浮点加法不满足结合律，而编译器和并行硬件都必须重排加法——所以「同一份代码给出同一个答案」不是默认，是一个你要付性能去买的属性。

最后一章不再引入新东西。它把前十九章压成一张自查表（合上书之后能带走的那部分）、一张「什么时候该用什么」的决策表，以及一份继续往下走的路线图。还会回答一个从第 1 章就悬着的问题：这套东西真正的用途，其实不在浮点数上。