math-functions 0.2.1.0 → 0.3.0.0

raw patch · 11 files changed

+493/−237 lines, 11 filesdep +data-default-classdep ~vector-th-unboxPVP ok

version bump matches the API change (PVP)

Dependencies added: data-default-class

Dependency ranges changed: vector-th-unbox

API changes (from Hackage documentation)

Files

• Numeric/RootFinding.hs +336/−93
• Numeric/SpecFunctions/Internal.hs +6/−5
• Numeric/Sum.hs +34/−2
• README.markdown +17/−3
• benchmark/Summation.hs +0/−15
• benchmark/bench.hs +0/−89
• changelog.md +11/−1
• math-functions.cabal +28/−24
• tests/Tests/RootFinding.hs +44/−0
• tests/Tests/SpecFunctions.hs +13/−3
• tests/tests.hs +4/−2

@@ -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<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.

@@ -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

@@ -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'

@@ -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.

@@ -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-    ]

@@ -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 ]-  ]

@@ -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

@@ -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 @@     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

@@ -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))

@@ -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

@@ -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                    ]