diff --git a/Numeric/MathFunctions/Constants.hs b/Numeric/MathFunctions/Constants.hs --- a/Numeric/MathFunctions/Constants.hs +++ b/Numeric/MathFunctions/Constants.hs @@ -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 #-} diff --git a/Numeric/Polynomial.hs b/Numeric/Polynomial.hs new file mode 100644 --- /dev/null +++ b/Numeric/Polynomial.hs @@ -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 diff --git a/Numeric/Polynomial/Chebyshev.hs b/Numeric/Polynomial/Chebyshev.hs --- a/Numeric/Polynomial/Chebyshev.hs +++ b/Numeric/Polynomial/Chebyshev.hs @@ -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. diff --git a/Numeric/SpecFunctions.hs b/Numeric/SpecFunctions.hs --- a/Numeric/SpecFunctions.hs +++ b/Numeric/SpecFunctions.hs @@ -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. +-- +-- +-- * 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 +-- +-- -- * Lanczos, C. (1964) A precision approximation of the gamma -- function. /SIAM Journal on Numerical Analysis B/ -- 1:86–96. @@ -576,10 +825,6 @@ -- /Journal of the Royal Statistical Society, Series C (Applied Statistics)/ -- 38(2):397–402. -- --- * Shea, B. (1988) Algorithm AS 239: Chi-squared and incomplete --- gamma integral. /Applied Statistics/ --- 37(3):466–473. --- -- * 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 -- -- --- * 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 --- +-- * Shea, B. (1988) Algorithm AS 239: Chi-squared and incomplete +-- gamma integral. /Applied Statistics/ +-- 37(3):466–473. 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.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 diff --git a/tests/Tests/Chebyshev.hs b/tests/Tests/Chebyshev.hs --- a/tests/Tests/Chebyshev.hs +++ b/tests/Tests/Chebyshev.hs @@ -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 diff --git a/tests/Tests/SpecFunctions.hs b/tests/Tests/SpecFunctions.hs --- a/tests/Tests/SpecFunctions.hs +++ b/tests/Tests/SpecFunctions.hs @@ -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 ---------------------------------------------------------------- diff --git a/tests/Tests/SpecFunctions/Tables.hs b/tests/Tests/SpecFunctions/Tables.hs --- a/tests/Tests/SpecFunctions/Tables.hs +++ b/tests/Tests/SpecFunctions/Tables.hs @@ -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) ] diff --git a/tests/Tests/SpecFunctions/gen.py b/tests/Tests/SpecFunctions/gen.py --- a/tests/Tests/SpecFunctions/gen.py +++ b/tests/Tests/SpecFunctions/gen.py @@ -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) ] + ] + )