diff --git a/Numeric/SpecFunctions.hs b/Numeric/SpecFunctions.hs --- a/Numeric/SpecFunctions.hs +++ b/Numeric/SpecFunctions.hs @@ -47,10 +47,6 @@ ) where import Numeric.SpecFunctions.Internal -#if MIN_VERSION_base(4,9,0) -import GHC.Float (log1p, expm1) -#endif - -- $log1p -- diff --git a/Numeric/SpecFunctions/Compat.hs b/Numeric/SpecFunctions/Compat.hs new file mode 100644 --- /dev/null +++ b/Numeric/SpecFunctions/Compat.hs @@ -0,0 +1,156 @@ +{-# LANGUAGE CPP #-} +{-# LANGUAGE ForeignFunctionInterface #-} +-- | +-- Functions which have different implementations on different platforms +module Numeric.SpecFunctions.Compat ( + erf + , erfc + , log1p + , expm1 + ) where + +import qualified Data.Vector.Unboxed as U +import Numeric.MathFunctions.Constants +import Numeric.Polynomial.Chebyshev (chebyshev,chebyshevBroucke) + +-- GHC.Float provides log1p and expm1 since base-4.9.0 (GHC8.0). GHCJS +-- doesn't +#define USE_GHC_LOG1P_EXP1M (MIN_VERSION_base(4,9,0) && !defined(__GHCJS__)) +#if USE_GHC_LOG1P_EXP1M +import GHC.Float (log1p,expm1) +#endif + + +---------------------------------------------------------------- +-- erf & erfc +-- +-- We provide pure haskell implementation for GHCJS and accesible on +-- GHC via flag +---------------------------------------------------------------- + +#if USE_SYSTEM_ERF && !defined(__GHCJS__) + +erf :: Double -> Double +erf = c_erf +{-# INLINE erf #-} + +erfc :: Double -> Double +erfc = c_erfc +{-# INLINE erfc #-} + +foreign import ccall unsafe "erf" c_erf :: Double -> Double +foreign import ccall unsafe "erfc" c_erfc :: Double -> Double + +#else + +erf :: Double -> Double +erf x | x < 0 = (-1) + erfcCheb (-x) + | otherwise = 1 - erfcCheb x + +erfc :: Double -> Double +erfc x | x < 0 = 2 - erfcCheb (-x) + | otherwise = erfcCheb x + +-- Adapted from Numerical Recipes §6.2.2 +erfcCheb :: Double -> Double +erfcCheb z + = t * exp( -z * z + chebyshev ty erfcCoef ) + where + -- We're using approximation: + -- + -- erfc(z) ≈ t·exp(-z² + P(t)) + -- t = 2 / (2 + z) + t = 2 / (2 + z) + ty = 2 * t - 1 + +erfcCoef :: U.Vector Double +{-# NOINLINE erfcCoef #-} +erfcCoef = U.fromList + [ -0.6513268598908546 , 6.4196979235649026e-1 , 1.9476473204185836e-2 + , -9.561514786808631e-3 , -9.46595344482036e-4 , 3.66839497852761e-4 + , 4.2523324806907e-5 , -2.0278578112534e-5 , -1.624290004647e-6 + , 1.303655835580e-6 , 1.5626441722e-8 , -8.5238095915e-8 + , 6.529054439e-9 , 5.059343495e-9 , -9.91364156e-10 + , -2.27365122e-10 , 9.6467911e-11 , 2.394038e-12 + , -6.886027e-12 , 8.94487e-13 , 3.13092e-13 + , -1.12708e-13 , 3.81e-16 , 7.106e-15 + , -1.523e-15 , -9.4e-17 , 1.21e-16 + , -2.8e-17 + ] + +#endif + + +---------------------------------------------------------------- +-- expm1 +-- +-- We use version provided by GHC is available otherwise we can either +-- get from libc or if everything else fails use one from library +---------------------------------------------------------------- + +#if !USE_GHC_LOG1P_EXP1M +-- | 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 unsafe "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 + + +---------------------------------------------------------------- +-- log1p +-- +-- Basically same as exm1 +---------------------------------------------------------------- + +#if !USE_GHC_LOG1P_EXP1M +-- | 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. +log1p :: Double -> Double +log1p x + | x == 0 = 0 + | x == -1 = m_neg_inf + | x < -1 = m_NaN + | x' < m_epsilon * 0.5 = x + | (x >= 0 && x < 1e-8) || (x >= -1e-9 && x < 0) + = x * (1 - x * 0.5) + | x' < 0.375 = x * (1 - x * chebyshevBroucke (x / 0.375) coeffs) + | otherwise = log (1 + x) + where + x' = abs x + coeffs = U.fromList [ + 0.10378693562743769800686267719098e+1, + -0.13364301504908918098766041553133e+0, + 0.19408249135520563357926199374750e-1, + -0.30107551127535777690376537776592e-2, + 0.48694614797154850090456366509137e-3, + -0.81054881893175356066809943008622e-4, + 0.13778847799559524782938251496059e-4, + -0.23802210894358970251369992914935e-5, + 0.41640416213865183476391859901989e-6, + -0.73595828378075994984266837031998e-7, + 0.13117611876241674949152294345011e-7, + -0.23546709317742425136696092330175e-8, + 0.42522773276034997775638052962567e-9, + -0.77190894134840796826108107493300e-10, + 0.14075746481359069909215356472191e-10, + -0.25769072058024680627537078627584e-11, + 0.47342406666294421849154395005938e-12, + -0.87249012674742641745301263292675e-13, + 0.16124614902740551465739833119115e-13, + -0.29875652015665773006710792416815e-14, + 0.55480701209082887983041321697279e-15, + -0.10324619158271569595141333961932e-15 + ] +#endif 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 CPP, BangPatterns, ScopedTypeVariables, ForeignFunctionInterface #-} +{-# LANGUAGE BangPatterns, ScopedTypeVariables #-} -- | -- Module : Numeric.SpecFunctions.Internal -- Copyright : (c) 2009, 2011, 2012 Bryan O'Sullivan @@ -9,11 +9,13 @@ -- Portability : portable -- -- Internal module with implementation of special functions. -module Numeric.SpecFunctions.Internal where +module Numeric.SpecFunctions.Internal + ( module Numeric.SpecFunctions.Internal + , Compat.log1p + , Compat.expm1 + ) where -#if !MIN_VERSION_base(4,9,0) import Control.Applicative -#endif import Data.Bits ((.&.), (.|.), shiftR) import Data.Int (Int64) import Data.Word (Word) @@ -21,17 +23,14 @@ 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, NewtonParam(..), Tolerance(..)) import Numeric.Series import Numeric.MathFunctions.Constants - - +import Numeric.SpecFunctions.Compat (log1p) +import qualified Numeric.SpecFunctions.Compat as Compat ---------------------------------------------------------------- -- Error function @@ -53,8 +52,8 @@ -- \end{aligned} -- \] erf :: Double -> Double +erf = Compat.erf {-# INLINE erf #-} -erf = c_erf -- | Complementary error function. -- @@ -72,13 +71,9 @@ -- \end{aligned} -- \] erfc :: Double -> Double +erfc = Compat.erfc {-# INLINE erfc #-} -erfc = c_erfc -foreign import ccall "erf" c_erf :: Double -> Double -foreign import ccall "erfc" c_erfc :: Double -> Double - - -- | Inverse of 'erf'. invErf :: Double -- ^ /p/ ∈ [-1,1] -> Double @@ -721,65 +716,6 @@ ---------------------------------------------------------------- -- 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. -log1p :: Double -> Double -log1p x - | x == 0 = 0 - | x == -1 = m_neg_inf - | x < -1 = m_NaN - | x' < m_epsilon * 0.5 = x - | (x >= 0 && x < 1e-8) || (x >= -1e-9 && x < 0) - = x * (1 - x * 0.5) - | x' < 0.375 = x * (1 - x * chebyshevBroucke (x / 0.375) coeffs) - | otherwise = log (1 + x) - where - x' = abs x - coeffs = U.fromList [ - 0.10378693562743769800686267719098e+1, - -0.13364301504908918098766041553133e+0, - 0.19408249135520563357926199374750e-1, - -0.30107551127535777690376537776592e-2, - 0.48694614797154850090456366509137e-3, - -0.81054881893175356066809943008622e-4, - 0.13778847799559524782938251496059e-4, - -0.23802210894358970251369992914935e-5, - 0.41640416213865183476391859901989e-6, - -0.73595828378075994984266837031998e-7, - 0.13117611876241674949152294345011e-7, - -0.23546709317742425136696092330175e-8, - 0.42522773276034997775638052962567e-9, - -0.77190894134840796826108107493300e-10, - 0.14075746481359069909215356472191e-10, - -0.25769072058024680627537078627584e-11, - 0.47342406666294421849154395005938e-12, - -0.87249012674742641745301263292675e-13, - 0.16124614902740551465739833119115e-13, - -0.29875652015665773006710792416815e-14, - 0.55480701209082887983041321697279e-15, - -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 diff --git a/Numeric/Sum.hs b/Numeric/Sum.hs --- a/Numeric/Sum.hs +++ b/Numeric/Sum.hs @@ -78,7 +78,7 @@ -- | Sum a collection of values. -- -- Example: - -- @foo = 'sum' 'kbn' [1,2,3]@ + -- @foo = 'Numeric.Sum.sum' 'kbn' [1,2,3]@ sum :: (F.Foldable f) => (s -> Double) -> f Double -> Double sum f = f . F.foldl' add zero {-# INLINE sum #-} @@ -255,7 +255,7 @@ -- where s' = 'add' s x -- @ -- --- In most instances, you can simply use the much more general 'sum' +-- In most instances, you can simply use the much more general 'Numeric.Sum.sum' -- function instead of writing a summation function by hand. -- -- @ @@ -263,7 +263,7 @@ -- import Prelude hiding (sum) -- -- -- betterSumList :: [Double] -> Double --- betterSumList xs = 'sum' 'kbn' xs +-- betterSumList xs = 'Numeric.Sum.sum' 'kbn' xs -- @ -- Note well the use of 'seq' in the example above to force the diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -1,3 +1,9 @@ +## Changes in 0.3.2.0 + + * GHCJS is now supported + + * Flag `system-expm1` is set to true by default. Only affects GHC<8.0 + ## Changes in 0.3.1.0 * Exported data types for iteration steps in root finding 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.3.1.0 +version: 0.3.2.0 cabal-version: >= 1.10 license: BSD2 license-file: LICENSE @@ -19,6 +19,17 @@ for real functions, polynomial summation and Chebyshev polynomials. +tested-with: + GHC ==7.4.2 + || ==7.6.3 + || ==7.8.4 + || ==7.10.3 + || ==8.0.2 + || ==8.2.2 + || ==8.4.4 + || ==8.6.5 + , GHCJS ==8.4 + extra-source-files: changelog.md README.markdown @@ -49,6 +60,8 @@ , vector-th-unbox >= 0.2.1.6 if flag(system-expm1) || !os(windows) cpp-options: -DUSE_SYSTEM_EXPM1 + if flag(system-erf) && !impl(ghcjs) + cpp-options: -DUSE_SYSTEM_ERF exposed-modules: Numeric.MathFunctions.Constants Numeric.MathFunctions.Comparison @@ -61,10 +74,22 @@ Numeric.Sum other-modules: Numeric.SpecFunctions.Internal + Numeric.SpecFunctions.Compat flag system-expm1 - description: Use expm1 provided by system. Only have effect on windows - default: False + description: Use expm1 provided by system. For GHC newer that + 8.0, GHCJS, and on Windows has no effect. GHC>=8.0 + provides expm1 so it's used. On GHCJS and on Windows + we don't have C implementation so bundled one is + used instead. + default: True + manual: True + +flag system-erf + description: Use erf and erfc provided by system. On GHCJS + version provided by library is used regardless of + flag for that lack of libc. + default: True manual: True test-suite tests diff --git a/tests/Tests/SpecFunctions.hs b/tests/Tests/SpecFunctions.hs --- a/tests/Tests/SpecFunctions.hs +++ b/tests/Tests/SpecFunctions.hs @@ -4,6 +4,7 @@ tests ) where +import Control.Monad import qualified Data.Vector as V import Data.Vector ((!)) @@ -16,7 +17,7 @@ import Tests.Helpers import Tests.SpecFunctions.Tables import Numeric.SpecFunctions -import Numeric.MathFunctions.Comparison (within,relativeError) +import Numeric.MathFunctions.Comparison (within,relativeError,ulpDistance) import Numeric.MathFunctions.Constants (m_epsilon,m_tiny) tests :: Test @@ -28,11 +29,29 @@ , testProperty "0 <= I[B] <= 1" $ incompleteBetaInRange , testProperty "invIncompleteGamma = gamma^-1" $ invIGammaIsInverse -- XXX FIXME DISABLED due to failures - -- , testProperty "invIncompleteBeta = B^-1" $ invIBetaIsInverse + , testProperty "invIncompleteBeta = B^-1" $ invIBetaIsInverse , testProperty "gamma - increases" $ \(abs -> s) (abs -> x) (abs -> y) -> s > 0 ==> monotonicallyIncreases (incompleteGamma s) x y , testProperty "invErfc = erfc^-1" $ invErfcIsInverse , testProperty "invErf = erf^-1" $ invErfIsInverse + -- Tests for erfc mostly are to test implementation bundled with + -- library. libc's one is accurate within 1 ulp + , testCase "erfc table" $ forM_ tableErfc $ \(x,exact) -> do + let val = erfc x + assertBool (unlines [ " x = " ++ show x + , " expected = " ++ show exact + , " got = " ++ show val + , " ulps diff = " ++ show (ulpDistance exact val) + ]) + (within 64 exact val) + , testCase "erf table" $ forM_ tableErf $ \(x,exact) -> do + let val = erf x + assertBool (unlines [ " x = " ++ show x + , " expected = " ++ show exact + , " got = " ++ show val + , " ulps diff = " ++ show (ulpDistance exact val) + ]) + (within 24 exact val) -- Unit tests , testAssertion "Factorial is expected to be precise at 1e-15 level" $ and [ eq 1e-15 (factorial (fromIntegral n :: Int)) @@ -224,3 +243,81 @@ -- Truncate double to [0,1] range01 :: Double -> Double range01 = abs . (snd :: (Integer, Double) -> Double) . properFraction + + +-- Table of values for erfc. +-- +-- Values are computed using python's mpmath up to 30 significant +-- digits +tableErfc :: [(Double,Double)] +tableErfc = + [ (0.000000, 1.0) + , (0.020000, 0.977435425308155055306039814223) + , (0.040000, 0.954888893854875246972188637445) + , (0.060000, 0.932378405606691560417070009221) + , (0.080000, 0.909921874158981837409467036376) + , (0.100000, 0.887537083981715101595287748986) + , (0.200000, 0.777297410789521533823546968791) + , (0.300000, 0.671373240540872583810382014682) + , (0.400000, 0.571607644953331523545890372692) + , (0.500000, 0.479500122186953462317253346108) + , (0.600000, 0.396143909152074094917693241426) + , (0.700000, 0.322198806162581557723141845649) + , (0.800000, 0.257899035292339487410212644387) + , (0.900000, 0.20309178757716786033533383966) + , (1.000000, 0.157299207050285130658779364917) + , (1.100000, 0.119794930425918270342740490744) + , (1.200000, 0.0896860217703646316340682061529) + , (1.300000, 0.0659920550593475541498146384224) + , (1.400000, 0.047714880237351203600376783391) + , (1.500000, 0.0338948535246892729330237383541) + , (1.600000, 0.0236516166553559844782198079153) + , (1.700000, 0.0162095414092254391586870541911) + , (1.800000, 0.0109094983642692838537604396016) + , (1.900000, 0.0072095707647425327627840328679) + , (2.000000, 0.00467773498104726583793074363275) + , (2.0009765625, 0.00465759175242884900812001805563) + , (2.100000, 0.00297946665633298428569058244218) + , (2.200000, 0.00186284629798188985855863885328) + , (2.300000, 0.00114317659735665247591992820336) + , (2.400000, 0.000688513896645078885549974809715) + , (2.500000, 0.000406952017444958939564215739975) + , (3.000000, 0.0000220904969985854413727761295823) + , (3.500000, 0.000000743098372341412745523683756096) + , (11.000000, 1.44086613794369468033980970286e-54) + , (23.000000, 4.44126594808805724407488442895e-232) + ] +tableErf :: [(Double,Double)] +tableErf = + [ (0.000000, 0.0) + , (0.020000, 0.0225645746918449446939601857765) + , (0.040000, 0.0451111061451247530278113625549) + , (0.060000, 0.0676215943933084395829299907792) + , (0.080000, 0.0900781258410181625905329636245) + , (0.100000, 0.112462916018284898404712251014) + , (0.200000, 0.222702589210478466176453031209) + , (0.300000, 0.328626759459127416189617985318) + , (0.400000, 0.428392355046668476454109627308) + , (0.500000, 0.520499877813046537682746653892) + , (0.600000, 0.603856090847925905082306758574) + , (0.700000, 0.677801193837418442276858154351) + , (0.800000, 0.742100964707660512589787355613) + , (0.900000, 0.79690821242283213966466616034) + , (1.000000, 0.842700792949714869341220635083) + , (1.100000, 0.880205069574081729657259509256) + , (1.200000, 0.910313978229635368365931793847) + , (1.300000, 0.934007944940652445850185361578) + , (1.400000, 0.952285119762648796399623216609) + , (1.500000, 0.966105146475310727066976261646) + , (1.600000, 0.976348383344644015521780192085) + , (1.700000, 0.983790458590774560841312945809) + , (1.800000, 0.989090501635730716146239560398) + , (1.900000, 0.992790429235257467237215967132) + , (2.000000, 0.995322265018952734162069256367) + , (2.100000, 0.997020533343667015714309417558) + , (2.200000, 0.998137153702018110141441361147) + , (2.300000, 0.998856823402643347524080071797) + , (2.400000, 0.99931148610335492111445002519) + , (2.500000, 0.99959304798255504106043578426) + , (3.000000, 0.99997790950300141455862722387) + ]