Skip to content

simd: mul_add rounds once on every backend (wasm emulation, exact scalar FMA, no libm) - #374

Merged
AdaWorldAPI merged 11 commits into
masterfrom
fix/wasm-fused-mul-add-20261024
Oct 24, 2026
Merged

AdaWorldAPI merged 11 commits into
masterfrom
fix/wasm-fused-mul-add-20261024

Conversation

@AdaWorldAPI

@AdaWorldAPI AdaWorldAPI commented Oct 24, 2026 •

Copy link
Copy Markdown
Owner

Summary

mul_add on the float wrappers now rounds once on every backend. This makes it IEEE 754 fusedMultiplyAdd: one rounding of the exact a * b + c, bit-identical across all eight realizations. Backend selection stays compile-time only (simd.rs → simd_{arch}.rs). No dependency is added.

Before this PR:

  • wasm F32x16/F64x8 computed mul then add, which rounds twice. Under relaxed-simd they used relaxed_madd, which fuses on some engines and not others.
  • Targets without an FMA instruction (scalar backend, x86-64-v2, AVX without FMA3) called f32::mul_add / f64::mul_add. There these lower to a call into the toolchain's compiler_builtins fma/fmaf, which picks its implementation at run time by CPUID (fma_with_fma / fma_with_fma4 / fma_fallback). That is run-time dispatch inside the polyfill, which the polyfill exists to rule out. It also put VEX FMA instructions into the v2 and AVX binaries, which simd-matrix's instruction scans correctly flagged.
  • AVX2 F32x16/F64x8 ran mul_add as a scalar loop even when FMA3 was available.

The mul_add contract (reference for further SIMD work)

Result. a.mul_add(b, c) is the correctly rounded a * b + c, with round-to-nearest-even and one rounding. Subnormals are included, and NaN in gives NaN out. A backend that has the instruction must use it (FMA is mandatory where the CPU has it). A backend that lacks it must emulate it exactly, never as mul + add, and never through libm.

backend (selected in simd.rs) F32x16 / F32x8 F64x8 / F64x4 scalar mul_add_f32/f64
AVX-512 (simd_avx512.rs) _mm512_fmadd_ps _mm512_fmadd_pd mul_add → vfmadd
AVX2 + FMA3 (simd_avx2.rs, F32x8/F64x4 in simd_avx512.rs) _mm256_fmadd_ps on each half _mm256_fmadd_pd on each half mul_add → vfmadd
AVX without FMA3, x86-64-v2, scalar (simd_scalar.rs) exact emulation per lane (simd_fma::fma_f32) exact emulation per lane (simd_fma::fma_f64) the same
NEON (simd_neon.rs) vfmaq_f32 vfmaq_f64 mul_add → fmadd
wasm simd128 (simd_wasm.rs) f64 product + round-to-odd sum, vectorized Dekker + round-to-odd, vectorized; out-of-range lanes go to simd_fma::fma_f64 simd_fma
nightly core::simd (simd_nightly/) StdFloat::mul_add StdFloat::mul_add —

Where it lives.

  • src/simd_fma.rs (new, pub(crate)) holds the one ExactFma trait. With target_feature = "fma" on x86, and on aarch64, it is mul_add, i.e. the hardware instruction. Otherwise it is the exact integer emulation:
    • fma_f32: the product is exact in f64. TwoSum gives the error of p + c, the sum is rounded to odd (one ulp toward the error when inexact and even), and the single rounding to f32 is then correct (Boldo & Melquiond, Emulation of FMA and correctly rounded sums: proved algorithms using rounding to odd, IEEE Trans. Computers 57(4), 2008).
    • fma_f64: the significands multiply exactly in a u128, c is aligned into the same window with a sticky bit, and the sum is rounded to nearest-even once, including into the subnormal range and overflow.
  • crate::simd::mul_add_f32 / mul_add_f64 (in simd_ops.rs, re-exported from simd.rs) are the scalar entry points. Scalar tails such as add_mul_f32/f64 and the gemm_f64_tiled_fma tail now use them instead of hand-rolling f32::mul_add.
  • wasm (simd_wasm.rs) keeps the vector path from the first commit:
    • F32x16: promote, exact f64 product, add_round_to_odd, demote.
    • F64x8: Dekker split (2^27 + 1), TwoSum, round-to-odd sum of the error terms, one final rounding.
    • That F64x8 path is exact only while the split cannot overflow or underflow (2^-460 <= |a|, |b| <= 2^460, |c| <= 2^960, finite, non-zero). A v128 pair with any lane outside that range goes lane by lane through simd_fma::fma_f64, no longer through f64::mul_add.

