summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorAlexeyKhudyakov <>2018-07-11 16:23:00 (GMT)
committerhdiff <hdiff@hdiff.luite.com>2018-07-11 16:23:00 (GMT)
commit50b0ab2a976531d8b631f75dfb66ac432266c376 (patch)
tree7761379d04bf13c1fd03f56722a7d9819313003d
parentd93234603620c56692c8d21f405021bc43b273bb (diff)
version 0.3.0.00.3.0.0
-rw-r--r--Numeric/RootFinding.hs429
-rw-r--r--Numeric/SpecFunctions/Internal.hs11
-rw-r--r--Numeric/Sum.hs36
-rw-r--r--README.markdown20
-rw-r--r--benchmark/Summation.hs15
-rw-r--r--benchmark/bench.hs89
-rw-r--r--changelog.md12
-rw-r--r--math-functions.cabal52
-rw-r--r--tests/Tests/RootFinding.hs44
-rw-r--r--tests/Tests/SpecFunctions.hs16
-rw-r--r--tests/tests.hs6
11 files changed, 493 insertions, 237 deletions
diff --git a/Numeric/RootFinding.hs b/Numeric/RootFinding.hs
index 5985427..0f90f39 100644
--- a/Numeric/RootFinding.hs
+++ b/Numeric/RootFinding.hs
@@ -1,33 +1,63 @@
-{-# LANGUAGE BangPatterns, DeriveDataTypeable, DeriveGeneric, CPP #-}
+{-# LANGUAGE BangPatterns #-}
+{-# LANGUAGE CPP #-}
+{-# LANGUAGE DeriveDataTypeable #-}
+{-# LANGUAGE DeriveFoldable #-}
+{-# LANGUAGE DeriveGeneric #-}
+{-# LANGUAGE DeriveTraversable #-}
+{-# LANGUAGE TypeFamilies #-}
-- |
-- Module : Numeric.RootFinding
--- Copyright : (c) 2011 Bryan O'Sullivan
+-- Copyright : (c) 2011 Bryan O'Sullivan, 2018 Alexey Khudyakov
-- License : BSD3
--
-- Maintainer : bos@serpentine.com
-- Stability : experimental
-- Portability : portable
--
--- Haskell functions for finding the roots of real functions of real arguments.
+-- Haskell functions for finding the roots of real functions of real
+-- arguments. These algorithms are iterative so we provide both
+-- function returning root (or failure to find root) and list of
+-- iterations.
module Numeric.RootFinding
- (
+ ( -- * Data types
Root(..)
, fromRoot
+ , Tolerance(..)
+ , withinTolerance
+ , IterationStep(..)
+ , findRoot
+ -- * Ridders algorithm
+ , RiddersParam(..)
, ridders
+ , riddersIterations
+ -- * Newton-Raphson algorithm
+ , NewtonParam(..)
, newtonRaphson
+ , newtonRaphsonIterations
-- * References
-- $references
) where
import Control.Applicative (Alternative(..), Applicative(..))
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)
+import Numeric.MathFunctions.Comparison (within,eqRelErr)
+import Numeric.MathFunctions.Constants (m_epsilon)
+
+----------------------------------------------------------------
+-- Data types
+----------------------------------------------------------------
+
-- | The result of searching for a root of a mathematical function.
data Root a = NotBracketed
-- ^ The function does not have opposite signs when
@@ -35,139 +65,352 @@ data Root a = NotBracketed
| SearchFailed
-- ^ The search failed to converge to within the given
-- error tolerance after the given number of iterations.
- | Root a
+ | Root !a
-- ^ A root was successfully found.
- deriving (Eq, Read, Show, Typeable, Data
+ deriving (Eq, Read, Show, Typeable, Data, Foldable, Traversable
#if __GLASGOW_HASKELL__ > 704
, Generic
#endif
)
+instance (NFData a) => NFData (Root a) where
+ rnf NotBracketed = ()
+ rnf SearchFailed = ()
+ rnf (Root a) = rnf a
instance Functor Root where
fmap _ NotBracketed = NotBracketed
fmap _ SearchFailed = SearchFailed
fmap f (Root a) = Root (f a)
+instance Applicative Root where
+ pure = return
+ (<*>) = ap
+
instance Monad Root where
NotBracketed >>= _ = NotBracketed
SearchFailed >>= _ = SearchFailed
- Root a >>= m = m a
-
+ Root a >>= f = f a
return = Root
instance MonadPlus Root where
- mzero = SearchFailed
-
- r@(Root _) `mplus` _ = r
- _ `mplus` p = p
-
-instance Applicative Root where
- pure = Root
- (<*>) = ap
+ mzero = empty
+ mplus = (<|>)
instance Alternative Root where
- empty = SearchFailed
-
- r@(Root _) <|> _ = r
- _ <|> p = p
+ empty = NotBracketed
+ r@Root{} <|> _ = r
+ _ <|> r@Root{} = r
+ NotBracketed <|> r = r
+ r <|> NotBracketed = r
+ _ <|> r = r
-- | Returns either the result of a search for a root, or the default
-- value if the search failed.
-fromRoot :: a -- ^ Default value.
- -> Root a -- ^ Result of search for a root.
+fromRoot :: a -- ^ Default value.
+ -> Root a -- ^ Result of search for a root.
-> a
fromRoot _ (Root a) = a
fromRoot a _ = a
--- | Use the method of Ridders to compute a root of a function.
+-- | Error tolerance for finding root. It describes when root finding
+-- algorithm should stop trying to improve approximation.
+data Tolerance
+ = RelTol !Double
+ -- ^ Relative error tolerance. Given @RelTol ε@ two values are
+ -- considered approximately equal if
+ -- \[ |a - b| / |\operatorname{max}(a,b)} < \vareps \]
+ | AbsTol !Double
+ -- ^ Absolute error tolerance. Given @AbsTol δ@ two values are
+ -- 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
+ )
+
+-- | Check that two values are approximately equal. In addition to
+-- specification values are considered equal if they're within 1ulp
+-- of precision. No further improvement could be done anyway.
+withinTolerance :: Tolerance -> Double -> Double -> Bool
+withinTolerance _ a b
+ | within 1 a b = True
+withinTolerance (RelTol eps) a b = eqRelErr eps a b
+withinTolerance (AbsTol tol) a b = abs (a - b) < tol
+
+-- | Type class for checking whether iteration converged already.
+class IterationStep a where
+ -- | Return @Just root@ is current iteration converged within
+ -- required error tolerance. Returns @Nothing@ otherwise.
+ matchRoot :: Tolerance -> a -> Maybe (Root Double)
+
+-- | Find root in lazy list of iterations.
+findRoot :: IterationStep a
+ => Int -- ^ Maximum
+ -> Tolerance -- ^ Error tolerance
+ -> [a]
+ -> Root Double
+findRoot maxN tol = go 0
+ where
+ go !i _ | i >= maxN = SearchFailed
+ go !_ [] = SearchFailed
+ go i (x:xs) = case matchRoot tol x of
+ Just r -> r
+ Nothing -> go (i+1) xs
+{-# INLINABLE findRoot #-}
+{-# SPECIALIZE findRoot :: Int -> Tolerance -> [RiddersStep] -> Root Double #-}
+{-# SPECIALIZE findRoot :: Int -> Tolerance -> [NewtonStep] -> Root Double #-}
+
+
+----------------------------------------------------------------
+-- Attaching information to roots
+----------------------------------------------------------------
+
+-- | Parameters for 'ridders' root finding
+data RiddersParam = RiddersParam
+ { riddersMaxIter :: !Int
+ -- ^ Maximum number of iterations.
+ , riddersTol :: !Tolerance
+ -- ^ Error tolerance for root approximation.
+ }
+ deriving (Eq, Read, Show, Typeable, Data
+#if __GLASGOW_HASKELL__ > 704
+ , Generic
+#endif
+ )
+
+instance Default RiddersParam where
+ def = RiddersParam
+ { riddersMaxIter = 100
+ , riddersTol = RelTol (4 * m_epsilon)
+ }
+
+-- | Single Ridders step. It's a bracket of root
+data RiddersStep
+ = RiddersStep !Double !Double
+ -- ^ Ridders step. Parameters are bracket for the root
+ | RiddersBisect !Double !Double
+ -- ^ Bisection step. It's fallback which is taken when Ridders
+ -- update takes us out of bracket
+ | RiddersRoot !Double
+ -- ^ Root found
+ | RiddersNoBracket
+ -- ^ Root is not bracketed
+ deriving (Eq, Read, Show, Typeable, Data
+#if __GLASGOW_HASKELL__ > 704
+ , Generic
+#endif
+ )
+
+instance NFData RiddersStep where
+ rnf x = x `seq` ()
+
+instance IterationStep RiddersStep where
+ matchRoot tol r = case r of
+ RiddersRoot x -> Just $ Root x
+ RiddersNoBracket -> Just NotBracketed
+ RiddersStep a b
+ | withinTolerance tol a b -> Just $ Root ((a + b) / 2)
+ | otherwise -> Nothing
+ RiddersBisect a b
+ | withinTolerance tol a b -> Just $ Root ((a + b) / 2)
+ | otherwise -> Nothing
+
+
+-- | Use the method of Ridders[Ridders1979] to compute a root of a
+-- function. It doesn't require derivative and provide quadratic
+-- convergence (number of significant digits grows quadratically
+-- with number of iterations).
--
--- The function must have opposite signs when evaluated at the lower
--- and upper bounds of the search (i.e. the root must be bracketed).
-ridders :: Double -- ^ Absolute error tolerance.
- -> (Double,Double) -- ^ Lower and upper bounds for the search.
- -> (Double -> Double) -- ^ Function to find the roots of.
- -> Root Double
-ridders tol (lo,hi) f
- | flo == 0 = Root lo
- | fhi == 0 = Root hi
- | flo*fhi > 0 = NotBracketed -- root is not bracketed
- | otherwise = go lo flo hi fhi 0
+-- The function must have opposite signs when evaluated at the lower
+-- and upper bounds of the search (i.e. the root must be
+-- bracketed). If there's more that one root in the bracket
+-- iteration will converge to some root in the bracket.
+ridders
+ :: RiddersParam -- ^ Parameters for algorithms. @def@
+ -- provides reasonable defaults
+ -> (Double,Double) -- ^ Bracket for root
+ -> (Double -> Double) -- ^ Function to find roots
+ -> Root Double
+ridders p bracket fun
+ = findRoot (riddersMaxIter p) (riddersTol p)
+ $ riddersIterations bracket fun
+
+-- | List of iterations for Ridders methods. See 'ridders' for
+-- documentation of parameters
+riddersIterations :: (Double,Double) -> (Double -> Double) -> [RiddersStep]
+riddersIterations (lo,hi) f
+ | flo == 0 = [RiddersRoot lo]
+ | fhi == 0 = [RiddersRoot hi]
+ -- root is not bracketed
+ | flo*fhi > 0 = [RiddersNoBracket]
+ -- Ensure that a<b in iterations
+ | lo < hi = RiddersStep lo hi : go lo flo hi fhi
+ | otherwise = RiddersStep lo hi : go hi fhi lo flo
where
- go !a !fa !b !fb !i
- -- Root is bracketed within 1 ulp. No improvement could be made
- | within 1 a b = Root a
- -- Root is found. Check that f(m) == 0 is nessesary to ensure
- -- that root is never passed to 'go'
- | fm == 0 = Root m
- | fn == 0 = Root n
- | d < tol = Root n
- -- Too many iterations performed. Fail
- | i >= (100 :: Int) = SearchFailed
- -- Ridder's approximation coincide with one of old
- -- bounds. Revert to bisection
- | n == a || n == b = case () of
- _| fm*fa < 0 -> go a fa m fm (i+1)
- | otherwise -> go m fm b fb (i+1)
- -- Proceed as usual
- | fn*fm < 0 = go n fn m fm (i+1)
- | fn*fa < 0 = go a fa n fn (i+1)
- | otherwise = go n fn b fb (i+1)
+ flo = f lo
+ fhi = f hi
+ --
+ go !a !fa !b !fb
+ | fm == 0 = [RiddersRoot m]
+ | fn == 0 = [RiddersRoot n]
+ -- Ridder's approximation coincide with one of old bounds or
+ -- went out of (a,b) range due to numerical problems. Revert
+ -- to bisection
+ | n <= a || n >= b = case () of
+ _| fm*fa < 0 -> recBisect a fa m fm
+ | otherwise -> recBisect m fm b fb
+ | fn*fm < 0 = recRidders n fn m fm
+ | fn*fa < 0 = recRidders a fa n fn
+ | otherwise = recRidders n fn b fb
where
- d = abs (b - a)
- dm = (b - a) * 0.5
- !m = a + dm
- !fm = f m
- !dn = signum (fb - fa) * dm * fm / sqrt(fm*fm - fa*fb)
- !n = m - signum dn * min (abs dn) (abs dm - 0.5 * tol)
- !fn = f n
- !flo = f lo
- !fhi = f hi
+ recBisect x fx y fy = RiddersBisect x y : go x fx y fy
+ recRidders x fx y fy = RiddersStep x y : go x fx y fy
+ --
+ dm = (b - a) * 0.5
+ -- Mean point
+ m = (a + b) / 2
+ fm = f m
+ -- Ridders update
+ n = m - signum (fb - fa) * dm * fm / sqrt(fm*fm - fa*fb)
+ fn = f n
+
+
+
+----------------------------------------------------------------
+-- Newton-Raphson algorithm
+----------------------------------------------------------------
+
+-- | Parameters for 'ridders' root finding
+data NewtonParam = NewtonParam
+ { newtonMaxIter :: !Int
+ -- ^ Maximum number of iterations.
+ , newtonTol :: !Tolerance
+ -- ^ Error tolerance for root approximation.
+ }
+ deriving (Eq, Read, Show, Typeable, Data
+#if __GLASGOW_HASKELL__ > 704
+ , Generic
+#endif
+ )
+
+instance Default NewtonParam where
+ def = NewtonParam
+ { newtonMaxIter = 50
+ , newtonTol = RelTol (4 * m_epsilon)
+ }
+
+-- | Steps for Newton iterations
+data NewtonStep
+ = NewtonStep !Double !Double
+ -- ^ Normal Newton-Raphson update. Parameters are: old guess, new guess
+ | NewtonBisection !Double !Double
+ -- ^ Bisection fallback when Newton-Raphson iteration doesn't
+ -- work. Parameters are bracket on root
+ | NewtonRoot !Double
+ -- ^ Root is found
+ | NewtonNoBracket
+ -- ^ Root is not bracketed
+ deriving (Eq, Read, Show, Typeable, Data
+#if __GLASGOW_HASKELL__ > 704
+ , Generic
+#endif
+ )
+instance NFData NewtonStep where
+ rnf x = x `seq` ()
+
+instance IterationStep NewtonStep where
+ matchRoot tol r = case r of
+ NewtonRoot x -> Just (Root x)
+ NewtonNoBracket -> Just NotBracketed
+ NewtonStep x x'
+ | withinTolerance tol x x' -> Just (Root x')
+ | otherwise -> Nothing
+ NewtonBisection a b
+ | withinTolerance tol a b -> Just (Root ((a + b) / 2))
+ | otherwise -> Nothing
+ {-# INLINE matchRoot #-}
-- | Solve equation using Newton-Raphson iterations.
--
--- This method require both initial guess and bounds for root. If
--- Newton step takes us out of bounds on root function reverts to
--- bisection.
+-- This method require both initial guess and bounds for root. If
+-- Newton step takes us out of bounds on root function reverts to
+-- bisection.
newtonRaphson
- :: Double
- -- ^ Required precision
- -> (Double,Double,Double)
- -- ^ (lower bound, initial guess, upper bound). Iterations will no
- -- go outside of the interval
- -> (Double -> (Double,Double))
- -- ^ Function to finds roots. It returns pair of function value and
- -- its derivative
+ :: NewtonParam -- ^ Parameters for algorithm. @def@
+ -- provide reasonable defaults.
+ -> (Double,Double,Double) -- ^ Triple of @(low bound, initial
+ -- guess, upper bound)@. If initial
+ -- guess if out of bracket middle
+ -- of bracket is taken as
+ -- approximation
+ -> (Double -> (Double,Double)) -- ^ Function to find root of. It
+ -- returns pair of function value and
+ -- its first derivative
-> Root Double
-newtonRaphson !prec (!low,!guess,!hi) function
- = go low guess hi
+newtonRaphson p guess fun
+ = findRoot (newtonMaxIter p) (newtonTol p)
+ $ newtonRaphsonIterations guess fun
+
+-- | List of iteration for Newton-Raphson algorithm. See documentation
+-- for 'newtonRaphson' for meaning of parameters.
+newtonRaphsonIterations :: (Double,Double,Double) -> (Double -> (Double,Double)) -> [NewtonStep]
+newtonRaphsonIterations (lo,guess,hi) function
+ | flo == 0 = [NewtonRoot lo]
+ | fhi == 0 = [NewtonRoot hi]
+ | flo*fhi > 0 = [NewtonNoBracket]
+ -- Ensure that function value on low bound is negative
+ | flo > 0 = go hi guess' lo
+ | otherwise = go lo guess hi
where
- go !xMin !x !xMax
- | f == 0 = Root x
- | abs (dx / x) < prec = Root x
- | otherwise = go xMin' x' xMax'
+ (flo,_) = function lo
+ (fhi,_) = function hi
+ -- Ensure that initial guess is within bracket
+ guess'
+ | guess >= lo && guess <= hi = guess
+ | guess >= hi && guess <= lo = guess
+ | otherwise = (lo + hi) / 2
+ -- Newton iterations. Invariant:
+ -- > f xA < 0
+ -- > f xB > 0
+ go xA x xB
+ | f == 0 = [NewtonRoot x]
+ | f' == 0 = bisectionStep
+ -- Accept Newton step since it stays within bracket.
+ | (x' - xA) * (x' - xB) < 0 = newtonStep
+ -- Otherwise bracket root and pick new approximation as
+ -- midpoint.
+ | otherwise = bisectionStep
where
+ -- Calculate Newton step
(f,f') = function x
- -- Calculate Newton-Raphson step
- delta | f' == 0 = error "handle f'==0"
- | otherwise = f / f'
- -- Calculate new approximation and actual change of approximation
- (dx,x') | z <= xMin = let d = 0.5*(x - xMin) in (d, x - d)
- | z >= xMax = let d = 0.5*(x - xMax) in (d, x - d)
- | otherwise = (delta, z)
- where z = x - delta
- -- Update root bracket
- xMin' | dx < 0 = x
- | otherwise = xMin
- xMax' | dx > 0 = x
- | otherwise = xMax
+ x' = x - f / f'
+ -- Newton step
+ newtonStep
+ | f > 0 = NewtonStep x x' : go xA x' x
+ | otherwise = NewtonStep x x' : go x x' xB
+ -- Fallback bisection step
+ bisectionStep
+ | f > 0 = NewtonBisection xA x : go xA ((xA + x) / 2) x
+ | otherwise = NewtonBisection x xB : go x ((x + xB) / 2) xB
+----------------------------------------------------------------
+-- Internal functions
+----------------------------------------------------------------
+
-- $references
--
-- * Ridders, C.F.J. (1979) A new algorithm for computing a single
-- root of a real continuous function.
-- /IEEE Transactions on Circuits and Systems/ 26:979&#8211;980.
+--
+-- * Press W.H.; Teukolsky S.A.; Vetterling W.T.; Flannery B.P.
+-- (2007). \"Section 9.2.1. Ridders' Method\". /Numerical Recipes: The
+-- Art of Scientific Computing (3rd ed.)./ New York: Cambridge
+-- University Press. ISBN 978-0-521-88068-8.
diff --git a/Numeric/SpecFunctions/Internal.hs b/Numeric/SpecFunctions/Internal.hs
index a45ac06..386c8cc 100644
--- a/Numeric/SpecFunctions/Internal.hs
+++ b/Numeric/SpecFunctions/Internal.hs
@@ -14,9 +14,10 @@ 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 Data.Bits ((.&.), (.|.), shiftR)
+import Data.Int (Int64)
+import Data.Word (Word)
+import Data.Default.Class
import qualified Data.Vector.Unboxed as U
import Data.Vector.Unboxed ((!))
import Text.Printf
@@ -26,7 +27,7 @@ import GHC.Float (log1p,expm1)
import Numeric.Polynomial.Chebyshev (chebyshevBroucke)
import Numeric.Polynomial (evaluatePolynomialL,evaluateEvenPolynomialL,evaluateOddPolynomialL)
-import Numeric.RootFinding (Root(..), newtonRaphson)
+import Numeric.RootFinding (Root(..), newtonRaphson, NewtonParam(..), Tolerance(..))
import Numeric.Series
import Numeric.MathFunctions.Constants
@@ -652,7 +653,7 @@ invIncBetaGuess beta a b p
func x = ( u + log x + mu*log(1 - x)
, 1/x - mu/(1-x)
)
- Root x0 = newtonRaphson 1e-8 (lower, x_guess, upper) func
+ Root x0 = newtonRaphson def{newtonTol=RelTol 1e-8} (lower, x_guess, upper) func
in x0
-- For large a and b approximation from AS109 (Carter
-- approximation). It's reasonably good in this region
diff --git a/Numeric/Sum.hs b/Numeric/Sum.hs
index b1b3bdf..d890622 100644
--- a/Numeric/Sum.hs
+++ b/Numeric/Sum.hs
@@ -1,5 +1,5 @@
{-# LANGUAGE BangPatterns, DeriveDataTypeable, FlexibleContexts,
- MultiParamTypeClasses, TemplateHaskell, TypeFamilies #-}
+ MultiParamTypeClasses, TemplateHaskell, TypeFamilies, CPP #-}
{-# OPTIONS_GHC -fno-warn-name-shadowing #-}
-- |
-- Module : Numeric.Sum
@@ -53,7 +53,11 @@ import Control.Arrow ((***))
import Control.DeepSeq (NFData(..))
import Data.Bits (shiftR)
import Data.Data (Typeable, Data)
-import Data.Vector.Generic (Vector(..), foldl')
+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.Unboxed.Deriving (derivingUnbox)
-- Needed for GHC 7.2 & 7.4 to derive Unbox instances
import Data.Vector.Generic.Mutable (MVector(..))
@@ -106,6 +110,15 @@ instance Summation KahanSum where
instance NFData KahanSum where
rnf !_ = ()
+instance Monoid KahanSum where
+ mempty = zero
+ s `mappend` KahanSum s' _ = add s s'
+
+#if MIN_VERSION_base(4,9,0)
+instance Semigroup KahanSum where
+ (<>) = mappend
+#endif
+
kahanAdd :: KahanSum -> Double -> KahanSum
kahanAdd (KahanSum sum c) x = KahanSum sum' c'
where sum' = sum + y
@@ -134,6 +147,15 @@ instance Summation KBNSum where
instance NFData KBNSum where
rnf !_ = ()
+instance Monoid KBNSum where
+ mempty = zero
+ s `mappend` KBNSum s' c' = add (add s s') c'
+
+#if MIN_VERSION_base(4,9,0)
+instance Semigroup KBNSum where
+ (<>) = mappend
+#endif
+
kbnAdd :: KBNSum -> Double -> KBNSum
kbnAdd (KBNSum sum c) x = KBNSum sum' c'
where c' | abs sum >= abs x = c + ((sum - sum') + x)
@@ -169,6 +191,16 @@ instance Summation KB2Sum where
instance NFData KB2Sum where
rnf !_ = ()
+instance Monoid KB2Sum where
+ mempty = zero
+ s `mappend` KB2Sum s' c' cc' = add (add (add s s') c') cc'
+
+#if MIN_VERSION_base(4,9,0)
+instance Semigroup KB2Sum where
+ (<>) = mappend
+#endif
+
+
kb2Add :: KB2Sum -> Double -> KB2Sum
kb2Add (KB2Sum sum c cc) x = KB2Sum sum' c' cc'
where sum' = sum + x
diff --git a/README.markdown b/README.markdown
index ec09286..3ab920d 100644
--- a/README.markdown
+++ b/README.markdown
@@ -1,7 +1,21 @@
-# math-functions: efficient, special purpose mathematical functions
+# math-functions: collection of tools for numeric computations
-This package provides a number of special-purpose mathematical
-functions used in statistical and numerical computing.
+[![Build Status](https://travis-ci.org/Shimuuar/math-functions.png?branch=master)](https://travis-ci.org/Shimuuar/math-functions)
+[![Build status](https://ci.appveyor.com/api/projects/status/6xexxj9g6rnbg2q4/branch/master?svg=true)](https://ci.appveyor.com/project/Shimuuar/math-functions/branch/master)
+
+This package provides collection of various tools for numeric
+computations. Namely:
+
+ - Number pure haskell implementations of special function which are used in
+ statistical and numerical computing.
+
+ - Compensated summation (Kahan summation) which allows to
+
+ - Root finding for functions of single real variable
+
+ - Series summation
+
+ - Functions for comparing IEEE754 numbers
Where possible, we give citations and computational complexity
estimates for the algorithms used.
diff --git a/benchmark/Summation.hs b/benchmark/Summation.hs
deleted file mode 100644
index 0b9dbad..0000000
--- a/benchmark/Summation.hs
+++ /dev/null
@@ -1,15 +0,0 @@
-import Criterion.Main
-import Numeric.Sum as Sum
-import System.Random.MWC
-import qualified Data.Vector.Unboxed as U
-
-main = do
- gen <- createSystemRandom
- v <- uniformVector gen 10000000 :: IO (U.Vector Double)
- defaultMain [
- bench "naive" $ whnf U.sum v
- , bench "pairwise" $ whnf pairwiseSum v
- , bench "kahan" $ whnf (sumVector kahan) v
- , bench "kbn" $ whnf (sumVector kbn) v
- , bench "kb2" $ whnf (sumVector kb2) v
- ]
diff --git a/benchmark/bench.hs b/benchmark/bench.hs
deleted file mode 100644
index 2fb1abb..0000000
--- a/benchmark/bench.hs
+++ /dev/null
@@ -1,89 +0,0 @@
-import Criterion.Main
-import qualified Data.Vector.Unboxed as U
-import Numeric.SpecFunctions
-import Numeric.Polynomial
-import Text.Printf
-
--- Uniformly sample logGamma performance between 10^-6 to 10^6
-benchmarkLogGamma logG =
- [ bench (printf "%.3g" x) $ nf logG x
- | x <- [ m * 10**n | n <- [ -8 .. 8 ]
- , m <- [ 10**(i / tics) | i <- [0 .. tics-1] ]
- ]
- ]
- where tics = 3
-{-# INLINE benchmarkLogGamma #-}
-
-
--- Power of polynomial to be evaluated (In other words length of coefficients vector)
-coef_size :: [Int]
-coef_size = [ 1,2,3,4,5,6,7,8,9
- , 10, 30
- , 100, 300
- , 1000, 3000
- , 10000, 30000
- ]
-{-# INLINE coef_size #-}
-
--- Precalculated coefficients
-coef_list :: [U.Vector Double]
-coef_list = [ U.replicate n 1.2 | n <- coef_size]
-{-# NOINLINE coef_list #-}
-
-
-
-main :: IO ()
-main = defaultMain
- [ bgroup "logGamma" $
- benchmarkLogGamma logGamma
- , bgroup "logGammaL" $
- benchmarkLogGamma logGammaL
- , bgroup "incompleteGamma" $
- [ bench (show p) $ nf (incompleteGamma p) p
- | p <- [ 0.1
- , 1, 3
- , 10, 30
- , 100, 300
- , 999, 1000
- ]
- ]
- , bgroup "factorial"
- [ bench (show n) $ nf factorial n
- | n <- [ 0, 1, 3, 6, 9, 11, 15
- , 20, 30, 40, 50, 60, 70, 80, 90, 100
- ]
- ]
- , bgroup "incompleteBeta"
- [ bench (show (p,q,x)) $ nf (incompleteBeta p q) x
- | (p,q,x) <- [ (10, 10, 0.5)
- , (101, 101, 0.5)
- , (1010, 1010, 0.5)
- , (10100, 10100, 0.5)
- , (100100, 100100, 0.5)
- , (1001000, 1001000, 0.5)
- , (10010000,10010000,0.5)
- ]
- ]
- , bgroup "log1p"
- [ bench (show x) $ nf log1p x
- | x <- [ -0.9
- , -0.5
- , -0.1
- , 0.1
- , 0.5
- , 1
- , 10
- , 100
- ]
- ]
- , bgroup "sinc" $
- bench "sin" (nf sin (0.55 :: Double))
- : [ bench (show x) $ nf sinc x
- | x <- [0, 1e-6, 1e-3, 0.5]
- ]
- , bgroup "poly"
- $ [ bench ("vector_"++show (U.length coefs)) $ nf (\x -> evaluatePolynomial x coefs) (1 :: Double)
- | coefs <- coef_list ]
- ++ [ bench ("unpacked_"++show n) $ nf (\x -> evaluatePolynomialL x (map fromIntegral [1..n])) (1 :: Double)
- | n <- coef_size ]
- ]
diff --git a/changelog.md b/changelog.md
index b9cf4e4..06c0bad 100644
--- a/changelog.md
+++ b/changelog.md
@@ -1,3 +1,13 @@
+## Changes in 0.3.0.0
+
+ * `Semigroup` and `Monoid` instances added for data types from `Numeric.Sum`
+
+ * API for finding roots of real functions reworked. 1) All algorithm
+ parameters are now tweakable. 2) Functions for getting list of iterations
+ added.
+
+ * `Foldable` and `Traversable` instances for `Root` were added.
+
## Changes in 0.2.1.0
* `log1p` and `expm1` are simply reexported from `GHC.Float`. They're methods
@@ -20,7 +30,7 @@
* New much more precise implementation for `incompleteGamma`
- * Dependency on `erf` pacakge dropped. `erf` and `erfc` just do direct calls
+ * Dependency on `erf` package dropped. `erf` and `erfc` just do direct calls
to C.
* `Numeric.SpecFunctions.expm1` added
diff --git a/math-functions.cabal b/math-functions.cabal
index c5ad286..082b031 100644
--- a/math-functions.cabal
+++ b/math-functions.cabal
@@ -1,25 +1,27 @@
name: math-functions
-version: 0.2.1.0
+version: 0.3.0.0
cabal-version: >= 1.10
license: BSD3
license-file: LICENSE
author: Bryan O'Sullivan <bos@serpentine.com>,
- Aleksey Khudyakov <alexey.skladnoy@gmail.com>
-maintainer: 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
category: Math, Numeric
build-type: Simple
-synopsis: Special functions and Chebyshev polynomials
+synopsis: Collection of tools for numeric computations
description:
- This library provides implementations of special mathematical
- functions and Chebyshev polynomials. These functions are often
- useful in statistical and numerical computing.
+
+ This library contain collection of various utilities for numerical
+ computing. So far there're special mathematical functions,
+ compensated summation algorithm, summation of series, root finding
+ for real functions, polynomial summation and Chebyshev
+ polynomials.
extra-source-files:
changelog.md
README.markdown
- benchmark/*.hs
tests/*.hs
tests/Tests/*.hs
tests/Tests/SpecFunctions/gen.py
@@ -39,11 +41,12 @@ library
DeriveGeneric
ghc-options: -Wall -O2
- build-depends: base >=4.5 && <5
+ build-depends: base >= 4.5 && < 5
, deepseq
- , vector >= 0.7
+ , data-default-class >= 0.1.2.0
+ , vector >= 0.7
, primitive
- , vector-th-unbox
+ , vector-th-unbox >= 0.2.1.6
if flag(system-expm1) || !os(windows)
cpp-options: -DUSE_SYSTEM_EXPM1
exposed-modules:
@@ -79,22 +82,23 @@ test-suite tests
Tests.Helpers
Tests.Chebyshev
Tests.Comparison
+ Tests.RootFinding
Tests.SpecFunctions
Tests.SpecFunctions.Tables
Tests.Sum
- build-depends:
- math-functions,
- base >=4.5 && <5,
- deepseq,
- primitive,
- vector >= 0.7,
- vector-th-unbox,
- erf,
- HUnit >= 1.2,
- QuickCheck >= 2.7,
- test-framework,
- test-framework-hunit,
- test-framework-quickcheck2
+ build-depends: base >=4.5 && <5
+ , math-functions
+ , data-default-class >= 0.1.2.0
+ , deepseq
+ , primitive
+ , vector >= 0.7
+ , vector-th-unbox
+ , erf
+ , HUnit >= 1.2
+ , QuickCheck >= 2.7
+ , test-framework
+ , test-framework-hunit
+ , test-framework-quickcheck2
source-repository head
type: git
diff --git a/tests/Tests/RootFinding.hs b/tests/Tests/RootFinding.hs
new file mode 100644
index 0000000..745130f
--- /dev/null
+++ b/tests/Tests/RootFinding.hs
@@ -0,0 +1,44 @@
+-- |
+module Tests.RootFinding ( tests ) where
+
+import Data.Default.Class
+import Test.Framework
+import Test.Framework.Providers.HUnit
+
+import Numeric.RootFinding
+import Tests.Helpers
+
+
+tests :: Test
+tests = testGroup "Root finding"
+ [ testGroup "Ridders"
+ [ testAssertion "sin x - 0.525 [exact]" $ testRiddersSin0_525 (AbsTol 0)
+ , testAssertion "sin x - 0.525 [abs 1e-12]" $ testRiddersSin0_525 (AbsTol 1e-12)
+ , testAssertion "sin x - 0.525 [abs 1e-6]" $ testRiddersSin0_525 (AbsTol 1e-6)
+ , testAssertion "sin x - 0.525 [rel 1e-12]" $ testRiddersSin0_525 (RelTol 1e-12)
+ , testAssertion "sin x - 0.525 [rel 1e-6]" $ testRiddersSin0_525 (RelTol 1e-6)
+ ]
+ , testGroup "Newton-Raphson"
+ [ testAssertion "sin x - 0.525 [rel 1e-12]" $ testNewtonSin0_525 (RelTol 1e-12)
+ , testAssertion "sin x - 0.525 [rel 1e-6]" $ testNewtonSin0_525 (RelTol 1e-6)
+ , testAssertion "sin x - 0.525 [abs 1e-12]" $ testNewtonSin0_525 (AbsTol 1e-12)
+ , testAssertion "sin x - 0.525 [abs 1e-6]" $ testNewtonSin0_525 (AbsTol 1e-6)
+ , testAssertion "1/x - 0.5 [0]" $
+ let Root r = newtonRaphson def{newtonTol=RelTol 0} (1,1000,1000)
+ (\x -> (1/x - 0.5, -1/(x*x)))
+ in r == 2
+ ]
+ ]
+ where
+ -- Exact root for equation: sin x - 0.525 = 0
+ exactRoot = 0.5527151130967832
+ --
+ testRiddersSin0_525 tol
+ = withinTolerance tol r exactRoot
+ where
+ Root r = ridders def{riddersTol = tol} (0, pi/2) (\x -> sin x - 0.525)
+ --
+ testNewtonSin0_525 tol
+ = withinTolerance tol r exactRoot
+ where
+ Root r = newtonRaphson def{newtonTol=tol} (0, pi/4, pi/2) (\x -> (sin x - 0.525, cos x))
diff --git a/tests/Tests/SpecFunctions.hs b/tests/Tests/SpecFunctions.hs
index c5d6ff3..c4da33d 100644
--- a/tests/Tests/SpecFunctions.hs
+++ b/tests/Tests/SpecFunctions.hs
@@ -10,6 +10,8 @@ import Data.Vector ((!))
import Test.QuickCheck hiding (choose,within)
import Test.Framework
import Test.Framework.Providers.QuickCheck2
+import Test.Framework.Providers.HUnit
+import Test.HUnit (assertBool)
import Tests.Helpers
import Tests.SpecFunctions.Tables
@@ -65,8 +67,16 @@ tests = testGroup "Special functions"
-- Relative precision is lost when digamma(x) ≈ 0
, testAssertion "digamma is expected to be precise at 1e-12"
$ and [ eq 1e-12 r (digamma x) | (x,r) <- tableDigamma ]
- , testAssertion "incompleteBeta is expected to be precise at 32*m_epsilon level"
- $ and [ eq (32 * m_epsilon) (incompleteBeta p q x) ib | (p,q,x,ib) <- tableIncompleteBeta ]
+ --
+ , let deviations = [ ( "p=",p, "q=",q, "x=",x
+ , "ib=",ib, "ib'=",ib'
+ , "err=",relativeError ib ib' / m_epsilon)
+ | (p,q,x,ib) <- tableIncompleteBeta
+ , let ib' = incompleteBeta p q x
+ , not $ eq (64 * m_epsilon) ib' ib
+ ]
+ in testCase "incompleteBeta is expected to be precise at 32*m_epsilon level"
+ $ assertBool (unlines (map show deviations)) (null deviations)
, testAssertion "incompleteBeta with p > 3000 and q > 3000"
$ and [ eq 1e-11 (incompleteBeta p q x) ib | (x,p,q,ib) <-
[ (0.495, 3001, 3001, 0.2192546757957825068677527085659175689142653854877723)
@@ -82,7 +92,7 @@ tests = testGroup "Special functions"
$ and [ let n' = fromIntegral n
k' = fromIntegral k
in within 2 (logChoose n' k') (log $ choose n' k')
- | n <- [0..1000], k <- [0..n]]
+ | n <- [0::Int .. 1000], k <- [0 .. n]]
----------------------------------------------------------------
-- Self tests
, testProperty "Self-test: 0 <= range01 <= 1" $ \x -> let f = range01 x in f <= 1 && f >= 0
diff --git a/tests/tests.hs b/tests/tests.hs
index 86ee021..b8aec23 100644
--- a/tests/tests.hs
+++ b/tests/tests.hs
@@ -1,8 +1,9 @@
import Test.Framework (defaultMain)
-import qualified Tests.SpecFunctions
import qualified Tests.Chebyshev
-import qualified Tests.Sum
import qualified Tests.Comparison
+import qualified Tests.RootFinding
+import qualified Tests.SpecFunctions
+import qualified Tests.Sum
main :: IO ()
main = defaultMain [ Tests.SpecFunctions.tests
@@ -10,4 +11,5 @@ main = defaultMain [ Tests.SpecFunctions.tests
-- , Tests.Chebyshev.tests
, Tests.Sum.tests
, Tests.Comparison.tests
+ , Tests.RootFinding.tests
]