Skip to content

Commit

Permalink
Use expm1 defined locally only on windows
Browse files Browse the repository at this point in the history
  • Loading branch information
Shimuuar committed Dec 1, 2016
1 parent 004cf1e commit 32d66a5
Show file tree
Hide file tree
Showing 3 changed files with 23 additions and 12 deletions.
3 changes: 2 additions & 1 deletion Numeric/SpecFunctions.hs
Expand Up @@ -51,10 +51,11 @@ module Numeric.SpecFunctions (
) where

import Numeric.SpecFunctions.Internal
#if MIN_VERSION_base(4,9,0)
#if MIN_VERSION_base(44,9,0)
import GHC.Float (log1p, expm1)
#endif


-- $references
--
-- * Bernardo, J. (1976) Algorithm AS 103: Psi (digamma)
Expand Down
15 changes: 9 additions & 6 deletions Numeric/SpecFunctions/Internal.hs
Expand Up @@ -11,7 +11,7 @@
-- Internal module with implementation of special functions.
module Numeric.SpecFunctions.Internal where

#if !MIN_VERSION_base(4,9,0)
#if !MIN_VERSION_base(44,9,0)
import Control.Applicative
#endif
import Data.Bits ((.&.), (.|.), shiftR)
Expand All @@ -20,7 +20,7 @@ import Data.Word (Word)
import qualified Data.Vector.Unboxed as U
import Data.Vector.Unboxed ((!))
import Text.Printf
#if MIN_VERSION_base(4,9,0)
#if MIN_VERSION_base(44,9,0)
import GHC.Float (log1p,expm1)
#endif

Expand Down Expand Up @@ -722,7 +722,7 @@ sinc x
----------------------------------------------------------------

-- GHC.Float provides log1p and expm1 since 4.9.0
#if !MIN_VERSION_base(4,9,0)
#if !MIN_VERSION_base(44,9,0)
-- | Compute the natural logarithm of 1 + @x@. This is accurate even
-- for values of @x@ near zero, where use of @log(1+x)@ would lose
-- precision.
Expand Down Expand Up @@ -764,16 +764,19 @@ log1p x
]
-- | Compute @exp x - 1@ without loss of accuracy for x near zero.
expm1 :: Double -> Double
#ifdef USE_SYSTEM_GHC
expm1 = c_expm1

foreign import ccall "expm1" c_expm1 :: Double -> Double
#else
-- NOTE: this is simplest implementation and not terribly efficient.
expm1 x
| x < (-37.42994775023705) = -1
| x > m_max_log = m_pos_inf
| abs x > 0.5 = exp x - 1
| otherwise = sumSeries $ liftA2 (*) (scanSequence (*) x (pure x))
(1 / scanSequence (*) 1 (enumSequenceFrom 2))

-- foreign import ccall "expm1" c_expm1 :: Double -> Double
-- MIN_VERSION_base(4,9,0)
#endif
#endif

-- | Compute log(1+x)-x:
Expand Down
17 changes: 12 additions & 5 deletions math-functions.cabal
Expand Up @@ -39,11 +39,13 @@ library
DeriveGeneric

ghc-options: -Wall -O2
build-depends: base >=4.5 && <5,
deepseq,
vector >= 0.7,
primitive,
vector-th-unbox
build-depends: base >=4.5 && <5
, deepseq
, vector >= 0.7
, primitive
, vector-th-unbox
if flag(system-expm1) || !os(windows)
cpp-options: -DUSE_SYSTEM_EXPM1
exposed-modules:
Numeric.MathFunctions.Constants
Numeric.MathFunctions.Comparison
Expand All @@ -57,6 +59,11 @@ library
other-modules:
Numeric.SpecFunctions.Internal

flag system-expm1
description: Use expm1 provided by system. Only have effect on windows
default: False
manual: True

test-suite tests
default-language: Haskell2010
other-extensions: ViewPatterns
Expand Down

0 comments on commit 32d66a5

Please sign in to comment.