Reverted. An intermediate commit whitelisted the libm fma symbols in the simd-matrix instruction scans (a03450e, ebefc2f). That hid the run-time dispatch instead of removing it, and is reverted (b21226b, c6b285d). The workflow is identical to master.

Verification

  • Parity group 0x11xx in crates/simd-masking-parity:

    • It compares each backend's mul_add bit-exactly against an independent reference, oracle_fma_f32 / oracle_fma_f64. The reference is a 4608-bit integer place-value sum and contains no mul_add.
    • Fixtures: random bit patterns over every exponent, near-cancelling c, exact f32 and f64 ties that a tiny c must break, exact f64 midpoints that must round to the even neighbour, and f64 lanes outside the wasm fast range.
    • Anti-vacuity: the unfused result must differ on at least 256 lanes per type.
    • oracle_tests checks the reference itself against hardware FMA.
  • All 8 arms on b41760d, 18 groups bit-identical: native v4, native v3, avx-qemu (Sandy Bridge), v2-qemu (Nehalem), neon-qemu, wasm (simd128), wasm-scalar, nightly core::simd. Each arm prints its own CPU flags, and those were read.

  • The unmodified simd-matrix scans on the parity binaries: v2 (%ymm|%zmm|%k[0-7]|VEX xmm) gives 0 hits. avx (AVX2 ymm integer ops, FMA3, F16C) gives 0 hits. objdump -t finds no fma_with_fma* symbol in either binary.

  • simd_fma unit tests (on v4 against hardware mul_add; also run on v2): 200,000 f32 and 200,000 f64 random cases; more than 1000 of them must differ from unfused, so the test has teeth. Special values are checked too: every triple drawn from {±0, ±1, 1.5, 3, MIN_POSITIVE, ±smallest subnormal, ±MAX, ±inf, NaN}, which covers signed-zero sums, overflow and subnormal results.

  • Disable runs, each red:

    removed fails at
    wasm f32 round-to-odd (plain f64 add) 0x112e
    wasm f64 round-to-odd in the error-term sum 0x1150
    the whole wasm f32 path (unfused mul + add) 0x1110
    simd_fma::fma_f32 round-to-odd, on the v2 arm 0x112e, plus the simd_fma unit test
    simd_fma::fma_f64 ties-to-even → ties-away, on the v2 arm 0x1160, plus both simd_fma unit tests

    The last row first came back green in the parity group (the unit test alone caught it). Every f64 tie round added a tiny c that broke the tie, so the final rounding never saw an exact midpoint. Round 4 (b41760d) now cancels the product's middle terms, so a*b + c is exactly halfway between two f64 neighbours. The even neighbour is below on some lanes and above on others, in both signs, over 2^±900. The D2 row was re-run after that change and still fails at 0x1150. The two f32 wasm rows were run before b41760d, which does not touch the f32 rounds. All disable runs were done before the rebase onto simd: saturating narrows down the 512-bit ladder (16x32 -> 32x16 -> 64x8) #373. The rebase only merged the narrow ladder next to this code; the mul_add paths are unchanged.

  • Clippy --all-targets -D warnings on v4, v3 and v2: clean. No-std (--no-default-features, also with approx,serde) under -D warnings: clean. simd_fma is gated on std like every module that uses it (a119298). The CI no-std jobs had failed on b27d401 because of this.

  • Codegen on b41760d (probe crate, objdump): on v4 mul_add_f32/f64 is one vfmadd213ss/sd, and F32x16/F64x8::mul_add is one vfmadd on zmm. On v3 the scalar entry points are one vfmadd213ss/sd, and the vectors are two vfmadd213ps/pd on ymm (the halves). On v2 there is no FMA instruction and no libm call: the f32 emulation is inlined, and F64x8 calls simd_fma::fma_f64 per lane. With constant arguments, v4/v3 and the v2 f32 path fold to the constant; the v2 f64 path does not (see below).

  • Scalar call sites: the 9 scalar f32/f64::mul_add calls in src/hpc (randomized_signature, signature_pde, renderer) and the two lane loops of BF16x16/F16x16::fma now call simd::mul_add_f32/f64 (efdc7d3). Their module tests pass on v2, where the emulation is the active path (68 tests).

