1-- |
2-- Module    : Numeric.SpecFunctions.Extra
3-- Copyright : (c) 2009, 2011 Bryan O'Sullivan
4-- License   : BSD3
5--
6-- Maintainer  : bos@serpentine.com
7-- Stability   : experimental
8-- Portability : portable
9--
10-- Less common mathematical functions.
11module Numeric.SpecFunctions.Extra (
12    bd0
13  , chooseExact
14  , logChooseFast
15  , logGammaAS245
16  , logGammaCorrection
17  ) where
18
19import Numeric.MathFunctions.Constants (m_NaN,m_pos_inf)
20import Numeric.SpecFunctions.Internal  (chooseExact,logChooseFast,logGammaCorrection)
21
22-- | Evaluate the deviance term @x log(x/np) + np - x@.
23bd0 :: Double                   -- ^ @x@
24    -> Double                   -- ^ @np@
25    -> Double
26bd0 x np
27  | isInfinite x || isInfinite np || np == 0 = m_NaN
28  | abs x_np >= 0.1*(x+np)                   = x * log (x/np) - x_np
29  | otherwise                                = loop 1 (ej0*vv) s0
30  where
31    x_np = x - np
32    v    = x_np / (x+np)
33    s0   = x_np * v
34    ej0  = 2*x*v
35    vv   = v*v
36    loop j ej s = case s + ej/(2*j+1) of
37                    s' | s' == s   -> s'  -- FIXME: Comparing Doubles for equality!
38                       | otherwise -> loop (j+1) (ej*vv) s'
39
40
41
42-- | Compute the logarithm of the gamma function Γ(/x/).  Uses
43-- Algorithm AS 245 by Macleod.
44--
45-- Gives an accuracy of 10-12 significant decimal digits, except
46-- for small regions around /x/ = 1 and /x/ = 2, where the function
47-- goes to zero.  For greater accuracy, use 'logGammaL'.
48--
49-- Returns ∞ if the input is outside of the range (0 < /x/ ≤ 1e305).
50logGammaAS245 :: Double -> Double
51-- Adapted from http://people.sc.fsu.edu/~burkardt/f_src/asa245/asa245.html
52logGammaAS245 x
53    | x <= 0    = m_pos_inf
54    -- Handle positive infinity. logGamma overflows before 1e308 so
55    -- it's safe
56    | x > 1e308 = m_pos_inf
57    -- Normal cases
58    | x < 1.5   = a + c *
59                  ((((r1_4 * b + r1_3) * b + r1_2) * b + r1_1) * b + r1_0) /
60                  ((((b + r1_8) * b + r1_7) * b + r1_6) * b + r1_5)
61    | x < 4     = (x - 2) *
62                  ((((r2_4 * x + r2_3) * x + r2_2) * x + r2_1) * x + r2_0) /
63                  ((((x + r2_8) * x + r2_7) * x + r2_6) * x + r2_5)
64    | x < 12    = ((((r3_4 * x + r3_3) * x + r3_2) * x + r3_1) * x + r3_0) /
65                  ((((x + r3_8) * x + r3_7) * x + r3_6) * x + r3_5)
66    | x > 3e6   = k
67    | otherwise = k + x1 *
68                  ((r4_2 * x2 + r4_1) * x2 + r4_0) /
69                  ((x2 + r4_4) * x2 + r4_3)
70  where
71    (a , b , c)
72        | x < 0.5   = (-y , x + 1 , x)
73        | otherwise = (0  , x     , x - 1)
74
75    y      = log x
76    k      = x * (y-1) - 0.5 * y + alr2pi
77    alr2pi = 0.918938533204673
78
79    x1 = 1 / x
80    x2 = x1 * x1
81
82    r1_0 =  -2.66685511495;   r1_1 =  -24.4387534237;    r1_2 = -21.9698958928
83    r1_3 =  11.1667541262;    r1_4 =    3.13060547623;   r1_5 =   0.607771387771
84    r1_6 =  11.9400905721;    r1_7 =   31.4690115749;    r1_8 =  15.2346874070
85
86    r2_0 = -78.3359299449;    r2_1 = -142.046296688;     r2_2 = 137.519416416
87    r2_3 =  78.6994924154;    r2_4 =    4.16438922228;   r2_5 =  47.0668766060
88    r2_6 = 313.399215894;     r2_7 =  263.505074721;     r2_8 =  43.3400022514
89
90    r3_0 =  -2.12159572323e5; r3_1 =    2.30661510616e5; r3_2 =   2.74647644705e4
91    r3_3 =  -4.02621119975e4; r3_4 =   -2.29660729780e3; r3_5 =  -1.16328495004e5
92    r3_6 =  -1.46025937511e5; r3_7 =   -2.42357409629e4; r3_8 =  -5.70691009324e2
93
94    r4_0 = 0.279195317918525;  r4_1 = 0.4917317610505968;
95    r4_2 = 0.0692910599291889; r4_3 = 3.350343815022304
96    r4_4 = 6.012459259764103
97