math-functions 0.1.1.2 → 0.1.3.0
raw patch · 9 files changed
+673/−138 lines, 9 files
Files
- Numeric/MathFunctions/Constants.hs +36/−18
- Numeric/Polynomial.hs +57/−0
- Numeric/Polynomial/Chebyshev.hs +4/−2
- Numeric/SpecFunctions.hs +315/−73
- math-functions.cabal +13/−1
- tests/Tests/Chebyshev.hs +39/−15
- tests/Tests/SpecFunctions.hs +67/−9
- tests/Tests/SpecFunctions/Tables.hs +122/−19
- tests/Tests/SpecFunctions/gen.py +20/−1
Numeric/MathFunctions/Constants.hs view
@@ -11,20 +11,27 @@ module Numeric.MathFunctions.Constants (+ -- * IEE754 constants m_epsilon , m_huge , m_tiny+ , m_max_exp+ , m_pos_inf+ , m_neg_inf+ , m_NaN+ -- * Mathematical constants , m_1_sqrt_2 , m_2_sqrt_pi , m_ln_sqrt_2_pi- , m_max_exp , m_sqrt_2 , m_sqrt_2_pi- , m_pos_inf- , m_neg_inf- , m_NaN+ , m_eulerMascheroni ) where +----------------------------------------------------------------+-- IEE754 constants+----------------------------------------------------------------+ -- | A very large number. m_huge :: Double m_huge = 1.7976931348623157e308@@ -39,6 +46,27 @@ m_max_exp :: Int m_max_exp = 1024 +-- | Positive infinity.+m_pos_inf :: Double+m_pos_inf = 1/0+{-# INLINE m_pos_inf #-}++-- | Negative infinity.+m_neg_inf :: Double+m_neg_inf = -1/0+{-# INLINE m_neg_inf #-}++-- | Not a number.+m_NaN :: Double+m_NaN = 0/0+{-# INLINE m_NaN #-}++++----------------------------------------------------------------+-- Mathematical constants+----------------------------------------------------------------+ -- | @sqrt 2@ m_sqrt_2 :: Double m_sqrt_2 = 1.4142135623730950488016887242096980785696718753769480731766@@ -69,17 +97,7 @@ m_ln_sqrt_2_pi = 0.9189385332046727417803297364056176398613974736377834128171 {-# INLINE m_ln_sqrt_2_pi #-} --- | Positive infinity.-m_pos_inf :: Double-m_pos_inf = 1/0-{-# INLINE m_pos_inf #-}---- | Negative infinity.-m_neg_inf :: Double-m_neg_inf = -1/0-{-# INLINE m_neg_inf #-}---- | Not a number.-m_NaN :: Double-m_NaN = 0/0-{-# INLINE m_NaN #-}+-- | Euler–Mascheroni constant (γ = 0.57721...)+m_eulerMascheroni :: Double+m_eulerMascheroni = 0.5772156649015328606065121+{-# INLINE m_eulerMascheroni #-}
+ Numeric/Polynomial.hs view
@@ -0,0 +1,57 @@+-- |+-- Module : Numeric.Polynomial+-- Copyright : (c) 2012 Aleksey Khudyakov+-- License : BSD3+--+-- Maintainer : bos@serpentine.com+-- Stability : experimental+-- Portability : portable+--+-- Function for evaluating polynomials using Horher's method.+module Numeric.Polynomial (+ evaluatePolynomial+ , evaluateEvenPolynomial+ , evaluateOddPolynomial+ ) where++import qualified Data.Vector.Generic as G+import Data.Vector.Generic (Vector)+++-- | Evaluate polynomial using Horner's method. Coefficients starts+-- from lowest. In pseudocode:+--+-- > evaluateOddPolynomial x [1,2,3] = 1 + 2*x + 3*x^2+evaluatePolynomial :: (Vector v a, Num a)+ => a -- ^ /x/+ -> v a -- ^ Coefficients+ -> a+{-# INLINE evaluatePolynomial #-}+evaluatePolynomial x coefs+ = G.foldr (\a r -> a + r*x) 0 coefs++-- | Evaluate polynomial with only even powers using Horner's method.+-- Coefficients starts from lowest. In pseudocode:+--+-- > evaluateOddPolynomial x [1,2,3] = 1 + 2*x^2 + 3*x^4+evaluateEvenPolynomial :: (Vector v a, Num a)+ => a -- ^ /x/+ -> v a -- ^ Coefficients+ -> a+{-# INLINE evaluateEvenPolynomial #-}+evaluateEvenPolynomial x coefs+ = G.foldr (\a r -> a + r*x2) 0 coefs+ where x2 = x * x++-- | Evaluate polynomial with only odd powers using Horner's method.+-- Coefficients starts from lowest. In pseudocode:+--+-- > evaluateOddPolynomial x [1,2,3] = 1*x + 2*x^3 + 3*x^5+evaluateOddPolynomial :: (Vector v a, Num a)+ => a -- ^ /x/+ -> v a -- ^ Coefficients+ -> a+{-# INLINE evaluateOddPolynomial #-}+evaluateOddPolynomial x coefs+ = x * G.foldr (\a r -> a + r*x2) 0 coefs+ where x2 = x * x
Numeric/Polynomial/Chebyshev.hs view
@@ -48,8 +48,10 @@ data B = B {-# UNPACK #-} !Double {-# UNPACK #-} !Double {-# UNPACK #-} !Double -- | Evaluate a Chebyshev polynomial of the first kind. Uses Broucke's--- ECHEB algorithm, and his convention for coefficient handling, and so--- gives different results than 'chebyshev' for the same inputs.+-- ECHEB algorithm, and his convention for coefficient handling. It+-- treat 0th coefficient different so+--+-- > chebyshev x [a0,a1,a2...] == chebyshevBroucke [2*a0,a1,a2...] chebyshevBroucke :: (G.Vector v Double) => Double -- ^ Parameter of each function. -> v Double -- ^ Coefficients of each polynomial term, in increasing order.
Numeric/SpecFunctions.hs view
@@ -10,11 +10,17 @@ -- -- Special functions and factorials. module Numeric.SpecFunctions (+ -- * Error function+ erf+ , erfc+ , invErf+ , invErfc -- * Gamma function- logGamma+ , logGamma , logGammaL , incompleteGamma , invIncompleteGamma+ , digamma -- * Beta function , logBeta , incompleteBeta@@ -35,34 +41,91 @@ import Data.Bits ((.&.), (.|.), shiftR) import Data.Int (Int64)-import Data.Word (Word64)-import Data.Number.Erf (erfc)+import qualified Data.Number.Erf as Erf (erfc,erf) import qualified Data.Vector.Unboxed as U import Numeric.Polynomial.Chebyshev (chebyshevBroucke)-import Numeric.MathFunctions.Constants (m_epsilon, m_sqrt_2_pi, m_ln_sqrt_2_pi, - m_NaN, m_neg_inf, m_pos_inf, m_sqrt_2)+import Numeric.Polynomial (evaluateEvenPolynomial)+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+ )+import Text.Printf +----------------------------------------------------------------+-- Error function+---------------------------------------------------------------- +-- | Error function.+--+-- > erf -∞ = -1+-- > erf 0 = 0+-- > erf +∞ = 1+erf :: Double -> Double+{-# INLINE erf #-}+erf = Erf.erf++-- | Complementary error function.+--+-- > erfc -∞ = 2+-- > erfc 0 = 1+-- > errc +∞ = 0+erfc :: Double -> Double+{-# INLINE erfc #-}+erfc = Erf.erfc+++-- | Inverse of 'erf'.+invErf :: Double -- ^ /p/ ∈ [-1,1]+ -> Double+invErf p = invErfc (1 - p)++-- | Inverse of 'erfc'.+invErfc :: Double -- ^ /p/ ∈ [0,2]+ -> Double+invErfc p+ | p == 2 = m_neg_inf+ | p == 0 = m_pos_inf+ | p >0 && p < 2 = if p <= 1 then r else -r+ | otherwise = modErr $ "invErfc: p must be in [0,2] got " ++ show p+ where+ pp = if p <= 1 then p else 2 - p+ t = sqrt $ -2 * log( 0.5 * pp)+ -- Initial guess+ x0 = -0.70711 * ((2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t)+ r = loop 0 x0+ --+ loop :: Int -> Double -> Double+ loop !j !x+ | j >= 2 = x+ | otherwise = let err = erfc x - pp+ x' = x + err / (1.12837916709551257 * exp(-x * x) - x * err) -- // Halley+ in loop (j+1) x'+++ ---------------------------------------------------------------- -- Gamma function ---------------------------------------------------------------- -- Adapted from http://people.sc.fsu.edu/~burkardt/f_src/asa245/asa245.html --- | Compute the logarithm of the gamma function Γ(/x/). Uses+-- | Compute the logarithm of the gamma function Γ(/x/). Uses -- Algorithm AS 245 by Macleod. ----- Gives an accuracy of 10–12 significant decimal digits, except+-- Gives an accuracy of 10-12 significant decimal digits, except -- for small regions around /x/ = 1 and /x/ = 2, where the function -- goes to zero. For greater accuracy, use 'logGammaL'. ----- Returns ∞ if the input is outside of the range (0 < /x/--- ≤ 1e305).+-- Returns ∞ if the input is outside of the range (0 < /x/ ≤ 1e305). logGamma :: Double -> Double logGamma x | x <= 0 = m_pos_inf+ -- Handle positive infinity. logGamma overflows before 1e308 so+ -- it's safe+ | x > 1e308 = m_pos_inf+ -- Normal cases | x < 1.5 = a + c * ((((r1_4 * b + r1_3) * b + r1_2) * b + r1_1) * b + r1_0) / ((((b + r1_8) * b + r1_7) * b + r1_6) * b + r1_5)@@ -119,6 +182,8 @@ logGammaL :: Double -> Double logGammaL x | x <= 0 = m_pos_inf+ -- Lanroz approximation loses precision for small arguments+ | x <= 1e-3 = logGamma x | otherwise = fini . U.foldl' go (L 0 (x+7)) $ a where fini (L l _) = log (l+a0) + log m_sqrt_2_pi - x65 + (x-0.5) * log x65 go (L l t) k = L (l + k / t) (t-1)@@ -162,16 +227,19 @@ -- | Compute the normalized lower incomplete gamma function--- γ(/s/,/x/). Normalization means that--- γ(/s/,∞)=1. Uses Algorithm AS 239 by Shea.-incompleteGamma :: Double -- ^ /s/- -> Double -- ^ /x/+-- γ(/s/,/x/). Normalization means that+-- γ(/s/,∞)=1. Uses Algorithm AS 239 by Shea.+incompleteGamma :: Double -- ^ /s/ ∈ (0,∞)+ -> Double -- ^ /x/ ∈ (0,∞) -> Double incompleteGamma p x | isNaN p || isNaN x = m_NaN | x < 0 || p <= 0 = m_pos_inf | x == 0 = 0- | p >= 1000 = norm (3 * sqrt p * ((x/p) ** (1/3) + 1/(9*p) - 1))+ -- For very large `p' normal approximation gives <1e-10 error+ | p >= 2e5 = norm (3 * sqrt p * ((x/p) ** (1/3) + 1/(9*p) - 1))+ | p >= 500 = approx+ -- Dubious approximation | x >= 1e8 = 1 | x <= 1 || x < p = let a = p * log x - x - logGamma (p + 1) g = a + log (pearson p 1 1)@@ -179,7 +247,26 @@ | otherwise = let g = p * log x - x - logGamma p + log cf in if g > limit then 1 - exp g else 1 where+ -- CDF for standard normal distributions norm a = 0.5 * erfc (- a / m_sqrt_2)+ -- For large values of `p' we use 18-point Gauss-Legendre+ -- integration.+ approx+ | ans > 0 = 1 - ans+ | otherwise = -ans+ where+ -- Set upper limit for integration+ xu | x > p1 = (p1 + 11.5*sqrtP1) `max` (x + 6*sqrtP1)+ | otherwise = max 0 $ (p1 - 7.5*sqrtP1) `min` (x - 5*sqrtP1)+ s = U.sum $ U.zipWith go coefY coefW+ go y w = let t = x + (xu - x)*y+ in w * exp( -(t-p1) + p1*(log t - lnP1) )+ ans = s * (xu - x) * exp( p1 * (lnP1 - 1) - logGamma p)+ --+ p1 = p - 1+ lnP1 = log p1+ sqrtP1 = sqrt p1+ -- pearson !a !c !g | c' <= tolerance = g' | otherwise = pearson a' c' g'@@ -216,15 +303,14 @@ -- approximately holds: -- -- > invIncompleteGamma s . incompleteGamma s = id------ For @invIncompleteGamma s p@ /s/ must be positive and /p/ must be--- in [0,1] range.-invIncompleteGamma :: Double -> Double -> Double+invIncompleteGamma :: Double -- ^ /s/ ∈ (0,∞)+ -> Double -- ^ /p/ ∈ [0,1]+ -> Double invIncompleteGamma a p- | a <= 0 = - error $ "Statistics.Math.invIncompleteGamma: a must be positive. Got: " ++ show a- | p < 0 || p > 1 = - error $ "Statistics.Math.invIncompleteGamma: p must be in [0,1] range. Got: " ++ show p+ | a <= 0 =+ modErr $ printf "invIncompleteGamma: a must be positive. a=%g p=%g" a p+ | p < 0 || p > 1 =+ modErr $ printf "invIncompleteGamma: p must be in [0,1] range. a=%g p=%g" a p | p == 0 = 0 | p == 1 = 1 / 0 | otherwise = loop 0 guess@@ -253,7 +339,7 @@ | otherwise = x - dx -- Calculate inital guess for root guess- -- + -- | a > 1 = let t = sqrt $ -2 * log(if p < 0.5 then p else 1 - p) x1 = (2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t@@ -299,7 +385,8 @@ c = logGammaCorrection q - logGammaCorrection pq -- | Regularized incomplete beta function. Uses algorithm AS63 by--- Majumder and Bhattachrjee.+-- Majumder and Bhattachrjee and quadrature approximation for large+-- /p/ and /q/. incompleteBeta :: Double -- ^ /p/ > 0 -> Double -- ^ /q/ > 0 -> Double -- ^ /x/, must lie in [0,1] range@@ -308,22 +395,54 @@ -- | Regularized incomplete beta function. Same as 'incompleteBeta' -- but also takes logarithm of beta function as parameter.-incompleteBeta_ :: Double -- ^ logarithm of beta function+incompleteBeta_ :: Double -- ^ logarithm of beta function for given /p/ and /q/ -> Double -- ^ /p/ > 0 -> Double -- ^ /q/ > 0 -> Double -- ^ /x/, must lie in [0,1] range -> Double incompleteBeta_ beta p q x- | p <= 0 || q <= 0 = error "p <= 0 || q <= 0"- | x < 0 || x > 1 || isNaN x = error "x out of [0,1] range"+ | p <= 0 || q <= 0 =+ modErr $ printf "incompleteBeta_: p <= 0 || q <= 0. p=%g q=%g x=%g" p q x+ | x < 0 || x > 1 || isNaN x =+ modErr $ printf "incompletBeta_: x out of [0,1] range. p=%g q=%g x=%g" p q x | x == 0 || x == 1 = x | p >= (p+q) * x = incompleteBetaWorker beta p q x | otherwise = 1 - incompleteBetaWorker beta q p (1 - x) ++-- Approximation of incomplete beta by quandrature.+--+-- Note that x =< p/(p+q)+incompleteBetaApprox :: Double -> Double -> Double -> Double -> Double+incompleteBetaApprox beta p q x+ | ans > 0 = 1 - ans+ | otherwise = -ans+ where+ -- Constants+ p1 = p - 1+ q1 = q - 1+ mu = p / (p + q)+ lnmu = log mu+ lnmuc = log (1 - mu)+ -- Upper limit for integration+ xu = max 0 $ min (mu - 10*t) (x - 5*t)+ where+ t = sqrt $ p*q / ( (p+q) * (p+q) * (p + q + 1) )+ -- Calculate incomplete beta by quadrature+ go y w = let t = x + (xu - x) * y+ in w * exp( p1 * (log t - lnmu) + q1 * (log(1-t) - lnmuc) )+ s = U.sum $ U.zipWith go coefY coefW+ ans = s * (xu - x) * exp( p1 * lnmu + q1 * lnmuc - beta )++ -- Worker for incomplete beta function. It is separate function to -- avoid confusion with parameter during parameter swapping incompleteBetaWorker :: Double -> Double -> Double -> Double -> Double-incompleteBetaWorker beta p q x = loop (p+q) (truncate $ q + cx * (p+q)) 1 1 1+incompleteBetaWorker beta p q x+ -- For very large p and q this method becomes very slow so another+ -- method is used.+ | p > 3000 && q > 3000 = incompleteBetaApprox beta p q x+ | otherwise = loop (p+q) (truncate $ q + cx * (p+q)) 1 1 1 where -- Constants eps = 1e-15@@ -346,60 +465,109 @@ -- | Compute inverse of regularized incomplete beta function. Uses--- initial approximation from AS109 and Halley method to solve equation.-invIncompleteBeta :: Double -- ^ /p/- -> Double -- ^ /q/- -> Double -- ^ /a/+-- initial approximation from AS109, AS64 and Halley method to solve+-- equation.+invIncompleteBeta :: Double -- ^ /p/ > 0+ -> Double -- ^ /q/ > 0+ -> Double -- ^ /a/ ∈ [0,1] -> Double invIncompleteBeta p q a- | p <= 0 || q <= 0 = error "p <= 0 || q <= 0"- | a < 0 || a > 1 = error "bad a"+ | p <= 0 || q <= 0 =+ modErr $ printf "invIncompleteBeta p <= 0 || q <= 0. p=%g q=%g a=%g" p q a+ | a < 0 || a > 1 =+ modErr $ printf "invIncompleteBeta x must be in [0,1]. p=%g q=%g a=%g" p q a | a == 0 || a == 1 = a | a > 0.5 = 1 - invIncompleteBetaWorker (logBeta p q) q p (1 - a)- | otherwise = invIncompleteBetaWorker (logBeta p q) p q a+ | otherwise = invIncompleteBetaWorker (logBeta p q) p q a + invIncompleteBetaWorker :: Double -> Double -> Double -> Double -> Double-invIncompleteBetaWorker beta p q a = loop (0::Int) guess+-- NOTE: p <= 0.5.+invIncompleteBetaWorker beta a b p = loop (0::Int) guess where- p1 = p - 1- q1 = q - 1+ a1 = a - 1+ b1 = b - 1 -- Solve equation using Halley method loop !i !x+ -- We cannot continue at this point so we simply return `x' | x == 0 || x == 1 = x- | i >= 10 = x- | abs dx <= 16 * m_epsilon * x = x+ -- When derivative becomes infinite we cannot continue+ -- iterations. It cat 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+ -- solution have good precision already.+ | i >= 1000 = x+ -- Solution converges+ | abs dx <= 16 * m_epsilon * x = x' | otherwise = loop (i+1) x' where- f = incompleteBeta_ beta p q x - a- f' = exp $ p1 * log x + q1 * log (1 - x) - beta+ -- Calculate Halley step.+ f = incompleteBeta_ beta a b x - p+ f' = exp $ a1 * log x + b1 * log (1 - x) - beta u = f / f'- dx = u / (1 - 0.5 * min 1 (u * (p1 / x - q1 / (1 - x))))+ dx = u / (1 - 0.5 * min 1 (u * (a1 / x - b1 / (1 - x))))+ -- Next approximation. If Halley step leas us out of [0,1]+ -- range we revert to bisection. x' | z < 0 = x / 2 | z > 1 = (x + 1) / 2 | otherwise = z where z = x - dx- -- Calculate initial guess- guess - | p > 1 && q > 1 = - let rr = (y*y - 3) / 6- ss = 1 / (2*p - 1)- tt = 1 / (2*q - 1)- hh = 2 / (ss + tt)- ww = y * sqrt(hh + rr) / hh - (tt - ss) * (rr + 5/6 - 2 / (3 * hh))- in p / (p + q * exp(2 * ww))- | t' <= 0 = 1 - exp( (log((1 - a) * q) + beta) / q )- | t'' <= 1 = exp( (log(a * p) + beta) / p )- | otherwise = 1 - 2 / (t'' + 1)+ -- 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]+ -- 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.+ | a > 1 && b > 1 =+ let r = (y*y - 3) / 6+ s = 1 / (2*a - 1)+ t = 1 / (2*b - 1)+ h = 2 / (s + t)+ w = y * sqrt(h + r) / h - (t - s) * (r + 5/6 - 2 / (3 * h))+ in a / (a + b * exp(2 * w))+ -- Otherwise we revert to approximation from AS64 derived from+ -- [AS64 2] when it's applicable.+ --+ -- It slightly reduces average number of iterations when `a' and+ -- `b' have different magnitudes.+ | chi2 > 0 && ratio > 1 = 1 - 2 / (ratio + 1)+ -- If all else fails we use approximation from "Numerical+ -- Recipes". It's very similar to approximations [AS64 4,5] but+ -- it never goes out of [0,1] interval.+ | otherwise = case () of+ _| p < t / w -> (a * p * w) ** (1/a)+ | otherwise -> 1 - (b * (1 - p) * w) ** (1/b)+ where+ lna = log $ a / (a+b)+ lnb = log $ b / (a+b)+ t = exp( a * lna ) / a+ u = exp( b * lnb ) / b+ w = t + u where- r = sqrt ( - log ( a * a ) )+ -- Formula [2]+ ratio = (4*a + 2*b - 2) / chi2+ -- Quantile of chi-squared distribution. Formula [3].+ chi2 = 2 * b * (1 - t + y * sqrt t) ** 3+ where+ t = 1 / (9 * b)+ -- `y' is Hasting's approximation of p'th quantile of standard+ -- normal distribution. y = r - ( 2.30753 + 0.27061 * r )- / ( 1.0 + ( 0.99229 + 0.04481 * r ) * r )- t = 1 / (9 * q)- t' = 2 * q * (1 - t + y * sqrt t) ** 3- t'' = (4*p + 2*q - 2) / t'+ / ( 1.0 + ( 0.99229 + 0.04481 * r ) * r )+ where+ r = sqrt $ - 2 * log p + ---------------------------------------------------------------- -- Logarithm ----------------------------------------------------------------@@ -448,7 +616,7 @@ -- | /O(log n)/ Compute the logarithm in base 2 of the given value. log2 :: Int -> Int log2 v0- | v0 <= 0 = error "Statistics.Math.log2: invalid input"+ | v0 <= 0 = modErr $ "log2: negative input, got " ++ show v0 | otherwise = go 5 0 v0 where go !i !r !v | i == -1 = r@@ -464,12 +632,12 @@ -- Factorial ---------------------------------------------------------------- --- | Compute the factorial function /n/!. Returns ∞ if the+-- | Compute the factorial function /n/!. Returns +∞ if the -- input is above 170 (above which the result cannot be represented by -- a 64-bit 'Double'). factorial :: Int -> Double factorial n- | n < 0 = error "Statistics.Math.factorial: negative input"+ | n < 0 = error "Numeric.SpecFunctions.factorial: negative input" | n <= 1 = 1 | n <= 170 = U.product $ U.map fromIntegral $ U.enumFromTo 2 n | otherwise = m_pos_inf@@ -559,11 +727,92 @@ max64 = fromIntegral (maxBound :: Int64) round64 x = round x :: Int64 +-- | Compute ψ0(/x/), the first logarithmic derivative of the gamma+-- function. Uses Algorithm AS 103 by Bernardo, based on Minka's C+-- implementation.+digamma :: Double -> Double+digamma x+ | isNaN x || isInfinite x = m_NaN+ -- FIXME:+ -- This is ugly. We are testing here that number is in fact+ -- integer. It's somewhat tricky question to answer. When ε for+ -- given number becomes 1 or greater every number is represents+ -- an integer. We also must make sure that excess precision+ -- won't bite us.+ | x <= 0 && fromIntegral (truncate x :: Int64) == x = m_neg_inf+ -- Jeffery's reflection formula+ | x < 0 = digamma (1 - x) + pi / tan (negate pi * x)+ | x <= 1e-6 = - γ - 1/x + trigamma1 * x+ | 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+ ]+ where+ γ = m_eulerMascheroni+ c = 12+ -- Reduce to digamma (x + n) where (x + n) >= c+ (r, x') = reduce 0 x+ where+ reduce !s y+ | y < c = reduce (s - 1 / y) (y + 1)+ | otherwise = (s, y) +----------------------------------------------------------------+-- Constants+----------------------------------------------------------------++-- Coefficients for 18-point Gauss-Legendre integration. They are+-- used in implementation of incomplete gamma and beta functions.+coefW,coefY :: U.Vector Double+coefW = U.fromList [ 0.0055657196642445571, 0.012915947284065419, 0.020181515297735382+ , 0.027298621498568734, 0.034213810770299537, 0.040875750923643261+ , 0.047235083490265582, 0.053244713977759692, 0.058860144245324798+ , 0.064039797355015485, 0.068745323835736408, 0.072941885005653087+ , 0.076598410645870640, 0.079687828912071670, 0.082187266704339706+ , 0.084078218979661945, 0.085346685739338721, 0.085983275670394821+ ]+coefY = U.fromList [ 0.0021695375159141994, 0.011413521097787704, 0.027972308950302116+ , 0.051727015600492421, 0.082502225484340941, 0.12007019910960293+ , 0.16415283300752470, 0.21442376986779355, 0.27051082840644336+ , 0.33199876341447887, 0.39843234186401943, 0.46931971407375483+ , 0.54413605556657973, 0.62232745288031077, 0.70331500465597174+ , 0.78649910768313447, 0.87126389619061517, 0.95698180152629142+ ]+{-# NOINLINE coefW #-}+{-# NOINLINE coefY #-}++trigamma1 :: Double+trigamma1 = 1.6449340668482264365 -- pi**2 / 6++modErr :: String -> a+modErr msg = error $ "Numeric.SpecFunctions." ++ msg+++ -- $references --+-- * Bernardo, J. (1976) Algorithm AS 103: Psi (digamma)+-- function. /Journal of the Royal Statistical Society. Series C+-- (Applied Statistics)/ 25(3):315-317.+-- <http://www.jstor.org/stable/2347257>+--+-- * Cran, G.W., Martin, K.J., Thomas, G.E. (1977) Remark AS R19+-- and Algorithm AS 109: A Remark on Algorithms: AS 63: The+-- Incomplete Beta Integral AS 64: Inverse of the Incomplete Beta+-- Function Ratio. /Journal of the Royal Statistical Society. Series+-- C (Applied Statistics)/ Vol. 26, No. 1 (1977), pp. 111-114+-- <http://www.jstor.org/pss/2346887>+-- -- * Lanczos, C. (1964) A precision approximation of the gamma -- function. /SIAM Journal on Numerical Analysis B/ -- 1:86–96. <http://www.jstor.org/stable/2949767>@@ -576,10 +825,6 @@ -- /Journal of the Royal Statistical Society, Series C (Applied Statistics)/ -- 38(2):397–402. <http://www.jstor.org/stable/2348078> ----- * Shea, B. (1988) Algorithm AS 239: Chi-squared and incomplete--- gamma integral. /Applied Statistics/--- 37(3):466–473. <http://www.jstor.org/stable/2347328>--- -- * Majumder, K.L., Bhattacharjee, G.P. (1973) Algorithm AS 63: The -- Incomplete Beta Integral. /Journal of the Royal Statistical -- Society. Series C (Applied Statistics)/ Vol. 22, No. 3 (1973),@@ -591,9 +836,6 @@ -- Vol. 22, No. 3 (1973), pp. 411-414 -- <http://www.jstor.org/pss/2346798> ----- * Cran, G.W., Martin, K.J., Thomas, G.E. (1977) Remark AS R19--- and Algorithm AS 109: A Remark on Algorithms: AS 63: The--- Incomplete Beta Integral AS 64: Inverse of the Incomplete Beta--- Function Ratio. /Journal of the Royal Statistical Society. Series--- C (Applied Statistics)/ Vol. 26, No. 1 (1977), pp. 111-114--- <http://www.jstor.org/pss/2346887>+-- * Shea, B. (1988) Algorithm AS 239: Chi-squared and incomplete+-- gamma integral. /Applied Statistics/+-- 37(3):466–473. <http://www.jstor.org/stable/2347328>
math-functions.cabal view
@@ -1,5 +1,5 @@ name: math-functions-version: 0.1.1.2+version: 0.1.3.0 cabal-version: >= 1.8 license: BSD3 license-file: LICENSE@@ -14,6 +14,16 @@ 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@@ -21,12 +31,14 @@ tests/Tests/SpecFunctions/gen.py library+ ghc-options: -Wall build-depends: base >=3 && <5, vector >= 0.7, erf >= 2 exposed-modules: Numeric.SpecFunctions Numeric.SpecFunctions.Extra+ Numeric.Polynomial Numeric.Polynomial.Chebyshev Numeric.MathFunctions.Constants
tests/Tests/Chebyshev.hs view
@@ -5,7 +5,7 @@ import Data.Vector.Unboxed (fromList) import Test.Framework import Test.Framework.Providers.QuickCheck2-import Test.QuickCheck (Arbitrary(..))+import Test.QuickCheck (Arbitrary(..),printTestCase,Property) import Tests.Helpers import Numeric.Polynomial.Chebyshev@@ -14,31 +14,55 @@ tests :: Test tests = testGroup "Chebyshev polynomials" [ testProperty "Chebyshev 0" $ \a0 (Ch x) ->- (ch0 x * a0) ≈ (chebyshev x $ fromList [a0])+ testCheb [a0] x , testProperty "Chebyshev 1" $ \a0 a1 (Ch x) ->- (a0*ch0 x + a1*ch1 x) ≈ (chebyshev x $ fromList [a0,a1])+ testCheb [a0,a1] x , testProperty "Chebyshev 2" $ \a0 a1 a2 (Ch x) ->- (a0*ch0 x + a1*ch1 x + a2*ch2 x) ≈ (chebyshev x $ fromList [a0,a1,a2])+ testCheb [a0,a1,a2] x , testProperty "Chebyshev 3" $ \a0 a1 a2 a3 (Ch x) ->- (a0*ch0 x + a1*ch1 x + a2*ch2 x + a3*ch3 x) ≈ (chebyshev x $ fromList [a0,a1,a2,a3])+ testCheb [a0,a1,a2,a3] x , testProperty "Chebyshev 4" $ \a0 a1 a2 a3 a4 (Ch x) ->- (a0*ch0 x + a1*ch1 x + a2*ch2 x + a3*ch3 x + a4*ch4 x) ≈ (chebyshev x $ fromList [a0,a1,a2,a3,a4])+ testCheb [a0,a1,a2,a3,a4] x+ , testProperty "Broucke" $ testBroucke ]- where (≈) = eq 1e-12+ where +testBroucke _ [] = True+testBroucke (Ch x) (c:cs) = let c1 = chebyshev x (fromList $ c : cs)+ cb = chebyshevBroucke x (fromList $ c*2 : cs)+ in eq 1e-15 c1 cb --- Chebyshev polynomials of low order-ch0,ch1,ch2,ch3,ch4 :: Double -> Double-ch0 _ = 1-ch1 x = x-ch2 x = 2*x^2 - 1-ch3 x = 4*x^3 - 3*x-ch4 x = 8*x^4 - 8*x^2 + 1+testCheb :: [Double] -> Double -> Property+testCheb as x+ = printTestCase (">>> Exact = " ++ show exact)+ $ printTestCase (">>> Numeric = " ++ show num )+ $ printTestCase (">>> rel.err.= " ++ show err )+ $ eq 1e-12 num exact+ where+ exact = evalCheb as x+ num = chebyshev x (fromList as)+ err = abs (num - exact) / abs exact +evalCheb :: [Double] -> Double -> Double+evalCheb as x+ = realToFrac+ $ sum+ $ zipWith (*) (map realToFrac as)+ $ map ($ realToFrac x) cheb +-- Chebyshev polynomials of low order+cheb :: [Rational -> Rational]+cheb =+ [ \_ -> 1+ , \x -> x+ , \x -> 2*x^2 - 1+ , \x -> 4*x^3 - 3*x+ , \x -> 8*x^4 - 8*x^2 + 1+ ]+ -- Double in the [-1 .. 1] range newtype Ch = Ch Double deriving Show instance Arbitrary Ch where arbitrary = do x <- arbitrary- return $ Ch $ 2 * (snd . properFraction) x - 1+ return $ Ch $ 2 * (abs . snd . properFraction) x - 1
tests/Tests/SpecFunctions.hs view
@@ -18,15 +18,17 @@ tests :: Test tests = testGroup "Special functions"- [ testProperty "Γ(x+1) = x·Γ(x) logGamma" $ gammaReccurence logGamma 3e-8- , testProperty "Γ(x+1) = x·Γ(x) logGammaL" $ gammaReccurence logGammaL 2e-13- , testProperty "γ(1,x) = 1 - exp(-x)" $ incompleteGammaAt1Check- , testProperty "0 <= γ <= 1" $ incompleteGammaInRange- , testProperty "γ - increases" $+ [ testProperty "Gamma(x+1) = x*Gamma(x) [logGamma]" $ gammaReccurence logGamma 3e-8+ , testProperty "Gamma(x+1) = x*Gamma(x) [logGammaL]" $ gammaReccurence logGammaL 2e-13+ , testProperty "gamma(1,x) = 1 - exp(-x)" $ incompleteGammaAt1Check+ , testProperty "0 <= gamma <= 1" $ incompleteGammaInRange+ , testProperty "gamma - increases" $ \s x y -> s > 0 && x > 0 && y > 0 ==> monotonicallyIncreases (incompleteGamma s) x y- , testProperty "invIncompleteGamma = γ^-1" $ invIGammaIsInverse+ , testProperty "invIncompleteGamma = gamma^-1" $ invIGammaIsInverse , testProperty "0 <= I[B] <= 1" $ incompleteBetaInRange , testProperty "invIncompleteBeta = B^-1" $ invIBetaIsInverse+ , testProperty "invErfc = erfc^-1" $ invErfcIsInverse+ , testProperty "invErf = erf^-1" $ invErfIsInverse -- Unit tests , testAssertion "Factorial is expected to be precise at 1e-15 level" $ and [ eq 1e-15 (factorial (fromIntegral n))@@ -46,16 +48,33 @@ $ and [ eq 1e-15 (logGammaL (fromIntegral n)) (logFactorial (n-1)) | n <- [3..10000]]- , testAssertion "logGammaL is expected to be precise at 1e-9 level [fractional points]"+ -- FIXME: Too low!+ , testAssertion "logGammaL is expected to be precise at 1e-10 level [fractional points]" $ and [ eq 1e-10 (logGammaL x) lg | (x,lg) <- tableLogGamma ]+ -- FIXME: loss of precision when logBeta p q ≈ 0.+ -- Relative error doesn't work properly in this case. , testAssertion "logBeta is expected to be precise at 1e-6 level" $ and [ eq 1e-6 (logBeta p q) (logGammaL p + logGammaL q - logGammaL (p+q)) | p <- [0.1,0.2 .. 0.9] ++ [2 .. 20]- , q <- [0.1,0.2 .. 0.9] ++ [2 .. 20]]- -- FIXME: Why 1e-8? Is it due to poor precision of logBeta?+ , q <- [0.1,0.2 .. 0.9] ++ [2 .. 20]+ ]+ , testAssertion "digamma is expected to be precise at 1e-14 [integers]"+ $ digammaTestIntegers 1e-14+ -- 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 ]+ -- FIXME: Why 1e-8? Is it due to poor precision of logBeta? , testAssertion "incompleteBeta is expected to be precise at 1e-8 level" $ and [ eq 1e-8 (incompleteBeta p q x) ib | (p,q,x,ib) <- tableIncompleteBeta ]+ , 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)+ , (0.501, 3001, 3001, 0.5615652382981522803424365187631195161665429270531389)+ , (0.531, 3500, 3200, 0.9209758089734407825580172472327758548870610822321278)+ , (0.501, 13500, 13200, 0.0656209987264794057358373443387716674955276089622780)+ ]+ ] , testAssertion "choose is expected to precise at 1e-12 level" $ and [ eq 1e-12 (choose (fromIntegral n) (fromIntegral k)) (fromIntegral $ choose' n k) | n <- [0..300], k <- [0..n]]@@ -103,6 +122,30 @@ x = invIncompleteGamma a p p' = incompleteGamma a x +-- invErfc is inverse of erfc+invErfcIsInverse :: Double -> Property+invErfcIsInverse ((*2) . range01 -> p)+ = printTestCase ("p = " ++ show p )+ $ printTestCase ("x = " ++ show x )+ $ printTestCase ("p' = " ++ show p')+ $ abs (p - p') <= 1e-14+ where+ x = invErfc p+ p' = erfc x++-- invErf is inverse of erf+invErfIsInverse :: Double -> Property+invErfIsInverse a+ = printTestCase ("p = " ++ show p )+ $ printTestCase ("x = " ++ show x )+ $ printTestCase ("p' = " ++ show p')+ $ abs (p - p') <= 1e-14+ where+ x = invErf p+ p' = erf x+ p | a < 0 = - range01 a+ | otherwise = range01 a+ -- B(s,x) is in [0,1] range incompleteBetaInRange :: Double -> Double -> Double -> Property incompleteBetaInRange (abs -> p) (abs -> q) (range01 -> x) =@@ -123,6 +166,21 @@ x' = incompleteBeta p q a a = invIncompleteBeta p q x +-- Table for digamma function:+--+-- Uses equality ψ(n) = H_{n-1} - γ where+-- H_{n} = Σ 1/k, k = [1 .. n] - harmonic number+-- γ = 0.57721566490153286060 - Euler-Mascheroni number+digammaTestIntegers :: Double -> Bool+digammaTestIntegers eps+ = all (uncurry $ eq eps) $ take 3000 digammaInt+ where+ ok approx exact = approx+ -- Harmonic numbers starting from 0+ harmN = scanl (\a n -> a + 1/n) 0 [1::Rational .. ]+ gam = 0.57721566490153286060+ -- Digamma values+ digammaInt = zipWith (\i h -> (digamma i, realToFrac h - gam)) [1..] harmN ----------------------------------------------------------------
tests/Tests/SpecFunctions/Tables.hs view
@@ -2,22 +2,22 @@ tableLogGamma :: [(Double,Double)] tableLogGamma =- [(0.000001250000000, 13.592366285131769033)+ [(0.000001250000000, 13.592366285131767256) , (0.000068200000000, 9.5930266308318756785)- , (0.000246000000000, 8.3100370767447966358)- , (0.000880000000000, 7.03508133735248542)- , (0.003120000000000, 5.768129358365567505)- , (0.026700000000000, 3.6082588918892977148)+ , (0.000246000000000, 8.3100370767447948595)+ , (0.000880000000000, 7.0350813373524845318)+ , (0.003120000000000, 5.7681293583655666168)+ , (0.026700000000000, 3.6082588918892972707) , (0.077700000000000, 2.5148371858768232556)- , (0.234000000000000, 1.3579557559432759994)- , (0.860000000000000, 0.098146578027685615897)+ , (0.234000000000000, 1.3579557559432757774)+ , (0.860000000000000, 0.098146578027685588141) , (1.340000000000000, -0.11404757557207759189)- , (1.890000000000000, -0.0425116422978701336)- , (2.450000000000000, 0.25014296569217625565)+ , (1.890000000000000, -0.042511642297870140539)+ , (2.450000000000000, 0.25014296569217620014) , (3.650000000000000, 1.3701041997380685178)- , (4.560000000000000, 2.5375143317949580002)+ , (4.560000000000000, 2.5375143317949575561) , (6.660000000000000, 5.9515377269550207018)- , (8.250000000000000, 9.0331869196051233217)+ , (8.250000000000000, 9.0331869196051215454) , (11.300000000000001, 15.814180681373947834) , (25.600000000000001, 56.711261598328121636) , (50.399999999999999, 146.12815158702164808)@@ -26,22 +26,125 @@ , (853.399999999999977, 4903.9359135978220365) , (2923.300000000000182, 20402.93198938705973) , (8764.299999999999272, 70798.268343590112636)- , (12630.000000000000000, 106641.77264982508495)+ , (12630.000000000000000, 106641.7726498250704) , (34500.000000000000000, 325976.34838781820145) , (82340.000000000000000, 849629.79603036714252)- , (234800.000000000000000, 2668846.4390507959761)- , (834300.000000000000000, 10540830.912557534873)+ , (234800.000000000000000, 2668846.4390507955104)+ , (834300.000000000000000, 10540830.912557533011) , (1230000.000000000000000, 16017699.322315014899) ] tableIncompleteBeta :: [(Double,Double,Double,Double)] tableIncompleteBeta =- [(2.000000000000000, 3.000000000000000, 0.030000000000000, 0.0051864299999999996862)+ [(2.000000000000000, 3.000000000000000, 0.030000000000000, 0.0051864299999999988189) , (2.000000000000000, 3.000000000000000, 0.230000000000000, 0.22845923000000001313)- , (2.000000000000000, 3.000000000000000, 0.760000000000000, 0.95465728000000005249)- , (4.000000000000000, 2.300000000000000, 0.890000000000000, 0.93829812158347802864)+ , (2.000000000000000, 3.000000000000000, 0.760000000000000, 0.95465727999999994147)+ , (4.000000000000000, 2.300000000000000, 0.890000000000000, 0.93829812158347791762) , (1.000000000000000, 1.000000000000000, 0.550000000000000, 0.55000000000000004441)- , (0.300000000000000, 12.199999999999999, 0.110000000000000, 0.95063000053947077639)+ , (0.300000000000000, 12.199999999999999, 0.110000000000000, 0.95063000053947066537) , (13.100000000000000, 9.800000000000001, 0.120000000000000, 1.3483109941962659385e-07) , (13.100000000000000, 9.800000000000001, 0.420000000000000, 0.071321857831804780226)- , (13.100000000000000, 9.800000000000001, 0.920000000000000, 0.99999578339197081611)+ , (13.100000000000000, 9.800000000000001, 0.920000000000000, 0.99999578339197070509)+ ]+tableDigamma :: [(Double,Double)]+tableDigamma =+ [(10.0261172557341425, 2.2544954834170942704)+ , (0.9070101446062873, -0.74152778337908598072)+ , (3.4679213262860156, 1.0925031389314479036)+ , (28.5703089405901878, 3.3347652650101657912)+ , (5.9700184459319399, 1.7006665338476731897)+ , (20.5303177686997920, 2.9973508205248808878)+ , (5.6622605630542511, 1.6429280447671743559)+ , (4.4741465342999014, 1.3824198603491071324)+ , (21.4416006516504787, 3.0418326144933285349)+ , (47.6946291432301663, 3.8542988022858128971)+ , (11.2357450115053670, 2.37393979612347783)+ , (0.3352840110772935, -3.1124447967622668187)+ , (2.5037441860153118, 0.70499097759044615508)+ , (0.5241560861477529, -1.8489960634174653631)+ , (0.1972018552655726, -5.3635382066874592866)+ , (0.8289440927562556, -0.90024805153750442344)+ , (2.0717397641759350, 0.4680412969073853291)+ , (9.1173553049782452, 2.1543380160183831507)+ , (1.1815938184339669, -0.31262126373727594508)+ , (7.3600347508772019, 1.9265946441432049152)+ , (19.7457045917841398, 2.9574003365402390386)+ , (4.1956416643620571, 1.3101672771843546617)+ , (7.3868205159465790, 1.9304848277860633399)+ , (1.2786090750546355, -0.19373178842778399078)+ , (10.6498308581562604, 2.3178608134278069208)+ , (10.6750266252851169, 2.3203381265132185796)+ , (10.6883248506773985, 2.3216431742802625671)+ , (14.3373372205836365, 2.6275879484098640937)+ , (3.3932538441985769, 1.0672611106295626371)+ , (11.4168205413938768, 2.3906538776946248959)+ , (3.2500957742991048, 1.0170253699094919941)+ , (2.7573211981404855, 0.82209952378707851217)+ , (21.8943170241258827, 3.063216323919045081)+ , (16.7950471612825254, 2.7910180230044043803)+ , (9.2578640399661225, 2.1704940538770385317)+ , (5.3213868642873896, 1.5748408574979930741)+ , (9.4381079039564071, 2.1908443398518979706)+ , (13.1568457441413429, 2.538458049596743038)+ , (10.6478950333943825, 2.3176702242110884811)+ , (6.4894496431749733, 1.7911554320176725774)+ , (20.3998669454332315, 2.9908182167188113176)+ , (3.6989463639934752, 1.1668268193484248041)+ , (3.4716258279958572, 1.093739186127963281)+ , (24.7013029455164919, 3.1864775907749920414)+ , (1.1608524325026863, -0.33982067949719851896)+ , (1.9482800424522431, 0.3888762195060542215)+ , (30.4956621109554185, 3.4010990755913685923)+ , (16.3105956379859052, 2.7608468922073350349)+ , (10.6908820268137070, 2.3218939328714371939)+ , (3.4369121607821915, 1.082096765647714065)+ , (2.2914619096171260, 0.5953971130541900747)+ , (24.1273989930028883, 3.1624816269998849982)+ , (14.9455957898231535, 2.6705890837495616097)+ , (32.2002179941400826, 3.4563650137673369578)+ , (1.7232417075599473, 0.22682264125689588496)+ , (9.9662376350778192, 2.248195612105357899)+ , (10.9702870318273966, 2.348920912357223223)+ , (18.8934063317711676, 2.912115343761407793)+ , (8.6720493874148570, 2.1013420151521415846)+ , (20.4905634096258815, 2.9953645521238549954)+ , (1.4654265058258678, 0.0036653372399428492921)+ , (15.4401781010745509, 2.7042406258657996077)+ , (13.6688064138713390, 2.5780909087521290957)+ , (2.4073661551765566, 0.65668881914974130964)+ , (0.8108729056729371, -0.94026521559981879328)+ , (29.5024809785193902, 3.367430902728568487)+ , (7.5321882978878660, 1.9513375601887514854)+ , (3.3716588961200955, 1.0598414578703589939)+ , (2.9310065630306474, 0.89516303667430119351)+ , (7.2023118361897769, 1.9033764996201536501)+ , (3.1362387322050900, 0.97520764792577085966)+ , (6.5709053027851487, 1.8046329737306385788)+ , (3.7348491113356177, 1.1779005641199544741)+ , (1.2328105814385013, -0.24823346907893503732)+ , (7.9098387372709587, 2.0035651569967258823)+ , (2.8590898311999715, 0.86554629114604864082)+ , (2.1964374279534344, 0.54225028515290207842)+ , (3.8933394033155189, 1.2253803767351847398)+ , (10.7410508007627694, 2.3268008547643748152)+ , (2.4921048837305193, 0.69927782909414781809)+ , (2.2101710538553756, 0.55010424351998354897)+ , (14.0357118427322334, 2.6055587167248708269)+ , (4.1320729121597584, 1.2929216807716104043)+ , (0.2766365979680845, -3.8108738889017752527)+ , (27.9448247140513644, 3.3122329205038494315)+ , (9.3081256750537182, 2.1762105230057038341)+ , (1.4222181352589696, -0.038843893649701873028)+ , (1.5107587188614726, 0.046499571962236106726)+ , (3.3467578222470555, 1.0512176183500512305)+ , (12.2373583939228876, 2.4630788434421742039)+ , (0.9385094944630431, -0.68317598609698348966)+ , (5.8655552400886410, 1.6814385243672138603)+ , (17.1377048621110468, 2.8118219246156086477)+ , (4.0502102843199079, 1.2702685434611069581)+ , (2.2041235084734976, 0.54665320805956585382)+ , (0.9498749870396368, -0.66283138696545962354)+ , (5.5020466797149687, 1.6115010556650317675)+ , (1.8741725410778542, 0.33826100356492333487)+ , (14.1730624058772161, 2.6156503142962224118)+ , (1.0704026637921555, -0.46701211139417769802) ]
tests/Tests/SpecFunctions/gen.py view
@@ -3,13 +3,20 @@ """ from mpmath import *+import random +# Set very-very large precision+mp.dps = 100+# Set fixed seed in order to get repeatable results+random.seed( 279570842 )+ def printListLiteral(lines) : print " [" + "\n , ".join(lines) + "\n ]" + ################################################################ # Generate header-print "module Tests.Math.Tables where"+print "module Tests.SpecFunctions.Tables where" print ################################################################@@ -49,3 +56,15 @@ [ '(%.15f, %.15f, %.15f, %.20g)' % (p,q,x, betainc(p,q,0,x, regularized=True)) for (p,q,x) in incompleteBetaArg ])+++################################################################+## Generate table for digamma++print "tableDigamma :: [(Double,Double)]"+print "tableDigamma ="+printListLiteral(+ [ '(%.16f, %.20g)' % (x, digamma(x)) for x in+ [ random.expovariate(0.1) for i in xrange(100) ]+ ]+ )