Migrate away from nonfunctional fenv stubs

Many routines have some form of handling for rounding mode and floating
point exceptions, which are implemented via a combination of stubs and
`force_eval!` use. This is suboptimal, however, because:

1. Rust does not interact with the floating point environment, so most
   of this code does nothing.
2. The parts of the code that are not dead are not testable.
3. `force_eval!` blocks optimizations, which is unnecessary because we
   do not rely on its side effects.

We cannot ensure correct rounding and exception handling in all cases
without some form of arithmetic operations that are aware of this
behavior. However, the cases where rounding mode is explicitly handled
or exceptions are explicitly raised are testable. Make this possible
here for functions that depend on `math::fenv` by moving the
implementation to a nonpublic function that takes a `Round` and returns
a `Status`.

Link: https://github.com/rust-lang/libm/issues/480

Author: tgross35

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<f64> {
     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 {

src/math/fenv.rs

@@ -1,20 +1,26 @@
 /* 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
 /// is still specialized only for `f64`. Computes `(x * y) + z`.
 #[cfg_attr(all(test, assert_no_panic), no_panic::no_panic)]
 pub fn fma<F>(x: F, y: F, z: F) -> F
 where
-    F: Float + FmaHelper,
+    F: Float,
+    F: CastFrom<F::SignedInt>,
+    F: CastFrom<i8>,
+    F::Int: HInt,
+    u32: CastInto<F::Int>,
+{
+    fma_round(x, y, z, Round::Nearest).val
+}
+
+pub fn fma_round<F>(x: F, y: F, z: F, _round: Round) -> FpResult<F>
+where
+    F: Float,
     F: CastFrom<F::SignedInt>,
     F: CastFrom<i8>,
     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<F, B>(x: F, y: F, z: F) -> F
+where
+    F: Float + HFloat<D = B>,
+    B: Float + DFloat<H = F>,
+    B::Int: CastInto<i32>,
+    i32: CastFrom<i32>,
+{
+    fma_wide_round(x, y, z, Round::Nearest).val
+}
+
+pub fn fma_wide_round<F, B>(x: F, y: F, z: F, round: Round) -> FpResult<F>
 where
     F: Float + HFloat<D = B>,
     B: Float + DFloat<H = F>,
@@ -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::<B>::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<F: Float> Norm<F> {
     }
 }
 
-/// 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<F>()
     where
-        F: Float + FmaHelper,
+        F: Float,
         F: CastFrom<F::SignedInt>,
         F: CastFrom<i8>,
         F::Int: HInt,
@@ -401,6 +387,29 @@ mod tests {
     #[test]
     fn spec_test_f64() {
         spec_test::<f64>();
+
+        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]