Add mul_add_precise() which guarantees single-rounding precision

[Shnatsel]

This is useful for building algorithms with guaranteed precision. For example, most of the literature on getting accurate SIMD trigonometry requires precise fused multiply-add.

Lowers into:

• accurate FMA on WASM even with relaxed SIMD
• Custom scalar FMA based on a formally verified algorithm from a 2008 paper to work around a bug in Rust libm and musl libc
• SSE4.2 f32 emulates it on f64 SIMD vectors using the same algorithm from the same paper

The paper is "Emulation of FMA and correctly-rounded sums: proved algorithms using rounding to odd" by Sylvie Boldo and Guillaume Melquiond.

The SSE4.2 version is 5x faster than the libm scalar implementation, and 2.5x faster when the rare fixup branch always needs to be taken. On branches: Nehalem misprediction recovery is 17 cycles which amounts to a 40% throughput penalty, so this should never be worse than scalar even in the worst case for branch misprediction. Tremont timings are hard to find but it gets a very beefy predictor and its throughput is awful, so it shouldn't be any worse. Zen4 benchmarks confirm that worst case misprediction is roughly in line with scalar. A branchles formulation with more work in the hot path would be possible but with misprediction being rare, it isn't worth it in practice.

Also fixes mul_sub() on NEON which used to produce zero with an incorrect sign when all inputs are 1s with different signs. The fixed assembly has identical performance to the previous one on llvm-mca in both throughput and latency. It's always two instructions because NEON doesn't have mul_sub with the semantics we need, this PR only changes which two instructions are selected.

It is unfortunately not possible to remove the dedicated mul_sub operation and just write mul_add(a, b, -c) because LLVM does not optimize it reliably enough, see #18

[Shnatsel]

Pedantry for posterity:

Commit 33d9f98 says Tremont (2021 Atom) improves in latency from 20 to 12, but that's a mistake: it improves in reciprocal throughput from 20 to 12, and in latency from 21 to 16.

But this is already a very rarely taken branch so it doesn't really matter either way.

[awxkee]

It seems emulation currently fails a bit https://gist.github.com/awxkee/afd3ee135602056a56806db9ddfcb9c9. For pair (0x19ffe002, 0x1a001001) it fails almost half of results

[Shnatsel]

Oh! Good catch, thanks!

I've been wondering why the exhaustive test missed it, but in retrospect it wasn't exhaustive at all, I've been feeding the same value into all three lanes 🤦

I'll rework this and also add far more in-depth random testing.

[Shnatsel]

@awxkee I think you've just found a bug in libm 👀

I'll investigate in more detail and report, but basically its software FMA uses the same ideas as this one and it seems to fail on your inputs in forced-software mode (hardware FMA disabled).

[awxkee]

I strongly recommend checking glibc's or LLVM's libm, since the Rust implementation is ported from musl, and musl isn't the most robust one (and LLVM is the unreadable one). The typical FMA code seems to work great.

Code

#[inline]
pub fn soft_fma(x: f32, y: f32, z: f32) -> f32 {
    let xy = f64::from(x) * f64::from(y);
    let z = f64::from(z);
    let result = xy + z;
    let mut u = result.to_bits();
    if u & 0x0fff_ffff != 0 {
        return result as f32;
    }

    if u & 0x1000_0000 == 0 && (u >> 52) & 0x7ff > 1023 - 126 {
        return result as f32;
    }

    if result - xy == z && result - z == xy {
        return result as f32;
    }

    let neg = u >> 63 != 0;
    let err = if neg == (z > xy) {
        xy - result + z
    } else {
        z - result + xy
    };
    if neg == (err < 0.0) {
        u += 1;
    } else {
        u -= 1;
    }
    f64::from_bits(u) as f32
}

[Shnatsel]

The musl implementation was also submitted to glibc: https://sourceware.org/pipermail/libc-alpha/2025-November/172449.html

So I don't even know what to use as the source of truth anymore.

WASM CI also fails so maybe our wasmtime is built with musl?

[awxkee]

Rust have an unique 'rounding branch', it enters it on those values and fails. Therefore it seems it's a porting problem :)

WASM always use Rust standard libm as far as I understand. WASM itself doesn't have libm at all.

[Shnatsel]

The wasmtime failure in CI is even more bizarre:

wasm trap: wasm unreachable instruction executed

[Shnatsel]

wasm trap: wasm unreachable instruction executed

Okay it's the same libm bug, it's just that the wasmtime runner doesn't report the assertion failure properly.

[Shnatsel]

The same subnormals rounding bug is present in musl libc: https://git.musl-libc.org/cgit/musl/tree/src/math/fmaf.c?id=f21a96538f78fa8e2040831b4209b35f2fb581da

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

this only handles normal floats but not subnormals.

This paper has a fix - it provides the correct algorithm with a correctness proof in Coq: https://guillaume.melquiond.fr/doc/08-tc.pdf

@awxkee wanna report it upstream since it's your finding?

[awxkee]

Nah, feel free to report

[Shnatsel]

I've reported the bug to musl: https://www.openwall.com/lists/musl/2026/08/10/1

And opened a PR for Rust's libm.

And worked around it using the formally proven FMA formulation from the paper, using your implementation as a test oracle.

[LaurenzV]

Not gonna pretend I understand much of the implementation, but if you think this is important to have in the library, LGTM!

[LaurenzV]

Not super urgent, but I think this (and the FloatExt trait) are duplicated in multiple files. Maybe there should be a util.rs file that contains those?

[Shnatsel]

In this case not really because there are 2 instances of the plain function and 1 instance of it wrapped in kernel! which we need to keep separate. The kernel!-wrapped one is there for runtime detection of SSE2 on i586 targets which get x87 f32 arithmetic otherwise which is not IEEE 754 compliant and breaks the function's precision contract.

But beyond this helper - yes, there's probably something that could be moved into shared code.

[awxkee]

Nitpicking: Knuth's summation and double-double representation, triple-double etc don't track a residual of the operation, they track an error term or "rounding error". s = a + b has an error, and Knuth's 2Sum captures it exactly and unconditionally. While it can be considered as residual, an "error" is much more common in the literature.

[Shnatsel]

That's clearer. Thanks!