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 da16cd8e2874..48d19c50d19a 100644 --- a/library/compiler-builtins/libm/crates/libm-macros/src/shared.rs +++ b/library/compiler-builtins/libm/crates/libm-macros/src/shared.rs @@ -106,6 +106,13 @@ const ALL_OPERATIONS_NESTED: &[(FloatTy, Signature, Option, &[&str])] None, &["fma"], ), + ( + // `(f128, f128, f128) -> f128` + FloatTy::F128, + Signature { args: &[Ty::F128, Ty::F128, Ty::F128], returns: &[Ty::F128] }, + None, + &["fmaf128"], + ), ( // `(f32) -> i32` FloatTy::F32, 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 53ecb5a37c4b..c41cef24e54d 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 66486a56a2a9..6e8a334795a2 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 9720f68e9612..302d5c391846 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 { @@ -250,7 +253,7 @@ fn fma_cases() -> Vec> { 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> { vec![] } +#[cfg(f128_enabled)] +fn fmaf128_cases() -> Vec> { + 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> { 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 ab77d541c812..f4a9ff7ffd5d 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 596f91fe1efa..20aa96b6aba3 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 927cb25afc47..7fa77e832b17 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 e5d6f374ad85..0f845a1c4d0d 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 243862075ff0..5742ed585f58 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 c92eaf9e23eb..1c9c5e3bc337 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 0768839c7799..68f1fb362666 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 000000000000..50f7360deb45 --- /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 b0e2117ea098..ac53acadfe1a 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::::ONE; let zero = IntTy::::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::::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::::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::::from((nz.m << (sbits - d)) != zero); + // Shift s.t. `nz.m` fits into `zlo` + let sticky = IntTy::::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::::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::::from(borrow))); @@ -99,129 +120,226 @@ where rhi = rhi.signed().wrapping_neg().unsigned() - IntTy::::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::::from((rlo << d) != zero); + rhi = (rhi << d) | (rlo >> (sbits - d)); + // Update sticky + rhi |= IntTy::::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); - 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; + /* Account for subnormal and rounding */ + + // Unbiased exponent for the maximum value of `r` + let max_pow = F::BITS - 1 + F::EXP_BIAS; + + 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(); + r = F::cast_from_lossy(i); + + // 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::::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::::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 { - /// 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 Norm { + /// 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() + where + F: Float + FmaHelper, + F: CastFrom, + F: CastFrom, + F::Int: HInt, + u32: CastInto, + { + 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::(); + } + + #[test] + #[cfg(f128_enabled)] + fn spec_test_f128() { + spec_test::(); + } +} diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs index 7ad808cf7522..677ed8d6e931 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;