diff options
| author | Artem Pelenitsyn <a.pelenitsyn@gmail.com> | 2020-06-04 05:09:48 +0200 |
|---|---|---|
| committer | Ben Gamari <ben@smart-cactus.org> | 2020-07-08 12:54:35 -0400 |
| commit | 14401143ce3df29cd8fe82b23b3e8cba33b7cb78 (patch) | |
| tree | 8c3f78014c4cd967e3e8d5b82c79413d7be27492 | |
| parent | 407b98f710a16a677abe1a0c2e5160ff2da31d57 (diff) | |
| download | haskell-14401143ce3df29cd8fe82b23b3e8cba33b7cb78.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.
(cherry picked from commit af5e3a885ddd09dd5f550552c535af3661ff3dbf)
| -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 75f6f8bd68..7c6ffad09a 100644 --- a/libraries/base/GHC/Float.hs +++ b/libraries/base/GHC/Float.hs @@ -141,6 +141,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 @@ -398,9 +406,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) @@ -539,9 +545,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) |