fmaf128: fix exponent calculation for subnormals

When `fmaf128` was introduced in [1], it included a bug where `self`
gets returned rather than the expected minimum positive value. Resolve
this and add a regression test.

[1]: #494

Author: tgross35

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
 }

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 {