Add `fmaf128`

Resolve all remaining `f64`-specific items in the generic version of
`fma`, then expose `fmaf128`.

14 files changed, 237 insertions, 67 deletions

diff --git a/library/compiler-builtins/libm/crates/libm-macros/src/shared.rs b/library/compiler-builtins/libm/crates/libm-macros/src/shared.rs
index da16cd8e287..48d19c50d19 100644
--- a/library/compiler-builtins/libm/crates/libm-macros/src/shared.rs
+++ b/library/compiler-builtins/libm/crates/libm-macros/src/shared.rs
@@ -107,6 +107,13 @@ const ALL_OPERATIONS_NESTED: &[(FloatTy, Signature, Option<Signature>, &[&str])]
         &["fma"],
     ),
     (
+        // `(f128, f128, f128) -> f128`
+        FloatTy::F128,
+        Signature { args: &[Ty::F128, Ty::F128, Ty::F128], returns: &[Ty::F128] },
+        None,
+        &["fmaf128"],
+    ),
+    (
         // `(f32) -> i32`
         FloatTy::F32,
         Signature { args: &[Ty::F32], returns: &[Ty::I32] },
diff --git a/library/compiler-builtins/libm/crates/libm-test/benches/icount.rs b/library/compiler-builtins/libm/crates/libm-test/benches/icount.rs
index 53ecb5a37c4..c41cef24e54 100644
--- a/library/compiler-builtins/libm/crates/libm-test/benches/icount.rs
+++ b/library/compiler-builtins/libm/crates/libm-test/benches/icount.rs
@@ -108,6 +108,7 @@ main!(
     icount_bench_floorf16_group,
     icount_bench_floorf_group,
     icount_bench_fma_group,
+    icount_bench_fmaf128_group,
     icount_bench_fmaf_group,
     icount_bench_fmax_group,
     icount_bench_fmaxf128_group,
diff --git a/library/compiler-builtins/libm/crates/libm-test/benches/random.rs b/library/compiler-builtins/libm/crates/libm-test/benches/random.rs
index 66486a56a2a..6e8a334795a 100644
--- a/library/compiler-builtins/libm/crates/libm-test/benches/random.rs
+++ b/library/compiler-builtins/libm/crates/libm-test/benches/random.rs
@@ -127,6 +127,7 @@ libm_macros::for_each_function! {
         | fdimf16
         | floorf128
         | floorf16
+        | fmaf128
         | fmaxf128
         | fmaxf16
         | fminf128
diff --git a/library/compiler-builtins/libm/crates/libm-test/src/gen/case_list.rs b/library/compiler-builtins/libm/crates/libm-test/src/gen/case_list.rs
index 9720f68e961..302d5c39184 100644
--- a/library/compiler-builtins/libm/crates/libm-test/src/gen/case_list.rs
+++ b/library/compiler-builtins/libm/crates/libm-test/src/gen/case_list.rs
@@ -6,6 +6,9 @@
 //!
 //! This is useful for adding regression tests or expected failures.
 
+#[cfg(f128_enabled)]
+use libm::hf128;
+
 use crate::{CheckBasis, CheckCtx, GeneratorKind, MathOp, op};
 
 pub struct TestCase<Op: MathOp> {
@@ -250,7 +253,7 @@ fn fma_cases() -> Vec<TestCase<op::fma::Routine>> {
     TestCase::append_pairs(
         &mut v,
         &[
-            // Previously failure with incorrect sign
+            // Previous failure with incorrect sign
             ((5e-324, -5e-324, 0.0), Some(-0.0)),
         ],
     );
@@ -261,6 +264,24 @@ fn fmaf_cases() -> Vec<TestCase<op::fmaf::Routine>> {
     vec![]
 }
 
+#[cfg(f128_enabled)]
+fn fmaf128_cases() -> Vec<TestCase<op::fmaf128::Routine>> {
+    let mut v = vec![];
+    TestCase::append_pairs(
+        &mut v,
+        &[(
+            // Tricky rounding case that previously failed in extensive tests
+            (
+                hf128!("-0x1.1966cc01966cc01966cc01966f06p-25"),
+                hf128!("-0x1.669933fe69933fe69933fe6997c9p-16358"),
+                hf128!("-0x0.000000000000000000000000048ap-16382"),
+            ),
+            Some(hf128!("0x0.c5171470a3ff5e0f68d751491b18p-16382")),
+        )],
+    );
+    v
+}
+
 fn fmax_cases() -> Vec<TestCase<op::fmax::Routine>> {
     vec![]
 }
diff --git a/library/compiler-builtins/libm/crates/libm-test/src/mpfloat.rs b/library/compiler-builtins/libm/crates/libm-test/src/mpfloat.rs
index ab77d541c81..f4a9ff7ffd5 100644
--- a/library/compiler-builtins/libm/crates/libm-test/src/mpfloat.rs
+++ b/library/compiler-builtins/libm/crates/libm-test/src/mpfloat.rs
@@ -196,7 +196,7 @@ libm_macros::for_each_function! {
         expm1 | expm1f => exp_m1,
         fabs | fabsf => abs,
         fdim | fdimf | fdimf16 | fdimf128  => positive_diff,
-        fma | fmaf => mul_add,
+        fma | fmaf | fmaf128 => mul_add,
         fmax | fmaxf | fmaxf16 | fmaxf128 => max,
         fmin | fminf | fminf16 | fminf128 => min,
         lgamma | lgammaf => ln_gamma,
diff --git a/library/compiler-builtins/libm/crates/libm-test/src/precision.rs b/library/compiler-builtins/libm/crates/libm-test/src/precision.rs
index 596f91fe1ef..20aa96b6aba 100644
--- a/library/compiler-builtins/libm/crates/libm-test/src/precision.rs
+++ b/library/compiler-builtins/libm/crates/libm-test/src/precision.rs
@@ -560,3 +560,5 @@ impl MaybeOverride<(f128, i32)> for SpecialCase {}
 
 impl MaybeOverride<(f32, f32, f32)> for SpecialCase {}
 impl MaybeOverride<(f64, f64, f64)> for SpecialCase {}
+#[cfg(f128_enabled)]
+impl MaybeOverride<(f128, f128, f128)> for SpecialCase {}
diff --git a/library/compiler-builtins/libm/crates/libm-test/tests/compare_built_musl.rs b/library/compiler-builtins/libm/crates/libm-test/tests/compare_built_musl.rs
index 927cb25afc4..7fa77e832b1 100644
--- a/library/compiler-builtins/libm/crates/libm-test/tests/compare_built_musl.rs
+++ b/library/compiler-builtins/libm/crates/libm-test/tests/compare_built_musl.rs
@@ -99,6 +99,7 @@ libm_macros::for_each_function! {
         fdimf16,
         floorf128,
         floorf16,
+        fmaf128,
         fmaxf128,
         fmaxf16,
         fminf128,
diff --git a/library/compiler-builtins/libm/crates/util/src/main.rs b/library/compiler-builtins/libm/crates/util/src/main.rs
index e5d6f374ad8..0f845a1c4d0 100644
--- a/library/compiler-builtins/libm/crates/util/src/main.rs
+++ b/library/compiler-builtins/libm/crates/util/src/main.rs
@@ -96,6 +96,7 @@ fn do_eval(basis: &str, op: &str, inputs: &[&str]) {
             | fdimf16
             | floorf128
             | floorf16
+            | fmaf128
             | fmaxf128
             | fmaxf16
             | fminf128
diff --git a/library/compiler-builtins/libm/etc/function-definitions.json b/library/compiler-builtins/libm/etc/function-definitions.json
index 243862075ff..5742ed585f5 100644
--- a/library/compiler-builtins/libm/etc/function-definitions.json
+++ b/library/compiler-builtins/libm/etc/function-definitions.json
@@ -356,6 +356,13 @@
         ],
         "type": "f32"
     },
+    "fmaf128": {
+        "sources": [
+            "src/math/fmaf128.rs",
+            "src/math/generic/fma.rs"
+        ],
+        "type": "f128"
+    },
     "fmax": {
         "sources": [
             "src/math/fmax.rs",
diff --git a/library/compiler-builtins/libm/etc/function-list.txt b/library/compiler-builtins/libm/etc/function-list.txt
index c92eaf9e23e..1c9c5e3bc33 100644
--- a/library/compiler-builtins/libm/etc/function-list.txt
+++ b/library/compiler-builtins/libm/etc/function-list.txt
@@ -53,6 +53,7 @@ floorf128
 floorf16
 fma
 fmaf
+fmaf128
 fmax
 fmaxf
 fmaxf128
diff --git a/library/compiler-builtins/libm/src/libm_helper.rs b/library/compiler-builtins/libm/src/libm_helper.rs
index 0768839c779..68f1fb36266 100644
--- a/library/compiler-builtins/libm/src/libm_helper.rs
+++ b/library/compiler-builtins/libm/src/libm_helper.rs
@@ -208,6 +208,7 @@ libm_helper! {
         (fn fabs(x: f128) -> (f128);                => fabsf128);
         (fn fdim(x: f128, y: f128) -> (f128);       => fdimf128);
         (fn floor(x: f128) -> (f128);               => floorf128);
+        (fn fmaf128(x: f128, y: f128, z: f128) -> (f128); => fmaf128);
         (fn fmax(x: f128, y: f128) -> (f128);       => fmaxf128);
         (fn fmin(x: f128, y: f128) -> (f128);       => fminf128);
         (fn fmod(x: f128, y: f128) -> (f128);       => fmodf128);
diff --git a/library/compiler-builtins/libm/src/math/fmaf128.rs b/library/compiler-builtins/libm/src/math/fmaf128.rs
new file mode 100644
index 00000000000..50f7360deb4
--- /dev/null
+++ b/library/compiler-builtins/libm/src/math/fmaf128.rs
@@ -0,0 +1,7 @@
+/// Fused multiply add (f128)
+///
+/// Computes `(x*y)+z`, rounded as one ternary operation (i.e. calculated with infinite precision).
+#[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
+pub fn fmaf128(x: f128, y: f128, z: f128) -> f128 {
+    return super::generic::fma(x, y, z);
+}
diff --git a/library/compiler-builtins/libm/src/math/generic/fma.rs b/library/compiler-builtins/libm/src/math/generic/fma.rs
index b0e2117ea09..ac53acadfe1 100644
--- a/library/compiler-builtins/libm/src/math/generic/fma.rs
+++ b/library/compiler-builtins/libm/src/math/generic/fma.rs
@@ -1,10 +1,11 @@
+/* SPDX-License-Identifier: MIT */
+/* origin: musl src/math/fma.c. Ported to generic Rust algorithm in 2025, TG. */
+
 use core::{f32, f64};
 
 use super::super::support::{DInt, HInt, IntTy};
 use super::super::{CastFrom, CastInto, Float, Int, MinInt};
 
-const ZEROINFNAN: i32 = 0x7ff - 0x3ff - 52 - 1;
-
 /// Fused multiply-add that works when there is not a larger float size available. Currently this
 /// is still specialized only for `f64`. Computes `(x * y) + z`.
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
@@ -18,79 +19,99 @@ where
 {
     let one = IntTy::<F>::ONE;
     let zero = IntTy::<F>::ZERO;
-    let magic = F::from_parts(false, F::BITS - 1 + F::EXP_BIAS, zero);
 
-    /* normalize so top 10bits and last bit are 0 */
+    // Normalize such that the top of the mantissa is zero and we have a guard bit.
     let nx = Norm::from_float(x);
     let ny = Norm::from_float(y);
     let nz = Norm::from_float(z);
 
-    if nx.e >= ZEROINFNAN || ny.e >= ZEROINFNAN {
+    if nx.is_zero_nan_inf() || ny.is_zero_nan_inf() {
+        // Value will overflow, defer to non-fused operations.
         return x * y + z;
     }
-    if nz.e >= ZEROINFNAN {
-        if nz.e > ZEROINFNAN {
-            /* z==0 */
+
+    if nz.is_zero_nan_inf() {
+        if nz.is_zero() {
+            // Empty add component means we only need to multiply.
             return x * y;
         }
+        // `z` is NaN or infinity, which sets the result.
         return z;
     }
 
-    /* mul: r = x*y */
+    // multiply: r = x * y
     let zhi: F::Int;
     let zlo: F::Int;
     let (mut rlo, mut rhi) = nx.m.widen_mul(ny.m).lo_hi();
 
-    /* either top 20 or 21 bits of rhi and last 2 bits of rlo are 0 */
-
-    /* align exponents */
+    // Exponent result of multiplication
     let mut e: i32 = nx.e + ny.e;
+    // Needed shift to align `z` to the multiplication result
     let mut d: i32 = nz.e - e;
     let sbits = F::BITS as i32;
 
-    /* shift bits z<<=kz, r>>=kr, so kz+kr == d, set e = e+kr (== ez-kz) */
+    // Scale `z`. Shift `z <<= kz`, `r >>= kr`, so `kz+kr == d`, set `e = e+kr` (== ez-kz)
     if d > 0 {
+        // The magnitude of `z` is larger than `x * y`
         if d < sbits {
+            // Maximum shift of one `F::BITS` means shifted `z` will fit into `2 * F::BITS`. Shift
+            // it into `(zhi, zlo)`. No exponent adjustment necessary.
             zlo = nz.m << d;
             zhi = nz.m >> (sbits - d);
         } else {
+            // Shift larger than `sbits`, `z` only needs the top half `zhi`. Place it there (acts
+            // as a shift by `sbits`).
             zlo = zero;
             zhi = nz.m;
-            e = nz.e - sbits;
             d -= sbits;
+
+            // `z`'s exponent is large enough that it now needs to be taken into account.
+            e = nz.e - sbits;
+
             if d == 0 {
+                // Exactly `sbits`, nothing to do
             } else if d < sbits {
-                rlo = (rhi << (sbits - d))
-                    | (rlo >> d)
-                    | IntTy::<F>::from((rlo << (sbits - d)) != zero);
+                // Remaining shift fits within `sbits`. Leave `z` in place, shift `x * y`
+                rlo = (rhi << (sbits - d)) | (rlo >> d);
+                // Set the sticky bit
+                rlo |= IntTy::<F>::from((rlo << (sbits - d)) != zero);
                 rhi = rhi >> d;
             } else {
+                // `z`'s magnitude is enough that `x * y` is irrelevant. It was nonzero, so set
+                // the sticky bit.
                 rlo = one;
                 rhi = zero;
             }
         }
     } else {
+        // `z`'s magnitude once shifted fits entirely within `zlo`
         zhi = zero;
         d = -d;
         if d == 0 {
+            // No shift needed
             zlo = nz.m;
         } else if d < sbits {
-            zlo = (nz.m >> d) | IntTy::<F>::from((nz.m << (sbits - d)) != zero);
+            // Shift s.t. `nz.m` fits into `zlo`
+            let sticky = IntTy::<F>::from((nz.m << (sbits - d)) != zero);
+            zlo = (nz.m >> d) | sticky;
         } else {
+            // Would be entirely shifted out, only set the sticky bit
             zlo = one;
         }
     }
 
-    /* add */
+    /* addition */
+
     let mut neg = nx.neg ^ ny.neg;
     let samesign: bool = !neg ^ nz.neg;
-    let mut nonzero: i32 = 1;
+    let mut rhi_nonzero = true;
+
     if samesign {
-        /* r += z */
+        // r += z
         rlo = rlo.wrapping_add(zlo);
         rhi += zhi + IntTy::<F>::from(rlo < zlo);
     } else {
-        /* r -= z */
+        // r -= z
         let (res, borrow) = rlo.overflowing_sub(zlo);
         rlo = res;
         rhi = rhi.wrapping_sub(zhi.wrapping_add(IntTy::<F>::from(borrow)));
@@ -99,129 +120,226 @@ where
             rhi = rhi.signed().wrapping_neg().unsigned() - IntTy::<F>::from(rlo != zero);
             neg = !neg;
         }
-        nonzero = (rhi != zero) as i32;
+        rhi_nonzero = rhi != zero;
     }
 
-    /* set rhi to top 63bit of the result (last bit is sticky) */
-    if nonzero != 0 {
+    /* Construct result */
+
+    // Shift result into `rhi`, left-aligned. Last bit is sticky
+    if rhi_nonzero {
+        // `d` > 0, need to shift both `rhi` and `rlo` into result
         e += sbits;
         d = rhi.leading_zeros() as i32 - 1;
-        /* note: d > 0 */
-        rhi = (rhi << d) | (rlo >> (sbits - d)) | IntTy::<F>::from((rlo << d) != zero);
+        rhi = (rhi << d) | (rlo >> (sbits - d));
+        // Update sticky
+        rhi |= IntTy::<F>::from((rlo << d) != zero);
     } else if rlo != zero {
+        // `rhi` is zero, `rlo` is the entire result and needs to be shifted
         d = rlo.leading_zeros() as i32 - 1;
         if d < 0 {
+            // Shift and set sticky
             rhi = (rlo >> 1) | (rlo & one);
         } else {
             rhi = rlo << d;
         }
     } else {
-        /* exact +-0 */
+        // exact +/- 0.0
         return x * y + z;
     }
     e -= d;
 
-    /* convert to double */
-    let mut i: F::SignedInt = rhi.signed(); /* i is in [1<<62,(1<<63)-1] */
+    // Use int->float conversion to populate the significand.
+    // i is in [1 << (BITS - 2), (1 << (BITS - 1)) - 1]
+    let mut i: F::SignedInt = rhi.signed();
+
     if neg {
         i = -i;
     }
 
-    let mut r: F = F::cast_from_lossy(i); /* |r| is in [0x1p62,0x1p63] */
+    // `|r|` is in `[0x1p62,0x1p63]` for `f64`
+    let mut r: F = F::cast_from_lossy(i);
+
+    /* Account for subnormal and rounding */
+
+    // Unbiased exponent for the maximum value of `r`
+    let max_pow = F::BITS - 1 + F::EXP_BIAS;
 
-    if e < -(F::EXP_BIAS as i32 - 1) - (sbits - 2) {
-        /* result is subnormal before rounding */
-        if e == -(F::EXP_BIAS as i32 - 1) - (sbits - 1) {
-            let mut c: F = magic;
+    if e < -(max_pow as i32 - 2) {
+        // Result is subnormal before rounding
+        if e == -(max_pow as i32 - 1) {
+            let mut c = F::from_parts(false, max_pow, zero);
             if neg {
                 c = -c;
             }
+
             if r == c {
-                /* min normal after rounding, underflow depends
-                 * on arch behaviour which can be imitated by
-                 * a double to float conversion */
-                return r.raise_underflow();
+                // Min normal after rounding,
+                return r.raise_underflow_ret_self();
             }
-            /* one bit is lost when scaled, add another top bit to
-             * only round once at conversion if it is inexact */
-            if (rhi << F::SIG_BITS) != zero {
-                let iu: F::Int = (rhi >> 1) | (rhi & one) | (one << 62);
+
+            if (rhi << (F::SIG_BITS + 1)) != zero {
+                // Account for truncated bits. One bit will be lost in the `scalbn` call, add
+                // another top bit to avoid double rounding if inexact.
+                let iu: F::Int = (rhi >> 1) | (rhi & one) | (one << (F::BITS - 2));
                 i = iu.signed();
+
                 if neg {
                     i = -i;
                 }
+
                 r = F::cast_from_lossy(i);
-                r = F::cast_from(2i8) * r - c; /* remove top bit */
 
-                /* raise underflow portably, such that it
-                 * cannot be optimized away */
-                r += r.raise_underflow2();
+                // Remove the top bit
+                r = F::cast_from(2i8) * r - c;
+                r += r.raise_underflow_ret_zero();
             }
         } else {
-            /* only round once when scaled */
-            d = 10;
-            i = (((rhi >> d) | IntTy::<F>::from(rhi << (F::BITS as i32 - d) != zero)) << d)
-                .signed();
+            // Only round once when scaled
+            d = F::EXP_BITS as i32 - 1;
+            let sticky = IntTy::<F>::from(rhi << (F::BITS as i32 - d) != zero);
+            i = (((rhi >> d) | sticky) << d).signed();
+
             if neg {
                 i = -i;
             }
-            r = F::cast_from(i);
+
+            r = F::cast_from_lossy(i);
         }
     }
 
+    // Use our exponent to scale the final value.
     super::scalbn(r, e)
 }
 
 /// Representation of `F` that has handled subnormals.
+#[derive(Clone, Copy, Debug)]
 struct Norm<F: Float> {
-    /// Normalized significand with one guard bit.
+    /// Normalized significand with one guard bit, unsigned.
     m: F::Int,
-    /// Unbiased exponent, normalized.
+    /// Exponent of the mantissa such that `m * 2^e = x`. Accounts for the shift in the mantissa
+    /// and the guard bit; that is, 1.0 will normalize as `m = 1 << 53` and `e = -53`.
     e: i32,
     neg: bool,
 }
 
 impl<F: Float> Norm<F> {
+    /// Unbias the exponent and account for the mantissa's precision, including the guard bit.
+    const EXP_UNBIAS: u32 = F::EXP_BIAS + F::SIG_BITS + 1;
+
+    /// Values greater than this had a saturated exponent (infinity or NaN), OR were zero and we
+    /// adjusted the exponent such that it exceeds this threashold.
+    const ZERO_INF_NAN: u32 = F::EXP_SAT - Self::EXP_UNBIAS;
+
     fn from_float(x: F) -> Self {
         let mut ix = x.to_bits();
         let mut e = x.exp() as i32;
         let neg = x.is_sign_negative();
         if e == 0 {
             // Normalize subnormals by multiplication
-            let magic = F::from_parts(false, F::BITS - 1 + F::EXP_BIAS, F::Int::ZERO);
-            let scaled = x * magic;
+            let scale_i = F::BITS - 1;
+            let scale_f = F::from_parts(false, scale_i + F::EXP_BIAS, F::Int::ZERO);
+            let scaled = x * scale_f;
             ix = scaled.to_bits();
             e = scaled.exp() as i32;
-            e = if e != 0 { e - (F::BITS as i32 - 1) } else { 0x800 };
+            e = if e == 0 {
+                // If the exponent is still zero, the input was zero. Artifically set this value
+                // such that the final `e` will exceed `ZERO_INF_NAN`.
+                1 << F::EXP_BITS
+            } else {
+                // Otherwise, account for the scaling we just did.
+                e - scale_i as i32
+            };
         }
 
-        e -= F::EXP_BIAS as i32 + 52 + 1;
+        e -= Self::EXP_UNBIAS as i32;
 
+        // Absolute  value, set the implicit bit, and shift to create a guard bit
         ix &= F::SIG_MASK;
         ix |= F::IMPLICIT_BIT;
-        ix <<= 1; // add a guard bit
+        ix <<= 1;
 
         Self { m: ix, e, neg }
     }
+
+    /// True if the value was zero, infinity, or NaN.
+    fn is_zero_nan_inf(self) -> bool {
+        self.e >= Self::ZERO_INF_NAN as i32
+    }
+
+    /// The only value we have
+    fn is_zero(self) -> bool {
+        // The only exponent that strictly exceeds this value is our sentinel value for zero.
+        self.e > Self::ZERO_INF_NAN as i32
+    }
 }
 
 /// Type-specific helpers that are not needed outside of fma.
 pub trait FmaHelper {
-    fn raise_underflow(self) -> Self;
-    fn raise_underflow2(self) -> Self;
+    fn raise_underflow_ret_self(self) -> Self;
+    fn raise_underflow_ret_zero(self) -> Self;
 }
 
 impl FmaHelper for f64 {
-    fn raise_underflow(self) -> Self {
-        let x0_ffffff8p_63 = f64::from_bits(0x3bfffffff0000000); // 0x0.ffffff8p-63
-        let fltmin: f32 = (x0_ffffff8p_63 * f32::MIN_POSITIVE as f64 * self) as f32;
+    fn raise_underflow_ret_self(self) -> Self {
+        /* min normal after rounding, underflow depends
+         * on arch behaviour which can be imitated by
+         * a double to float conversion */
+        let fltmin: f32 = (hf64!("0x0.ffffff8p-63") * f32::MIN_POSITIVE as f64 * self) as f32;
         f64::MIN_POSITIVE / f32::MIN_POSITIVE as f64 * fltmin as f64
     }
 
-    fn raise_underflow2(self) -> Self {
+    fn raise_underflow_ret_zero(self) -> Self {
         /* raise underflow portably, such that it
          * cannot be optimized away */
         let tiny: f64 = f64::MIN_POSITIVE / f32::MIN_POSITIVE as f64 * self;
         (tiny * tiny) * (self - self)
     }
 }
+
+#[cfg(f128_enabled)]
+impl FmaHelper for f128 {
+    fn raise_underflow_ret_self(self) -> Self {
+        self
+    }
+
+    fn raise_underflow_ret_zero(self) -> Self {
+        f128::ZERO
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    fn spec_test<F>()
+    where
+        F: Float + FmaHelper,
+        F: CastFrom<F::SignedInt>,
+        F: CastFrom<i8>,
+        F::Int: HInt,
+        u32: CastInto<F::Int>,
+    {
+        let x = F::from_bits(F::Int::ONE);
+        let y = F::from_bits(F::Int::ONE);
+        let z = F::ZERO;
+
+        // 754-2020 says "When the exact result of (a × b) + c is non-zero yet the result of
+        // fusedMultiplyAdd is zero because of rounding, the zero result takes the sign of the
+        // exact result"
+        assert_biteq!(fma(x, y, z), F::ZERO);
+        assert_biteq!(fma(x, -y, z), F::NEG_ZERO);
+        assert_biteq!(fma(-x, y, z), F::NEG_ZERO);
+        assert_biteq!(fma(-x, -y, z), F::ZERO);
+    }
+
+    #[test]
+    fn spec_test_f64() {
+        spec_test::<f64>();
+    }
+
+    #[test]
+    #[cfg(f128_enabled)]
+    fn spec_test_f128() {
+        spec_test::<f128>();
+    }
+}
diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs
index 7ad808cf752..677ed8d6e93 100644
--- a/library/compiler-builtins/libm/src/math/mod.rs
+++ b/library/compiler-builtins/libm/src/math/mod.rs
@@ -385,6 +385,7 @@ cfg_if! {
         mod fabsf128;
         mod fdimf128;
         mod floorf128;
+        mod fmaf128;
         mod fmaxf128;
         mod fminf128;
         mod fmodf128;
@@ -402,6 +403,7 @@ cfg_if! {
         pub use self::fabsf128::fabsf128;
         pub use self::fdimf128::fdimf128;
         pub use self::floorf128::floorf128;
+        pub use self::fmaf128::fmaf128;
         pub use self::fmaxf128::fmaxf128;
         pub use self::fminf128::fminf128;
         pub use self::fmodf128::fmodf128;