math-functions 0.1.3.0 → 0.1.4.0
raw patch · 5 files changed
+183/−42 lines, 5 filesnew-uploader
Files
- ChangeLog +16/−0
- Numeric/Polynomial.hs +36/−7
- Numeric/SpecFunctions.hs +25/−24
- math-functions.cabal +4/−11
- tests/view.hs +102/−0
+ ChangeLog view
@@ -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.
Numeric/Polynomial.hs view
@@ -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 #-}
Numeric/SpecFunctions.hs view
@@ -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
math-functions.cabal view
@@ -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 <alexey.skladnoy@gmail.com> maintainer: Bryan O'Sullivan <bos@serpentine.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@@ -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
+ tests/view.hs view
@@ -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)+