1{-# LANGUAGE BangPatterns, ScopedTypeVariables #-} 2{-# OPTIONS_HADDOCK hide #-} 3-- | 4-- Module : Numeric.SpecFunctions.Internal 5-- Copyright : (c) 2009, 2011, 2012 Bryan O'Sullivan 6-- License : BSD3 7-- 8-- Maintainer : bos@serpentine.com 9-- Stability : experimental 10-- Portability : portable 11-- 12-- Internal module with implementation of special functions. 13module Numeric.SpecFunctions.Internal 14 ( module Numeric.SpecFunctions.Internal 15 , Compat.log1p 16 , Compat.expm1 17 ) where 18 19import Control.Applicative 20import Data.Bits ((.&.), (.|.), shiftR) 21import Data.Int (Int64) 22import Data.Word (Word) 23import Data.Default.Class 24import qualified Data.Vector.Unboxed as U 25import Data.Vector.Unboxed ((!)) 26import Text.Printf 27 28import Numeric.Polynomial.Chebyshev (chebyshevBroucke) 29import Numeric.Polynomial (evaluatePolynomial, evaluatePolynomialL, evaluateEvenPolynomialL 30 ,evaluateOddPolynomialL) 31import Numeric.RootFinding (Root(..), newtonRaphson, NewtonParam(..), Tolerance(..)) 32import Numeric.Series 33import Numeric.MathFunctions.Constants 34import Numeric.SpecFunctions.Compat (log1p) 35import qualified Numeric.SpecFunctions.Compat as Compat 36 37---------------------------------------------------------------- 38-- Error function 39---------------------------------------------------------------- 40 41-- | Error function. 42-- 43-- \[ 44-- \operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} \exp(-t^2) dt 45-- \] 46-- 47-- Function limits are: 48-- 49-- \[ 50-- \begin{aligned} 51-- &\operatorname{erf}(-\infty) &=& -1 \\ 52-- &\operatorname{erf}(0) &=& \phantom{-}\,0 \\ 53-- &\operatorname{erf}(+\infty) &=& \phantom{-}\,1 \\ 54-- \end{aligned} 55-- \] 56erf :: Double -> Double 57erf = Compat.erf 58{-# INLINE erf #-} 59 60-- | Complementary error function. 61-- 62-- \[ 63-- \operatorname{erfc}(x) = 1 - \operatorname{erf}(x) 64-- \] 65-- 66-- Function limits are: 67-- 68-- \[ 69-- \begin{aligned} 70-- &\operatorname{erf}(-\infty) &=&\, 2 \\ 71-- &\operatorname{erf}(0) &=&\, 1 \\ 72-- &\operatorname{erf}(+\infty) &=&\, 0 \\ 73-- \end{aligned} 74-- \] 75erfc :: Double -> Double 76erfc = Compat.erfc 77{-# INLINE erfc #-} 78 79-- | Inverse of 'erf'. 80invErf :: Double -- ^ /p/ ∈ [-1,1] 81 -> Double 82invErf p 83 | p == 1 = m_pos_inf 84 | p == -1 = m_neg_inf 85 | p < 1 && p > -1 = if p > 0 then r else -r 86 | otherwise = error "invErf: p must in [-1,1] range" 87 where 88 -- We solve equation same as in invErfc. We're able to ruse same 89 -- Halley step by solving equation: 90 -- > pp - erf x = 0 91 -- instead of 92 -- > erf x - pp = 0 93 pp = abs p 94 r = step $ step $ guessInvErfc $ 1 - pp 95 step x = invErfcHalleyStep (pp - erf x) x 96 97-- | Inverse of 'erfc'. 98invErfc :: Double -- ^ /p/ ∈ [0,2] 99 -> Double 100invErfc p 101 | p == 2 = m_neg_inf 102 | p == 0 = m_pos_inf 103 | p >0 && p < 2 = if p <= 1 then r else -r 104 | otherwise = modErr $ "invErfc: p must be in [0,2] got " ++ show p 105 where 106 pp | p <= 1 = p 107 | otherwise = 2 - p 108 -- We perform 2 Halley steps in order to get to solution 109 r = step $ step $ guessInvErfc pp 110 step x = invErfcHalleyStep (erfc x - pp) x 111 112-- Initial guess for invErfc & invErf 113guessInvErfc :: Double -> Double 114guessInvErfc p 115 = -0.70711 * ((2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t) 116 where 117 t = sqrt $ -2 * log( 0.5 * p) 118 119-- Halley step for solving invErfc 120invErfcHalleyStep :: Double -> Double -> Double 121invErfcHalleyStep err x 122 = x + err / (1.12837916709551257 * exp(-x * x) - x * err) 123 124---------------------------------------------------------------- 125-- Gamma function 126---------------------------------------------------------------- 127 128data L = L {-# UNPACK #-} !Double {-# UNPACK #-} !Double 129 130-- | Compute the logarithm of the gamma function, Γ(/x/). 131-- 132-- \[ 133-- \Gamma(x) = \int_0^{\infty}t^{x-1}e^{-t}\,dt = (x - 1)! 134-- \] 135-- 136-- This implementation uses Lanczos approximation. It gives 14 or more 137-- significant decimal digits, except around /x/ = 1 and /x/ = 2, 138-- where the function goes to zero. 139-- 140-- Returns ∞ if the input is outside of the range (0 < /x/ 141-- ≤ 1e305). 142logGamma :: Double -> Double 143logGamma z 144 | z <= 0 = m_pos_inf 145 -- For very small values z we can just use Laurent expansion 146 | z < m_sqrt_eps = log (1/z - m_eulerMascheroni) 147 -- For z<1 we use recurrence. Γ(z+1) = z·Γ(z) Note that in order to 148 -- avoid precision loss we have to compute parameter to 149 -- approximations here: 150 -- 151 -- > (z + 1) - 1 = z 152 -- > (z + 1) - 2 = z - 1 153 -- 154 -- Simple passing (z + 1) to piecewise approxiations and computing 155 -- difference leads to bad loss of precision near 1. 156 -- This is reason lgamma1_15 & lgamma15_2 have three parameters 157 | z < 0.5 = lgamma1_15 z (z - 1) - log z 158 | z < 1 = lgamma15_2 z (z - 1) - log z 159 -- Piecewise polynomial approximations 160 | z <= 1.5 = lgamma1_15 (z - 1) (z - 2) 161 | z < 2 = lgamma15_2 (z - 1) (z - 2) 162 | z < 15 = lgammaSmall z 163 -- Otherwise we switch to Lanczos approximation 164 | otherwise = lanczosApprox z 165 166 167-- | Synonym for 'logGamma'. Retained for compatibility 168logGammaL :: Double -> Double 169logGammaL = logGamma 170{-# DEPRECATED logGammaL "Use logGamma instead" #-} 171 172 173 174-- Polynomial expansion used in interval (1,1.5] 175-- 176-- > logΓ(z) = (z-1)(z-2)(Y + R(z-1)) 177lgamma1_15 :: Double -> Double -> Double 178lgamma1_15 zm1 zm2 179 = r * y + r * ( evaluatePolynomial zm1 tableLogGamma_1_15P 180 / evaluatePolynomial zm1 tableLogGamma_1_15Q 181 ) 182 where 183 r = zm1 * zm2 184 y = 0.52815341949462890625 185 186tableLogGamma_1_15P,tableLogGamma_1_15Q :: U.Vector Double 187tableLogGamma_1_15P = U.fromList 188 [ 0.490622454069039543534e-1 189 , -0.969117530159521214579e-1 190 , -0.414983358359495381969e0 191 , -0.406567124211938417342e0 192 , -0.158413586390692192217e0 193 , -0.240149820648571559892e-1 194 , -0.100346687696279557415e-2 195 ] 196{-# NOINLINE tableLogGamma_1_15P #-} 197tableLogGamma_1_15Q = U.fromList 198 [ 1 199 , 0.302349829846463038743e1 200 , 0.348739585360723852576e1 201 , 0.191415588274426679201e1 202 , 0.507137738614363510846e0 203 , 0.577039722690451849648e-1 204 , 0.195768102601107189171e-2 205 ] 206{-# NOINLINE tableLogGamma_1_15Q #-} 207 208 209 210-- Polynomial expansion used in interval (1.5,2) 211-- 212-- > logΓ(z) = (2-z)(1-z)(Y + R(2-z)) 213lgamma15_2 :: Double -> Double -> Double 214lgamma15_2 zm1 zm2 215 = r * y + r * ( evaluatePolynomial (-zm2) tableLogGamma_15_2P 216 / evaluatePolynomial (-zm2) tableLogGamma_15_2Q 217 ) 218 where 219 r = zm1 * zm2 220 y = 0.452017307281494140625 221 222tableLogGamma_15_2P,tableLogGamma_15_2Q :: U.Vector Double 223tableLogGamma_15_2P = U.fromList 224 [ -0.292329721830270012337e-1 225 , 0.144216267757192309184e0 226 , -0.142440390738631274135e0 227 , 0.542809694055053558157e-1 228 , -0.850535976868336437746e-2 229 , 0.431171342679297331241e-3 230 ] 231{-# NOINLINE tableLogGamma_15_2P #-} 232tableLogGamma_15_2Q = U.fromList 233 [ 1 234 , -0.150169356054485044494e1 235 , 0.846973248876495016101e0 236 , -0.220095151814995745555e0 237 , 0.25582797155975869989e-1 238 , -0.100666795539143372762e-2 239 , -0.827193521891290553639e-6 240 ] 241{-# NOINLINE tableLogGamma_15_2Q #-} 242 243 244 245-- Polynomial expansion used in interval (2,3) 246-- 247-- > logΓ(z) = (z - 2)(z + 1)(Y + R(z-2)) 248lgamma2_3 :: Double -> Double 249lgamma2_3 z 250 = r * y + r * ( evaluatePolynomial zm2 tableLogGamma_2_3P 251 / evaluatePolynomial zm2 tableLogGamma_2_3Q 252 ) 253 where 254 r = zm2 * (z + 1) 255 zm2 = z - 2 256 y = 0.158963680267333984375e0 257 258 259tableLogGamma_2_3P,tableLogGamma_2_3Q :: U.Vector Double 260tableLogGamma_2_3P = U.fromList 261 [ -0.180355685678449379109e-1 262 , 0.25126649619989678683e-1 263 , 0.494103151567532234274e-1 264 , 0.172491608709613993966e-1 265 , -0.259453563205438108893e-3 266 , -0.541009869215204396339e-3 267 , -0.324588649825948492091e-4 268 ] 269{-# NOINLINE tableLogGamma_2_3P #-} 270tableLogGamma_2_3Q = U.fromList 271 [ 1 272 , 0.196202987197795200688e1 273 , 0.148019669424231326694e1 274 , 0.541391432071720958364e0 275 , 0.988504251128010129477e-1 276 , 0.82130967464889339326e-2 277 , 0.224936291922115757597e-3 278 , -0.223352763208617092964e-6 279 ] 280{-# NOINLINE tableLogGamma_2_3Q #-} 281 282 283 284-- For small z we can just use Gamma function recurrence and reduce 285-- problem to interval [2,3] and use polynomial approximation 286-- there. Surpringly it gives very good precision 287lgammaSmall :: Double -> Double 288lgammaSmall = go 0 289 where 290 go acc z | z < 3 = acc + lgamma2_3 z 291 | otherwise = go (acc + log zm1) zm1 292 where 293 zm1 = z - 1 294 295 296-- Lanczos approximation for gamma function. 297-- 298-- > Γ(z) = sqrt(2π)(z + g - 0.5)^(z - 0.5)·exp{-(z + g - 0.5)}·A_g(z) 299-- 300-- Coeffients are taken from boost. Constants are absorbed into 301-- polynomial's coefficients. 302lanczosApprox :: Double -> Double 303lanczosApprox z 304 = (log (z + g - 0.5) - 1) * (z - 0.5) 305 + log (evalRatio tableLanczos z) 306 where 307 g = 6.024680040776729583740234375 308 309tableLanczos :: U.Vector (Double,Double) 310{-# NOINLINE tableLanczos #-} 311tableLanczos = U.fromList 312 [ (56906521.91347156388090791033559122686859 , 0) 313 , (103794043.1163445451906271053616070238554 , 39916800) 314 , (86363131.28813859145546927288977868422342 , 120543840) 315 , (43338889.32467613834773723740590533316085 , 150917976) 316 , (14605578.08768506808414169982791359218571 , 105258076) 317 , (3481712.15498064590882071018964774556468 , 45995730) 318 , (601859.6171681098786670226533699352302507 , 13339535) 319 , (75999.29304014542649875303443598909137092 , 2637558) 320 , (6955.999602515376140356310115515198987526 , 357423) 321 , (449.9445569063168119446858607650988409623 , 32670) 322 , (19.51992788247617482847860966235652136208 , 1925) 323 , (0.5098416655656676188125178644804694509993 , 66) 324 , (0.006061842346248906525783753964555936883222 , 1) 325 ] 326 327-- Evaluate rational function. Polynomials in both numerator and 328-- denominator must have same order. Function seems to be too specific 329-- so it's not exposed 330-- 331-- Special care taken in order to avoid overflow for large values of x 332evalRatio :: U.Vector (Double,Double) -> Double -> Double 333evalRatio coef x 334 | x > 1 = fini $ U.foldl' stepL (L 0 0) coef 335 | otherwise = fini $ U.foldr' stepR (L 0 0) coef 336 where 337 fini (L num den) = num / den 338 stepR (a,b) (L num den) = L (num * x + a) (den * x + b) 339 stepL (L num den) (a,b) = L (num * rx + a) (den * rx + b) 340 rx = recip x 341 342 343 344-- | 345-- Compute the log gamma correction factor for Stirling 346-- approximation for @x@ ≥ 10. This correction factor is 347-- suitable for an alternate (but less numerically accurate) 348-- definition of 'logGamma': 349-- 350-- \[ 351-- \log\Gamma(x) = \frac{1}{2}\log(2\pi) + (x-\frac{1}{2})\log x - x + \operatorname{logGammaCorrection}(x) 352-- \] 353logGammaCorrection :: Double -> Double 354logGammaCorrection x 355 | x < 10 = m_NaN 356 | x < big = chebyshevBroucke (t * t * 2 - 1) coeffs / x 357 | otherwise = 1 / (x * 12) 358 where 359 big = 94906265.62425156 360 t = 10 / x 361 coeffs = U.fromList [ 362 0.1666389480451863247205729650822e+0, 363 -0.1384948176067563840732986059135e-4, 364 0.9810825646924729426157171547487e-8, 365 -0.1809129475572494194263306266719e-10, 366 0.6221098041892605227126015543416e-13, 367 -0.3399615005417721944303330599666e-15, 368 0.2683181998482698748957538846666e-17 369 ] 370 371 372 373-- | Compute the normalized lower incomplete gamma function 374-- γ(/z/,/x/). Normalization means that γ(/z/,∞)=1 375-- 376-- \[ 377-- \gamma(z,x) = \frac{1}{\Gamma(z)}\int_0^{x}t^{z-1}e^{-t}\,dt 378-- \] 379-- 380-- Uses Algorithm AS 239 by Shea. 381incompleteGamma :: Double -- ^ /z/ ∈ (0,∞) 382 -> Double -- ^ /x/ ∈ (0,∞) 383 -> Double 384-- Notation used: 385-- + P(a,x) - regularized lower incomplete gamma 386-- + Q(a,x) - regularized upper incomplete gamma 387incompleteGamma a x 388 | a <= 0 || x < 0 = error 389 $ "incompleteGamma: Domain error z=" ++ show a ++ " x=" ++ show x 390 | x == 0 = 0 391 | x == m_pos_inf = 1 392 -- For very small x we use following expansion for P: 393 -- 394 -- See http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/ 395 | x < sqrt m_epsilon && a > 1 396 = x**a / a / exp (logGamma a) * (1 - a*x / (a + 1)) 397 | x < 0.5 = case () of 398 _| (-0.4)/log x < a -> taylorSeriesP 399 | otherwise -> taylorSeriesComplQ 400 | x < 1.1 = case () of 401 _| 0.75*x < a -> taylorSeriesP 402 | otherwise -> taylorSeriesComplQ 403 | a > 20 && useTemme = uniformExpansion 404 | x - (1 / (3 * x)) < a = taylorSeriesP 405 | otherwise = contFraction 406 where 407 mu = (x - a) / a 408 useTemme = (a > 200 && 20/a > mu*mu) 409 || (abs mu < 0.4) 410 -- Gautschi's algorithm. 411 -- 412 -- Evaluate series for P(a,x). See [Temme1994] Eq. 5.5 and [NOTE: 413 -- incompleteGamma.taylorP] 414 factorP 415 | a < 10 = x ** a 416 / (exp x * exp (logGamma (a + 1))) 417 | a < 1182.5 = (x * exp 1 / a) ** a 418 / exp x 419 / sqrt (2*pi*a) 420 / exp (logGammaCorrection a) 421 | otherwise = (x * exp 1 / a * exp (-x/a)) ** a 422 / sqrt (2*pi*a) 423 / exp (logGammaCorrection a) 424 taylorSeriesP 425 = sumPowerSeries x (scanSequence (/) 1 $ enumSequenceFrom (a+1)) 426 * factorP 427 -- Series for 1-Q(a,x). See [Temme1994] Eq. 5.5 428 taylorSeriesComplQ 429 = sumPowerSeries (-x) (scanSequence (/) 1 (enumSequenceFrom 1) / enumSequenceFrom a) 430 * x**a / exp(logGamma a) 431 -- Legendre continued fractions 432 contFraction = 1 - ( exp ( log x * a - x - logGamma a ) 433 / evalContFractionB frac 434 ) 435 where 436 frac = (\k -> (k*(a-k), x - a + 2*k + 1)) <$> enumSequenceFrom 0 437 -- Evaluation based on uniform expansions. See [Temme1994] 5.2 438 uniformExpansion = 439 let -- Coefficients f_m in paper 440 fm :: U.Vector Double 441 fm = U.fromList [ 1.00000000000000000000e+00 442 ,-3.33333333333333370341e-01 443 , 8.33333333333333287074e-02 444 ,-1.48148148148148153802e-02 445 , 1.15740740740740734316e-03 446 , 3.52733686067019369930e-04 447 ,-1.78755144032921825352e-04 448 , 3.91926317852243766954e-05 449 ,-2.18544851067999240532e-06 450 ,-1.85406221071515996597e-06 451 , 8.29671134095308545622e-07 452 ,-1.76659527368260808474e-07 453 , 6.70785354340149841119e-09 454 , 1.02618097842403069078e-08 455 ,-4.38203601845335376897e-09 456 , 9.14769958223679020897e-10 457 ,-2.55141939949462514346e-11 458 ,-5.83077213255042560744e-11 459 , 2.43619480206674150369e-11 460 ,-5.02766928011417632057e-12 461 , 1.10043920319561347525e-13 462 , 3.37176326240098513631e-13 463 ] 464 y = - log1pmx mu 465 eta = sqrt (2 * y) * signum mu 466 -- Evaluate S_α (Eq. 5.9) 467 loop !_ !_ u 0 = u 468 loop bm1 bm0 u i = let t = (fm ! i) + (fromIntegral i + 1)*bm1 / a 469 u' = eta * u + t 470 in loop bm0 t u' (i-1) 471 s_a = let n = U.length fm 472 in loop (fm ! (n-1)) (fm ! (n-2)) 0 (n-3) 473 / exp (logGammaCorrection a) 474 in 1/2 * erfc(-eta*sqrt(a/2)) - exp(-(a*y)) / sqrt (2*pi*a) * s_a 475 476 477 478-- Adapted from Numerical Recipes §6.2.1 479 480-- | Inverse incomplete gamma function. It's approximately inverse of 481-- 'incompleteGamma' for the same /z/. So following equality 482-- approximately holds: 483-- 484-- > invIncompleteGamma z . incompleteGamma z ≈ id 485invIncompleteGamma :: Double -- ^ /z/ ∈ (0,∞) 486 -> Double -- ^ /p/ ∈ [0,1] 487 -> Double 488invIncompleteGamma a p 489 | a <= 0 = 490 modErr $ printf "invIncompleteGamma: a must be positive. a=%g p=%g" a p 491 | p < 0 || p > 1 = 492 modErr $ printf "invIncompleteGamma: p must be in [0,1] range. a=%g p=%g" a p 493 | p == 0 = 0 494 | p == 1 = 1 / 0 495 | otherwise = loop 0 guess 496 where 497 -- Solve equation γ(a,x) = p using Halley method 498 loop :: Int -> Double -> Double 499 loop i x 500 | i >= 12 = x' 501 -- For small s derivative becomes approximately 1/x*exp(-x) and 502 -- skyrockets for small x. If it happens correct answer is 0. 503 | isInfinite f' = 0 504 | abs dx < eps * x' = x' 505 | otherwise = loop (i + 1) x' 506 where 507 -- Value of γ(a,x) - p 508 f = incompleteGamma a x - p 509 -- dγ(a,x)/dx 510 f' | a > 1 = afac * exp( -(x - a1) + a1 * (log x - lna1)) 511 | otherwise = exp( -x + a1 * log x - gln) 512 u = f / f' 513 -- Halley correction to Newton-Rapson step 514 corr = u * (a1 / x - 1) 515 dx = u / (1 - 0.5 * min 1.0 corr) 516 -- New approximation to x 517 x' | x < dx = 0.5 * x -- Do not go below 0 518 | otherwise = x - dx 519 -- Calculate inital guess for root 520 guess 521 -- 522 | a > 1 = 523 let t = sqrt $ -2 * log(if p < 0.5 then p else 1 - p) 524 x1 = (2.30753 + t * 0.27061) / (1 + t * (0.99229 + t * 0.04481)) - t 525 x2 = if p < 0.5 then -x1 else x1 526 in max 1e-3 (a * (1 - 1/(9*a) - x2 / (3 * sqrt a)) ** 3) 527 -- For a <= 1 use following approximations: 528 -- γ(a,1) ≈ 0.253a + 0.12a² 529 -- 530 -- γ(a,x) ≈ γ(a,1)·x^a x < 1 531 -- γ(a,x) ≈ γ(a,1) + (1 - γ(a,1))(1 - exp(1 - x)) x >= 1 532 | otherwise = 533 let t = 1 - a * (0.253 + a*0.12) 534 in if p < t 535 then (p / t) ** (1 / a) 536 else 1 - log( 1 - (p-t) / (1-t)) 537 -- Constants 538 a1 = a - 1 539 lna1 = log a1 540 afac = exp( a1 * (lna1 - 1) - gln ) 541 gln = logGamma a 542 eps = 1e-8 543 544 545 546---------------------------------------------------------------- 547-- Beta function 548---------------------------------------------------------------- 549 550-- | Compute the natural logarithm of the beta function. 551-- 552-- \[ 553-- B(a,b) = \int_0^1 t^{a-1}(1-t)^{b-1}\,dt = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)} 554-- \] 555logBeta 556 :: Double -- ^ /a/ > 0 557 -> Double -- ^ /b/ > 0 558 -> Double 559logBeta a b 560 | p < 0 = m_NaN 561 | p == 0 = m_pos_inf 562 | p >= 10 = allStirling 563 | q >= 10 = twoStirling 564 -- This order of summands marginally improves precision 565 | otherwise = logGamma p + (logGamma q - logGamma pq) 566 where 567 p = min a b 568 q = max a b 569 ppq = p / pq 570 pq = p + q 571 -- When both parameters are large than 10 we can use Stirling 572 -- approximation with correction. It's more precise than sum of 573 -- logarithms of gamma functions 574 allStirling 575 = log q * (-0.5) 576 + m_ln_sqrt_2_pi 577 + logGammaCorrection p 578 + (logGammaCorrection q - logGammaCorrection pq) 579 + (p - 0.5) * log ppq 580 + q * log1p(-ppq) 581 -- Otherwise only two of three gamma functions use Stirling 582 -- approximation 583 twoStirling 584 = logGamma p 585 + (logGammaCorrection q - logGammaCorrection pq) 586 + p 587 - p * log pq 588 + (q - 0.5) * log1p(-ppq) 589 590 591-- | Regularized incomplete beta function. 592-- 593-- \[ 594-- I(x;a,b) = \frac{1}{B(a,b)} \int_0^x t^{a-1}(1-t)^{b-1}\,dt 595-- \] 596-- 597-- Uses algorithm AS63 by Majumder and Bhattachrjee and quadrature 598-- approximation for large /p/ and /q/. 599incompleteBeta :: Double -- ^ /a/ > 0 600 -> Double -- ^ /b/ > 0 601 -> Double -- ^ /x/, must lie in [0,1] range 602 -> Double 603incompleteBeta p q = incompleteBeta_ (logBeta p q) p q 604 605-- | Regularized incomplete beta function. Same as 'incompleteBeta' 606-- but also takes logarithm of beta function as parameter. 607incompleteBeta_ :: Double -- ^ logarithm of beta function for given /p/ and /q/ 608 -> Double -- ^ /a/ > 0 609 -> Double -- ^ /b/ > 0 610 -> Double -- ^ /x/, must lie in [0,1] range 611 -> Double 612incompleteBeta_ beta p q x 613 | p <= 0 || q <= 0 = 614 modErr $ printf "incompleteBeta_: p <= 0 || q <= 0. p=%g q=%g x=%g" p q x 615 | x < 0 || x > 1 || isNaN x = 616 modErr $ printf "incompleteBeta_: x out of [0,1] range. p=%g q=%g x=%g" p q x 617 | x == 0 || x == 1 = x 618 | p >= (p+q) * x = incompleteBetaWorker beta p q x 619 | otherwise = 1 - incompleteBetaWorker beta q p (1 - x) 620 621 622-- Approximation of incomplete beta by quandrature. 623-- 624-- Note that x =< p/(p+q) 625incompleteBetaApprox :: Double -> Double -> Double -> Double -> Double 626incompleteBetaApprox beta p q x 627 | ans > 0 = 1 - ans 628 | otherwise = -ans 629 where 630 -- Constants 631 p1 = p - 1 632 q1 = q - 1 633 mu = p / (p + q) 634 lnmu = log mu 635 lnmuc = log1p (-mu) 636 -- Upper limit for integration 637 xu = max 0 $ min (mu - 10*t) (x - 5*t) 638 where 639 t = sqrt $ p*q / ( (p+q) * (p+q) * (p + q + 1) ) 640 -- Calculate incomplete beta by quadrature 641 go y w = let t = x + (xu - x) * y 642 in w * exp( p1 * (log t - lnmu) + q1 * (log(1-t) - lnmuc) ) 643 s = U.sum $ U.zipWith go coefY coefW 644 ans = s * (xu - x) * exp( p1 * lnmu + q1 * lnmuc - beta ) 645 646 647-- Worker for incomplete beta function. It is separate function to 648-- avoid confusion with parameter during parameter swapping 649incompleteBetaWorker :: Double -> Double -> Double -> Double -> Double 650incompleteBetaWorker beta p q x 651 -- For very large p and q this method becomes very slow so another 652 -- method is used. 653 | p > 3000 && q > 3000 = incompleteBetaApprox beta p q x 654 | otherwise = loop (p+q) (truncate $ q + cx * (p+q)) 1 1 1 655 where 656 -- Constants 657 eps = 1e-15 658 cx = 1 - x 659 -- Common multiplies for expansion. Accurate calculation is a bit 660 -- tricky. Performing calculation in log-domain leads to slight 661 -- loss of precision for small x, while using ** prone to 662 -- underflows. 663 -- 664 -- If either beta function of x**p·(1-x)**(q-1) underflows we 665 -- switch to log domain. It could waste work but there's no easy 666 -- switch criterion. 667 factor 668 | beta < m_min_log || prod < m_tiny = exp( p * log x + (q - 1) * log cx - beta) 669 | otherwise = prod / exp beta 670 where 671 prod = x**p * cx**(q - 1) 672 -- Soper's expansion of incomplete beta function 673 loop !psq (ns :: Int) ai term betain 674 | done = betain' * factor / p 675 | otherwise = loop psq' (ns - 1) (ai + 1) term' betain' 676 where 677 -- New values 678 term' = term * fact / (p + ai) 679 betain' = betain + term' 680 fact | ns > 0 = (q - ai) * x/cx 681 | ns == 0 = (q - ai) * x 682 | otherwise = psq * x 683 -- Iterations are complete 684 done = db <= eps && db <= eps*betain' where db = abs term' 685 psq' = if ns < 0 then psq + 1 else psq 686 687 688 689-- | Compute inverse of regularized incomplete beta function. Uses 690-- initial approximation from AS109, AS64 and Halley method to solve 691-- equation. 692invIncompleteBeta :: Double -- ^ /a/ > 0 693 -> Double -- ^ /b/ > 0 694 -> Double -- ^ /x/ ∈ [0,1] 695 -> Double 696invIncompleteBeta p q a 697 | p <= 0 || q <= 0 = 698 modErr $ printf "invIncompleteBeta p <= 0 || q <= 0. p=%g q=%g a=%g" p q a 699 | a < 0 || a > 1 = 700 modErr $ printf "invIncompleteBeta x must be in [0,1]. p=%g q=%g a=%g" p q a 701 | a == 0 || a == 1 = a 702 | otherwise = invIncompleteBetaWorker (logBeta p q) p q a 703 704 705invIncompleteBetaWorker :: Double -> Double -> Double -> Double -> Double 706invIncompleteBetaWorker beta a b p = loop (0::Int) (invIncBetaGuess beta a b p) 707 where 708 a1 = a - 1 709 b1 = b - 1 710 -- Solve equation using Halley method 711 loop !i !x 712 -- We cannot continue at this point so we simply return `x' 713 | x == 0 || x == 1 = x 714 -- When derivative becomes infinite we cannot continue 715 -- iterations. It can only happen in vicinity of 0 or 1. It's 716 -- hardly possible to get good answer in such circumstances but 717 -- `x' is already reasonable. 718 | isInfinite f' = x 719 -- Iterations limit reached. Most of the time solution will 720 -- converge to answer because of discreteness of Double. But 721 -- solution have good precision already. 722 | i >= 10 = x 723 -- Solution converges 724 | abs dx <= 16 * m_epsilon * x = x' 725 | otherwise = loop (i+1) x' 726 where 727 -- Calculate Halley step. 728 f = incompleteBeta_ beta a b x - p 729 f' = exp $ a1 * log x + b1 * log1p (-x) - beta 730 u = f / f' 731 -- We bound Halley correction to Newton-Raphson to (-1,1) range 732 corr | d > 1 = 1 733 | d < -1 = -1 734 | isNaN d = 0 735 | otherwise = d 736 where 737 d = u * (a1 / x - b1 / (1 - x)) 738 dx = u / (1 - 0.5 * corr) 739 -- Next approximation. If Halley step leads us out of [0,1] 740 -- range we revert to bisection. 741 x' | z < 0 = x / 2 742 | z > 1 = (x + 1) / 2 743 | otherwise = z 744 where z = x - dx 745 746 747-- Calculate initial guess for inverse incomplete beta function. 748invIncBetaGuess :: Double -> Double -> Double -> Double -> Double 749-- Calculate initial guess. for solving equation for inverse incomplete beta. 750-- It's really hodgepodge of different approximations accumulated over years. 751-- 752-- Equations are referred to by name of paper and number e.g. [AS64 2] 753-- In AS64 papers equations are not numbered so they are refered to by 754-- number of appearance starting from definition of incomplete beta. 755invIncBetaGuess beta a b p 756 -- If both a and b are less than 1 incomplete beta have inflection 757 -- point. 758 -- 759 -- > x = (1 - a) / (2 - a - b) 760 -- 761 -- We approximate incomplete beta by neglecting one of factors under 762 -- integral and then rescaling result of integration into [0,1] 763 -- range. 764 | a < 1 && b < 1 = 765 let x_infl = (1 - a) / (2 - a - b) 766 p_infl = incompleteBeta a b x_infl 767 x | p < p_infl = let xg = (a * p * exp beta) ** (1/a) in xg / (1+xg) 768 | otherwise = let xg = (b * (1-p) * exp beta) ** (1/b) in 1 - xg/(1+xg) 769 in x 770 -- If both a and b larger or equal that 1 but not too big we use 771 -- same approximation as above but calculate it a bit differently 772 | a+b <= 6 && a>1 && b>1 = 773 let x_infl = (a - 1) / (a + b - 2) 774 p_infl = incompleteBeta a b x_infl 775 x | p < p_infl = exp ((log(p * a) + beta) / a) 776 | otherwise = 1 - exp((log((1-p) * b) + beta) / b) 777 in x 778 -- For small a and not too big b we use approximation from boost. 779 | b < 5 && a <= 1 = 780 let x | p**(1/a) < 0.5 = (p ** (a * exp beta)) ** (1/a) 781 | otherwise = 1 - (1 - p ** (b * exp beta))**(1/b) 782 in x 783 -- When a>>b and both are large approximation from [Temme1992], 784 -- section 4 "the incomplete gamma function case" used. In this 785 -- region it greatly improves over other approximation (AS109, AS64, 786 -- "Numerical Recipes") 787 -- 788 -- FIXME: It could be used when b>>a too but it require inverse of 789 -- upper incomplete gamma to be precise enough. In current 790 -- implementation it loses precision in horrible way (40 791 -- order of magnitude off for sufficiently small p) 792 | a+b > 5 && a/b > 4 = 793 let -- Calculate initial approximation to eta using eq 4.1 794 eta0 = invIncompleteGamma b (1-p) / a 795 mu = b / a -- Eq. 4.3 796 -- A lot of helpers for calculation of 797 w = sqrt(1 + mu) -- Eq. 4.9 798 w_2 = w * w 799 w_3 = w_2 * w 800 w_4 = w_2 * w_2 801 w_5 = w_3 * w_2 802 w_6 = w_3 * w_3 803 w_7 = w_4 * w_3 804 w_8 = w_4 * w_4 805 w_9 = w_5 * w_4 806 w_10 = w_5 * w_5 807 d = eta0 - mu 808 d_2 = d * d 809 d_3 = d_2 * d 810 d_4 = d_2 * d_2 811 w1 = w + 1 812 w1_2 = w1 * w1 813 w1_3 = w1 * w1_2 814 w1_4 = w1_2 * w1_2 815 -- Evaluation of eq 4.10 816 e1 = (w + 2) * (w - 1) / (3 * w) 817 + (w_3 + 9 * w_2 + 21 * w + 5) * d 818 / (36 * w_2 * w1) 819 - (w_4 - 13 * w_3 + 69 * w_2 + 167 * w + 46) * d_2 820 / (1620 * w1_2 * w_3) 821 - (7 * w_5 + 21 * w_4 + 70 * w_3 + 26 * w_2 - 93 * w - 31) * d_3 822 / (6480 * w1_3 * w_4) 823 - (75 * w_6 + 202 * w_5 + 188 * w_4 - 888 * w_3 - 1345 * w_2 + 118 * w + 138) * d_4 824 / (272160 * w1_4 * w_5) 825 e2 = (28 * w_4 + 131 * w_3 + 402 * w_2 + 581 * w + 208) * (w - 1) 826 / (1620 * w1 * w_3) 827 - (35 * w_6 - 154 * w_5 - 623 * w_4 - 1636 * w_3 - 3983 * w_2 - 3514 * w - 925) * d 828 / (12960 * w1_2 * w_4) 829 - ( 2132 * w_7 + 7915 * w_6 + 16821 * w_5 + 35066 * w_4 + 87490 * w_3 830 + 141183 * w_2 + 95993 * w + 21640 831 ) * d_2 832 / (816480 * w_5 * w1_3) 833 - ( 11053 * w_8 + 53308 * w_7 + 117010 * w_6 + 163924 * w_5 + 116188 * w_4 834 - 258428 * w_3 - 677042 * w_2 - 481940 * w - 105497 835 ) * d_3 836 / (14696640 * w1_4 * w_6) 837 e3 = -( (3592 * w_7 + 8375 * w_6 - 1323 * w_5 - 29198 * w_4 - 89578 * w_3 838 - 154413 * w_2 - 116063 * w - 29632 839 ) * (w - 1) 840 ) 841 / (816480 * w_5 * w1_2) 842 - ( 442043 * w_9 + 2054169 * w_8 + 3803094 * w_7 + 3470754 * w_6 + 2141568 * w_5 843 - 2393568 * w_4 - 19904934 * w_3 - 34714674 * w_2 - 23128299 * w - 5253353 844 ) * d 845 / (146966400 * w_6 * w1_3) 846 - ( 116932 * w_10 + 819281 * w_9 + 2378172 * w_8 + 4341330 * w_7 + 6806004 * w_6 847 + 10622748 * w_5 + 18739500 * w_4 + 30651894 * w_3 + 30869976 * w_2 848 + 15431867 * w + 2919016 849 ) * d_2 850 / (146966400 * w1_4 * w_7) 851 eta = evaluatePolynomialL (1/a) [eta0, e1, e2, e3] 852 -- Now we solve eq 4.2 to recover x using Newton iterations 853 u = eta - mu * log eta + (1 + mu) * log(1 + mu) - mu 854 cross = 1 / (1 + mu); 855 lower = if eta < mu then cross else 0 856 upper = if eta < mu then 1 else cross 857 x_guess = (lower + upper) / 2 858 func x = ( u + log x + mu*log(1 - x) 859 , 1/x - mu/(1-x) 860 ) 861 Root x0 = newtonRaphson def{newtonTol=RelTol 1e-8} (lower, x_guess, upper) func 862 in x0 863 -- For large a and b approximation from AS109 (Carter 864 -- approximation). It's reasonably good in this region 865 | a > 1 && b > 1 = 866 let r = (y*y - 3) / 6 867 s = 1 / (2*a - 1) 868 t = 1 / (2*b - 1) 869 h = 2 / (s + t) 870 w = y * sqrt(h + r) / h - (t - s) * (r + 5/6 - 2 / (3 * h)) 871 in a / (a + b * exp(2 * w)) 872 -- Otherwise we revert to approximation from AS64 derived from 873 -- [AS64 2] when it's applicable. 874 -- 875 -- It slightly reduces average number of iterations when `a' and 876 -- `b' have different magnitudes. 877 | chi2 > 0 && ratio > 1 = 1 - 2 / (ratio + 1) 878 -- If all else fails we use approximation from "Numerical 879 -- Recipes". It's very similar to approximations [AS64 4,5] but 880 -- it never goes out of [0,1] interval. 881 | otherwise = case () of 882 _| p < t / w -> (a * p * w) ** (1/a) 883 | otherwise -> 1 - (b * (1 - p) * w) ** (1/b) 884 where 885 lna = log $ a / (a+b) 886 lnb = log $ b / (a+b) 887 t = exp( a * lna ) / a 888 u = exp( b * lnb ) / b 889 w = t + u 890 where 891 -- Formula [AS64 2] 892 ratio = (4*a + 2*b - 2) / chi2 893 -- Quantile of chi-squared distribution. Formula [AS64 3]. 894 chi2 = 2 * b * (1 - t + y * sqrt t) ** 3 895 where 896 t = 1 / (9 * b) 897 -- `y' is Hasting's approximation of p'th quantile of standard 898 -- normal distribution. 899 y = r - ( 2.30753 + 0.27061 * r ) 900 / ( 1.0 + ( 0.99229 + 0.04481 * r ) * r ) 901 where 902 r = sqrt $ - 2 * log p 903 904 905 906---------------------------------------------------------------- 907-- Sinc function 908---------------------------------------------------------------- 909 910-- | Compute sinc function @sin(x)\/x@ 911sinc :: Double -> Double 912sinc x 913 | ax < eps_0 = 1 914 | ax < eps_2 = 1 - x2/6 915 | ax < eps_4 = 1 - x2/6 + x2*x2/120 916 | otherwise = sin x / x 917 where 918 ax = abs x 919 x2 = x*x 920 -- For explanation of choice see `doc/sinc.hs' 921 eps_0 = 1.8250120749944284e-8 -- sqrt (6ε/4) 922 eps_2 = 1.4284346431400855e-4 -- (30ε)**(1/4) / 2 923 eps_4 = 4.043633626430947e-3 -- (1206ε)**(1/6) / 2 924 925 926---------------------------------------------------------------- 927-- Logarithm 928---------------------------------------------------------------- 929 930-- | Compute log(1+x)-x: 931log1pmx :: Double -> Double 932log1pmx x 933 | x < -1 = error "Domain error" 934 | x == -1 = m_neg_inf 935 | ax > 0.95 = log(1 + x) - x 936 | ax < m_epsilon = -(x * x) /2 937 | otherwise = - x * x * sumPowerSeries (-x) (recip <$> enumSequenceFrom 2) 938 where 939 ax = abs x 940 941-- | /O(log n)/ Compute the logarithm in base 2 of the given value. 942log2 :: Int -> Int 943log2 v0 944 | v0 <= 0 = modErr $ "log2: nonpositive input, got " ++ show v0 945 | otherwise = go 5 0 v0 946 where 947 go !i !r !v | i == -1 = r 948 | v .&. b i /= 0 = let si = U.unsafeIndex sv i 949 in go (i-1) (r .|. si) (v `shiftR` si) 950 | otherwise = go (i-1) r v 951 b = U.unsafeIndex bv 952 !bv = U.fromList [ 0x02, 0x0c, 0xf0, 0xff00 953 , fromIntegral (0xffff0000 :: Word) 954 , fromIntegral (0xffffffff00000000 :: Word)] 955 !sv = U.fromList [1,2,4,8,16,32] 956 957 958---------------------------------------------------------------- 959-- Factorial 960---------------------------------------------------------------- 961 962-- | Compute the factorial function /n/!. Returns +∞ if the input is 963-- above 170 (above which the result cannot be represented by a 964-- 64-bit 'Double'). 965factorial :: Int -> Double 966factorial n 967 | n < 0 = error "Numeric.SpecFunctions.factorial: negative input" 968 | n > 170 = m_pos_inf 969 | otherwise = U.unsafeIndex factorialTable n 970 971-- | Compute the natural logarithm of the factorial function. Gives 972-- 16 decimal digits of precision. 973logFactorial :: Integral a => a -> Double 974logFactorial n 975 | n < 0 = error "Numeric.SpecFunctions.logFactorial: negative input" 976 -- For smaller inputs we just look up table 977 | n <= 170 = log $ U.unsafeIndex factorialTable (fromIntegral n) 978 -- Otherwise we use asymptotic Stirling's series. Number of terms 979 -- necessary depends on the argument. 980 | n < 1500 = stirling + rx * ((1/12) - (1/360)*rx*rx) 981 | otherwise = stirling + (1/12)*rx 982 where 983 stirling = (x - 0.5) * log x - x + m_ln_sqrt_2_pi 984 x = fromIntegral n + 1 985 rx = 1 / x 986{-# SPECIALIZE logFactorial :: Int -> Double #-} 987 988 989-- | Calculate the error term of the Stirling approximation. This is 990-- only defined for non-negative values. 991-- 992-- \[ 993-- \operatorname{stirlingError}(n) = \log(n!) - \log(\sqrt{2\pi n}\frac{n}{e}^n) 994-- \] 995stirlingError :: Double -> Double 996stirlingError n 997 | n <= 15.0 = case properFraction (n+n) of 998 (i,0) -> sfe `U.unsafeIndex` i 999 _ -> logGamma (n+1.0) - (n+0.5) * log n + n - 1000 m_ln_sqrt_2_pi 1001 | n > 500 = evaluateOddPolynomialL (1/n) [s0,-s1] 1002 | n > 80 = evaluateOddPolynomialL (1/n) [s0,-s1,s2] 1003 | n > 35 = evaluateOddPolynomialL (1/n) [s0,-s1,s2,-s3] 1004 | otherwise = evaluateOddPolynomialL (1/n) [s0,-s1,s2,-s3,s4] 1005 where 1006 s0 = 0.083333333333333333333 -- 1/12 1007 s1 = 0.00277777777777777777778 -- 1/360 1008 s2 = 0.00079365079365079365079365 -- 1/1260 1009 s3 = 0.000595238095238095238095238 -- 1/1680 1010 s4 = 0.0008417508417508417508417508 -- 1/1188 1011 sfe = U.fromList [ 0.0, 1012 0.1534264097200273452913848, 0.0810614667953272582196702, 1013 0.0548141210519176538961390, 0.0413406959554092940938221, 1014 0.03316287351993628748511048, 0.02767792568499833914878929, 1015 0.02374616365629749597132920, 0.02079067210376509311152277, 1016 0.01848845053267318523077934, 0.01664469118982119216319487, 1017 0.01513497322191737887351255, 0.01387612882307074799874573, 1018 0.01281046524292022692424986, 0.01189670994589177009505572, 1019 0.01110455975820691732662991, 0.010411265261972096497478567, 1020 0.009799416126158803298389475, 0.009255462182712732917728637, 1021 0.008768700134139385462952823, 0.008330563433362871256469318, 1022 0.007934114564314020547248100, 0.007573675487951840794972024, 1023 0.007244554301320383179543912, 0.006942840107209529865664152, 1024 0.006665247032707682442354394, 0.006408994188004207068439631, 1025 0.006171712263039457647532867, 0.005951370112758847735624416, 1026 0.005746216513010115682023589, 0.005554733551962801371038690 ] 1027 1028 1029---------------------------------------------------------------- 1030-- Combinatorics 1031---------------------------------------------------------------- 1032 1033-- | 1034-- Quickly compute the natural logarithm of /n/ @`choose`@ /k/, with 1035-- no checking. 1036-- 1037-- Less numerically stable: 1038-- 1039-- > exp $ lg (n+1) - lg (k+1) - lg (n-k+1) 1040-- > where lg = logGamma . fromIntegral 1041logChooseFast :: Double -> Double -> Double 1042logChooseFast n k = -log (n + 1) - logBeta (n - k + 1) (k + 1) 1043 1044-- | Calculate binomial coefficient using exact formula 1045chooseExact :: Int -> Int -> Double 1046n `chooseExact` k 1047 = U.foldl' go 1 $ U.enumFromTo 1 k 1048 where 1049 go a i = a * (nk + j) / j 1050 where j = fromIntegral i :: Double 1051 nk = fromIntegral (n - k) 1052 1053-- | Compute logarithm of the binomial coefficient. 1054logChoose :: Int -> Int -> Double 1055n `logChoose` k 1056 | k > n = (-1) / 0 1057 -- For very large N exact algorithm overflows double so we 1058 -- switch to beta-function based one 1059 | k' < 50 && (n < 20000000) = log $ chooseExact n k' 1060 | otherwise = logChooseFast (fromIntegral n) (fromIntegral k) 1061 where 1062 k' = min k (n-k) 1063 1064-- | Compute the binomial coefficient /n/ @\``choose`\`@ /k/. For 1065-- values of /k/ > 50, this uses an approximation for performance 1066-- reasons. The approximation is accurate to 12 decimal places in the 1067-- worst case 1068-- 1069-- Example: 1070-- 1071-- > 7 `choose` 3 == 35 1072choose :: Int -> Int -> Double 1073n `choose` k 1074 | k > n = 0 1075 | k' < 50 = chooseExact n k' 1076 | approx < max64 = fromIntegral . round64 $ approx 1077 | otherwise = approx 1078 where 1079 k' = min k (n-k) 1080 approx = exp $ logChooseFast (fromIntegral n) (fromIntegral k') 1081 max64 = fromIntegral (maxBound :: Int64) 1082 round64 x = round x :: Int64 1083 1084-- | Compute ψ(/x/), the first logarithmic derivative of the gamma 1085-- function. 1086-- 1087-- \[ 1088-- \psi(x) = \frac{d}{dx} \ln \left(\Gamma(x)\right) = \frac{\Gamma'(x)}{\Gamma(x)} 1089-- \] 1090-- 1091-- Uses Algorithm AS 103 by Bernardo, based on Minka's C implementation. 1092digamma :: Double -> Double 1093digamma x 1094 | isNaN x || isInfinite x = m_NaN 1095 -- FIXME: 1096 -- This is ugly. We are testing here that number is in fact 1097 -- integer. It's somewhat tricky question to answer. When ε for 1098 -- given number becomes 1 or greater every number is represents 1099 -- an integer. We also must make sure that excess precision 1100 -- won't bite us. 1101 | x <= 0 && fromIntegral (truncate x :: Int64) == x = m_neg_inf 1102 -- Jeffery's reflection formula 1103 | x < 0 = digamma (1 - x) + pi / tan (negate pi * x) 1104 | x <= 1e-6 = - γ - 1/x + trigamma1 * x 1105 | x' < c = r 1106 -- De Moivre's expansion 1107 | otherwise = let s = 1/x' 1108 in evaluateEvenPolynomialL s 1109 [ r + log x' - 0.5 * s 1110 , - 1/12 1111 , 1/120 1112 , - 1/252 1113 , 1/240 1114 , - 1/132 1115 , 391/32760 1116 ] 1117 where 1118 γ = m_eulerMascheroni 1119 c = 12 1120 -- Reduce to digamma (x + n) where (x + n) >= c 1121 (r, x') = reduce 0 x 1122 where 1123 reduce !s y 1124 | y < c = reduce (s - 1 / y) (y + 1) 1125 | otherwise = (s, y) 1126 1127 1128 1129---------------------------------------------------------------- 1130-- Constants 1131---------------------------------------------------------------- 1132 1133-- Coefficients for 18-point Gauss-Legendre integration. They are 1134-- used in implementation of incomplete gamma and beta functions. 1135coefW,coefY :: U.Vector Double 1136coefW = U.fromList [ 0.0055657196642445571, 0.012915947284065419, 0.020181515297735382 1137 , 0.027298621498568734, 0.034213810770299537, 0.040875750923643261 1138 , 0.047235083490265582, 0.053244713977759692, 0.058860144245324798 1139 , 0.064039797355015485, 0.068745323835736408, 0.072941885005653087 1140 , 0.076598410645870640, 0.079687828912071670, 0.082187266704339706 1141 , 0.084078218979661945, 0.085346685739338721, 0.085983275670394821 1142 ] 1143coefY = U.fromList [ 0.0021695375159141994, 0.011413521097787704, 0.027972308950302116 1144 , 0.051727015600492421, 0.082502225484340941, 0.12007019910960293 1145 , 0.16415283300752470, 0.21442376986779355, 0.27051082840644336 1146 , 0.33199876341447887, 0.39843234186401943, 0.46931971407375483 1147 , 0.54413605556657973, 0.62232745288031077, 0.70331500465597174 1148 , 0.78649910768313447, 0.87126389619061517, 0.95698180152629142 1149 ] 1150{-# NOINLINE coefW #-} 1151{-# NOINLINE coefY #-} 1152 1153trigamma1 :: Double 1154trigamma1 = 1.6449340668482264365 -- pi**2 / 6 1155 1156modErr :: String -> a 1157modErr msg = error $ "Numeric.SpecFunctions." ++ msg 1158 1159factorialTable :: U.Vector Double 1160{-# NOINLINE factorialTable #-} 1161factorialTable = U.fromListN 171 1162 [ 1.0 1163 , 1.0 1164 , 2.0 1165 , 6.0 1166 , 24.0 1167 , 120.0 1168 , 720.0 1169 , 5040.0 1170 , 40320.0 1171 , 362880.0 1172 , 3628800.0 1173 , 3.99168e7 1174 , 4.790016e8 1175 , 6.2270208e9 1176 , 8.71782912e10 1177 , 1.307674368e12 1178 , 2.0922789888e13 1179 , 3.55687428096e14 1180 , 6.402373705728e15 1181 , 1.21645100408832e17 1182 , 2.43290200817664e18 1183 , 5.109094217170944e19 1184 , 1.1240007277776077e21 1185 , 2.5852016738884974e22 1186 , 6.204484017332394e23 1187 , 1.5511210043330984e25 1188 , 4.032914611266056e26 1189 , 1.0888869450418352e28 1190 , 3.0488834461171384e29 1191 , 8.841761993739702e30 1192 , 2.6525285981219103e32 1193 , 8.222838654177922e33 1194 , 2.631308369336935e35 1195 , 8.683317618811886e36 1196 , 2.9523279903960412e38 1197 , 1.0333147966386144e40 1198 , 3.719933267899012e41 1199 , 1.3763753091226343e43 1200 , 5.23022617466601e44 1201 , 2.0397882081197442e46 1202 , 8.159152832478977e47 1203 , 3.3452526613163803e49 1204 , 1.4050061177528798e51 1205 , 6.041526306337383e52 1206 , 2.6582715747884485e54 1207 , 1.1962222086548019e56 1208 , 5.5026221598120885e57 1209 , 2.5862324151116818e59 1210 , 1.2413915592536073e61 1211 , 6.082818640342675e62 1212 , 3.0414093201713376e64 1213 , 1.5511187532873822e66 1214 , 8.065817517094388e67 1215 , 4.2748832840600255e69 1216 , 2.308436973392414e71 1217 , 1.2696403353658275e73 1218 , 7.109985878048634e74 1219 , 4.0526919504877214e76 1220 , 2.3505613312828785e78 1221 , 1.386831185456898e80 1222 , 8.32098711274139e81 1223 , 5.075802138772247e83 1224 , 3.146997326038793e85 1225 , 1.9826083154044399e87 1226 , 1.2688693218588415e89 1227 , 8.24765059208247e90 1228 , 5.44344939077443e92 1229 , 3.647111091818868e94 1230 , 2.4800355424368305e96 1231 , 1.711224524281413e98 1232 , 1.197857166996989e100 1233 , 8.504785885678623e101 1234 , 6.1234458376886085e103 1235 , 4.470115461512684e105 1236 , 3.307885441519386e107 1237 , 2.4809140811395396e109 1238 , 1.88549470166605e111 1239 , 1.4518309202828586e113 1240 , 1.1324281178206297e115 1241 , 8.946182130782974e116 1242 , 7.15694570462638e118 1243 , 5.797126020747368e120 1244 , 4.753643337012841e122 1245 , 3.9455239697206583e124 1246 , 3.314240134565353e126 1247 , 2.81710411438055e128 1248 , 2.422709538367273e130 1249 , 2.1077572983795275e132 1250 , 1.8548264225739844e134 1251 , 1.650795516090846e136 1252 , 1.4857159644817613e138 1253 , 1.352001527678403e140 1254 , 1.2438414054641305e142 1255 , 1.1567725070816416e144 1256 , 1.087366156656743e146 1257 , 1.0329978488239058e148 1258 , 9.916779348709496e149 1259 , 9.619275968248211e151 1260 , 9.426890448883246e153 1261 , 9.332621544394413e155 1262 , 9.332621544394415e157 1263 , 9.425947759838358e159 1264 , 9.614466715035125e161 1265 , 9.902900716486179e163 1266 , 1.0299016745145626e166 1267 , 1.0813967582402908e168 1268 , 1.1462805637347082e170 1269 , 1.2265202031961378e172 1270 , 1.3246418194518288e174 1271 , 1.4438595832024934e176 1272 , 1.5882455415227428e178 1273 , 1.7629525510902446e180 1274 , 1.974506857221074e182 1275 , 2.2311927486598134e184 1276 , 2.543559733472187e186 1277 , 2.9250936934930154e188 1278 , 3.393108684451898e190 1279 , 3.9699371608087206e192 1280 , 4.68452584975429e194 1281 , 5.574585761207606e196 1282 , 6.689502913449126e198 1283 , 8.094298525273443e200 1284 , 9.875044200833601e202 1285 , 1.214630436702533e205 1286 , 1.5061417415111406e207 1287 , 1.8826771768889257e209 1288 , 2.372173242880047e211 1289 , 3.0126600184576594e213 1290 , 3.856204823625804e215 1291 , 4.974504222477286e217 1292 , 6.466855489220473e219 1293 , 8.471580690878819e221 1294 , 1.1182486511960041e224 1295 , 1.4872707060906857e226 1296 , 1.9929427461615188e228 1297 , 2.6904727073180504e230 1298 , 3.6590428819525483e232 1299 , 5.012888748274991e234 1300 , 6.917786472619488e236 1301 , 9.615723196941088e238 1302 , 1.3462012475717523e241 1303 , 1.898143759076171e243 1304 , 2.6953641378881624e245 1305 , 3.8543707171800725e247 1306 , 5.5502938327393044e249 1307 , 8.047926057471992e251 1308 , 1.1749972043909107e254 1309 , 1.7272458904546386e256 1310 , 2.5563239178728654e258 1311 , 3.808922637630569e260 1312 , 5.713383956445854e262 1313 , 8.62720977423324e264 1314 , 1.3113358856834524e267 1315 , 2.0063439050956823e269 1316 , 3.0897696138473508e271 1317 , 4.789142901463393e273 1318 , 7.471062926282894e275 1319 , 1.1729568794264143e278 1320 , 1.8532718694937346e280 1321 , 2.946702272495038e282 1322 , 4.714723635992061e284 1323 , 7.590705053947218e286 1324 , 1.2296942187394494e289 1325 , 2.0044015765453023e291 1326 , 3.287218585534296e293 1327 , 5.423910666131589e295 1328 , 9.003691705778436e297 1329 , 1.5036165148649988e300 1330 , 2.526075744973198e302 1331 , 4.269068009004705e304 1332 , 7.257415615307998e306 1333 ] 1334 1335 1336-- [NOTE: incompleteGamma.taylorP] 1337-- 1338-- Incompltete gamma uses several algorithms for different parts of 1339-- parameter space. Most troublesome is P(a,x) Taylor series 1340-- [Temme1994,Eq.5.5] which requires to evaluate rather nasty 1341-- expression: 1342-- 1343-- x^a x^a 1344-- ------------- = ------------- 1345-- exp(x)·Γ(a+1) exp(x)·a·Γ(a) 1346-- 1347-- Conditions: 1348-- | 0.5<x<1.1 = x < 4/3*a 1349-- | otherwise = x < a 1350-- 1351-- For small `a` computation could be performed directly. However for 1352-- largish values of `a` it's possible some of factor in the 1353-- expression overflow. Values below take into account ranges for 1354-- Taylor P approximation: 1355-- 1356-- · a > 155 - x^a could overflow 1357-- · a > 1182.5 - exp(x) could overflow 1358-- 1359-- Usual way to avoid overflow problem is to perform calculations in 1360-- the log domain. It however doesn't work very well in this case 1361-- since we encounter catastrophic cancellations and could easily lose 1362-- up to 6(!) digits for large `a`. 1363-- 1364-- So we take another approach and use Stirling approximation with 1365-- correction (logGammaCorrection). 1366-- 1367-- x^a / x·e \^a 1 1368-- ≈ ------------------------- = | --- | · ---------------- 1369-- exp(x)·sqrt(2πa)·(a/e)^a) \ a / exp(x)·sqrt(2πa) 1370-- 1371-- We're using this approach as soon as logGammaCorrection starts 1372-- working (a>10) because we don't have implementation for gamma 1373-- function and exp(logGamma z) results in errors for large a. 1374-- 1375-- Once we get into region when exp(x) could overflow we rewrite 1376-- expression above once more: 1377-- 1378-- / x·e \^a 1 1379-- | --- · e^(-x/a) | · --------- 1380-- \ a / sqrt(2πa) 1381-- 1382-- This approach doesn't work very well but it's still big improvement 1383-- over calculations in the log domain. 1384