math-functions 0.3.4.2 → 0.3.4.3
raw patch · 13 files changed
+117/−138 lines, 13 filesdep +tasty-benchdep −gaugedep ~basePVP ok
version bump matches the API change (PVP)
Dependencies added: tasty-bench
Dependencies removed: gauge
Dependency ranges changed: base
API changes (from Hackage documentation)
Files
- Numeric/Polynomial.hs +1/−1
- Numeric/Polynomial/Chebyshev.hs +1/−1
- Numeric/RootFinding.hs +7/−34
- Numeric/SpecFunctions.hs +0/−1
- Numeric/SpecFunctions/Compat.hs +20/−30
- Numeric/SpecFunctions/Internal.hs +6/−7
- Numeric/Sum.hs +3/−12
- bench/bench.hs +2/−1
- changelog.md +7/−0
- math-functions.cabal +17/−21
- tests/Tests/SpecFunctions.hs +11/−2
- tests/Tests/Sum.hs +41/−27
- tests/tables/generate.py +1/−1
Numeric/Polynomial.hs view
@@ -69,7 +69,7 @@ -- $lists -- -- When all coefficients are known statically it's more convenient to--- pass coefficient in a list instad of vector. Functions below+-- pass coefficient in a list instead of vector. Functions below -- provide just that functionality. If list is known statically it -- will be inlined anyway.
Numeric/Polynomial/Chebyshev.hs view
@@ -71,7 +71,7 @@ -- -- * Broucke, R. (1973) Algorithm 446: Ten subroutines for the -- manipulation of Chebyshev series. /Communications of the ACM/--- 16(4):254–256. <http://doi.acm.org/10.1145/362003.362037>+-- 16(4):254–256. <http://doi.acm.org/10.1145/362003.362037> -- -- * Clenshaw, C.W. (1962) Chebyshev series for mathematical -- functions. /National Physical Laboratory Mathematical Tables 5/,
Numeric/RootFinding.hs view
@@ -1,5 +1,4 @@ {-# LANGUAGE BangPatterns #-}-{-# LANGUAGE CPP #-} {-# LANGUAGE DeriveDataTypeable #-} {-# LANGUAGE DeriveFoldable #-} {-# LANGUAGE DeriveGeneric #-}@@ -40,17 +39,12 @@ -- $references ) where -import Control.Applicative (Alternative(..), Applicative(..))+import Control.Applicative (Alternative(..)) import Control.Monad (MonadPlus(..), ap) import Control.DeepSeq (NFData(..)) import Data.Data (Data, Typeable)-import Data.Monoid (Monoid(..))-import Data.Foldable (Foldable)-import Data.Traversable (Traversable) import Data.Default.Class-#if __GLASGOW_HASKELL__ > 704 import GHC.Generics (Generic)-#endif import Numeric.MathFunctions.Comparison (within,eqRelErr) import Numeric.MathFunctions.Constants (m_epsilon) @@ -70,9 +64,7 @@ | Root !a -- ^ A root was successfully found. deriving (Eq, Read, Show, Typeable, Data, Foldable, Traversable-#if __GLASGOW_HASKELL__ > 704 , Generic-#endif ) instance (NFData a) => NFData (Root a) where@@ -128,11 +120,7 @@ -- considered approximately equal if \[ |a - b| < \delta \]. -- Note that @AbsTol 0@ could be used to require to find -- approximation within machine precision.- deriving (Eq, Read, Show, Typeable, Data-#if __GLASGOW_HASKELL__ > 704- , Generic-#endif- )+ deriving (Eq, Read, Show, Typeable, Data, Generic) -- | Check that two values are approximately equal. In addition to -- specification values are considered equal if they're within 1ulp@@ -179,11 +167,7 @@ -- ^ Error tolerance for root approximation. Default is relative -- error 4·ε, where ε is machine precision. }- deriving (Eq, Read, Show, Typeable, Data-#if __GLASGOW_HASKELL__ > 704- , Generic-#endif- )+ deriving (Eq, Read, Show, Typeable, Data, Generic) instance Default RiddersParam where def = RiddersParam@@ -202,11 +186,7 @@ -- ^ Root found | RiddersNoBracket -- ^ Root is not bracketed- deriving (Eq, Read, Show, Typeable, Data-#if __GLASGOW_HASKELL__ > 704- , Generic-#endif- )+ deriving (Eq, Read, Show, Typeable, Data, Generic) instance NFData RiddersStep where rnf x = x `seq` ()@@ -295,11 +275,7 @@ -- ^ Error tolerance for root approximation. Default is relative -- error 4·ε, where ε is machine precision }- deriving (Eq, Read, Show, Typeable, Data-#if __GLASGOW_HASKELL__ > 704- , Generic-#endif- )+ deriving (Eq, Read, Show, Typeable, Data, Generic) instance Default NewtonParam where def = NewtonParam@@ -318,11 +294,8 @@ -- ^ Root is found | NewtonNoBracket -- ^ Root is not bracketed- deriving (Eq, Read, Show, Typeable, Data-#if __GLASGOW_HASKELL__ > 704- , Generic-#endif- )+ deriving (Eq, Read, Show, Typeable, Data, Generic)+ instance NFData NewtonStep where rnf x = x `seq` ()
Numeric/SpecFunctions.hs view
@@ -1,4 +1,3 @@-{-# LANGUAGE CPP #-} -- | -- Module : Numeric.SpecFunctions -- Copyright : (c) 2009, 2011, 2012 Bryan O'Sullivan
Numeric/SpecFunctions/Compat.hs view
@@ -9,17 +9,23 @@ , expm1 ) where -import Control.Applicative+#if !defined(USE_SYSTEM_ERF) || !defined(USE_SYSTEM_EXPM1) import qualified Data.Vector.Unboxed as U-import Numeric.MathFunctions.Constants-import Numeric.Polynomial.Chebyshev (chebyshev,chebyshevBroucke)+#endif++#if !defined(USE_SYSTEM_ERF)+import Numeric.Polynomial.Chebyshev (chebyshev) import Numeric.Polynomial (evaluateOddPolynomial)-import Numeric.Series+#endif --- 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+#if !defined(USE_SYSTEM_EXPM1)+import Control.Applicative (liftA2)+import Numeric.Polynomial.Chebyshev (chebyshevBroucke)+import Numeric.Series (scanSequence,sumSeries,enumSequenceFrom)+import Numeric.MathFunctions.Constants+#endif++#if defined(USE_SYSTEM_EXPM1) import GHC.Float (log1p,expm1) #endif @@ -27,11 +33,11 @@ ---------------------------------------------------------------- -- erf & erfc ----- We provide pure haskell implementation for GHCJS and accesible on+-- We provide pure haskell implementation for GHCJS and accessible on -- GHC via flag ---------------------------------------------------------------- -#if USE_SYSTEM_ERF && !defined(__GHCJS__)+#if defined(USE_SYSTEM_ERF) erf :: Double -> Double erf = c_erf@@ -103,20 +109,15 @@ ------------------------------------------------------------------- expm1+-- expm1 & log1p ----- We use version provided by GHC is available otherwise we can either--- get from libc or if everything else fails use one from library+-- We use one provided by base of for GHCJS use hand-coded one ---------------------------------------------------------------- -#if !USE_GHC_LOG1P_EXP1M+#if !defined(USE_SYSTEM_EXPM1)+ -- | Compute @exp x - 1@ without loss of accuracy for x near zero. expm1 :: Double -> Double-#if USE_SYSTEM_EXPM1 && !defined(__GHCJS__)-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@@ -124,17 +125,6 @@ | 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.
Numeric/SpecFunctions/Internal.hs view
@@ -16,7 +16,6 @@ , Compat.expm1 ) where -import Control.Applicative import Data.Bits ((.&.), (.|.), shiftR) import Data.Int (Int64) import Data.Word (Word)@@ -151,7 +150,7 @@ -- > (z + 1) - 1 = z -- > (z + 1) - 2 = z - 1 --- -- Simple passing (z + 1) to piecewise approxiations and computing+ -- Simple passing (z + 1) to piecewise approximations and computing -- difference leads to bad loss of precision near 1. -- This is reason lgamma1_15 & lgamma15_2 have three parameters | z < 0.5 = lgamma1_15 z (z - 1) - log z@@ -283,7 +282,7 @@ -- For small z we can just use Gamma function recurrence and reduce -- problem to interval [2,3] and use polynomial approximation--- there. Surpringly it gives very good precision+-- there. Surprisingly it gives very good precision lgammaSmall :: Double -> Double lgammaSmall = go 0 where@@ -297,7 +296,7 @@ -- -- > Γ(z) = sqrt(2π)(z + g - 0.5)^(z - 0.5)·exp{-(z + g - 0.5)}·A_g(z) ----- Coeffients are taken from boost. Constants are absorbed into+-- Coefficients are taken from boost. Constants are absorbed into -- polynomial's coefficients. lanczosApprox :: Double -> Double lanczosApprox z@@ -516,7 +515,7 @@ -- New approximation to x x' | x < dx = 0.5 * x -- Do not go below 0 | otherwise = x - dx- -- Calculate inital guess for root+ -- Calculate initial guess for root guess -- | a > 1 =@@ -619,7 +618,7 @@ | otherwise = 1 - incompleteBetaWorker beta q p (1 - x) --- Approximation of incomplete beta by quandrature.+-- Approximation of incomplete beta by quadrature. -- -- Note that x =< p/(p+q) incompleteBetaApprox :: Double -> Double -> Double -> Double -> Double@@ -750,7 +749,7 @@ -- It's really hodgepodge of different approximations accumulated over years. -- -- Equations are referred to by name of paper and number e.g. [AS64 2]--- In AS64 papers equations are not numbered so they are refered to by+-- In AS64 papers equations are not numbered so they are referred to by -- number of appearance starting from definition of incomplete beta. invIncBetaGuess beta a b p -- If both a and b are less than 1 incomplete beta have inflection
Numeric/Sum.hs view
@@ -1,5 +1,5 @@ {-# LANGUAGE BangPatterns, DeriveDataTypeable, FlexibleContexts,- MultiParamTypeClasses, TypeFamilies, CPP #-}+ MultiParamTypeClasses, TypeFamilies #-} {-# OPTIONS_GHC -fno-warn-name-shadowing #-} -- | -- Module : Numeric.Sum@@ -53,11 +53,8 @@ import Control.DeepSeq (NFData(..)) import Data.Bits (shiftR) import Data.Data (Typeable, Data)-import Data.Monoid (Monoid(..))-#if MIN_VERSION_base(4,9,0) import Data.Semigroup (Semigroup(..))-#endif-import Data.Vector.Generic (Vector(..), foldl')+import Data.Vector.Generic (Vector(..)) -- Needed for GHC 7.2 & 7.4 to derive Unbox instances import Control.Monad (liftM) import Data.Vector.Generic.Mutable (MVector(..))@@ -159,11 +156,9 @@ mempty = zero s `mappend` KahanSum s' _ = add s s' -#if MIN_VERSION_base(4,9,0) -- | @since 0.3.0.0 instance Semigroup KahanSum where (<>) = mappend-#endif kahanAdd :: KahanSum -> Double -> KahanSum kahanAdd (KahanSum sum c) x = KahanSum sum' c'@@ -241,11 +236,9 @@ mempty = zero s `mappend` KBNSum s' c' = add (add s s') c' -#if MIN_VERSION_base(4,9,0) -- | @since 0.3.0.0 instance Semigroup KBNSum where (<>) = mappend-#endif kbnAdd :: KBNSum -> Double -> KBNSum kbnAdd (KBNSum sum c) x = KBNSum sum' c'@@ -329,11 +322,9 @@ mempty = zero s `mappend` KB2Sum s' c' cc' = add (add (add s s') c') cc' -#if MIN_VERSION_base(4,9,0) -- | @since 0.3.0.0 instance Semigroup KB2Sum where (<>) = mappend-#endif kb2Add :: KB2Sum -> Double -> KB2Sum@@ -352,7 +343,7 @@ -- | /O(n)/ Sum a vector of values. sumVector :: (Vector v Double, Summation s) => (s -> Double) -> v Double -> Double-sumVector f = f . foldl' add zero+sumVector f = f . G.foldl' add zero {-# INLINE sumVector #-} -- | /O(n)/ Sum a vector of values using pairwise summation.
bench/bench.hs view
@@ -1,8 +1,8 @@ {-# LANGUAGE NumDecimals #-}-import Gauge.Main import Data.Default.Class import qualified Data.Vector.Unboxed as U import Text.Printf+import Test.Tasty.Bench import System.Random (randomIO) import qualified Numeric.Sum as Sum@@ -13,6 +13,7 @@ -- Uniformly sample logGamma performance between 10^-6 to 10^6+benchmarkLogGamma :: (Double -> Double) -> [Benchmark] benchmarkLogGamma logG = [ bench (printf "%.3g" x) $ nf logG x | x <- [ m * 10**n | n <- [ -8 .. 8 ]
changelog.md view
@@ -1,3 +1,10 @@+## Changes in 0.3.4.3+ + * Support for `QuickCheck >= 2.14`. Test no longer fail++ * Support for GHC<8.0 dropped+ + ## Changes in 0.3.4.2 * Fixed crash in `invIncompleteBeta` (#68) for some inputs initial approximation
math-functions.cabal view
@@ -1,13 +1,13 @@ name: math-functions-version: 0.3.4.2+version: 0.3.4.3 cabal-version: >= 1.10 license: BSD2 license-file: LICENSE author: Bryan O'Sullivan <bos@serpentine.com>, Alexey Khudyakov <alexey.skladnoy@gmail.com> maintainer: Alexey Khudyakov <alexey.skladnoy@gmail.com>-homepage: https://github.com/bos/math-functions-bug-reports: https://github.com/bos/math-functions/issues+homepage: https://github.com/haskell/math-functions+bug-reports: https://github.com/haskell/math-functions/issues category: Math, Numeric build-type: Simple synopsis: Collection of tools for numeric computations@@ -20,19 +20,18 @@ polynomials. tested-with:- GHC ==7.4.2- || ==7.6.3- || ==7.8.4- || ==7.10.3- || ==8.0.2+ GHC ==8.0.2 || ==8.2.2 || ==8.4.4 || ==8.6.5 || ==8.8.4- || ==8.10.2- || ==9.0.1- , GHCJS ==8.4+ || ==8.10.7+ || ==9.0.2+ || ==9.2.7+ || ==9.4.5+ || ==9.6.2 + extra-source-files: changelog.md README.markdown@@ -44,11 +43,8 @@ doc/sinc.hs flag system-expm1- 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.+ description: Use expm1 provided by GHC. On GHCJS we don't have one so we+ have to use hand-coded one. default: True manual: True @@ -72,14 +68,14 @@ DeriveGeneric ghc-options: -Wall -O2- build-depends: base >= 4.5 && < 5+ build-depends: base >= 4.9 && < 5 , deepseq , data-default-class >= 0.1.2.0 , vector >= 0.11 , primitive- if flag(system-expm1) && !os(windows)+ if flag(system-expm1) && !impl(ghcjs) cpp-options: -DUSE_SYSTEM_EXPM1- if flag(system-erf) && !impl(ghcjs)+ if flag(system-erf) && !impl(ghcjs) cpp-options: -DUSE_SYSTEM_ERF exposed-modules: Numeric.MathFunctions.Constants@@ -131,7 +127,7 @@ benchmark math-functions-bench type: exitcode-stdio-1.0- if impl(ghc <= 7.10 ) || impl(ghcjs)+ if impl(ghcjs) buildable: False default-language: Haskell2010 other-extensions:@@ -152,7 +148,7 @@ , data-default-class , vector , random- , gauge >=0.2.5+ , tasty-bench >=0.3.4 source-repository head type: git
tests/Tests/SpecFunctions.hs view
@@ -22,7 +22,7 @@ import Tests.SpecFunctions.Tables import Numeric.SpecFunctions import Numeric.SpecFunctions.Internal (factorialTable)-import Numeric.MathFunctions.Comparison (within,relativeError,ulpDistance)+import Numeric.MathFunctions.Comparison (within,ulpDistance) import Numeric.MathFunctions.Constants (m_epsilon,m_tiny) erfTol,erfcTol,erfcLargeTol :: Int@@ -36,6 +36,11 @@ erfTol = 2 erfcTol = 2 erfcLargeTol = 2+-- Windows' one is not very good too+#elif defined(mingw32_HOST_OS)+erfTol = 2+erfcTol = 2+erfcLargeTol = 4 #else erfTol = 1 erfcTol = 2@@ -60,7 +65,7 @@ tests :: TestTree tests = testGroup "Special functions" [ testGroup "erf"- [ -- implementation from numerical recipes loses presision for+ [ -- implementation from numerical recipes loses precision for -- large arguments testCase "erfc table" $ forTable "tests/tables/erfc.dat" $ \[x, exact] ->@@ -236,6 +241,10 @@ roundtrip_erfc_invErfc = (4,4) roundtrip_invErfc_erfc = (4,4) roundtrip_erf_invErf = (2,2)+#elif defined(mingw32_HOST_OS)+roundtrip_erfc_invErfc = (4,4)+roundtrip_invErfc_erfc = (4,4)+roundtrip_erf_invErf = (4,4) #else roundtrip_erfc_invErfc = (2,2) roundtrip_invErfc_erfc = (2,2)
tests/Tests/Sum.hs view
@@ -4,54 +4,68 @@ import Control.Applicative ((<$>)) import Numeric.Sum as Sum+import Numeric.MathFunctions.Comparison import Prelude hiding (sum) import Test.Tasty (TestTree, testGroup)-import Test.Tasty.QuickCheck (testProperty)+import Test.Tasty.QuickCheck import Test.QuickCheck (Arbitrary(..)) import qualified Prelude -t_sum :: ([Double] -> Double) -> [Double] -> Bool-t_sum f xs = f xs == trueSum xs+-- Test that summation result is same as exact sum. That should pass+-- if we're effectively working with quad precision+t_sum :: ([Double] -> Double) -> [Double] -> Property+t_sum f xs+ = counterexample ("APPROX = " ++ show approx)+ $ counterexample ("EXACT = " ++ show exact)+ $ counterexample ("DELTA = " ++ show (approx - exact))+ $ counterexample ("ULPS = " ++ show (ulpDistance approx exact))+ $ approx == exact+ where+ approx = f xs+ exact = trueSum xs -t_sum_error :: ([Double] -> Double) -> [Double] -> Bool-t_sum_error f xs = abs (ts - f xs) <= abs (ts - Prelude.sum xs)- where ts = trueSum xs+-- Test that summation has smaller error than naive summation or no+-- worse than given number of ulps. If we're close enough to exact+-- answer naive may get ahead+t_sum_error :: ([Double] -> Double) -> [Double] -> Property+t_sum_error f xs+ = counterexample ("APPROX = " ++ show approx)+ $ counterexample ("NAIVE = " ++ show naive)+ $ counterexample ("EXACT = " ++ show exact)+ $ counterexample ("A-EXACT = " ++ show (approx - exact))+ $ counterexample ("N-EXACT = " ++ show (naive - exact))+ $ counterexample ("ULPS[A] = " ++ show (ulpDistance approx exact))+ $ counterexample ("ULPS[N] = " ++ show (ulpDistance naive exact))+ $ abs (exact - approx) <= abs (exact - naive)+ where+ naive = Prelude.sum xs+ approx = f xs+ exact = trueSum xs -t_sum_shifted :: ([Double] -> Double) -> [Double] -> Bool+t_sum_shifted :: ([Double] -> Double) -> [Double] -> Property t_sum_shifted f = t_sum_error f . zipWith (+) badvec trueSum :: (Fractional b, Real a) => [a] -> b trueSum xs = fromRational . Prelude.sum . map toRational $ xs badvec :: [Double]-badvec = cycle [1,1e16,-1e16]+badvec = cycle [1, 1e14, -1e14] tests :: TestTree-tests = testGroup "Summation" [- testGroup "ID" [- -- plain summation loses precision quickly- -- testProperty "t_sum" $ t_sum (sum id)-- -- tautological tests:- -- testProperty "t_sum_error" $ t_sum_error (sum id)- -- testProperty "t_sum_shifted" $ t_sum_shifted (sum id)- ]- , testGroup "Kahan" [- -- tests that cannot pass:- -- testProprty "t_sum" $ t_sum (sum kahan)- -- testProperty "t_sum_error" $ t_sum_error (sum kahan)-- -- kahan summation only beats normal summation with large values+tests = testGroup "Summation"+ [ testGroup "Kahan" [+ -- Kahan summation only beats naive summation when truly+ -- catastrophic cancellation occurs testProperty "t_sum_shifted" $ t_sum_shifted (sum kahan) ] , testGroup "KBN" [- testProperty "t_sum" $ t_sum (sum kbn)- , testProperty "t_sum_error" $ t_sum_error (sum kbn)+ testProperty "t_sum" $ t_sum (sum kbn)+ , testProperty "t_sum_error" $ t_sum_error (sum kbn) , testProperty "t_sum_shifted" $ t_sum_shifted (sum kbn) ] , testGroup "KB2" [- testProperty "t_sum" $ t_sum (sum kb2)- , testProperty "t_sum_error" $ t_sum_error (sum kb2)+ testProperty "t_sum" $ t_sum (sum kb2)+ , testProperty "t_sum_error" $ t_sum_error (sum kb2) , testProperty "t_sum_shifted" $ t_sum_shifted (sum kb2) ] ]
tests/tables/generate.py view
@@ -46,7 +46,7 @@ def load_inputs_cartesian(path):- "Load inputs for several variables where we want to genrate all pair"+ "Load inputs for several variables where we want to generate all pair" with open(path) as f: for x in itertools.product(*tokenize_stream(skip_comments(f))): yield x