Implement phase 1 of expanded Floating
authorDan Doel <dan.doel@gmail.com>
Sun, 20 Dec 2015 14:19:52 +0000 (15:19 +0100)
committerBen Gamari <ben@smart-cactus.org>
Mon, 21 Dec 2015 12:01:36 +0000 (13:01 +0100)
- This part of the proposal is to add log1p, expm1, log1pexp and
  log1mexp to the Floating class, and export the full Floating class
  from Numeric

Reviewers: ekmett, #core_libraries_committee, bgamari, hvr, austin

Reviewed By: ekmett, #core_libraries_committee, bgamari

Subscribers: Phyx, RyanGlScott, ekmett, thomie

Differential Revision: https://phabricator.haskell.org/D1605

GHC Trac Issues: #11166

libraries/base/Data/Complex.hs
libraries/base/GHC/Float.hs
libraries/base/Numeric.hs
libraries/base/changelog.md
rts/RtsSymbols.c

index 31550d5..dd831bb 100644 (file)
@@ -37,6 +37,7 @@ module Data.Complex
         )  where
 
 import GHC.Generics (Generic, Generic1)
+import GHC.Float (Floating(..))
 import Data.Data (Data)
 import Foreign (Storable, castPtr, peek, poke, pokeElemOff, peekElemOff, sizeOf,
                 alignment)
@@ -195,6 +196,20 @@ instance  (RealFloat a) => Floating (Complex a) where
     acosh z        =  log (z + (z+1) * sqrt ((z-1)/(z+1)))
     atanh z        =  0.5 * log ((1.0+z) / (1.0-z))
 
