Commit 6457903e authored by dolio's avatar dolio Committed by Ben Gamari

Implement phase 1 of expanded Floating

- 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
parent edcf17bd
...@@ -37,6 +37,7 @@ module Data.Complex ...@@ -37,6 +37,7 @@ module Data.Complex
) where ) where
import GHC.Generics (Generic, Generic1) import GHC.Generics (Generic, Generic1)
import GHC.Float (Floating(..))
import Data.Data (Data) import Data.Data (Data)
import Foreign (Storable, castPtr, peek, poke, pokeElemOff, peekElemOff, sizeOf, import Foreign (Storable, castPtr, peek, poke, pokeElemOff, peekElemOff, sizeOf,
alignment) alignment)
...@@ -195,6 +196,20 @@ instance (RealFloat a) => Floating (Complex a) where ...@@ -195,6 +196,20 @@ instance (RealFloat a) => Floating (Complex a) where
acosh z = log (z + (z+1) * sqrt ((z-1)/(z+1))) acosh z = log (z + (z+1) * sqrt ((z-1)/(z+1)))
atanh z = 0.5 * log ((1.0+z) / (1.0-z)) 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 instance Storable a => Storable (Complex a) where
sizeOf a = 2 * sizeOf (realPart a) sizeOf a = 2 * sizeOf (realPart a)
alignment a = alignment (realPart a) alignment a = alignment (realPart a)
......
...@@ -4,6 +4,7 @@ ...@@ -4,6 +4,7 @@
, MagicHash , MagicHash
, UnboxedTuples , UnboxedTuples
#-} #-}
{-# LANGUAGE CApiFFI #-}
-- We believe we could deorphan this module, by moving lots of things -- We believe we could deorphan this module, by moving lots of things
-- around, but we haven't got there yet: -- around, but we haven't got there yet:
{-# OPTIONS_GHC -Wno-orphans #-} {-# OPTIONS_GHC -Wno-orphans #-}
...@@ -61,6 +62,46 @@ class (Fractional a) => Floating a where ...@@ -61,6 +62,46 @@ class (Fractional a) => Floating a where
sinh, cosh, tanh :: a -> a sinh, cosh, tanh :: a -> a
asinh, acosh, atanh :: 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 (**) #-}
{-# INLINE logBase #-} {-# INLINE logBase #-}
{-# INLINE sqrt #-} {-# INLINE sqrt #-}
...@@ -72,6 +113,15 @@ class (Fractional a) => Floating a where ...@@ -72,6 +113,15 @@ class (Fractional a) => Floating a where
tan x = sin x / cos x tan x = sin x / cos x
tanh x = sinh x / cosh 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 -- | Efficient, machine-independent access to the components of a
-- floating-point number. -- floating-point number.
class (RealFrac a, Floating a) => RealFloat a where class (RealFrac a, Floating a) => RealFloat a where
...@@ -307,6 +357,19 @@ instance Floating Float where ...@@ -307,6 +357,19 @@ instance Floating Float where
acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0))) 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)) 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 instance RealFloat Float where
floatRadix _ = FLT_RADIX -- from float.h floatRadix _ = FLT_RADIX -- from float.h
floatDigits _ = FLT_MANT_DIG -- ditto floatDigits _ = FLT_MANT_DIG -- ditto
...@@ -415,6 +478,19 @@ instance Floating Double where ...@@ -415,6 +478,19 @@ instance Floating Double where
acosh x = log (x + (x+1.0) * sqrt ((x-1.0)/(x+1.0))) 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)) 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 for Integer and Int
{-# RULES {-# RULES
"properFraction/Double->Integer" properFraction = properFractionDoubleInteger "properFraction/Double->Integer" properFraction = properFractionDoubleInteger
...@@ -1069,6 +1145,16 @@ foreign import ccall unsafe "isDoubleDenormalized" isDoubleDenormalized :: Doubl ...@@ -1069,6 +1145,16 @@ foreign import ccall unsafe "isDoubleDenormalized" isDoubleDenormalized :: Doubl
foreign import ccall unsafe "isDoubleNegativeZero" isDoubleNegativeZero :: Double -> Int foreign import ccall unsafe "isDoubleNegativeZero" isDoubleNegativeZero :: Double -> Int
foreign import ccall unsafe "isDoubleFinite" isDoubleFinite :: 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 -- Coercion rules
------------------------------------------------------------------------ ------------------------------------------------------------------------
......
...@@ -56,6 +56,7 @@ module Numeric ( ...@@ -56,6 +56,7 @@ module Numeric (
-- * Miscellaneous -- * Miscellaneous
fromRat, fromRat,
Floating(..)
) where ) where
......
...@@ -112,6 +112,9 @@ ...@@ -112,6 +112,9 @@
* Re-export `Const` from `Control.Applicative` for backwards compatibility. * 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* ## 4.8.2.0 *Oct 2015*
......
...@@ -157,6 +157,10 @@ ...@@ -157,6 +157,10 @@
SymI_HasProto(erfc) \ SymI_HasProto(erfc) \
SymI_HasProto(erff) \ SymI_HasProto(erff) \
SymI_HasProto(erfcf) \ SymI_HasProto(erfcf) \
SymI_HasProto(expm1) \
SymI_HasProto(expm1f) \
SymI_HasProto(log1p) \
SymI_HasProto(log1pf) \
SymI_HasProto(memcpy) \ SymI_HasProto(memcpy) \
SymI_HasProto(rts_InstallConsoleEvent) \ SymI_HasProto(rts_InstallConsoleEvent) \
SymI_HasProto(rts_ConsoleHandlerDone) \ SymI_HasProto(rts_ConsoleHandlerDone) \
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment