math-functions 0.3.1.0 → 0.3.2.0
raw patch · 7 files changed
+302/−86 lines, 7 filesPVP ok
version bump matches the API change (PVP)
API changes (from Hackage documentation)
Files
- Numeric/SpecFunctions.hs +0/−4
- Numeric/SpecFunctions/Compat.hs +156/−0
- Numeric/SpecFunctions/Internal.hs +10/−74
- Numeric/Sum.hs +3/−3
- changelog.md +6/−0
- math-functions.cabal +28/−3
- tests/Tests/SpecFunctions.hs +99/−2
Numeric/SpecFunctions.hs view
@@ -47,10 +47,6 @@ ) where import Numeric.SpecFunctions.Internal-#if MIN_VERSION_base(4,9,0)-import GHC.Float (log1p, expm1)-#endif- -- $log1p --
+ Numeric/SpecFunctions/Compat.hs view
@@ -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
Numeric/SpecFunctions/Internal.hs view
@@ -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
Numeric/Sum.hs view
@@ -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
changelog.md view
@@ -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
math-functions.cabal view
@@ -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
tests/Tests/SpecFunctions.hs view
@@ -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)+ ]