+    log1p x@(a :+ b)
+      | abs a < 0.5 && abs b < 0.5
+      , u <- 2*a + a*a + b*b = log1p (u/(1 + sqrt(u+1))) :+ atan2 (1 + a) b
+      | otherwise = log (1 + x)
+    {-# INLINE log1p #-}
+
+    expm1 x@(a :+ b)
+      | a*a + b*b < 1
+      , u <- expm1 a
+      , v <- sin (b/2)
+      , w <- -2*v*v = (u*w + u + w) :+ (u+1)*sin b
+      | otherwise = exp x - 1
+    {-# INLINE expm1 #-}
+
 instance Storable a => Storable (Complex a) where
     sizeOf a       = 2 * sizeOf (realPart a)
     alignment a    = alignment (realPart a)
index e6180af..ddf9cf0 100644 (file)
@@ -4,6 +4,7 @@
            , MagicHash
            , UnboxedTuples
   #-}
+{-# LANGUAGE CApiFFI #-}
 -- We believe we could deorphan this module, by moving lots of things
 -- around, but we haven't got there yet:
 {-# OPTIONS_GHC -Wno-orphans #-}
@@ -61,6 +62,46 @@ class  (Fractional a) => Floating a  where
     sinh, cosh, tanh    :: a -> a
     asinh, acosh, atanh :: a -> a
 
+    -- | @'log1p' x@ computes @'log' (1 + x)@, but provides more precise
+    -- results for small (absolute) values of @x@ if possible.
+    --
+    -- @since 4.9.0.0
+    log1p               :: a -> a
+
+    -- | @'expm1' x@ computes @'exp' x - 1@, but provides more precise
+    -- results for small (absolute) values of @x@ if possible.
+    --
+    -- @since 4.9.0.0
+    expm1               :: a -> a
+
+    -- | @'log1pexp' x@ computes @'log' (1 + 'exp' x)@, but provides more
+    -- precise results if possible.
+    --
+    -- Examples:
+    --
+    -- * if @x@ is a large negative number, @'log' (1 + 'exp' x)@ will be
+    --   imprecise for the reasons given in 'log1p'.
+    --
+    -- * if @'exp' x@ is close to @-1@, @'log' (1 + 'exp' x)@ will be
+    --   imprecise for the reasons given in 'expm1'.
+    --
+    -- @since 4.9.0.0
+    log1pexp            :: a -> a
+
+    -- | @'log1mexp' x@ computes @'log' (1 - 'exp' x)@, but provides more
+    -- precise results if possible.
+    --
+    -- Examples:
+    --
+    -- * if @x@ is a large negative number, @'log' (1 - 'exp' x)@ will be
+    --   imprecise for the reasons given in 'log1p'.
+    --
+    -- * if @'exp' x@ is close to @1@, @'log' (1 - 'exp' x)@ will be
+    --   imprecise for the reasons given in 'expm1'.
+    --
+    -- @since 4.9.0.0
+    log1mexp            :: a -> a
+
     {-# INLINE (**) #-}
     {-# INLINE logBase #-}
     {-# INLINE sqrt #-}
@@ -72,6 +113,15 @@ class  (Fractional a) => Floating a  where
     tan  x              =  sin  x / cos  x
     tanh x              =  sinh x / cosh x
 
+    {-# INLINE log1p #-}
+    {-# INLINE expm1 #-}
+    {-# INLINE log1pexp #-}
+    {-# INLINE log1mexp #-}
+    log1p x = log (1 + x)
+    expm1 x = exp x - 1
+    log1pexp x = log1p (exp x)
+    log1mexp x = log1p (negate (exp x))
+
 -- | Efficient, machine-independent access to the components of a
 -- floating-point number.
 class  (RealFrac a, Floating a) => RealFloat a  where
@@ -307,6 +357,19 @@ instance  Floating Float  where
     acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
     atanh x = 0.5 * log ((1.0+x) / (1.0-x))
 
+    log1p = log1pFloat
+    expm1 = expm1Float
+
+    log1mexp a
+      | a <= log 2 = log (negate (expm1Float a))
+      | otherwise  = log1pFloat (negate (exp a))
+    {-# INLINE log1mexp #-}
+    log1pexp a
+      | a <= 18   = log1pFloat (exp a)
+      | a <= 100  = a + exp (negate a)
+      | otherwise = a
+    {-# INLINE log1pexp #-}
+
 instance  RealFloat Float  where
     floatRadix _        =  FLT_RADIX        -- from float.h
     floatDigits _       =  FLT_MANT_DIG     -- ditto
@@ -415,6 +478,19 @@ instance  Floating Double  where
     acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0)))
     atanh x = 0.5 * log ((1.0+x) / (1.0-x))
 
+    log1p = log1pDouble
+    expm1 = expm1Double
+
+    log1mexp a
+      | a <= log 2 = log (negate (expm1Double a))
+      | otherwise  = log1pDouble (negate (exp a))
+    {-# INLINE log1mexp #-}
+    log1pexp a
+      | a <= 18   = log1pDouble (exp a)
+      | a <= 100  = a + exp (negate a)
+      | otherwise = a
+    {-# INLINE log1pexp #-}
+
 -- RULES for Integer and Int
 {-# RULES
 "properFraction/Double->Integer"    properFraction = properFractionDoubleInteger
@@ -1069,6 +1145,16 @@ foreign import ccall unsafe "isDoubleDenormalized" isDoubleDenormalized :: Doubl
 foreign import ccall unsafe "isDoubleNegativeZero" isDoubleNegativeZero :: Double -> Int
 foreign import ccall unsafe "isDoubleFinite" isDoubleFinite :: Double -> Int
 
+
+------------------------------------------------------------------------
+-- libm imports for extended floating
+------------------------------------------------------------------------
+foreign import capi unsafe "math.h log1p" log1pDouble :: Double -> Double
+foreign import capi unsafe "math.h expm1" expm1Double :: Double -> Double
+foreign import capi unsafe "math.h log1pf" log1pFloat :: Float -> Float
+foreign import capi unsafe "math.h expm1f" expm1Float :: Float -> Float
+
+
 ------------------------------------------------------------------------
 -- Coercion rules
 ------------------------------------------------------------------------
index 4e24bfe..51be3a1 100644 (file)
@@ -56,6 +56,7 @@ module Numeric (
         -- * Miscellaneous
 
         fromRat,
+        Floating(..)
 
         ) where
 
index 33a5114..fa57556 100644 (file)
 
   * Re-export `Const` from `Control.Applicative` for backwards compatibility.
 
+  * Expand `Floating` class to include operations that allow for better
+    precision: `log1p`, `expm1`, `log1pexp` and `log1mexp`. These are not
+    available from `Prelude`, but the full class is exported from `Numeric`.
 
 ## 4.8.2.0  *Oct 2015*
 
index 8413d31..4b0a1d5 100644 (file)
       SymI_HasProto(erfc)                                \
       SymI_HasProto(erff)                                \
       SymI_HasProto(erfcf)                               \
+      SymI_HasProto(expm1)                               \
+      SymI_HasProto(expm1f)                              \
+      SymI_HasProto(log1p)                               \
+      SymI_HasProto(log1pf)                              \
       SymI_HasProto(memcpy)                              \
       SymI_HasProto(rts_InstallConsoleEvent)             \
       SymI_HasProto(rts_ConsoleHandlerDone)              \