    -- |
    -- 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