diff options
author | Artem Pelenitsyn <a.pelenitsyn@gmail.com> | 2020-06-04 05:09:48 +0200 |
---|---|---|
committer | Marge Bot <ben+marge-bot@smart-cactus.org> | 2020-06-05 03:18:49 -0400 |
commit | af5e3a885ddd09dd5f550552c535af3661ff3dbf (patch) | |
tree | d7fba407eebfa84fd13c8d1f88868f1cf442fa1b | |
parent | 6a4098a4bb89b3d30cca26d82b82724913062536 (diff) | |
download | haskell-af5e3a885ddd09dd5f550552c535af3661ff3dbf.tar.gz |
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.
-rw-r--r-- | libraries/base/GHC/Float.hs | 16 |
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) |