diff --git a/Numeric/SpecFunctions.hs b/Numeric/SpecFunctions.hs --- a/Numeric/SpecFunctions.hs +++ b/Numeric/SpecFunctions.hs @@ -1,3 +1,4 @@ +{-# LANGUAGE CPP #-} -- | -- Module : Numeric.SpecFunctions -- Copyright : (c) 2009, 2011, 2012 Bryan O'Sullivan @@ -28,6 +29,7 @@ -- * Sinc , sinc -- * Logarithm + -- $log1p , log1p , log1pmx , log2 @@ -45,6 +47,17 @@ ) where import Numeric.SpecFunctions.Internal +#if MIN_VERSION_base(4,9,0) +import GHC.Float (log1p, expm1) +#endif + + +-- $log1p +-- +-- Base starting from @4.9.0@ (GHC 8.0) provides 'log1p' and 'expm1' +-- as method of class 'Floating'. In this case we simply reexport +-- these function. Otherwise we provide our own with more restrictive +-- signature @Double → Double@. -- $references -- diff --git a/Numeric/SpecFunctions/Internal.hs b/Numeric/SpecFunctions/Internal.hs --- a/Numeric/SpecFunctions/Internal.hs +++ b/Numeric/SpecFunctions/Internal.hs @@ -1,4 +1,4 @@ -{-# LANGUAGE BangPatterns, ScopedTypeVariables, ForeignFunctionInterface #-} +{-# LANGUAGE CPP, BangPatterns, ScopedTypeVariables, ForeignFunctionInterface #-} -- | -- Module : Numeric.SpecFunctions.Internal -- Copyright : (c) 2009, 2011, 2012 Bryan O'Sullivan @@ -11,24 +11,27 @@ -- Internal module with implementation of special functions. module Numeric.SpecFunctions.Internal where +#if !MIN_VERSION_base(4,9,0) import Control.Applicative +#endif import Data.Bits ((.&.), (.|.), shiftR) import Data.Int (Int64) 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) +import GHC.Float (log1p,expm1) +#endif import Numeric.Polynomial.Chebyshev (chebyshevBroucke) import Numeric.Polynomial (evaluatePolynomialL,evaluateEvenPolynomialL,evaluateOddPolynomialL) import Numeric.RootFinding (Root(..), newtonRaphson) -import Numeric.Series (sumPowerSeries,enumSequenceFrom,scanSequence,evalContFractionB) -import Numeric.MathFunctions.Constants ( m_epsilon, m_NaN, m_neg_inf, m_pos_inf - , m_sqrt_2_pi, m_ln_sqrt_2_pi, m_eulerMascheroni - , m_min_log, m_tiny - ) -import Text.Printf +import Numeric.Series +import Numeric.MathFunctions.Constants + ---------------------------------------------------------------- -- Error function ---------------------------------------------------------------- @@ -718,6 +721,8 @@ -- Logarithm ---------------------------------------------------------------- +-- GHC.Float provides log1p and expm1 since 4.9.0 +#if !MIN_VERSION_base(4,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. @@ -758,6 +763,23 @@ -0.10324619158271569595141333961932e-15 ] +-- | Compute @exp x - 1@ without loss of accuracy for x near zero. +expm1 :: Double -> Double +#ifdef USE_SYSTEM_EXPM1 +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)) +#endif +#endif + -- | Compute log(1+x)-x: log1pmx :: Double -> Double log1pmx x @@ -768,14 +790,6 @@ | otherwise = - x * x * sumPowerSeries (-x) (recip <$> enumSequenceFrom 2) where ax = abs x - --- | Compute @exp x - 1@ without loss of accuracy for x near zero. -expm1 :: Double -> Double -expm1 = c_expm1 - -foreign import ccall "expm1" c_expm1 :: Double -> Double - - -- | /O(log n)/ Compute the logarithm in base 2 of the given value. log2 :: Int -> Int diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -1,9 +1,17 @@ -Changes in 0.2.0.0 +## Changes in 0.2.1.0 + * `log1p` and `expm1` are simply reexported from `GHC.Float`. They're methods + of `Floating` type class. + + * On windows `expm1` is implemented in pure haskell for older GHCs. + + +## Changes in 0.2.0.0 + * Bug fixes and documentation tweaks -Changes in 0.2.0.0 +## Changes in 0.2.0.0 * `logGamma` now uses Lancsoz approximation and same as `logGammaL`. Old implementation of `logGamma` moved to `Numeric.SpecFunctions.Extra.logGammaAS245`. @@ -36,7 +44,7 @@ -Changes in 0.1.7.0 +## Changes in 0.1.7.0 * Module `statistics: Statistics.Function.Comparison` moved to `Numeric.MathFunctions.Comparison`. Old implementation if `within` compared @@ -49,7 +57,7 @@ * Precision of `logFactorial` is slightly improved. -Changes in 0.1.6.0 +## Changes in 0.1.6.0 * `logChoose` added for calculation of logarithm of binomial coefficient @@ -58,17 +66,17 @@ * `sinc` added -Changes in 0.1.5.3 +## Changes in 0.1.5.3 * Fix for test suite on 32bit platform -Changes in 0.1.5 +## Changes in 0.1.5 * Numeric.Sum: new module adds accurate floating point summation. -Changes in 0.1.4 +## Changes in 0.1.4 * logFactorial type is genberalized. It accepts any `Integral` type @@ -76,7 +84,7 @@ are store in lists added -Changes in 0.1.3 +## Changes in 0.1.3 * Error function and its inverse added. diff --git a/math-functions.cabal b/math-functions.cabal --- a/math-functions.cabal +++ b/math-functions.cabal @@ -1,5 +1,5 @@ name: math-functions -version: 0.2.0.2 +version: 0.2.1.0 cabal-version: >= 1.10 license: BSD3 license-file: LICENSE @@ -39,11 +39,13 @@ 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 @@ -56,6 +58,11 @@ Numeric.Sum 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