Merge pull request #497 from tgross35/edge-case-max-subnorm

fmaf128: fix exponent calculation for subnormals

Author: tgross35

.github/workflows/main.yaml

@@ -273,6 +273,9 @@ jobs:
             exit
           fi
 
+          # Run the non-extensive tests first to catch any easy failures
+          cargo t --profile release-checked -- "$CHANGED"
+
           LIBM_EXTENSIVE_TESTS="$CHANGED" cargo t \
             --features build-mpfr,unstable,force-soft-floats \
             --profile release-checked \

crates/libm-test/src/gen/case_list.rs

@@ -269,15 +269,26 @@ 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
+        &[
+            (
+                // 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")),
+            ),
             (
-                hf128!("-0x1.1966cc01966cc01966cc01966f06p-25"),
-                hf128!("-0x1.669933fe69933fe69933fe6997c9p-16358"),
-                hf128!("-0x0.000000000000000000000000048ap-16382"),
+                // Subnormal edge case that caused a failure
+                (
+                    hf128!("0x0.7ffffffffffffffffffffffffff7p-16382"),
+                    hf128!("0x1.ffffffffffffffffffffffffffffp-1"),
+                    hf128!("0x0.8000000000000000000000000009p-16382"),
+                ),
+                Some(hf128!("0x1.0000000000000000000000000000p-16382")),
             ),
-            Some(hf128!("0x0.c5171470a3ff5e0f68d751491b18p-16382")),
-        )],
+        ],
     );
     v
 }

crates/libm-test/src/gen/edge_cases.rs

@@ -1,9 +1,10 @@
 //! A generator that checks a handful of cases near infinities, zeros, asymptotes, and NaNs.
 
-use libm::support::{CastInto, Float, Int};
+use libm::support::{CastInto, Float, Int, MinInt};
 
 use crate::domain::get_domain;
 use crate::gen::KnownSize;
+use crate::op::OpITy;
 use crate::run_cfg::{check_near_count, check_point_count};
 use crate::{BaseName, CheckCtx, FloatExt, FloatTy, MathOp, test_log};
 
@@ -21,6 +22,7 @@ where
     Op: MathOp,
 {
     let mut ret = Vec::new();
+    let one = OpITy::<Op>::ONE;
     let values = &mut ret;
     let domain = get_domain::<_, i8>(ctx.fn_ident, argnum).unwrap_float();
     let domain_start = domain.range_start();
@@ -51,6 +53,22 @@ where
     values.push(Op::FTy::NAN);
     values.extend(Op::FTy::consts().iter());
 
+    // Check around the maximum subnormal value
+    let sub_max = Op::FTy::from_bits(Op::FTy::SIG_MASK);
+    count_up(sub_max, near_points, values);
+    count_down(sub_max, near_points, values);
+    count_up(-sub_max, near_points, values);
+    count_down(-sub_max, near_points, values);
+
+    // Check a few values around the subnormal range
+    for shift in (0..Op::FTy::SIG_BITS).step_by(Op::FTy::SIG_BITS as usize / 5) {
+        let v = Op::FTy::from_bits(one << shift);
+        count_up(v, 2, values);
+        count_down(v, 2, values);
+        count_up(-v, 2, values);
+        count_down(-v, 2, values);
+    }
+
     // Check around asymptotes
     if let Some(f) = domain.check_points {
         let iter = f();

crates/libm-test/src/run_cfg.rs

@@ -342,7 +342,7 @@ pub fn check_near_count(ctx: &CheckCtx) -> u64 {
             x => panic!("unexpected argument count {x}"),
         }
     } else {
-        10
+        8
     }
 }
 

src/math/generic/fma.rs

@@ -146,6 +146,7 @@ where
         // exact +/- 0.0
         return x * y + z;
     }
+
     e -= d;
 
     // Use int->float conversion to populate the significand.
@@ -174,7 +175,7 @@ where
 
             if r == c {
                 // Min normal after rounding,
-                return r.raise_underflow_ret_self();
+                return r.raise_underflow_as_min_positive();
             }
 
             if (rhi << (F::SIG_BITS + 1)) != zero {
@@ -275,12 +276,14 @@ impl<F: Float> Norm<F> {
 
 /// Type-specific helpers that are not needed outside of fma.
 pub trait FmaHelper {
-    fn raise_underflow_ret_self(self) -> Self;
+    /// 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_ret_self(self) -> Self {
+    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 */
@@ -298,8 +301,8 @@ impl FmaHelper for f64 {
 
 #[cfg(f128_enabled)]
 impl FmaHelper for f128 {
-    fn raise_underflow_ret_self(self) -> Self {
-        self
+    fn raise_underflow_as_min_positive(self) -> Self {
+        f128::MIN_POSITIVE.copysign(self)
     }
 
     fn raise_underflow_ret_zero(self) -> Self {