diff --git a/ChangeLog b/ChangeLog new file mode 100644 --- /dev/null +++ b/ChangeLog @@ -0,0 +1,16 @@ +Changes in 0.1.4 + + * logFactorial type is genberalized. It accepts any `Integral' type + + * Evaluation of polynomials using Horner's method where coefficients + are store in lists added + +Changes in 0.1.3 + + * Error function and its inverse added. + + * Digamma function added + + * Evaluation of polynomials using Horner's method added. + + * Crash bug in the inverse incomplete beta fixed. diff --git a/Numeric/Polynomial.hs b/Numeric/Polynomial.hs --- a/Numeric/Polynomial.hs +++ b/Numeric/Polynomial.hs @@ -9,12 +9,19 @@ -- -- Function for evaluating polynomials using Horher's method. module Numeric.Polynomial ( + -- * Polynomials evaluatePolynomial , evaluateEvenPolynomial , evaluateOddPolynomial + -- ** Lists + -- $list + , evaluatePolynomialL + , evaluateEvenPolynomialL + , evaluateOddPolynomialL ) where import qualified Data.Vector.Generic as G +import qualified Data.Vector as V import Data.Vector.Generic (Vector) @@ -27,8 +34,9 @@ -> v a -- ^ Coefficients -> a {-# INLINE evaluatePolynomial #-} -evaluatePolynomial x coefs - = G.foldr (\a r -> a + r*x) 0 coefs +evaluatePolynomial x v + | G.null v = 0 + | otherwise = G.foldr1 (\a r -> a + r*x) v -- | Evaluate polynomial with only even powers using Horner's method. -- Coefficients starts from lowest. In pseudocode: @@ -39,10 +47,10 @@ -> v a -- ^ Coefficients -> a {-# INLINE evaluateEvenPolynomial #-} -evaluateEvenPolynomial x coefs - = G.foldr (\a r -> a + r*x2) 0 coefs - where x2 = x * x +evaluateEvenPolynomial x + = evaluatePolynomial (x*x) + -- | Evaluate polynomial with only odd powers using Horner's method. -- Coefficients starts from lowest. In pseudocode: -- @@ -53,5 +61,26 @@ -> a {-# INLINE evaluateOddPolynomial #-} evaluateOddPolynomial x coefs - = x * G.foldr (\a r -> a + r*x2) 0 coefs - where x2 = x * x + = x * evaluatePolynomial (x*x) coefs + + + + +-- $lists +-- +-- When all coefficients are known statically it's more convenient to +-- pass coefficient in a list instad of vector. Functions below +-- provide just that functionality. If list is known statically it +-- will be inlined anyway. + +evaluatePolynomialL :: (Num a) => a -> [a] -> a +evaluatePolynomialL x = evaluatePolynomial x . V.fromList +{-# INLINE evaluatePolynomialL #-} + +evaluateEvenPolynomialL :: (Num a) => a -> [a] -> a +evaluateEvenPolynomialL x = evaluateEvenPolynomial x . V.fromList +{-# INLINE evaluateEvenPolynomialL #-} + +evaluateOddPolynomialL :: (Num a) => a -> [a] -> a +evaluateOddPolynomialL x = evaluateOddPolynomial x . V.fromList +{-# INLINE evaluateOddPolynomialL #-} diff --git a/Numeric/SpecFunctions.hs b/Numeric/SpecFunctions.hs --- a/Numeric/SpecFunctions.hs +++ b/Numeric/SpecFunctions.hs @@ -45,7 +45,7 @@ import qualified Data.Vector.Unboxed as U import Numeric.Polynomial.Chebyshev (chebyshevBroucke) -import Numeric.Polynomial (evaluateEvenPolynomial) +import Numeric.Polynomial (evaluateEvenPolynomialL,evaluateOddPolynomialL) import Numeric.MathFunctions.Constants ( m_epsilon, m_NaN, m_neg_inf, m_pos_inf , m_sqrt_2_pi, m_ln_sqrt_2_pi, m_sqrt_2 , m_eulerMascheroni @@ -492,14 +492,14 @@ -- We cannot continue at this point so we simply return `x' | x == 0 || x == 1 = x -- When derivative becomes infinite we cannot continue - -- iterations. It cat only happen in vicinity of 0 or 1. It's + -- iterations. It can only happen in vicinity of 0 or 1. It's -- hardly possible to get good answer in such circumstances but -- `x' is already reasonable. | isInfinite f' = x -- Iterations limit reached. Most of the time solution will - -- converge to answer because of discetenes of Double. But + -- converge to answer because of discreteness of Double. But -- solution have good precision already. - | i >= 1000 = x + | i >= 10 = x -- Solution converges | abs dx <= 16 * m_epsilon * x = x' | otherwise = loop (i+1) x' @@ -509,7 +509,7 @@ f' = exp $ a1 * log x + b1 * log (1 - x) - beta u = f / f' dx = u / (1 - 0.5 * min 1 (u * (a1 / x - b1 / (1 - x)))) - -- Next approximation. If Halley step leas us out of [0,1] + -- Next approximation. If Halley step leads us out of [0,1] -- range we revert to bisection. x' | z < 0 = x / 2 | z > 1 = (x + 1) / 2 @@ -518,14 +518,14 @@ -- Calculate initial guess. Approximations from AS64, AS109 and -- Numerical recipes are used. -- - -- Equations are refered to by name of paper and number e.g. [AS64 2] + -- Equations are referred to by name of paper and number e.g. [AS64 2] -- In AS64 papers equations are not numbered so they are refered -- to by number of appearance starting from definition of -- incomplete beta. guess -- In this region we use approximation from AS109 (Carter -- approximation). It's reasonably good (2 iterations on - -- average) and never crashes. + -- average) | a > 1 && b > 1 = let r = (y*y - 3) / 6 s = 1 / (2*a - 1) @@ -644,14 +644,16 @@ -- | Compute the natural logarithm of the factorial function. Gives -- 16 decimal digits of precision. -logFactorial :: Int -> Double +logFactorial :: Integral a => a -> Double logFactorial n - | n <= 14 = log (factorial n) + | n < 0 = error "Numeric.SpecFunctions.logFactorial: negative input" + | n <= 14 = log $ factorial $ fromIntegral n | otherwise = (x - 0.5) * log x - x + 9.1893853320467e-1 + z / x - where x = fromIntegral (n + 1) + where x = fromIntegral n + 1 y = 1 / (x * x) z = ((-(5.95238095238e-4 * y) + 7.936500793651e-4) * y - 2.7777777777778e-3) * y + 8.3333333333333e-2 +{-# SPECIALIZE logFactorial :: Int -> Double #-} -- | Calculate the error term of the Stirling approximation. This is -- only defined for non-negative values. @@ -663,12 +665,11 @@ (i,0) -> sfe `U.unsafeIndex` i _ -> logGamma (n+1.0) - (n+0.5) * log n + n - m_ln_sqrt_2_pi - | n > 500 = (s0-s1/nn)/n - | n > 80 = (s0-(s1-s2/nn)/nn)/n - | n > 35 = (s0-(s1-(s2-s3/nn)/nn)/nn)/n - | otherwise = (s0-(s1-(s2-(s3-s4/nn)/nn)/nn)/nn)/n + | n > 500 = evaluateOddPolynomialL (1/n) [s0,-s1] + | n > 80 = evaluateOddPolynomialL (1/n) [s0,-s1,s2] + | n > 35 = evaluateOddPolynomialL (1/n) [s0,-s1,s2,-s3] + | otherwise = evaluateOddPolynomialL (1/n) [s0,-s1,s2,-s3,s4] where - nn = n*n s0 = 0.083333333333333333333 -- 1/12 s1 = 0.00277777777777777777778 -- 1/360 s2 = 0.00079365079365079365079365 -- 1/1260 @@ -746,15 +747,15 @@ | x' < c = r -- De Moivre's expansion | otherwise = let s = 1/x' - in evaluateEvenPolynomial s $ - U.fromList [ r + log x' - 0.5 * s - , - 1/12 - , 1/120 - , - 1/252 - , 1/240 - , - 1/132 - , 391/32760 - ] + in evaluateEvenPolynomialL s + [ r + log x' - 0.5 * s + , - 1/12 + , 1/120 + , - 1/252 + , 1/240 + , - 1/132 + , 391/32760 + ] where γ = m_eulerMascheroni c = 12 diff --git a/math-functions.cabal b/math-functions.cabal --- a/math-functions.cabal +++ b/math-functions.cabal @@ -1,5 +1,5 @@ name: math-functions -version: 0.1.3.0 +version: 0.1.4.0 cabal-version: >= 1.8 license: BSD3 license-file: LICENSE @@ -7,6 +7,7 @@ Aleksey Khudyakov maintainer: Bryan O'Sullivan 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 @@ -14,21 +15,13 @@ This library provides implementations of special mathematical functions and Chebyshev polynomials. These functions are often useful in statistical and numerical computing. - . - Changes in 0.1.2 - . - * Error function and its inverse added. - . - * Digamma function added - . - * Evaluation of polynomials using Horner's method added. - . - * Crash bug in the inverse incomplete beta fixed. + extra-source-files: README.markdown tests/*.hs tests/Tests/*.hs tests/Tests/SpecFunctions/gen.py + ChangeLog library ghc-options: -Wall diff --git a/tests/view.hs b/tests/view.hs new file mode 100644 --- /dev/null +++ b/tests/view.hs @@ -0,0 +1,102 @@ +{-# LANGUAGE OverloadedStrings #-} +import Control.Applicative +import Control.Monad +import Numeric.SpecFunctions +import Numeric.MathFunctions.Constants +import CPython.Sugar +import CPython.MPMath +import qualified CPython as Py + +import HEP.ROOT.Plot + + +---------------------------------------------------------------- + + +viewBetaDelta = runPy $ do + addToPythonPath "." + m <- loadMPMath + mpmSetDps m 100 + xs <- forM pqBeta $ \(p,q) -> do x <- fromMPNum =<< mpmLog m =<< mpmBeta m (MPDouble p) (MPDouble q) + return (p,q, relErr x (logBeta p q)) + draws $ do + -- let xs = [ (p,q, logBeta p q `relErr` (logGammaL p + logGammaL q - logGammaL (q+p))) + -- | (p,q) <- pqBeta + -- ] + add $ Graph2D xs + + +pqBeta = [ (p,q) + | p <- logRange 50 0.3 0.6 + , q <- logRange 50 5 6 + ] + where + + + + +viewIBeta x = runPy $ do + addToPythonPath "." + m <- loadMPMath + mpmSetDps m 30 + -- + let n = 40 + let pq = (,) + <$> logRange n 100 1000 + <*> logRange n 100 1000 + -- + xs <- forM pq $ \(p,q) -> do + i <- fromMPNum =<< mpmIncompleteBeta m (MPDouble p) (MPDouble q) (MPDouble x) + return (p,q, incompleteBeta p q x `relErr` i) + -- + draws $ do + add $ Graph2D xs + + +go = runPy $ do + addToPythonPath "." + m <- loadMPMath + mpmSetDps m 16 + -- + print =<< fromMPNum =<< mpmIncompleteBeta m (MPDouble 10) (MPDouble 10) (MPDouble 0.4) + print $ incompleteBeta 10 10 0.4 + + + + +viewLancrox = runPy $ do + addToPythonPath "." + m <- loadMPMath + mpmSetDps m 50 + -- + let xs = logRange 10000 (1e-8) (1e-1) + pl <- forM xs $ \x -> do y0 <- fromMPNum =<< mpmLog m =<< mpmGamma m (MPDouble x) + return (x, y0) + draws $ do + add $ Graph $ [ (x, abs $ y `relErr` logGammaL x) | (x,y) <- pl ] + set $ lineColor RED + -- + add $ Graph $ [ (x, abs $ y `relErr` logGamma x) | (x,y) <- pl ] + set $ lineColor BLUE + -- + set $ xaxis $ logScale ON + -- set $ yaxis $ logScale ON + -- + add $ HLine m_epsilon + add $ HLine $ negate m_epsilon + + +---------------------------------------------------------------- + +relErr :: Double -> Double -> Double +relErr 0 0 = 0 +relErr x y = (x - y) / max (abs x) (abs y) + + + +logRange :: Int -> Double -> Double -> [Double] +logRange n a b + = [ a * r^i | i <- [0 .. n] ] + where + r = (b / a) ** (1 / fromIntegral n) +