Not covered

  • On targets without FMA, simd_fma::fma_f64 (u128 integer path) does not constant-fold. A call with constant arguments stays a runtime call there, and F64x8::mul_add makes 8 out-of-line calls per vector.
  • No timing was taken. The emulated paths are much slower than the instruction and exist only where there is no instruction.
  • The F16/BF16 fma in simd_half still rounds twice. That is a separate issue.
  • numeric/impl_numeric.rs (var/std) calls the generic num_traits::Float::mul_add. On a no-FMA build that still reaches compiler_builtins. It is generic upstream ndarray code and is not changed here.
  • The nightly arm's StdFloat::mul_add was only verified with FMA (v4). It was not checked on a build without FMA.
  • Rebased onto master after simd: saturating narrows down the 512-bit ladder (16x32 -> 32x16 -> 64x8) #373 (the narrow ladder) merged; the parity group list now carries both check_narrow_ladder and check_mul_add_single_rounding.

🤖 Generated with Claude Code

https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg

@coderabbitai

coderabbitai Bot commented Oct 24, 2026 •

Copy link
Copy Markdown
Contributor

Important

  • 🔍 Trigger review

This repository does not receive automatic reviews because it has fewer than 10 stars.

Run configuration
  • Source mode: Auto
  • Repository: Unavailable
  • Branch: Unavailable
  • File: No configuration file from this repository applied
  • Commit: Unavailable
  • Layers: Organization UI
  • Review profile: CHILL
  • Plan: Advanced
  • Run ID: a81e1747-d651-4734-808e-1fb3e3922077
  • Autofix · Keep fixing CodeRabbit findings and required CI, and resolving merge conflicts

Warning

Billing warning: we have not been able to collect payment for this subscription for more than 72 hours. Please update the payment method or pay any pending invoices in Billing to avoid service interruption.

Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out.

❤️ Share

Comment @coderabbitai help to get the list of available commands.

@AdaWorldAPI AdaWorldAPI changed the title simd(wasm): mul_add rounds once on F32x16 and F64x8 simd: mul_add rounds once on every backend (wasm emulation, exact scalar FMA, no libm) Oct 24, 2026
@AdaWorldAPI
AdaWorldAPI marked this pull request as ready for review October 24, 2026 17:41
@chatgpt-codex-connector

Copy link
Copy Markdown

Codex Review Summary

This comment shows the latest Codex review activity on this pull request.

Review Status Commit Review trigger
📝 Code Review 🔄 Running since 2026-10-24T17:41:53.906284Z fdaf898 Draft marked ready
ℹ️ About Codex in GitHub

Your team has set up Codex to review pull requests in this repo. Reviews are triggered when you

  • Open a pull request for review
  • Mark a draft as ready
  • Comment "@codex review" or "@codex security review".

Codex reacts with 👀 while any review is running, comments if it has suggestions, and reacts with 👍 once all reviews finish with no findings.

claude added 11 commits October 24, 2026 17:45
simd128 has no FMA, and relaxed-simd's relaxed_madd may or may not fuse
depending on the engine. wasm mul_add was mul-then-add (two roundings), or
relaxed_madd under relaxed-simd, so it differed from every other backend.

- F32x16: each lane in f64, where a * b is exact; the f64 sum is rounded to
  odd (two-sum error, then one ulp toward it when inexact and even), so the
  single rounding to f32 is correct, subnormals included (Boldo & Melquiond,
  IEEE TC 2008).
