diff --git a/library/compiler-builtins/libm/src/math/cbrt.rs b/library/compiler-builtins/libm/src/math/cbrt.rs index fbf81f77d2e3..8560d37abf32 100644 --- a/library/compiler-builtins/libm/src/math/cbrt.rs +++ b/library/compiler-builtins/libm/src/math/cbrt.rs @@ -5,12 +5,15 @@ */ use super::Float; -use super::fenv::Rounding; -use super::support::cold_path; +use super::support::{FpResult, Round, cold_path}; /// Compute the cube root of the argument. #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn cbrt(x: f64) -> f64 { + cbrt_round(x, Round::Nearest).val +} + +pub fn cbrt_round(x: f64, round: Round) -> FpResult { const ESCALE: [f64; 3] = [ 1.0, hf64!("0x1.428a2f98d728bp+0"), /* 2^(1/3) */ @@ -33,8 +36,6 @@ pub fn cbrt(x: f64) -> f64 { let off = [hf64!("0x1p-53"), 0.0, 0.0, 0.0]; - let rm = Rounding::get(); - /* rm=0 for rounding to nearest, and other values for directed roundings */ let hx: u64 = x.to_bits(); let mut mant: u64 = hx & f64::SIG_MASK; @@ -51,7 +52,7 @@ pub fn cbrt(x: f64) -> f64 { to that for x a signaling NaN, it correctly triggers the invalid exception. */ if e == f64::EXP_SAT || ix == 0 { - return x + x; + return FpResult::ok(x + x); } let nz = ix.leading_zeros() - 11; /* subnormal */ @@ -124,8 +125,8 @@ pub fn cbrt(x: f64) -> f64 { * from ulp(1); * for rounding to nearest, ady0 is tiny when dy is near from 1/2 ulp(1), * or from 3/2 ulp(1). */ - let mut ady0: f64 = (ady - off[rm as usize]).abs(); - let mut ady1: f64 = (ady - (hf64!("0x1p-52") + off[rm as usize])).abs(); + let mut ady0: f64 = (ady - off[round as usize]).abs(); + let mut ady1: f64 = (ady - (hf64!("0x1p-52") + off[round as usize])).abs(); if ady0 < hf64!("0x1p-75") || ady1 < hf64!("0x1p-75") { cold_path(); @@ -140,8 +141,8 @@ pub fn cbrt(x: f64) -> f64 { dy = (y1 - y) - dy; y1 = y; ady = dy.abs(); - ady0 = (ady - off[rm as usize]).abs(); - ady1 = (ady - (hf64!("0x1p-52") + off[rm as usize])).abs(); + ady0 = (ady - off[round as usize]).abs(); + ady1 = (ady - (hf64!("0x1p-52") + off[round as usize])).abs(); if ady0 < hf64!("0x1p-98") || ady1 < hf64!("0x1p-98") { cold_path(); @@ -157,7 +158,7 @@ pub fn cbrt(x: f64) -> f64 { y1 = hf64!("0x1.de87aa837820fp+0").copysign(zz); } - if rm != Rounding::Nearest { + if round != Round::Nearest { let wlist = [ (hf64!("0x1.3a9ccd7f022dbp+0"), hf64!("0x1.1236160ba9b93p+0")), // ~ 0x1.1236160ba9b930000000000001e7e8fap+0 (hf64!("0x1.7845d2faac6fep+0"), hf64!("0x1.23115e657e49cp+0")), // ~ 0x1.23115e657e49c0000000000001d7a799p+0 @@ -170,7 +171,7 @@ pub fn cbrt(x: f64) -> f64 { for (a, b) in wlist { if azz == a { - let tmp = if rm as u64 + sign == 2 { hf64!("0x1p-52") } else { 0.0 }; + let tmp = if round as u64 + sign == 2 { hf64!("0x1p-52") } else { 0.0 }; y1 = (b + tmp).copysign(zz); } } @@ -194,7 +195,7 @@ pub fn cbrt(x: f64) -> f64 { } } - f64::from_bits(cvt3) + FpResult::ok(f64::from_bits(cvt3)) } fn fmaf64(x: f64, y: f64, z: f64) -> f64 { diff --git a/library/compiler-builtins/libm/src/math/fenv.rs b/library/compiler-builtins/libm/src/math/fenv.rs deleted file mode 100644 index 328c9f3467c6..000000000000 --- a/library/compiler-builtins/libm/src/math/fenv.rs +++ /dev/null @@ -1,49 +0,0 @@ -// src: musl/src/fenv/fenv.c -/* Dummy functions for archs lacking fenv implementation */ - -pub(crate) const FE_UNDERFLOW: i32 = 0; -pub(crate) const FE_INEXACT: i32 = 0; - -pub(crate) const FE_TONEAREST: i32 = 0; -pub(crate) const FE_DOWNWARD: i32 = 1; -pub(crate) const FE_UPWARD: i32 = 2; -pub(crate) const FE_TOWARDZERO: i32 = 3; - -#[inline] -pub(crate) fn feclearexcept(_mask: i32) -> i32 { - 0 -} - -#[inline] -pub(crate) fn feraiseexcept(_mask: i32) -> i32 { - 0 -} - -#[inline] -pub(crate) fn fetestexcept(_mask: i32) -> i32 { - 0 -} - -#[inline] -pub(crate) fn fegetround() -> i32 { - FE_TONEAREST -} - -#[derive(Clone, Copy, Debug, PartialEq)] -pub(crate) enum Rounding { - Nearest = FE_TONEAREST as isize, - Downward = FE_DOWNWARD as isize, - Upward = FE_UPWARD as isize, - ToZero = FE_TOWARDZERO as isize, -} - -impl Rounding { - pub(crate) fn get() -> Self { - match fegetround() { - x if x == FE_DOWNWARD => Self::Downward, - x if x == FE_UPWARD => Self::Upward, - x if x == FE_TOWARDZERO => Self::ToZero, - _ => Self::Nearest, - } - } -} diff --git a/library/compiler-builtins/libm/src/math/generic/fma.rs b/library/compiler-builtins/libm/src/math/generic/fma.rs index a40d7aaaf5e0..821aee09028b 100644 --- a/library/compiler-builtins/libm/src/math/generic/fma.rs +++ b/library/compiler-builtins/libm/src/math/generic/fma.rs @@ -1,12 +1,7 @@ /* SPDX-License-Identifier: MIT */ /* origin: musl src/math/{fma,fmaf}.c. Ported to generic Rust algorithm in 2025, TG. */ -use core::{f32, f64}; - -use super::super::fenv::{ - FE_INEXACT, FE_TONEAREST, FE_UNDERFLOW, feclearexcept, fegetround, feraiseexcept, fetestexcept, -}; -use super::super::support::{DInt, HInt, IntTy}; +use super::super::support::{DInt, FpResult, HInt, IntTy, Round, Status}; use super::super::{CastFrom, CastInto, DFloat, Float, HFloat, Int, MinInt}; /// Fused multiply-add that works when there is not a larger float size available. Currently this @@ -14,7 +9,18 @@ use super::super::{CastFrom, CastInto, DFloat, Float, HFloat, Int, MinInt}; #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)] pub fn fma(x: F, y: F, z: F) -> F where - F: Float + FmaHelper, + F: Float, + F: CastFrom, + F: CastFrom, + F::Int: HInt, + u32: CastInto, +{ + fma_round(x, y, z, Round::Nearest).val +} + +pub fn fma_round(x: F, y: F, z: F, _round: Round) -> FpResult +where + F: Float, F: CastFrom, F: CastFrom, F::Int: HInt, @@ -30,16 +36,16 @@ where if nx.is_zero_nan_inf() || ny.is_zero_nan_inf() { // Value will overflow, defer to non-fused operations. - return x * y + z; + return FpResult::ok(x * y + z); } if nz.is_zero_nan_inf() { if nz.is_zero() { // Empty add component means we only need to multiply. - return x * y; + return FpResult::ok(x * y); } // `z` is NaN or infinity, which sets the result. - return z; + return FpResult::ok(z); } // multiply: r = x * y @@ -147,7 +153,7 @@ where } } else { // exact +/- 0.0 - return x * y + z; + return FpResult::ok(x * y + z); } e -= d; @@ -168,6 +174,8 @@ where // Unbiased exponent for the maximum value of `r` let max_pow = F::BITS - 1 + F::EXP_BIAS; + let mut status = Status::OK; + if e < -(max_pow as i32 - 2) { // Result is subnormal before rounding if e == -(max_pow as i32 - 1) { @@ -178,7 +186,9 @@ where if r == c { // Min normal after rounding, - return r.raise_underflow_as_min_positive(); + status.set_underflow(true); + r = F::MIN_POSITIVE_NORMAL.copysign(r); + return FpResult::new(r, status); } if (rhi << (F::SIG_BITS + 1)) != zero { @@ -195,7 +205,7 @@ where // Remove the top bit r = F::cast_from(2i8) * r - c; - r += r.raise_underflow_ret_zero(); + status.set_underflow(true); } } else { // Only round once when scaled @@ -212,12 +222,22 @@ where } // Use our exponent to scale the final value. - super::scalbn(r, e) + FpResult::new(super::scalbn(r, e), status) } /// Fma implementation when a hardware-backed larger float type is available. For `f32` and `f64`, /// `f64` has enough precision to represent the `f32` in its entirety, except for double rounding. pub fn fma_wide(x: F, y: F, z: F) -> F +where + F: Float + HFloat, + B: Float + DFloat, + B::Int: CastInto, + i32: CastFrom, +{ + fma_wide_round(x, y, z, Round::Nearest).val +} + +pub fn fma_wide_round(x: F, y: F, z: F, round: Round) -> FpResult where F: Float + HFloat, B: Float + DFloat, @@ -244,24 +264,26 @@ where // Or the result is exact || (result - xy == zb && result - zb == xy) // Or the mode is something other than round to nearest - || fegetround() != FE_TONEAREST + || round != Round::Nearest { let min_inexact_exp = (B::EXP_BIAS as i32 + F::EXP_MIN_SUBNORM) as u32; let max_inexact_exp = (B::EXP_BIAS as i32 + F::EXP_MIN) as u32; - if (min_inexact_exp..max_inexact_exp).contains(&re) && fetestexcept(FE_INEXACT) != 0 { - feclearexcept(FE_INEXACT); - // prevent `xy + vz` from being CSE'd with `xy + z` above - let vz: F = force_eval!(z); - result = xy + vz.widen(); - if fetestexcept(FE_INEXACT) != 0 { - feraiseexcept(FE_UNDERFLOW); + let mut status = Status::OK; + + if (min_inexact_exp..max_inexact_exp).contains(&re) && status.inexact() { + // This branch is never hit; requires previous operations to set a status + status.set_inexact(false); + + result = xy + z.widen(); + if status.inexact() { + status.set_underflow(true); } else { - feraiseexcept(FE_INEXACT); + status.set_inexact(true); } } - return result.narrow(); + return FpResult { val: result.narrow(), status }; } let neg = ui >> (B::BITS - 1) != IntTy::::ZERO; @@ -272,7 +294,7 @@ where ui -= one; } - B::from_bits(ui).narrow() + FpResult::ok(B::from_bits(ui).narrow()) } /// Representation of `F` that has handled subnormals. @@ -337,49 +359,13 @@ impl Norm { } } -/// Type-specific helpers that are not needed outside of fma. -pub trait FmaHelper { - /// Raise underflow and return the minimum positive normal value with the sign of `self`. - fn raise_underflow_as_min_positive(self) -> Self; - /// Raise underflow and return zero. - fn raise_underflow_ret_zero(self) -> Self; -} - -impl FmaHelper for f64 { - fn raise_underflow_as_min_positive(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_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_as_min_positive(self) -> Self { - f128::MIN_POSITIVE.copysign(self) - } - - fn raise_underflow_ret_zero(self) -> Self { - f128::ZERO - } -} - #[cfg(test)] mod tests { use super::*; fn spec_test() where - F: Float + FmaHelper, + F: Float, F: CastFrom, F: CastFrom, F::Int: HInt, @@ -401,6 +387,29 @@ mod tests { #[test] fn spec_test_f64() { spec_test::(); + + let expect_underflow = [ + ( + hf64!("0x1.0p-1070"), + hf64!("0x1.0p-1070"), + hf64!("0x1.ffffffffffffp-1023"), + hf64!("0x0.ffffffffffff8p-1022"), + ), + ( + // FIXME: we raise underflow but this should only be inexact (based on C and + // `rustc_apfloat`). + hf64!("0x1.0p-1070"), + hf64!("0x1.0p-1070"), + hf64!("-0x1.0p-1022"), + hf64!("-0x1.0p-1022"), + ), + ]; + + for (x, y, z, res) in expect_underflow { + let FpResult { val, status } = fma_round(x, y, z, Round::Nearest); + assert_biteq!(val, res); + assert_eq!(status, Status::UNDERFLOW); + } } #[test] diff --git a/library/compiler-builtins/libm/src/math/generic/sqrt.rs b/library/compiler-builtins/libm/src/math/generic/sqrt.rs index 90d6c01e9171..fdd612493b79 100644 --- a/library/compiler-builtins/libm/src/math/generic/sqrt.rs +++ b/library/compiler-builtins/libm/src/math/generic/sqrt.rs @@ -41,10 +41,23 @@ //! Goldschmidt has the advantage over Newton-Raphson that `sqrt(x)` and `1/sqrt(x)` are //! computed at the same time, i.e. there is no need to calculate `1/sqrt(x)` and invert it. -use super::super::support::{IntTy, cold_path, raise_invalid}; +use super::super::support::{FpResult, IntTy, Round, Status, cold_path}; use super::super::{CastFrom, CastInto, DInt, Float, HInt, Int, MinInt}; pub fn sqrt(x: F) -> F +where + F: Float + SqrtHelper, + F::Int: HInt, + F::Int: From, + F::Int: From, + F::Int: CastInto, + F::Int: CastInto, + u32: CastInto, +{ + sqrt_round(x, Round::Nearest).val +} + +pub fn sqrt_round(x: F, _round: Round) -> FpResult where F: Float + SqrtHelper, F::Int: HInt, @@ -78,17 +91,17 @@ where // +/-0 if ix << 1 == zero { - return x; + return FpResult::ok(x); } // Positive infinity if ix == F::EXP_MASK { - return x; + return FpResult::ok(x); } // NaN or negative if ix > F::EXP_MASK { - return raise_invalid(x); + return FpResult::new(F::NAN, Status::INVALID); } // Normalize subnormals by multiplying by 1.0 << SIG_BITS (e.g. 0x1p52 for doubles). @@ -215,7 +228,7 @@ where y = y + t; } - y + FpResult::ok(y) } /// Multiply at the wider integer size, returning the high half. @@ -329,7 +342,7 @@ impl SqrtHelper for f128 { /// A U0.16 representation of `1/sqrt(x)`. /// -// / The index is a 7-bit number consisting of a single exponent bit and 6 bits of significand. +/// The index is a 7-bit number consisting of a single exponent bit and 6 bits of significand. #[rustfmt::skip] static RSQRT_TAB: [u16; 128] = [ 0xb451, 0xb2f0, 0xb196, 0xb044, 0xaef9, 0xadb6, 0xac79, 0xab43, @@ -354,7 +367,7 @@ static RSQRT_TAB: [u16; 128] = [ mod tests { use super::*; - /// Test against edge cases from https://en.cppreference.com/w/cpp/numeric/math/sqrt + /// Test behavior specified in IEEE 754 `squareRoot`. fn spec_test() where F: Float + SqrtHelper, @@ -365,11 +378,22 @@ mod tests { F::Int: CastInto, u32: CastInto, { - // Not Asserted: FE_INVALID exception is raised if argument is negative. - assert!(sqrt(F::NEG_ONE).is_nan()); - assert!(sqrt(F::NAN).is_nan()); - for f in [F::ZERO, F::NEG_ZERO, F::INFINITY].iter().copied() { - assert_biteq!(sqrt(f), f); + // Values that should return a NaN and raise invalid + let nan = [F::NEG_INFINITY, F::NEG_ONE, F::NAN, F::MIN]; + + // Values that return unaltered + let roundtrip = [F::ZERO, F::NEG_ZERO, F::INFINITY]; + + for x in nan { + let FpResult { val, status } = sqrt_round(x, Round::Nearest); + assert!(val.is_nan()); + assert!(status == Status::INVALID); + } + + for x in roundtrip { + let FpResult { val, status } = sqrt_round(x, Round::Nearest); + assert_biteq!(val, x); + assert!(status == Status::OK); } } diff --git a/library/compiler-builtins/libm/src/math/mod.rs b/library/compiler-builtins/libm/src/math/mod.rs index e32045021dae..ae4a278f287a 100644 --- a/library/compiler-builtins/libm/src/math/mod.rs +++ b/library/compiler-builtins/libm/src/math/mod.rs @@ -94,7 +94,6 @@ cfg_if! { // Private modules mod arch; mod expo2; -mod fenv; mod k_cos; mod k_cosf; mod k_expo2; diff --git a/library/compiler-builtins/libm/src/math/support/env.rs b/library/compiler-builtins/libm/src/math/support/env.rs new file mode 100644 index 000000000000..7244381dacc6 --- /dev/null +++ b/library/compiler-builtins/libm/src/math/support/env.rs @@ -0,0 +1,118 @@ +//! Support for rounding directions and status flags as specified by IEEE 754. +//! +//! Rust does not support the floating point environment so rounding mode is passed as an argument +//! and status flags are returned as part of the result. There is currently not much support for +//! this; most existing ports from musl use a form of `force_eval!` to raise exceptions, but this +//! has no side effects in Rust. Further, correct behavior relies on elementary operations making +//! use of the correct rounding and raising relevant exceptions, which is not the case for Rust. +//! +//! This module exists so no functionality is lost when porting algorithms that respect floating +//! point environment, and so that some functionality may be tested (that which does not rely on +//! side effects from elementary operations). Full support would require wrappers around basic +//! operations, but there is no plan to add this at the current time. + +/// A value combined with a floating point status. +pub struct FpResult { + pub val: T, + #[cfg_attr(not(feature = "unstable-public-internals"), allow(dead_code))] + pub status: Status, +} + +impl FpResult { + pub fn new(val: T, status: Status) -> Self { + Self { val, status } + } + + /// Return `val` with `Status::OK`. + pub fn ok(val: T) -> Self { + Self { val, status: Status::OK } + } +} + +/// IEEE 754 rounding mode, excluding the optional `roundTiesToAway` version of nearest. +/// +/// Integer representation comes from what CORE-MATH uses for indexing. +#[cfg_attr(not(feature = "unstable-public-internals"), allow(dead_code))] +#[derive(Clone, Copy, Debug, PartialEq)] +pub enum Round { + /// IEEE 754 nearest, `roundTiesToEven`. + Nearest = 0, + /// IEEE 754 `roundTowardNegative`. + Negative = 1, + /// IEEE 754 `roundTowardPositive`. + Positive = 2, + /// IEEE 754 `roundTowardZero`. + Zero = 3, +} + +/// IEEE 754 exception status flags. +#[derive(Clone, Copy, Debug, PartialEq)] +pub struct Status(u8); + +impl Status { + /// Default status indicating no errors. + pub const OK: Self = Self(0); + + /// No definable result. + /// + /// Includes: + /// - Any ops on sNaN, with a few exceptions. + /// - `0 * inf`, `inf * 0`. + /// - `fma(0, inf, c)` or `fma(inf, 0, c)`, possibly excluding `c = qNaN`. + /// - `+inf + -inf` and similar (includes subtraction and fma). + /// - `0.0 / 0.0`, `inf / inf` + /// - `remainder(x, y)` if `y == 0.0` or `x == inf`, and neither is NaN. + /// - `sqrt(x)` with `x < 0.0`. + pub const INVALID: Self = Self(1); + + /// Division by zero. + /// + /// The default result for division is +/-inf based on operand sign. For `logB`, the default + /// result is -inf. + /// `x / y` when `x != 0.0` and `y == 0.0`, + + #[cfg_attr(not(feature = "unstable-public-internals"), allow(dead_code))] + pub const DIVIDE_BY_ZERO: Self = Self(1 << 2); + + /// The result exceeds the maximum finite value. + /// + /// The default result depends on rounding mode. `Nearest*` rounds to +/- infinity, sign based + /// on the intermediate result. `Zero` rounds to the signed maximum finite. `Positive` and + /// `Negative` round to signed maximum finite in one direction, signed infinity in the other. + #[cfg_attr(not(feature = "unstable-public-internals"), allow(dead_code))] + pub const OVERFLOW: Self = Self(1 << 3); + + /// The result is subnormal and lost precision. + pub const UNDERFLOW: Self = Self(1 << 4); + + /// The finite-precision result does not match that of infinite precision, and the reason + /// is not represented by one of the other flags. + pub const INEXACT: Self = Self(1 << 5); + + /// True if `UNDERFLOW` is set. + #[cfg_attr(not(feature = "unstable-public-internals"), allow(dead_code))] + pub fn underflow(self) -> bool { + self.0 & Self::UNDERFLOW.0 != 0 + } + + pub fn set_underflow(&mut self, val: bool) { + self.set_flag(val, Self::UNDERFLOW); + } + + /// True if `INEXACT` is set. + pub fn inexact(self) -> bool { + self.0 & Self::INEXACT.0 != 0 + } + + pub fn set_inexact(&mut self, val: bool) { + self.set_flag(val, Self::INEXACT); + } + + fn set_flag(&mut self, val: bool, mask: Self) { + if val { + self.0 |= mask.0; + } else { + self.0 &= !mask.0; + } + } +} diff --git a/library/compiler-builtins/libm/src/math/support/mod.rs b/library/compiler-builtins/libm/src/math/support/mod.rs index 28e9fd4132f4..ee3f2bbdf02d 100644 --- a/library/compiler-builtins/libm/src/math/support/mod.rs +++ b/library/compiler-builtins/libm/src/math/support/mod.rs @@ -1,12 +1,14 @@ #[macro_use] pub mod macros; mod big; +mod env; mod float_traits; pub mod hex_float; mod int_traits; #[allow(unused_imports)] pub use big::{i256, u256}; +pub use env::{FpResult, Round, Status}; #[allow(unused_imports)] pub use float_traits::{DFloat, Float, HFloat, IntTy}; pub(crate) use float_traits::{f32_from_bits, f64_from_bits}; @@ -25,8 +27,3 @@ pub fn cold_path() { #[cfg(intrinsics_enabled)] core::intrinsics::cold_path(); } - -/// Return `x`, first raising `FE_INVALID`. -pub fn raise_invalid(x: F) -> F { - (x - x) / (x - x) -}