diff --git a/Numeric/Polynomial.hs b/Numeric/Polynomial.hs --- a/Numeric/Polynomial.hs +++ b/Numeric/Polynomial.hs @@ -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. diff --git a/Numeric/Polynomial/Chebyshev.hs b/Numeric/Polynomial/Chebyshev.hs --- a/Numeric/Polynomial/Chebyshev.hs +++ b/Numeric/Polynomial/Chebyshev.hs @@ -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. +-- 16(4):254–256. -- -- * Clenshaw, C.W. (1962) Chebyshev series for mathematical -- functions. /National Physical Laboratory Mathematical Tables 5/, diff --git a/Numeric/RootFinding.hs b/Numeric/RootFinding.hs --- a/Numeric/RootFinding.hs +++ b/Numeric/RootFinding.hs @@ -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` () diff --git a/Numeric/SpecFunctions.hs b/Numeric/SpecFunctions.hs --- a/Numeric/SpecFunctions.hs +++ b/Numeric/SpecFunctions.hs @@ -1,4 +1,3 @@ -{-# LANGUAGE CPP #-} -- | -- Module : Numeric.SpecFunctions -- Copyright : (c) 2009, 2011, 2012 Bryan O'Sullivan diff --git a/Numeric/SpecFunctions/Compat.hs b/Numeric/SpecFunctions/Compat.hs --- a/Numeric/SpecFunctions/Compat.hs +++ b/Numeric/SpecFunctions/Compat.hs @@ -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. diff --git a/Numeric/SpecFunctions/Internal.hs b/Numeric/SpecFunctions/Internal.hs --- a/Numeric/SpecFunctions/Internal.hs +++ b/Numeric/SpecFunctions/Internal.hs @@ -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 diff --git a/Numeric/Sum.hs b/Numeric/Sum.hs --- a/Numeric/Sum.hs +++ b/Numeric/Sum.hs @@ -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. diff --git a/bench/bench.hs b/bench/bench.hs --- a/bench/bench.hs +++ b/bench/bench.hs @@ -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 ] diff --git a/changelog.md b/changelog.md --- a/changelog.md +++ b/changelog.md @@ -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 diff --git a/math-functions.cabal b/math-functions.cabal --- a/math-functions.cabal +++ b/math-functions.cabal @@ -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 , Alexey Khudyakov maintainer: Alexey Khudyakov -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 diff --git a/tests/Tests/SpecFunctions.hs b/tests/Tests/SpecFunctions.hs --- a/tests/Tests/SpecFunctions.hs +++ b/tests/Tests/SpecFunctions.hs @@ -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) diff --git a/tests/Tests/Sum.hs b/tests/Tests/Sum.hs --- a/tests/Tests/Sum.hs +++ b/tests/Tests/Sum.hs @@ -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) ] ] diff --git a/tests/tables/generate.py b/tests/tables/generate.py --- a/tests/tables/generate.py +++ b/tests/tables/generate.py @@ -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