- F64x8: Dekker exact product, two-sum with c, rounded-to-odd sum of the error
  terms, one final rounding (the same paper's emulated FMA). A v128 pair with
  a lane outside 2^-460 <= |a|, |b| <= 2^460, |c| <= 2^960 (or NaN, inf, zero)
  is computed lane by lane with f64::mul_add.

The relaxed-simd branches are gone. Parity group 0x11xx compares both against
the element mul_add bit for bit (NaN matches NaN): random bit patterns over
every exponent, near-cancelling c, ties that double rounding breaks the wrong
way, and lanes that take the f64 fallback; anti-vacuity requires the unfused
result to differ on at least 256 lanes per type.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
The mul_add parity group calls scalar f32/f64::mul_add as its reference,
which links compiler-builtins' libm fma/fmaf. Those carry fma_with_fma /
fma_with_fma4 variants picked at runtime by CPUID (the Nehalem qemu run
takes fma_fallback and passes), and their vfmadd* lines tripped the
"no AVX instruction" scan. The scan now skips exactly those four
compiler_builtins functions; any other AVX op still fails it (checked
against libc.so.6: 9620 hits remain).

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
The Sandy Bridge scan hit the same four compiler_builtins fma/fmaf
CPUID-dispatched variants as the v2 scan. Locally: 4 hits before, 0 after
on the avx-qemu parity binary; the filtered scan still reports 4481
AVX2/FMA/F16C lines in libc.so.6.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
…truction

On targets without a compile-time FMA instruction, f32/f64::mul_add is a
libm call that picks an FMA3, FMA4 or software path by CPUID at run time.
The polyfill selects its realization at compile time, so the per-lane
mul_add of the scalar, AVX-without-FMA and x86-64-v2 backends, and the wasm
F64x8 lanes outside the vector range, now use simd_fma::ExactFma:

- with target_feature = "fma" (x86) or on aarch64: mul_add, which lowers to
  the instruction;
- otherwise f32 through f64 with the sum rounded to odd, and f64 in integer
  arithmetic (exact sum in a 128-bit window with a sticky bit, one rounding).

The 0x11xx parity group now compares against its own exact reference (a
4608-bit integer sum rounded to nearest-even) instead of the element
mul_add, so the parity binary links no libm fma at all. Both the emulation
and the reference are checked against hardware FMA in unit tests
(400k + 300k cases plus all special-value triples).

This replaces the CI-scan exclusions of libm's CPUID-dispatched fma, which
are reverted.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
…imd_ops

- simd_avx2 F32x16/F64x8::mul_add were a per-lane scalar loop even on
  AVX2 + FMA. They now call the halves' F32x8/F64x4::mul_add:
  _mm256_fmadd_ps/_pd with target_feature = "fma", the exact scalar FMA per
  lane without it. F64x4 gains mul_add; F32x8's non-FMA branch uses the
  exact scalar FMA instead of f32::mul_add.
- simd_ops gains mul_add_f32 / mul_add_f64 (re-exported from ndarray::simd):
  the polyfill's single-rounding scalar FMA, so tails and consumers (wasm
  included) do not hand-roll one or fall back to f32::mul_add. The add_mul
  and gemm tails use it; the gemm doc no longer claims wasm is unfused.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
The no-std builds (tests/1.99.0 without std, nostd/thumbv6m, blas-msrv)
compiled simd_fma while every caller (simd, simd_ops, simd_wasm,
simd_avx512) is std-gated, so -D warnings failed on dead code.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
Nine scalar f32/f64::mul_add calls in hpc (randomized_signature,
signature_pde, renderer) and the two lane loops of BF16x16/F16x16::fma
called the std method, which on a target without the fma feature is a
call into compiler_builtins' run-time-dispatched fma. They now go
through the polyfill's exact scalar FMA. Results are unchanged: both
paths are the correctly rounded a*b+c. The BF16/F16 double rounding
(f32 FMA, then narrowing) is untouched.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
The f64 tie rounds always added a tiny c that broke the tie, so the
final rounding never saw an exact midpoint: a ties-away final rounding
in simd_fma::fma_f64 passed the whole group on the v2 arm (it was caught
only by simd_fma's own unit test). Round 4 now cancels the product's
middle terms so a*b + c is exactly halfway between two f64 neighbours,
with the even neighbour below on some lanes and above on others, in both
signs, over 2^-900..2^900.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01G22yT6htkcdyXsihxxXdrg
@AdaWorldAPI
AdaWorldAPI force-pushed the fix/wasm-fused-mul-add-20261024 branch from fdaf898 to b41760d Compare October 24, 2026 17:46
@AdaWorldAPI
AdaWorldAPI merged commit a78de7b into master Oct 24, 2026
29 checks passed
AdaWorldAPI pushed a commit that referenced this pull request Oct 27, 2026
… item closed

The 2026-10-11 open item and coresearch row X11 predate #374, which made wasm
mul_add round once and was never recorded here. New entry: fast domains, the
measured margin of the F64 guard, the edge fixtures and their disable runs,
V8 timings before and after the exact-product repair, and what is not claimed.
X11 regraded CLOSED.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01XAqKMPh2Seh2WLccCksAdY
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants