Implementing FMA and finding bugs in C and Rust standard libraries

This is a story of how I tried to compute a * b + c, and found that Rust and musl libc get it subtly wrong. Fused multiply-add (FMA) computes a * b + c with only one rounding error instead of two. It’s an important building block for e.g. trigonometric functions like sin(x) and cos(x) if you want to implement them accurately.

这是一个关于我如何尝试计算 a * b + c，却发现 Rust 和 musl libc 在实现上存在细微错误的故事。融合乘加（FMA）在计算 a * b + c 时仅产生一次舍入误差，而非两次。如果你想精确实现 sin(x) 和 cos(x) 等三角函数，它是非常重要的基础组件。

It’s a basic primitive that is usually implemented in hardware, but there is still some hardware out there that doesn’t have it. You’d think it would be something like cheap phones, but no, it’s Intel. Cheap ARM phones have it and it’s been required in 64-bit ARM since the very beginning, but Intel has been launching new parts without AVX2 or fused multiply-add as recently as 2021. (I’ve come to learn that whenever something is holding SIMD back, it’s usually Intel.)

这是一个通常由硬件实现的基本原语，但目前仍有一些硬件不支持它。你可能会认为这只会出现在廉价手机上，但事实并非如此，它是 Intel。廉价的 ARM 手机都具备该功能，且 64 位 ARM 从一开始就强制要求支持，然而直到 2021 年，Intel 仍在发布不带 AVX2 或融合乘加指令的新处理器。（我逐渐意识到，每当有什么东西阻碍了 SIMD 的发展，通常都是 Intel。）

15% of machines in the Firefox hardware survey don’t have AVX2 and hardware FMA that comes with it, so it has to be emulated for precise algorithms built on top of it to work correctly. On machines without hardware FMA, Rust’s std::simd gives up and runs scalar FMA on each f32 in [f32; 4] individually, which is slow. I wanted to do better in fearless_simd and provide an actually vectorized implementation.

在 Firefox 的硬件调查中，15% 的机器没有 AVX2 以及随之而来的硬件 FMA，因此为了让基于此构建的精确算法正常工作，必须对其进行模拟。在没有硬件 FMA 的机器上，Rust 的 std::simd 会放弃向量化，转而对 [f32; 4] 中的每个 f32 分别运行标量 FMA，这非常缓慢。我想在 fearless_simd 中做得更好，并提供一个真正的向量化实现。

Emulating FMA with SIMD

使用 SIMD 模拟 FMA

Since FMA works on three f32 values, the total number of possible inputs is 2 to the 96th power. It would take the world’s largest supercomputer only 500 years to try them all. We’ve come a long way! But I need working FMA later this year, so exhaustive verification isn’t really on the cards. The best I could do is some known values plus some random tests.

由于 FMA 处理三个 f32 值，可能的输入总数是 2 的 96 次方。即使是世界上最强大的超级计算机，也需要 500 年才能遍历所有情况。我们确实进步了！但我今年晚些时候就需要可用的 FMA，所以穷举验证是不可能的。我能做的最好方案就是测试一些已知值加上一些随机测试。

I followed the 2008 paper “Emulation of FMA and correctly-rounded sums: proved algorithms using rounding to odd” by Sylvie Boldo and Guillaume Melquiond, which has a formal proof of correctness in Coq. That way I don’t have to worry about trying to verify the algorithm myself. For f32, their algorithm is refreshingly simple: compute a * b + c in f64, then round it to f32. The only caveat is special handling of rounding errors in the conversion, in case the result falls exactly between two representable values.

我参考了 Sylvie Boldo 和 Guillaume Melquiond 在 2008 年发表的论文《Emulation of FMA and correctly-rounded sums: proved algorithms using rounding to odd》，该论文在 Coq 中提供了正确性的形式化证明。这样我就不必担心自己去验证算法了。对于 f32，他们的算法简单得令人耳目一新：在 f64 中计算 a * b + c，然后将其舍入到 f32。唯一的注意事项是在转换过程中对舍入误差进行特殊处理，以防结果恰好落在两个可表示值之间。

Translating the algorithm to SIMD was also straightforward: just do all that basic math per-lane. The special handling of rounding is needed very rarely - you hit it less than once in a million when processing values in the [-1, 1) range, so just put it under an if and it’ll be fine. Done!

将该算法转换为 SIMD 也非常直接：只需在每个通道（lane）上执行所有基本数学运算即可。舍入的特殊处理非常罕见——在处理 [-1, 1) 范围内的值时，触发概率不到百万分之一，所以只需加一个 if 判断就足够了。搞定！

Benchmarks look great: it’s 5x faster than std::simd, even bigger than the expected 4x speedup because Rust standard library has to check if FMA is available on the system in every call, and also has to worry about setting floating-point exception flags, both of which add overhead.

基准测试结果很棒：它比 std::simd 快 5 倍，甚至超过了预期的 4 倍加速，因为 Rust 标准库在每次调用时都必须检查系统是否支持 FMA，并且还要处理浮点异常标志的设置，这两者都会增加开销。

Follow the White Rabbit

追随白兔

Just a few hours later, a wild @awxkee appeared and posted some inputs on which my implementation diverged from the hardware results. Moments like these are why I love open source. I have no idea where he came from or how he even found this PR, he’s never contributed code to Fearless SIMD before. But there it was, a counter-example that broke my translation of a formally verified algorithm.

仅仅几个小时后，一位名为 @awxkee 的网友出现了，并发布了一些导致我的实现与硬件结果不一致的输入数据。这样的时刻就是我热爱开源的原因。我不知道他从哪里来，也不知道他是如何找到这个 PR 的，他以前从未为 Fearless SIMD 贡献过代码。但事实摆在眼前，这是一个打破了我对形式化验证算法翻译的“反例”。

Turns out I translated the paper into code incorrectly. I forgot to add special handling for values very close to zero, called subnormal values, which have a slightly different representation, and many of the usual floating-point “tricks” don’t work on them. The check that would apply the rounding fixup handled them incorrectly. (Aside: Subnormals are fascinating! Learn all about them here.)

原来是我将论文翻译成代码时出了错。我忘记了为非常接近零的值（称为非规格化数/subnormal values）添加特殊处理，它们的表示方式略有不同，许多常见的浮点“技巧”在它们身上并不适用。那个用于应用舍入修正的检查逻辑错误地处理了它们。（题外话：非规格化数非常迷人！点击此处了解详情。）

So I looked at the paper more carefully and fixed my implementation to match it more closely. Then I added tests covering the problematic inputs, random tests that try a million subnormals, and random tests that try a million values that should require fixup, for good measure.

于是我更仔细地研读了论文，并修正了我的实现以使其更贴合原意。然后，我添加了覆盖这些问题输入的测试，进行了百万次非规格化数的随机测试，以及百万次需要修正的数值的随机测试，以确保万无一失。

Wait… Why are tests still failing, but only on old systems? 等等……为什么测试仍然失败，而且只在旧系统上失败？

Down the rabbit hole

深入兔子洞

Fearless SIMD has CI configured to run tests in an emulator for every supported SIMD level, on top of running them on the host. This is the only way to test AVX-512 codepaths on CI. It would also fail if we tried to use instructions not available in hardware, although the Rust compiler already verifies this for us. It’s just good practice to test the codepaths on the hardware where they’d actually run.

Fearless SIMD 的 CI 配置为在模拟器中运行每个支持的 SIMD 级别的测试，同时也在宿主机上运行。这是在 CI 上测试 AVX-512 代码路径的唯一方法。如果我们尝试使用硬件不支持的指令，它也会失败，尽管 Rust 编译器已经为我们验证了这一点。在实际运行的硬件上测试代码路径是一种良好的实践。

The built-in f32::mul_add in the Rust standard library has the exact same bug! It fails to properly handle rounding for subnormals. As does std::simd. But the software implementation is only invoked for systems that don’t have FMA in hardware, so without the emulator we wouldn’t have noticed.

Rust 标准库中内置的 f32::mul_add 存在完全相同的 Bug！它无法正确处理非规格化数的舍入。std::simd 也是如此。但软件实现仅在没有硬件 FMA 的系统上被调用，所以如果没有模拟器，我们根本不会注意到。

So I report the bug to the Rust standard library, then transcribe the formally verified algorithm from the paper again (hopefully correctly) to insulate the scalar fallback inside fearless_simd from it. Tests pass.

于是我向 Rust 标准库报告了这个 Bug，然后再次（希望这次是正确的）从论文中转录了那个经过形式化验证的算法，以将 fearless_simd 内部的标量回退逻辑与标准库隔离开来。测试通过了。

Here’s the algorithm, it’s not that scary: 这就是那个算法，它并没有那么可怕：

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
}

Just the Rust standard library left to fix. 现在只剩下 Rust 标准库需要修复了。

How deep does this go?

这到底有多深？

Fixing the standard library is trickier because it not only computes the result but also sets the floating-point status flags. It also has this codepath not just for f32 using f64 but also for f64 using f128 on hardware where that’s available. So I write tests, for f32 and f64 both, then fix the code to the best of my understanding.

修复标准库更棘手，因为它不仅要计算结果，还要设置浮点状态标志。此外，它不仅有针对使用 f64 的 f32 的代码路径，还有在支持的硬件上使用 f128 的 f64 代码路径。所以我编写了针对 f32 和 f64 的测试，然后根据我的理解修复了代码。

Wait, what’s this at the top of the file? 等等，文件顶部这是什么？

/* origin: musl src/math/fmaf.c Ported to generic Rust algorithm in 2025, TG. */

I wonder if… 我在想……

if ((u.i & 0x1fffffff) != 0x10000000) /* not a halfway case */

Yep, that’s the same bug: this check ignores subnormals. So it’s… 没错，这就是同一个 Bug：这个检查忽略了非规格化数。所以它是……