-- |
-- Module    : Numeric.SpecFunctions.Extra
-- Copyright : (c) 2009, 2011 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Less common mathematical functions.
module Numeric.SpecFunctions.Extra (
    bd0
  , chooseExact
  , logChooseFast
  , logGammaAS245
  , logGammaCorrection
  ) where

import Numeric.MathFunctions.Constants (m_NaN,m_pos_inf)
import Numeric.SpecFunctions.Internal  (chooseExact,logChooseFast,logGammaCorrection)

-- | Evaluate the deviance term @x log(x/np) + np - x@.
bd0 :: Double                   -- ^ @x@
    -> Double                   -- ^ @np@
    -> Double 
bd0 :: Double -> Double -> Double
bd0 Double
x Double
np 
  | Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
np Bool -> Bool -> Bool
|| Double
np Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 = Double
m_NaN
  | Double -> Double
forall a. Num a => a -> a
abs Double
x_np Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
0.1Double -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
np)                   = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log (Double
xDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
np) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x_np
  | Bool
otherwise                                = Double -> Double -> Double -> Double
loop Double
1 (Double
ej0Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
vv) Double
s0
  where 
    x_np :: Double
x_np = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
np
    v :: Double
v    = Double
x_np Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
np)
    s0 :: Double
s0   = Double
x_np Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
v
    ej0 :: Double
ej0  = Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
v
    vv :: Double
vv   = Double
vDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
v
    loop :: Double -> Double -> Double -> Double
loop Double
j Double
ej Double
s = case Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
ejDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
jDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1) of
                    Double
s' | Double
s' Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
s   -> Double
s'  -- FIXME: Comparing Doubles for equality!
                       | Bool
otherwise -> Double -> Double -> Double -> Double
loop (Double
jDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1) (Double
ejDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
vv) Double
s'



-- | Compute the logarithm of the gamma function Γ(/x/).  Uses
-- Algorithm AS 245 by Macleod.
--
-- 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).
logGammaAS245 :: Double -> Double
-- Adapted from http://people.sc.fsu.edu/~burkardt/f_src/asa245/asa245.html
logGammaAS245 :: Double -> Double
logGammaAS245 Double
x
    | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0    = Double
m_pos_inf
    -- Handle positive infinity. logGamma overflows before 1e308 so
    -- it's safe
    | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1e308 = Double
m_pos_inf
    -- Normal cases
    | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1.5   = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
c Double -> Double -> Double
forall a. Num a => a -> a -> a
*
                  ((((Double
r1_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r1_3) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r1_2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r1_1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r1_0) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/
                  ((((Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r1_8) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r1_7) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r1_6) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r1_5)
    | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
4     = (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2) Double -> Double -> Double
forall a. Num a => a -> a -> a
*
                  ((((Double
r2_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r2_3) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r2_2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r2_1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r2_0) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/
                  ((((Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r2_8) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r2_7) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r2_6) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r2_5)
    | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
12    = ((((Double
r3_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r3_3) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r3_2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r3_1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r3_0) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/
                  ((((Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r3_8) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r3_7) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r3_6) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r3_5)
    | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
3e6   = Double
k
    | Bool
otherwise = Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x1 Double -> Double -> Double
forall a. Num a => a -> a -> a
*
                  ((Double
r4_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r4_1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r4_0) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/
                  ((Double
x2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r4_4) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r4_3)
  where
    (Double
a , Double
b , Double
c)
        | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.5   = (-Double
y , Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1 , Double
x)
        | Bool
otherwise = (Double
0  , Double
x     , Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)

    y :: Double
y      = Double -> Double
forall a. Floating a => a -> a
log Double
x
    k :: Double
k      = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
alr2pi
    alr2pi :: Double
alr2pi = Double
0.918938533204673

    x1 :: Double
x1 = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x
    x2 :: Double
x2 = Double
x1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x1

    r1_0 :: Double
r1_0 =  -Double
2.66685511495;   r1_1 :: Double
r1_1 =  -Double
24.4387534237;    r1_2 :: Double
r1_2 = -Double
21.9698958928
    r1_3 :: Double
r1_3 =  Double
11.1667541262;    r1_4 :: Double
r1_4 =    Double
3.13060547623;   r1_5 :: Double
r1_5 =   Double
0.607771387771
    r1_6 :: Double
r1_6 =  Double
11.9400905721;    r1_7 :: Double
r1_7 =   Double
31.4690115749;    r1_8 :: Double
r1_8 =  Double
15.2346874070

    r2_0 :: Double
r2_0 = -Double
78.3359299449;    r2_1 :: Double
r2_1 = -Double
142.046296688;     r2_2 :: Double
r2_2 = Double
137.519416416
    r2_3 :: Double
r2_3 =  Double
78.6994924154;    r2_4 :: Double
r2_4 =    Double
4.16438922228;   r2_5 :: Double
r2_5 =  Double
47.0668766060
    r2_6 :: Double
r2_6 = Double
313.399215894;     r2_7 :: Double
r2_7 =  Double
263.505074721;     r2_8 :: Double
r2_8 =  Double
43.3400022514

    r3_0 :: Double
r3_0 =  -Double
2.12159572323e5; r3_1 :: Double
r3_1 =    Double
2.30661510616e5; r3_2 :: Double
r3_2 =   Double
2.74647644705e4
    r3_3 :: Double
r3_3 =  -Double
4.02621119975e4; r3_4 :: Double
r3_4 =   -Double
2.29660729780e3; r3_5 :: Double
r3_5 =  -Double
1.16328495004e5
    r3_6 :: Double
r3_6 =  -Double
1.46025937511e5; r3_7 :: Double
r3_7 =   -Double
2.42357409629e4; r3_8 :: Double
r3_8 =  -Double
5.70691009324e2

    r4_0 :: Double
r4_0 = Double
0.279195317918525;  r4_1 :: Double
r4_1 = Double
0.4917317610505968;
    r4_2 :: Double
r4_2 = Double
0.0692910599291889; r4_3 :: Double
r4_3 = Double
3.350343815022304
    r4_4 :: Double
r4_4 = Double
6.012459259764103