diff --git a/Numeric/RootFinding.hs b/Numeric/RootFinding.hs --- 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 @@ | 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= (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–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 --- a/Numeric/SpecFunctions/Internal.hs +++ b/Numeric/SpecFunctions/Internal.hs @@ -14,9 +14,10 @@ #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 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 @@ 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 --- 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.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 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 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) @@ -168,6 +190,16 @@ 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' diff --git a/README.markdown b/README.markdown --- 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 --- 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 --- 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 --- 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 --- 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 , - Aleksey Khudyakov -maintainer: Bryan O'Sullivan + Alexey Khudyakov +maintainer: Alexey Khudyakov 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 @@ 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 @@ 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 --- /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 --- a/tests/Tests/SpecFunctions.hs +++ b/tests/Tests/SpecFunctions.hs @@ -10,6 +10,8 @@ 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 @@ -- 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 @@ $ 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 --- 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 @@ -- , Tests.Chebyshev.tests , Tests.Sum.tests , Tests.Comparison.tests + , Tests.RootFinding.tests ]