base: fix sign confusion in log1mexp implementation (fix #17125)

author: claude (https://gitlab.haskell.org/trac-claude)

The correct threshold for log1mexp is -(log 2) with the current specification
of log1mexp. This change improves accuracy for large negative inputs.

To avoid code duplication, a small helper function is added;
it isn't the default implementation in Floating because it needs Ord.

This patch does nothing to address that the Haskell specification is
different from that in common use in other languages.

1 files changed, 10 insertions, 6 deletions

diff --git a/libraries/base/GHC/Float.hs b/libraries/base/GHC/Float.hs
index 03e232c190..2f21daa57f 100644
--- a/libraries/base/GHC/Float.hs
+++ b/libraries/base/GHC/Float.hs
@@ -142,6 +142,14 @@ class  (Fractional a) => Floating a  where
     log1pexp x = log1p (exp x)
     log1mexp x = log1p (negate (exp x))
 
+-- | Default implementation for @'log1mexp'@ requiring @'Ord'@ to test
+-- against a threshold to decide which implementation variant to use.
+log1mexpOrd :: (Ord a, Floating a) => a -> a
+{-# INLINE log1mexpOrd #-}
+log1mexpOrd a
+    | a > -(log 2) = log (negate (expm1 a))
+    | otherwise  = log1p (negate (exp a))
+
 -- | Efficient, machine-independent access to the components of a
 -- floating-point number.
 class  (RealFrac a, Floating a) => RealFloat a  where
@@ -399,9 +407,7 @@ instance  Floating Float  where
     log1p = log1pFloat
     expm1 = expm1Float
 
-    log1mexp a
-      | a <= log 2 = log (negate (expm1Float a))
-      | otherwise  = log1pFloat (negate (exp a))
+    log1mexp x = log1mexpOrd x
     {-# INLINE log1mexp #-}
     log1pexp a
       | a <= 18   = log1pFloat (exp a)
@@ -540,9 +546,7 @@ instance  Floating Double  where
     log1p = log1pDouble
     expm1 = expm1Double
 
-    log1mexp a
-      | a <= log 2 = log (negate (expm1Double a))
-      | otherwise  = log1pDouble (negate (exp a))
+    log1mexp x = log1mexpOrd x
     {-# INLINE log1mexp #-}
     log1pexp a
       | a <= 18   = log1pDouble (exp a)