实现FMA并发现C和Rust标准库中的错误

# 实现FMA并发现C与Rust标准库中的漏洞 | Sergey "Shnatsel" Davidoff 来源：https://shnatsel.github.io/implementing-fma-finding-bugs-in-std/ 这是一个关于我尝试计算 `a * b + c`，却发现Rust和musl libc存在微妙错误的故事。 **融合乘加**（FMA）在计算 `a * b + c` 时仅产生一次舍入误差，而非两次。例如在实现`sin(x)`和`cos(x)`等三角函数时，若想获得精确结果，这是一个重要的基础构建块。 这种基本运算通常由硬件实现，但仍存在一些不具备该功能的硬件设备。你可能认为会是廉价手机之类的产品，但实际上是Intel。廉价的ARM手机已具备此功能，且自64位ARM架构诞生之初就要求支持FMA，但Intel直到2021年仍在推出不支持AVX2或融合乘加的新产品 (https://en.wikipedia.org/wiki/Tremont_(microarchitecture))。 （我逐渐了解到，每当某项功能阻碍了SIMD的发展，问题通常出在Intel身上。） 在Firefox硬件调查 (https://firefoxgraphics.github.io/telemetry/#view=system)中，15%的机器不具备AVX2及随之而来的硬件FMA，因此为了确保建立在其上的精确算法能正确运行，必须进行软件模拟。 在不具备硬件FMA的机器上，Rust的`std::simd`会放弃使用向量化指令，转而为`[f32; 4]`中的每个`f32`单独执行标量FMA，速度很慢。我想在fearless_simd (https://github.com/linebender/fearless_simd) 中做得更好，提供一个真正向量化的实现。 ## 使用SIMD模拟FMA https://shnatsel.github.io/implementing-fma-founding-bugs-in-std/#emulating-fma-with-simd 由于FMA处理三个`f32`值，总输入可能性为2的96次方。即使用世界上最大的超级计算机，也要耗时500年才能遍历所有情况。技术发展确实漫长！但今年晚些时候我就需要用到可用的FMA，因此 exhaustive verification 并不现实。我所能做的只是使用一些已知值加上随机测试。 我遵循了2008年论文 (https://guillaume.melquiond.fr/doc/08-tc.pdf)《使用舍入到奇数的FMA与正确舍入和模拟：经Coq证明的算法》（作者：Sylvie Boldo 和 Guillaume Melquiond），该论文在Coq中有形式化的正确性证明。这样我就不必亲自验证算法了。 对于`f32`，他们的算法简单明了：在`f64`中计算`a * b + c`，然后将结果舍入为`f32`。唯一的注意事项是转换过程中对舍入误差的特殊处理，以防结果恰好落在两个可表示数值的中间。 将该算法翻译为SIMD也很直接：只需在每个通道（lane）执行所有基础运算。舍入的特殊处理极少触发——在处理-1到1范围内的值时，出现概率低于百万分之一，因此只需放在`if`条件中处理即可。[已完成！](https://github.com/linebender/fearless_simd/pull/323) 基准测试结果出色：速度是`std::simd`的5倍，甚至超过了预期的4倍加速。这是因为Rust标准库在每次调用时都必须检查系统是否支持FMA，同时还需考虑设置浮点异常标志，这两者都增加了开销。 ## 追随白兔 https://shnatsel.github.io/implementing-fma-founding-bugs-in-std/#follow-the-white-rabbit 仅几小时后，一个名为 **@awxkee** (https://github.com/awxkee) 的用户突然出现，并发布了一些输入值 (https://github.com/linebender/fearless_simd/pull/323#issuecomment-5233682925)，在这些输入上我的实现与硬件结果产生了分歧。 正是这样的时刻让我热爱开源。我不知道他从何而来，也不清楚他是如何找到这个PR的，他从未向Fearless SIMD贡献过代码。但就在那里，一个反例打破了我对一个经过形式验证算法的实现。 原来是我将论文转化为代码时出错了。我忘记为非常接近零的值添加特殊处理，这些值称为**非规格化值（subnormal）**，它们的表示方式略有不同，许多常见的浮点“技巧”对它们不适用。处理舍入修正的检查对此类值的处理有误。 > **题外话：** 非规格化值很有趣！可以在此 (https://numerical-rust-cpu-81b2c3.pages.in2p3.fr/19-subnormal-entertainment.html) 了解更多。 因此我仔细重读了论文，并修正了实现，使其更贴近原文。接着，我添加了覆盖问题输入的测试、尝试一百万个非规格化值的随机测试，以及尝试一百万个应触发修正的值的随机测试，以防万一。 等等... 为什么测试仍然失败，但*仅*在旧系统上失败？ ## 兔子洞之旅 https://shnatsel.github.io/implementing-fma-founding-bugs-in-std/#down-the-rabbit-hole Fearless SIMD配置了CI，除了在主机上运行测试外，还为每个支持的SIMD级别在模拟器中运行测试。这是在CI中测试AVX-512代码路径的唯一方法。如果我们尝试使用硬件不支持的指令，它也会失败，尽管Rust编译器已经为我们验证了这一点 (https://shnatsel.github.io/safe-simd-in-rust-even-on-the-inside/)。在实际运行的硬件上测试代码路径仍然是良好的实践。 Rust标准库中的内置`f32::mul_add`**存在完全相同的错误！** 它未能正确处理非规格化值的舍入。`std::simd`也是如此。但软件实现仅在硬件没有FMA的系统上被调用，因此如果没有模拟器，我们不会注意到这个问题。 于是我向Rust标准库报告了该错误 (https://github.com/rust-lang/compiler-builtins/issues/1262)，然后再次（希望这次是正确的）从论文 (https://guillaume.melquiond.fr/doc/08-tc.pdf) 抄录经过形式验证的算法，以隔离`fearless_simd`内部的标量回退实现，使其不受影响。测试通过了。 以下是该算法，并不复杂： ``` fn scalar_mul_add_precise_f32(a: f32, b: f32, c: f32) -> f32 { let product = (a as f64) * (b as f64); let c = c as f64; let mut sum = product + c; if sum.is_finite() { let virtual_sum = sum - product; let rounding_error = (product - (sum - virtual_sum)) + (c - virtual_sum); let sum_bits = sum.to_bits(); if rounding_error != 0.0 && sum_bits & 1 == 0 { let corrected_bits = if sum.is_sign_negative() == rounding_error.is_sign_negative() { sum_bits.wrapping_add(1) } else { sum_bits.wrapping_sub(1) }; sum = f64::from_bits(corrected_bits); } } sum as f32 } ``` 只剩下标准库需要修复了。 ## 这问题有多深？ https://shnatsel.github.io/implementing-fma-founding-bugs-in-std/#how-deep-does-this-go 修复标准库更棘手，因为它不仅要计算结果，还要设置浮点状态标志 (https://en.wikipedia.org/wiki/FLAGS_register)。而且这段代码路径不仅适用于使用`f64`的`f32`，也适用于在支持的硬件上使用`f128`的`f64`。 因此我编写了`f32`和`f64`两者的测试，然后根据我的理解修复了代码 (https://github.com/rust-lang/compiler-builtins/pull/1270)。 等等，文件顶部这是什么？ ``` /* origin: musl src/math/fmaf.c Ported to generic Rust algorithm in 2025, TG. */ ``` 我在想... ``` if ((u.i & 0x1fffffff) != 0x10000000) /* not a halfway case */ ``` 没错，是同样的错误：这个检查忽略了非规格化值。 所以不仅Rust std有错误，**musl libc也有错误**。而且他们注明该实现源自FreeBSD (https://git.musl-libc.org/cgit/musl/tree/src/math/fmaf.c?id=f21a96538f78fa8e2040831b4209b35f2fb581da)，版权年份为2005-2011，因此谁知道这段有错误的代码在过去二十年里还被复制粘贴到了哪些地方。 Musl没有错误跟踪器，所以我们通过邮件列表报告 (https://www.openwall.com/lists/musl/2026/08/10/1)，并祈祷这不是一个有意的“优化”（邮件列表没有搜索功能，所以我无法检查），并且祈祷一旦它从最近消息页面消失后不会被永远遗忘。天哪，我真是爱死了这种30年前的开发流程！ ## 底线在哪里？ https://shnatsel.github.io/implementing-fma-founding-bugs-in-std/#where-is-the-bottom 哦嘿，我的标准库PR的CI检查完成了！ 等等，为什么它们失败了，但仅在32位ARM上？而且输出的完全是乱码？！ 不知何故，实验性的、仅限nightly的`f128`类型在32位ARM平台上被启用了，并且……完全失控了？我肯定不会去读那么多ARM汇编代码，但LLM可以轻松搞定 (https://github.com/rust-lang/compiler-builtins/pull/1270#issuecomment-5246545522)，然后……这是一个rustc的ABI bug。在这个类型上，在这个架构上，完全局限在我的PR中。 > **题外话：** 如果你好奇的话，使用LLM进行分析是Rust LLM策略 (https://blog.rust-lang.org/inside-rust/2026/08/05/rust-langrust-is-adopting-an-llm-policy/) 在特定条件下允许的。 因此现在也有关于此的错误报告 (https://github.com/rust-lang/compiler-builtins/issues/1271)，由一位能够验证LLM发现、并且比我更了解该领域的维护者提交。 ## 修复 https://shnatsel.github.io/implementing-fma-founding-bugs-in-std/#fixes `fearless_simd`很简单。已修复。你可以使用它来获得没有此错误的FMA。 Rust编译器中nightly专属的`f128`错误也已修复。 我为Rust标准库提交的PR (https://github.com/rust-lang/compiler-builtins/pull/1270) 在某些情况下破坏了浮点异常标志的设置。我想到用LLM搜索我可能引入的错误，它找到了一个反例。因此我投入了更多精力制作一个尽可能少改动代码的补丁。目前仍在等待审查。 令我惊讶的是，musl开发者认真对待了我的错误报告。基于其他开发者的反馈，一个最小化修复和一个对`fmaf()`函数的实质性重写方案很快被提出并迭代。在所有反馈被解决之后……什么也没发生。在我撰写本文时，它仍未被合并。 看来，无论你的开发工具有多年轻，代码审查都是瓶颈。 ## 这个错误重要吗？ https://shnatsel.github.io/implementing-fma-founding-bugs-in-std/#does-this-bug-matter 我不知道。 一方面，这只是不正确的舍入。结果在极少数情况下会偏离一个最低有效位。最初的反例与正确结果相差约百万分之十五。 另一方面，根据你如何使用结果，误差可能会被急剧放大。而那些具有精确误差界限的数学例程是用于*某些用途*的，因此无论用途是什么，得到的误差界限都是错误的。 或许最明显的是，在不同机器上运行的确定性模拟将会产生分歧。 ## FMA有错？在我的电脑里？ https://shnatsel.github.io/implementing-fma-founding-bugs-in-std/#buggy-fma-in-my-computer 在x86上，这仅在没有AVX2的廉价和/或非常老旧的Intel芯片，以及完全没有FMA的中国海光x86芯片上是个问题。64位ARM不受影响。 然而，这可能会在未来几年困扰32位ARM和各种嵌入式系统。如果那段有错误的FreeBSD实现被复制粘贴到许多不同平台的各种晦涩工具链中，我一点也不会感到惊讶。 如果FMA的准确性对你很重要，请检查`fmaf(a, b, c)`在以下输入比特模式下 ``` a = 0x97000800 b = 0x1cfff001 c = 0x00010002 ``` 计算结果是否为比特模式`0x00010001`（正确）而不是`0x00010002`（有错误）。