{-# LANGUAGE BangPatterns, ScopedTypeVariables #-}
{-# OPTIONS_HADDOCK hide #-}
-- |
-- Module    : Numeric.SpecFunctions.Internal
-- Copyright : (c) 2009, 2011, 2012 Bryan O'Sullivan
-- License   : BSD3
--
-- Maintainer  : bos@serpentine.com
-- Stability   : experimental
-- Portability : portable
--
-- Internal module with implementation of special functions.
module Numeric.SpecFunctions.Internal
    ( module Numeric.SpecFunctions.Internal
    , Compat.log1p
    , Compat.expm1
    ) where

import Data.Bits          ((.&.), (.|.), shiftR)
import Data.Int           (Int64)
import Data.Default.Class
import qualified Data.Vector.Unboxed as U
import           Data.Vector.Unboxed   ((!))
import Text.Printf

import Numeric.Polynomial.Chebyshev    (chebyshevBroucke)
import Numeric.Polynomial              (evaluatePolynomial, evaluatePolynomialL, evaluateEvenPolynomialL
                                       ,evaluateOddPolynomialL)
import Numeric.RootFinding             (Root(..), newtonRaphson, NewtonParam(..), Tolerance(..))
import Numeric.Series
import Numeric.MathFunctions.Constants
import Numeric.SpecFunctions.Compat (log1p)
import qualified Numeric.SpecFunctions.Compat as Compat

----------------------------------------------------------------
-- Error function
----------------------------------------------------------------

-- | Error function.
--
-- \[
-- \operatorname{erf}(x) = \frac{2}{\sqrt{\pi}} \int_{0}^{x} \exp(-t^2) dt
-- \]
--
-- Function limits are:
--
-- \[
-- \begin{aligned}
--  &\operatorname{erf}(-\infty) &=& -1 \\
--  &\operatorname{erf}(0)       &=& \phantom{-}\,0 \\
--  &\operatorname{erf}(+\infty) &=& \phantom{-}\,1 \\
-- \end{aligned}
-- \]
erf :: Double -> Double
erf :: Double -> Double
erf = Double -> Double
Compat.erf
{-# INLINE erf #-}

-- | Complementary error function.
--
-- \[
-- \operatorname{erfc}(x) = 1 - \operatorname{erf}(x)
-- \]
--
-- Function limits are:
--
-- \[
-- \begin{aligned}
--  &\operatorname{erf}(-\infty) &=&\, 2 \\
--  &\operatorname{erf}(0)       &=&\, 1 \\
--  &\operatorname{erf}(+\infty) &=&\, 0 \\
-- \end{aligned}
-- \]
erfc :: Double -> Double
erfc :: Double -> Double
erfc = Double -> Double
Compat.erfc
{-# INLINE erfc #-}

-- | Inverse of 'erf'.
invErf :: Double -- ^ /p/ ∈ [-1,1]
       -> Double
invErf :: Double -> Double
invErf Double
p
  | Double
p Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
==  Double
1         = Double
m_pos_inf
  | Double
p Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== -Double
1         = Double
m_neg_inf
  | Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1 Bool -> Bool -> Bool
&& Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> -Double
1 = if Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 then Double
r else -Double
r
  | Bool
otherwise       = String -> Double
forall a. HasCallStack => String -> a
error String
"invErf: p must in [-1,1] range"
  where
    -- We solve equation same as in invErfc. We're able to ruse same
    -- Halley step by solving equation:
    --   > pp - erf x = 0
    -- instead of
    --   > erf x - pp = 0
    pp :: Double
pp     = Double -> Double
forall a. Num a => a -> a
abs Double
p
    r :: Double
r      = Double -> Double
step (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double -> Double
step (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double -> Double
guessInvErfc (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
pp
    step :: Double -> Double
step Double
x = Double -> Double -> Double
invErfcHalleyStep (Double
pp Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
erf Double
x) Double
x

-- | Inverse of 'erfc'.
invErfc :: Double -- ^ /p/ ∈ [0,2]
        -> Double
invErfc :: Double -> Double
invErfc Double
p
  | Double
p Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
2        = Double
m_neg_inf
  | Double
p Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0        = Double
m_pos_inf
  | Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>Double
0 Bool -> Bool -> Bool
&& Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
2 = if Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1 then Double
r else -Double
r
  | Bool
otherwise     = String -> Double
forall a. String -> a
modErr (String -> Double) -> String -> Double
forall a b. (a -> b) -> a -> b
$ String
"invErfc: p must be in [0,2] got " String -> String -> String
forall a. [a] -> [a] -> [a]
++ Double -> String
forall a. Show a => a -> String
show Double
p
  where
    pp :: Double
pp | Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1    = Double
p
       | Bool
otherwise = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p
    -- We perform 2 Halley steps in order to get to solution
    r :: Double
r      = Double -> Double
step (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double -> Double
step (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double -> Double
guessInvErfc Double
pp
    step :: Double -> Double
step Double
x = Double -> Double -> Double
invErfcHalleyStep (Double -> Double
erfc Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
pp) Double
x

-- Initial guess for invErfc & invErf
guessInvErfc :: Double -> Double
guessInvErfc :: Double -> Double
guessInvErfc Double
p
  = -Double
0.70711 Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((Double
2.30753 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
0.27061) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
0.99229 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
0.04481)) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t)
  where
    t :: Double
t = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ -Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log( Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
p)

-- Halley step for solving invErfc
invErfcHalleyStep :: Double -> Double -> Double
invErfcHalleyStep :: Double -> Double -> Double
invErfcHalleyStep Double
err Double
x
  = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
err Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1.12837916709551257 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp(-Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
err)

----------------------------------------------------------------
-- Gamma function
----------------------------------------------------------------

data L = L {-# UNPACK #-} !Double {-# UNPACK #-} !Double

-- | Compute the logarithm of the gamma function, Γ(/x/).
--
-- \[
-- \Gamma(x) = \int_0^{\infty}t^{x-1}e^{-t}\,dt = (x - 1)!
-- \]
--
-- This implementation uses Lanczos approximation. It gives 14 or more
-- significant decimal digits, except around /x/ = 1 and /x/ = 2,
-- where the function goes to zero.
--
-- Returns &#8734; if the input is outside of the range (0 < /x/
-- &#8804; 1e305).
logGamma :: Double -> Double
logGamma :: Double -> Double
logGamma Double
z
  | Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0    = Double
m_pos_inf
  -- For very small values z we can just use Laurent expansion
  | Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
m_sqrt_eps = Double -> Double
forall a. Floating a => a -> a
log (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
m_eulerMascheroni)
  -- For z<1 we use recurrence. Γ(z+1) = z·Γ(z) Note that in order to
  -- avoid precision loss we have to compute parameter to
  -- approximations here:
  --
  -- > (z + 1) - 1 = z
  -- > (z + 1) - 2 = z - 1
  --
  -- Simple passing (z + 1) to piecewise approximations and computing
  -- difference leads to bad loss of precision near 1.
  -- This is reason lgamma1_15 & lgamma15_2 have three parameters
  | Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.5   = Double -> Double -> Double
lgamma1_15 Double
z (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
log Double
z
  | Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1     = Double -> Double -> Double
lgamma15_2 Double
z (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
log Double
z
  -- Piecewise polynomial approximations
  | Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1.5  = Double -> Double -> Double
lgamma1_15 (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2)
  | Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
2     = Double -> Double -> Double
lgamma15_2 (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2)
  | Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
15    = Double -> Double
lgammaSmall Double
z
  -- Otherwise we switch to Lanczos approximation
  | Bool
otherwise = Double -> Double
lanczosApprox Double
z


-- | Synonym for 'logGamma'. Retained for compatibility
logGammaL :: Double -> Double
logGammaL :: Double -> Double
logGammaL = Double -> Double
logGamma
{-# DEPRECATED logGammaL "Use logGamma instead" #-}



-- Polynomial expansion used in interval (1,1.5]
--
-- > logΓ(z) = (z-1)(z-2)(Y + R(z-1))
lgamma1_15 :: Double -> Double -> Double
lgamma1_15 :: Double -> Double -> Double
lgamma1_15 Double
zm1 Double
zm2
   = Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* ( Double -> Vector Double -> Double
forall (v :: * -> *) a. (Vector v a, Num a) => a -> v a -> a
evaluatePolynomial Double
zm1 Vector Double
tableLogGamma_1_15P
                 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Vector Double -> Double
forall (v :: * -> *) a. (Vector v a, Num a) => a -> v a -> a
evaluatePolynomial Double
zm1 Vector Double
tableLogGamma_1_15Q
                 )
   where
     r :: Double
r = Double
zm1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
zm2
     y :: Double
y = Double
0.52815341949462890625

tableLogGamma_1_15P,tableLogGamma_1_15Q :: U.Vector Double
tableLogGamma_1_15P :: Vector Double
tableLogGamma_1_15P = [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
U.fromList
  [  Double
0.490622454069039543534e-1
  , -Double
0.969117530159521214579e-1
  , -Double
0.414983358359495381969e0
  , -Double
0.406567124211938417342e0
  , -Double
0.158413586390692192217e0
  , -Double
0.240149820648571559892e-1
  , -Double
0.100346687696279557415e-2
  ]
{-# NOINLINE tableLogGamma_1_15P #-}
tableLogGamma_1_15Q :: Vector Double
tableLogGamma_1_15Q = [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
U.fromList
  [ Double
1
  , Double
0.302349829846463038743e1
  , Double
0.348739585360723852576e1
  , Double
0.191415588274426679201e1
  , Double
0.507137738614363510846e0
  , Double
0.577039722690451849648e-1
  , Double
0.195768102601107189171e-2
  ]
{-# NOINLINE tableLogGamma_1_15Q #-}



-- Polynomial expansion used in interval (1.5,2)
--
-- > logΓ(z) = (2-z)(1-z)(Y + R(2-z))
lgamma15_2 :: Double -> Double -> Double
lgamma15_2 :: Double -> Double -> Double
lgamma15_2 Double
zm1 Double
zm2
   = Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* ( Double -> Vector Double -> Double
forall (v :: * -> *) a. (Vector v a, Num a) => a -> v a -> a
evaluatePolynomial (-Double
zm2) Vector Double
tableLogGamma_15_2P
                 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Vector Double -> Double
forall (v :: * -> *) a. (Vector v a, Num a) => a -> v a -> a
evaluatePolynomial (-Double
zm2) Vector Double
tableLogGamma_15_2Q
                 )
   where
     r :: Double
r = Double
zm1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
zm2
     y :: Double
y = Double
0.452017307281494140625

tableLogGamma_15_2P,tableLogGamma_15_2Q :: U.Vector Double
tableLogGamma_15_2P :: Vector Double
tableLogGamma_15_2P = [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
U.fromList
  [ -Double
0.292329721830270012337e-1
  ,  Double
0.144216267757192309184e0
  , -Double
0.142440390738631274135e0
  ,  Double
0.542809694055053558157e-1
  , -Double
0.850535976868336437746e-2
  ,  Double
0.431171342679297331241e-3
  ]
{-# NOINLINE tableLogGamma_15_2P #-}
tableLogGamma_15_2Q :: Vector Double
tableLogGamma_15_2Q = [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
U.fromList
  [  Double
1
  , -Double
0.150169356054485044494e1
  ,  Double
0.846973248876495016101e0
  , -Double
0.220095151814995745555e0
  ,  Double
0.25582797155975869989e-1
  , -Double
0.100666795539143372762e-2
  , -Double
0.827193521891290553639e-6
  ]
{-# NOINLINE tableLogGamma_15_2Q #-}



-- Polynomial expansion used in interval (2,3)
--
-- > logΓ(z) = (z - 2)(z + 1)(Y + R(z-2))
lgamma2_3 :: Double -> Double
lgamma2_3 :: Double -> Double
lgamma2_3 Double
z
  = Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
* ( Double -> Vector Double -> Double
forall (v :: * -> *) a. (Vector v a, Num a) => a -> v a -> a
evaluatePolynomial Double
zm2 Vector Double
tableLogGamma_2_3P
                Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Vector Double -> Double
forall (v :: * -> *) a. (Vector v a, Num a) => a -> v a -> a
evaluatePolynomial Double
zm2 Vector Double
tableLogGamma_2_3Q
                )
  where
    r :: Double
r   = Double
zm2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1)
    zm2 :: Double
zm2 = Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2
    y :: Double
y   = Double
0.158963680267333984375e0


tableLogGamma_2_3P,tableLogGamma_2_3Q :: U.Vector Double
tableLogGamma_2_3P :: Vector Double
tableLogGamma_2_3P = [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
U.fromList
  [ -Double
0.180355685678449379109e-1
  ,  Double
0.25126649619989678683e-1
  ,  Double
0.494103151567532234274e-1
  ,  Double
0.172491608709613993966e-1
  , -Double
0.259453563205438108893e-3
  , -Double
0.541009869215204396339e-3
  , -Double
0.324588649825948492091e-4
  ]
{-# NOINLINE tableLogGamma_2_3P #-}
tableLogGamma_2_3Q :: Vector Double
tableLogGamma_2_3Q = [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
U.fromList
  [  Double
1
  ,  Double
0.196202987197795200688e1
  ,  Double
0.148019669424231326694e1
  ,  Double
0.541391432071720958364e0
  ,  Double
0.988504251128010129477e-1
  ,  Double
0.82130967464889339326e-2
  ,  Double
0.224936291922115757597e-3
  , -Double
0.223352763208617092964e-6
  ]
{-# NOINLINE tableLogGamma_2_3Q #-}



-- For small z we can just use Gamma function recurrence and reduce
-- problem to interval [2,3] and use polynomial approximation
-- there. Surprisingly it gives very good precision
lgammaSmall :: Double -> Double
lgammaSmall :: Double -> Double
lgammaSmall = Double -> Double -> Double
go Double
0
  where
    go :: Double -> Double -> Double
go Double
acc Double
z | Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
3     = Double
acc Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
lgamma2_3 Double
z
             | Bool
otherwise = Double -> Double -> Double
go (Double
acc Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
log Double
zm1) Double
zm1
             where
               zm1 :: Double
zm1 = Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1


-- Lanczos approximation for gamma function.
--
-- > Γ(z) = sqrt(2π)(z + g - 0.5)^(z - 0.5)·exp{-(z + g - 0.5)}·A_g(z)
--
-- Coefficients are taken from boost. Constants are absorbed into
-- polynomial's coefficients.
lanczosApprox :: Double -> Double
lanczosApprox :: Double -> Double
lanczosApprox Double
z
  = (Double -> Double
forall a. Floating a => a -> a
log (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
g Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5)
  Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
log (Vector (Double, Double) -> Double -> Double
evalRatio Vector (Double, Double)
tableLanczos Double
z)
  where
    g :: Double
g = Double
6.024680040776729583740234375

tableLanczos :: U.Vector (Double,Double)
{-# NOINLINE tableLanczos #-}
tableLanczos :: Vector (Double, Double)
tableLanczos = [(Double, Double)] -> Vector (Double, Double)
forall a. Unbox a => [a] -> Vector a
U.fromList
  [ (Double
56906521.91347156388090791033559122686859    , Double
0)
  , (Double
103794043.1163445451906271053616070238554    , Double
39916800)
  , (Double
86363131.28813859145546927288977868422342    , Double
120543840)
  , (Double
43338889.32467613834773723740590533316085    , Double
150917976)
  , (Double
14605578.08768506808414169982791359218571    , Double
105258076)
  , (Double
3481712.15498064590882071018964774556468     , Double
45995730)
  , (Double
601859.6171681098786670226533699352302507    , Double
13339535)
  , (Double
75999.29304014542649875303443598909137092    , Double
2637558)
  , (Double
6955.999602515376140356310115515198987526    , Double
357423)
  , (Double
449.9445569063168119446858607650988409623    , Double
32670)
  , (Double
19.51992788247617482847860966235652136208    , Double
1925)
  , (Double
0.5098416655656676188125178644804694509993   , Double
66)
  , (Double
0.006061842346248906525783753964555936883222 , Double
1)
  ]

-- Evaluate rational function. Polynomials in both numerator and
-- denominator must have same order. Function seems to be too specific
-- so it's not exposed
--
-- Special care taken in order to avoid overflow for large values of x
evalRatio :: U.Vector (Double,Double) -> Double -> Double
evalRatio :: Vector (Double, Double) -> Double -> Double
evalRatio Vector (Double, Double)
coef Double
x
  | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1     = L -> Double
fini (L -> Double) -> L -> Double
forall a b. (a -> b) -> a -> b
$ (L -> (Double, Double) -> L) -> L -> Vector (Double, Double) -> L
forall b a. Unbox b => (a -> b -> a) -> a -> Vector b -> a
U.foldl' L -> (Double, Double) -> L
stepL (Double -> Double -> L
L Double
0 Double
0) Vector (Double, Double)
coef
  | Bool
otherwise = L -> Double
fini (L -> Double) -> L -> Double
forall a b. (a -> b) -> a -> b
$ ((Double, Double) -> L -> L) -> L -> Vector (Double, Double) -> L
forall a b. Unbox a => (a -> b -> b) -> b -> Vector a -> b
U.foldr' (Double, Double) -> L -> L
stepR (Double -> Double -> L
L Double
0 Double
0) Vector (Double, Double)
coef
  where
    fini :: L -> Double
fini (L Double
num Double
den) = Double
num Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
den
    stepR :: (Double, Double) -> L -> L
stepR (Double
a,Double
b) (L Double
num Double
den) = Double -> Double -> L
L (Double
num Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x  Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
a) (Double
den Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x  Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b)
    stepL :: L -> (Double, Double) -> L
stepL (L Double
num Double
den) (Double
a,Double
b) = Double -> Double -> L
L (Double
num Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
rx Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
a) (Double
den Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
rx Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b)
    rx :: Double
rx = Double -> Double
forall a. Fractional a => a -> a
recip Double
x



-- |
-- Compute the log gamma correction factor for Stirling
-- approximation for @x@ &#8805; 10.  This correction factor is
-- suitable for an alternate (but less numerically accurate)
-- definition of 'logGamma':
--
-- \[
-- \log\Gamma(x) = \frac{1}{2}\log(2\pi) + (x-\frac{1}{2})\log x - x + \operatorname{logGammaCorrection}(x)
-- \]
logGammaCorrection :: Double -> Double
logGammaCorrection :: Double -> Double
logGammaCorrection Double
x
    | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
10    = Double
m_NaN
    | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
big   = Double -> Vector Double -> Double
forall (v :: * -> *).
Vector v Double =>
Double -> v Double -> Double
chebyshevBroucke (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Vector Double
coeffs Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x
    | Bool
otherwise = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
12)
  where
    big :: Double
big    = Double
94906265.62425156
    t :: Double
t      = Double
10 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x
    coeffs :: Vector Double
coeffs = [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
U.fromList [
               Double
0.1666389480451863247205729650822e+0,
              -Double
0.1384948176067563840732986059135e-4,
               Double
0.9810825646924729426157171547487e-8,
              -Double
0.1809129475572494194263306266719e-10,
               Double
0.6221098041892605227126015543416e-13,
              -Double
0.3399615005417721944303330599666e-15,
               Double
0.2683181998482698748957538846666e-17
             ]



-- | Compute the normalized lower incomplete gamma function
-- γ(/z/,/x/). Normalization means that γ(/z/,∞)=1
--
-- \[
-- \gamma(z,x) = \frac{1}{\Gamma(z)}\int_0^{x}t^{z-1}e^{-t}\,dt
-- \]
--
-- Uses Algorithm AS 239 by Shea.
incompleteGamma :: Double       -- ^ /z/ ∈ (0,∞)
                -> Double       -- ^ /x/ ∈ (0,∞)
                -> Double
-- Notation used:
--  + P(a,x) - regularized lower incomplete gamma
--  + Q(a,x) - regularized upper incomplete gamma
incompleteGamma :: Double -> Double -> Double
incompleteGamma Double
a Double
x
  | Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
|| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 = String -> Double
forall a. HasCallStack => String -> a
error
     (String -> Double) -> String -> Double
forall a b. (a -> b) -> a -> b
$ String
"incompleteGamma: Domain error z=" String -> String -> String
forall a. [a] -> [a] -> [a]
++ Double -> String
forall a. Show a => a -> String
show Double
a String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
" x=" String -> String -> String
forall a. [a] -> [a] -> [a]
++ Double -> String
forall a. Show a => a -> String
show Double
x
  | Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0          = Double
0
  | Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
m_pos_inf  = Double
1
  -- For very small x we use following expansion for P:
  --
  -- See http://functions.wolfram.com/GammaBetaErf/GammaRegularized/06/01/05/01/01/
  | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double -> Double
forall a. Floating a => a -> a
sqrt Double
m_epsilon Bool -> Bool -> Bool
&& Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1
    = Double
xDouble -> Double -> Double
forall a. Floating a => a -> a -> a
**Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
exp (Double -> Double
logGamma Double
a) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1))
  | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.5 = case () of
    ()
_| (-Double
0.4)Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double -> Double
forall a. Floating a => a -> a
log Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
a  -> Double
taylorSeriesP
     | Bool
otherwise         -> Double
taylorSeriesComplQ
  | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1.1 = case () of
    ()
_| Double
0.75Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
a        -> Double
taylorSeriesP
     | Bool
otherwise         -> Double
taylorSeriesComplQ
  | Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
20 Bool -> Bool -> Bool
&& Bool
useTemme    = Double
uniformExpansion
  | Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x)) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
a = Double
taylorSeriesP
  | Bool
otherwise             = Double
contFraction
  where
    mu :: Double
mu = (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a
    useTemme :: Bool
useTemme = (Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
200 Bool -> Bool -> Bool
&& Double
20Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
muDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
mu)
            Bool -> Bool -> Bool
|| (Double -> Double
forall a. Num a => a -> a
abs Double
mu Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.4)
    -- Gautschi's algorithm.
    --
    -- Evaluate series for P(a,x). See [Temme1994] Eq. 5.5 and [NOTE:
    -- incompleteGamma.taylorP]
    factorP :: Double
factorP
      | Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
10     = Double
x Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
a
                   Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double -> Double
forall a. Floating a => a -> a
exp Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp (Double -> Double
logGamma (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1)))
      | Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1182.5 = (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
a
                   Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
exp Double
x
                   Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a)
                   Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
exp (Double -> Double
logGammaCorrection Double
a)
      | Bool
otherwise  = (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp (-Double
xDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
a)) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
a
                   Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a)
                   Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
exp (Double -> Double
logGammaCorrection Double
a)
    taylorSeriesP :: Double
taylorSeriesP
      = Double -> Sequence Double -> Double
sumPowerSeries Double
x ((Double -> Double -> Double)
-> Double -> Sequence Double -> Sequence Double
forall b a. (b -> a -> b) -> b -> Sequence a -> Sequence b
scanSequence Double -> Double -> Double
forall a. Fractional a => a -> a -> a
(/) Double
1 (Sequence Double -> Sequence Double)
-> Sequence Double -> Sequence Double
forall a b. (a -> b) -> a -> b
$ Double -> Sequence Double
forall a. Num a => a -> Sequence a
enumSequenceFrom (Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1))
      Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
factorP
    -- Series for 1-Q(a,x). See [Temme1994] Eq. 5.5
    taylorSeriesComplQ :: Double
taylorSeriesComplQ
      = Double -> Sequence Double -> Double
sumPowerSeries (-Double
x) ((Double -> Double -> Double)
-> Double -> Sequence Double -> Sequence Double
forall b a. (b -> a -> b) -> b -> Sequence a -> Sequence b
scanSequence Double -> Double -> Double
forall a. Fractional a => a -> a -> a
(/) Double
1 (Double -> Sequence Double
forall a. Num a => a -> Sequence a
enumSequenceFrom Double
1) Sequence Double -> Sequence Double -> Sequence Double
forall a. Fractional a => a -> a -> a
/ Double -> Sequence Double
forall a. Num a => a -> Sequence a
enumSequenceFrom Double
a)
      Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
xDouble -> Double -> Double
forall a. Floating a => a -> a -> a
**Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
exp(Double -> Double
logGamma Double
a)
    -- Legendre continued fractions
    contFraction :: Double
contFraction = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- ( Double -> Double
forall a. Floating a => a -> a
exp ( Double -> Double
forall a. Floating a => a -> a
log Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
logGamma Double
a )
                       Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Sequence (Double, Double) -> Double
evalContFractionB Sequence (Double, Double)
frac
                       )
      where
        frac :: Sequence (Double, Double)
frac = (\Double
k -> (Double
kDouble -> Double -> Double
forall a. Num a => a -> a -> a
*(Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
k), Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1)) (Double -> (Double, Double))
-> Sequence Double -> Sequence (Double, Double)
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Double -> Sequence Double
forall a. Num a => a -> Sequence a
enumSequenceFrom Double
0
    -- Evaluation based on uniform expansions. See [Temme1994] 5.2
    uniformExpansion :: Double
uniformExpansion =
      let -- Coefficients f_m in paper
          fm :: U.Vector Double
          fm :: Vector Double
fm = [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
U.fromList [ Double
1.00000000000000000000e+00
                          ,-Double
3.33333333333333370341e-01
                          , Double
8.33333333333333287074e-02
                          ,-Double
1.48148148148148153802e-02
                          , Double
1.15740740740740734316e-03
                          , Double
3.52733686067019369930e-04
                          ,-Double
1.78755144032921825352e-04
                          , Double
3.91926317852243766954e-05
                          ,-Double
2.18544851067999240532e-06
                          ,-Double
1.85406221071515996597e-06
                          , Double
8.29671134095308545622e-07
                          ,-Double
1.76659527368260808474e-07
                          , Double
6.70785354340149841119e-09
                          , Double
1.02618097842403069078e-08
                          ,-Double
4.38203601845335376897e-09
                          , Double
9.14769958223679020897e-10
                          ,-Double
2.55141939949462514346e-11
                          ,-Double
5.83077213255042560744e-11
                          , Double
2.43619480206674150369e-11
                          ,-Double
5.02766928011417632057e-12
                          , Double
1.10043920319561347525e-13
                          , Double
3.37176326240098513631e-13
                          ]
          y :: Double
y   = - Double -> Double
log1pmx Double
mu
          eta :: Double
eta = Double -> Double
forall a. Floating a => a -> a
sqrt (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Num a => a -> a
signum Double
mu
          -- Evaluate S_α (Eq. 5.9)
          loop :: Double -> Double -> Double -> Int -> Double
loop !Double
_  !Double
_  Double
u Int
0 = Double
u
          loop Double
bm1 Double
bm0 Double
u Int
i = let t :: Double
t  = (Vector Double
fm Vector Double -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
! Int
i) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
bm1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a
                                 u' :: Double
u' = Double
eta Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t
                             in  Double -> Double -> Double -> Int -> Double
loop Double
bm0 Double
t Double
u' (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)
          s_a :: Double
s_a = let n :: Int
n = Vector Double -> Int
forall a. Unbox a => Vector a -> Int
U.length Vector Double
fm
                in  Double -> Double -> Double -> Int -> Double
loop (Vector Double
fm Vector Double -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
! (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1)) (Vector Double
fm Vector Double -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
! (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
2)) Double
0 (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
3)
                  Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
exp (Double -> Double
logGammaCorrection Double
a)
      in Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
erfc(-Double
etaDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double -> Double
forall a. Floating a => a -> a
sqrt(Double
aDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2)) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
exp(-(Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
y)) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
sqrt (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
piDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s_a



-- Adapted from Numerical Recipes §6.2.1

-- | Inverse incomplete gamma function. It's approximately inverse of
--   'incompleteGamma' for the same /z/. So following equality
--   approximately holds:
--
-- > invIncompleteGamma z . incompleteGamma z ≈ id
invIncompleteGamma :: Double    -- ^ /z/ ∈ (0,∞)
                   -> Double    -- ^ /p/ ∈ [0,1]
                   -> Double
invIncompleteGamma :: Double -> Double -> Double
invIncompleteGamma Double
a Double
p
  | Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0         =
      String -> Double
forall a. String -> a
modErr (String -> Double) -> String -> Double
forall a b. (a -> b) -> a -> b
$ String -> Double -> Double -> String
forall r. PrintfType r => String -> r
printf String
"invIncompleteGamma: a must be positive. a=%g p=%g" Double
a Double
p
  | Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 Bool -> Bool -> Bool
|| Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1 =
      String -> Double
forall a. String -> a
modErr (String -> Double) -> String -> Double
forall a b. (a -> b) -> a -> b
$ String -> Double -> Double -> String
forall r. PrintfType r => String -> r
printf String
"invIncompleteGamma: p must be in [0,1] range. a=%g p=%g" Double
a Double
p
  | Double
p Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0         = Double
0
  | Double
p Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
1         = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0
  | Bool
otherwise      = Int -> Double -> Double
loop Int
0 Double
guess
  where
    -- Solve equation γ(a,x) = p using Halley method
    loop :: Int -> Double -> Double
    loop :: Int -> Double -> Double
loop Int
i Double
x
      | Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= Int
12           = Double
x'
      -- For small s derivative becomes approximately 1/x*exp(-x) and
      -- skyrockets for small x. If it happens correct answer is 0.
      | Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
f'     = Double
0
      | Double -> Double
forall a. Num a => a -> a
abs Double
dx Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
eps Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x' = Double
x'
      | Bool
otherwise         = Int -> Double -> Double
loop (Int
i Int -> Int -> Int
forall a. Num a => a -> a -> a
+ Int
1) Double
x'
      where
        -- Value of γ(a,x) - p
        f :: Double
f    = Double -> Double -> Double
incompleteGamma Double
a Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p
        -- dγ(a,x)/dx
        f' :: Double
f'   | Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1     = Double
afac Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp( -(Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a1) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
a1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Floating a => a -> a
log Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
lna1))
             | Bool
otherwise = Double -> Double
forall a. Floating a => a -> a
exp( -Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
a1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
gln)
        u :: Double
u    = Double
f Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
f'
        -- Halley correction to Newton-Rapson step
        corr :: Double
corr = Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
a1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)
        dx :: Double
dx   = Double
u Double -> Double -> Double
forall a. Fractional 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 -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
1.0 Double
corr)
        -- New approximation to x
        x' :: Double
x'   | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
dx    = Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x -- Do not go below 0
             | Bool
otherwise = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
dx
    -- Calculate initial guess for root
    guess :: Double
guess
      --
      | Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1   =
         let t :: Double
t  = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ -Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log(if Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.5 then Double
p else Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p)
             x1 :: Double
x1 = (Double
2.30753 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
0.27061) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
0.99229 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
0.04481)) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t
             x2 :: Double
x2 = if Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.5 then -Double
x1 else Double
x1
         in Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
1e-3 (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
9Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt Double
a)) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
3)
      -- For a <= 1 use following approximations:
      --   γ(a,1) ≈ 0.253a + 0.12a²
      --
      --   γ(a,x) ≈ γ(a,1)·x^a                               x <  1
      --   γ(a,x) ≈ γ(a,1) + (1 - γ(a,1))(1 - exp(1 - x))    x >= 1
      | Bool
otherwise =
         let t :: Double
t = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
0.253 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
0.12)
         in if Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
t
            then (Double
p Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
t) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a)
            else Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
log( Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
pDouble -> Double -> Double
forall a. Num a => a -> a -> a
-Double
t) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
t))
    -- Constants
    a1 :: Double
a1   = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1
    lna1 :: Double
lna1 = Double -> Double
forall a. Floating a => a -> a
log Double
a1
    afac :: Double
afac = Double -> Double
forall a. Floating a => a -> a
exp( Double
a1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
lna1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
gln )
    gln :: Double
gln  = Double -> Double
logGamma Double
a
    eps :: Double
eps  = Double
1e-8



----------------------------------------------------------------
-- Beta function
----------------------------------------------------------------

-- | Compute the natural logarithm of the beta function.
--
-- \[
-- B(a,b) = \int_0^1 t^{a-1}(1-t)^{b-1}\,dt = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}
-- \]
logBeta
  :: Double                     -- ^ /a/ > 0
  -> Double                     -- ^ /b/ > 0
  -> Double
logBeta :: Double -> Double -> Double
logBeta Double
a Double
b
  | Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0     = Double
m_NaN
  | Double
p Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0    = Double
m_pos_inf
  | Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
10   = Double
allStirling
  | Double
q Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
10   = Double
twoStirling
  -- This order of summands marginally improves precision
  | Bool
otherwise = Double -> Double
logGamma Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double -> Double
logGamma Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
logGamma Double
pq)
  where
    p :: Double
p   = Double -> Double -> Double
forall a. Ord a => a -> a -> a
min Double
a Double
b
    q :: Double
q   = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
a Double
b
    ppq :: Double
ppq = Double
p Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
pq
    pq :: Double
pq  = Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
q
    -- When both parameters are large than 10 we can use Stirling
    -- approximation with correction. It's more precise than sum of
    -- logarithms of gamma functions
    allStirling :: Double
allStirling
      = Double -> Double
forall a. Floating a => a -> a
log Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-Double
0.5)
      Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
m_ln_sqrt_2_pi
      Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
logGammaCorrection Double
p
      Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double -> Double
logGammaCorrection Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
logGammaCorrection Double
pq)
      Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
ppq
      Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log1p(-Double
ppq)
    -- Otherwise only two of three gamma functions use Stirling
    -- approximation
    twoStirling :: Double
twoStirling
      = Double -> Double
logGamma Double
p
      Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double -> Double
logGammaCorrection Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
logGammaCorrection Double
pq)
      Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
p
      Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
pq
      Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log1p(-Double
ppq)


-- | Regularized incomplete beta function.
--
-- \[
-- I(x;a,b) = \frac{1}{B(a,b)} \int_0^x t^{a-1}(1-t)^{b-1}\,dt
-- \]
--
-- Uses algorithm AS63 by Majumder and Bhattachrjee and quadrature
-- approximation for large /p/ and /q/.
incompleteBeta :: Double -- ^ /a/ > 0
               -> Double -- ^ /b/ > 0
               -> Double -- ^ /x/, must lie in [0,1] range
               -> Double
incompleteBeta :: Double -> Double -> Double -> Double
incompleteBeta Double
p Double
q = Double -> Double -> Double -> Double -> Double
incompleteBeta_ (Double -> Double -> Double
logBeta Double
p Double
q) Double
p Double
q

-- | Regularized incomplete beta function. Same as 'incompleteBeta'
-- but also takes logarithm of beta function as parameter.
incompleteBeta_ :: Double -- ^ logarithm of beta function for given /p/ and /q/
                -> Double -- ^ /a/ > 0
                -> Double -- ^ /b/ > 0
                -> Double -- ^ /x/, must lie in [0,1] range
                -> Double
incompleteBeta_ :: Double -> Double -> Double -> Double -> Double
incompleteBeta_ Double
beta Double
p Double
q Double
x
  | Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
|| Double
q Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0            =
      String -> Double
forall a. String -> a
modErr (String -> Double) -> String -> Double
forall a b. (a -> b) -> a -> b
$ String -> Double -> Double -> Double -> String
forall r. PrintfType r => String -> r
printf String
"incompleteBeta_: p <= 0 || q <= 0. p=%g q=%g x=%g" Double
p Double
q Double
x
  | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<  Double
0 Bool -> Bool -> Bool
|| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>  Double
1 Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x =
      String -> Double
forall a. String -> a
modErr (String -> Double) -> String -> Double
forall a b. (a -> b) -> a -> b
$ String -> Double -> Double -> Double -> String
forall r. PrintfType r => String -> r
printf String
"incompleteBeta_: x out of [0,1] range. p=%g q=%g x=%g" Double
p Double
q Double
x
  | Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 Bool -> Bool -> Bool
|| Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
1            = Double
x
  | Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= (Double
pDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
q) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x   = Double -> Double -> Double -> Double -> Double
incompleteBetaWorker Double
beta Double
p Double
q Double
x
  | Bool
otherwise        = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double -> Double -> Double -> Double
incompleteBetaWorker Double
beta Double
q Double
p (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x)


-- Approximation of incomplete beta by quadrature.
--
-- Note that x =< p/(p+q)
incompleteBetaApprox :: Double -> Double -> Double -> Double -> Double
incompleteBetaApprox :: Double -> Double -> Double -> Double -> Double
incompleteBetaApprox Double
beta Double
p Double
q Double
x
  | Double
ans Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0   = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
ans
  | Bool
otherwise = -Double
ans
  where
    -- Constants
    p1 :: Double
p1    = Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1
    q1 :: Double
q1    = Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1
    mu :: Double
mu    = Double
p Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
q)
    lnmu :: Double
lnmu  = Double -> Double
forall a. Floating a => a -> a
log     Double
mu
    lnmuc :: Double
lnmuc = Double -> Double
forall a. Floating a => a -> a
log1p (-Double
mu)
    -- Upper limit for integration
    xu :: Double
xu = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max Double
0 (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double
forall a. Ord a => a -> a -> a
min (Double
mu Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
10Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t) (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
5Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t)
       where
         t :: Double
t = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
pDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
q Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ( (Double
pDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
q) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
pDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
q) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1) )
    -- Calculate incomplete beta by quadrature
    go :: Double -> Double -> Double
go Double
y Double
w = let t :: Double
t = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
xu Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
y
             in  Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp( Double
p1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Floating a => a -> a
log Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
lnmu) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
q1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double -> Double
forall a. Floating a => a -> a
log(Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
t) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
lnmuc) )
    s :: Double
s   = Vector Double -> Double
forall a. (Unbox a, Num a) => Vector a -> a
U.sum (Vector Double -> Double) -> Vector Double -> Double
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double)
-> Vector Double -> Vector Double -> Vector Double
forall a b c.
(Unbox a, Unbox b, Unbox c) =>
(a -> b -> c) -> Vector a -> Vector b -> Vector c
U.zipWith Double -> Double -> Double
go Vector Double
coefY Vector Double
coefW
    ans :: Double
ans = Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
xu Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp( Double
p1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
lnmu Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
q1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
lnmuc Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
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 :: Double -> Double -> Double -> Double -> Double
incompleteBetaWorker Double
beta Double
p Double
q Double
x
  -- For very large p and q this method becomes very slow so another
  -- method is used.
  | Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
3000 Bool -> Bool -> Bool
&& Double
q Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
3000 = Double -> Double -> Double -> Double -> Double
incompleteBetaApprox Double
beta Double
p Double
q Double
x
  | Bool
otherwise            = Double -> Int -> Double -> Double -> Double -> Double
loop (Double
pDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
q) (Double -> Int
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
truncate (Double -> Int) -> Double -> Int
forall a b. (a -> b) -> a -> b
$ Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
cx Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
pDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
q)) Double
1 Double
1 Double
1
  where
    -- Constants
    eps :: Double
eps = Double
1e-15
    cx :: Double
cx  = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x
    -- Common multiplies for expansion. Accurate calculation is a bit
    -- tricky. Performing calculation in log-domain leads to slight
    -- loss of precision for small x, while using ** prone to
    -- underflows.
    --
    -- If either beta function of x**p·(1-x)**(q-1) underflows we
    -- switch to log domain. It could waste work but there's no easy
    -- switch criterion.
    factor :: Double
factor
      | Double
beta Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
m_min_log Bool -> Bool -> Bool
|| Double
prod Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
m_tiny = Double -> Double
forall a. Floating a => a -> a
exp( Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
cx Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
beta)
      | Bool
otherwise                         = Double
prod Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
exp Double
beta
      where
        prod :: Double
prod =  Double
xDouble -> Double -> Double
forall a. Floating a => a -> a -> a
**Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
cxDouble -> Double -> Double
forall a. Floating a => a -> a -> a
**(Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)
    -- Soper's expansion of incomplete beta function
    loop :: Double -> Int -> Double -> Double -> Double -> Double
loop !Double
psq (Int
ns :: Int) Double
ai Double
term Double
betain
      | Bool
done      = Double
betain' Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
factor Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
p
      | Bool
otherwise = Double -> Int -> Double -> Double -> Double -> Double
loop Double
psq' (Int
ns Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
1) (Double
ai Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1) Double
term' Double
betain'
      where
        -- New values
        term' :: Double
term'   = Double
term Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
fact Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
ai)
        betain' :: Double
betain' = Double
betain Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
term'
        fact :: Double
fact | Int
ns Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>  Int
0   = (Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
ai) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
xDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
cx
             | Int
ns Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
0   = (Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
ai) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x
             | Bool
otherwise = Double
psq Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x
        -- Iterations are complete
        done :: Bool
done = Double
db Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
eps Bool -> Bool -> Bool
&& Double
db Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
epsDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
betain' where db :: Double
db = Double -> Double
forall a. Num a => a -> a
abs Double
term'
        psq' :: Double
psq' = if Int
ns Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 then Double
psq Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1 else Double
psq



-- | Compute inverse of regularized incomplete beta function. Uses
-- initial approximation from AS109, AS64 and Halley method to solve
-- equation.
invIncompleteBeta :: Double     -- ^ /a/ > 0
                  -> Double     -- ^ /b/ > 0
                  -> Double     -- ^ /x/ ∈ [0,1]
                  -> Double
invIncompleteBeta :: Double -> Double -> Double -> Double
invIncompleteBeta Double
p Double
q Double
a
  | Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
|| Double
q Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 =
      String -> Double
forall a. String -> a
modErr (String -> Double) -> String -> Double
forall a b. (a -> b) -> a -> b
$ String -> Double -> Double -> Double -> String
forall r. PrintfType r => String -> r
printf String
"invIncompleteBeta p <= 0 || q <= 0.  p=%g q=%g a=%g" Double
p Double
q Double
a
  | Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<  Double
0 Bool -> Bool -> Bool
|| Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>  Double
1 =
      String -> Double
forall a. String -> a
modErr (String -> Double) -> String -> Double
forall a b. (a -> b) -> a -> b
$ String -> Double -> Double -> Double -> String
forall r. PrintfType r => String -> r
printf String
"invIncompleteBeta x must be in [0,1].  p=%g q=%g a=%g" Double
p Double
q Double
a
  | Double
a Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 Bool -> Bool -> Bool
|| Double
a Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
1 = Double
a
  | Bool
otherwise        = Double -> Double -> Double -> Double -> Double
invIncompleteBetaWorker (Double -> Double -> Double
logBeta Double
p Double
q) Double
p Double
q  Double
a


invIncompleteBetaWorker :: Double -> Double -> Double -> Double -> Double
invIncompleteBetaWorker :: Double -> Double -> Double -> Double -> Double
invIncompleteBetaWorker Double
beta Double
a Double
b Double
p = Int -> Double -> Double
forall {t}. (Ord t, Num t) => t -> Double -> Double
loop (Int
0::Int) (Double -> Double -> Double -> Double -> Double
invIncBetaGuess Double
beta Double
a Double
b Double
p)
  where
    a1 :: Double
a1 = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1
    b1 :: Double
b1 = Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1
    -- Solve equation using Halley method
    loop :: t -> Double -> Double
loop !t
i !Double
x
      -- We cannot continue at this point so we simply return `x'
      | Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 Bool -> Bool -> Bool
|| Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
1             = Double
x
      -- When derivative becomes infinite we cannot continue
      -- iterations. It can only happen in vicinity of 0 or 1. It's
      -- hardly possible to get good answer in such circumstances but
      -- `x' is already reasonable.
      | Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
f'                = Double
x
      -- Iterations limit reached. Most of the time solution will
      -- converge to answer because of discreteness of Double. But
      -- solution have good precision already.
      | t
i t -> t -> Bool
forall a. Ord a => a -> a -> Bool
>= t
10                      = Double
x
      -- Solution converges
      | Double -> Double
forall a. Num a => a -> a
abs Double
dx Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
16 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
m_epsilon Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x = Double
x'
      | Bool
otherwise                    = t -> Double -> Double
loop (t
it -> t -> t
forall a. Num a => a -> a -> a
+t
1) Double
x'
      where
        -- Calculate Halley step.
        f :: Double
f   = Double -> Double -> Double -> Double -> Double
incompleteBeta_ Double
beta Double
a Double
b Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p
        f' :: Double
f'  = Double -> Double
forall a. Floating a => a -> a
exp (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
a1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log1p (-Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
beta
        u :: Double
u   = Double
f Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
f'
        -- We bound Halley correction to Newton-Raphson to (-1,1) range
        corr :: Double
corr | Double
d Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1     = Double
1
             | Double
d Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< -Double
1    = -Double
1
             | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
d   = Double
0
             | Bool
otherwise = Double
d
          where
            d :: Double
d = Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
a1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x))
        dx :: Double
dx  = Double
u Double -> Double -> Double
forall a. Fractional 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
corr)
        -- Next approximation. If Halley step leads us out of [0,1]
        -- range we revert to bisection.
        x' :: Double
x'  | Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0     = Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
            | Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1     = (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
            | Bool
otherwise = Double
z
            where z :: Double
z = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
dx


-- Calculate initial guess for inverse incomplete beta function.
invIncBetaGuess :: Double -> Double -> Double -> Double -> Double
-- Calculate initial guess. for solving equation for inverse incomplete beta.
-- It's really hodgepodge of different approximations accumulated over years.
--
-- Equations are referred to by name of paper and number e.g. [AS64 2]
-- In AS64 papers equations are not numbered so they are referred to by
-- number of appearance starting from definition of incomplete beta.
invIncBetaGuess :: Double -> Double -> Double -> Double -> Double
invIncBetaGuess Double
beta Double
a Double
b Double
p
  -- If both a and b are less than 1 incomplete beta have inflection
  -- point.
  --
  -- > x = (1 - a) / (2 - a - b)
  --
  -- We approximate incomplete beta by neglecting one of factors under
  -- integral and then rescaling result of integration into [0,1]
  -- range.
  | Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1 Bool -> Bool -> Bool
&& Double
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
1 =
    let x_infl :: Double
x_infl = (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
b)
        p_infl :: Double
p_infl = Double -> Double -> Double -> Double
incompleteBeta Double
a Double
b Double
x_infl
        x :: Double
x | Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
p_infl = let xg :: Double
xg = (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
p     Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp Double
beta) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
a) in Double
xg Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
xg)
          | Bool
otherwise  = let xg :: Double
xg = (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
p) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp Double
beta) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
b) in Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
xgDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
xg)
    in Double
x
  -- If both a and b larger or equal that 1 but not too big we use
  -- same approximation as above but calculate it a bit differently
  | Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
6 Bool -> Bool -> Bool
&& Double
aDouble -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>Double
1 Bool -> Bool -> Bool
&& Double
bDouble -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>Double
1 =
    let x_infl :: Double
x_infl = (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2)
        p_infl :: Double
p_infl = Double -> Double -> Double -> Double
incompleteBeta Double
a Double
b Double
x_infl
        x :: Double
x | Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
p_infl = Double -> Double
forall a. Floating a => a -> a
exp ((Double -> Double
forall a. Floating a => a -> a
log(Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
a) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
beta) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a)
          | Bool
otherwise  = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
exp((Double -> Double
forall a. Floating a => a -> a
log((Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
p) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
beta) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
b)
    in Double
x
  -- For small a and not too big b we use approximation from boost.
  | Double
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
5 Bool -> Bool -> Bool
&& Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1 =
    let x :: Double
x | Double
pDouble -> Double -> Double
forall a. Floating a => a -> a -> a
**(Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
a) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0.5 = (Double
p Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp Double
beta)) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
a)
          | Bool
otherwise      = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp Double
beta))Double -> Double -> Double
forall a. Floating a => a -> a -> a
**(Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
b)
    in Double
x
  -- When a>>b and both are large approximation from [Temme1992],
  -- section 4 "the incomplete gamma function case" used. In this
  -- region it greatly improves over other approximation (AS109, AS64,
  -- "Numerical Recipes")
  --
  -- FIXME: It could be used when b>>a too but it require inverse of
  --        upper incomplete gamma to be precise enough. In current
  --        implementation it loses precision in horrible way (40
  --        order of magnitude off for sufficiently small p)
  | Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
5 Bool -> Bool -> Bool
&&  Double
aDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
4 =
    let -- Calculate initial approximation to eta using eq 4.1
        eta0 :: Double
eta0 = Double -> Double -> Double
invIncompleteGamma Double
b (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
p) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a
        mu :: Double
mu   = Double
b Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a            -- Eq. 4.3
        -- A lot of helpers for calculation of
        w :: Double
w    = Double -> Double
forall a. Floating a => a -> a
sqrt(Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
mu)     -- Eq. 4.9
        w_2 :: Double
w_2  = Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w
        w_3 :: Double
w_3  = Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w
        w_4 :: Double
w_4  = Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2
        w_5 :: Double
w_5  = Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2
        w_6 :: Double
w_6  = Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3
        w_7 :: Double
w_7  = Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3
        w_8 :: Double
w_8  = Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4
        w_9 :: Double
w_9  = Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4
        w_10 :: Double
w_10 = Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5
        d :: Double
d    = Double
eta0 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
mu
        d_2 :: Double
d_2  = Double
d Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d
        d_3 :: Double
d_3  = Double
d_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d
        d_4 :: Double
d_4  = Double
d_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d_2
        w1 :: Double
w1   = Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1
        w1_2 :: Double
w1_2 = Double
w1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w1
        w1_3 :: Double
w1_3 = Double
w1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w1_2
        w1_4 :: Double
w1_4 = Double
w1_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w1_2
        -- Evaluation of eq 4.10
        e1 :: Double
e1 = (Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w)
           Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
9 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
21 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
5) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d
             Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
36 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w1)
           Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
13 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
69 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
167 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
46) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d_2
             Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1620 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w1_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3)
           Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
7 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
21 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
70 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
26 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
93 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
31) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d_3
             Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
6480 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w1_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4)
           Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
75 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_6 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
202 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
188 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
888 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1345 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
118 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
138) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d_4
             Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
272160 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w1_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5)
        e2 :: Double
e2 = (Double
28 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
131 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
402 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
581 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
208) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)
             Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1620 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3)
           Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
35 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_6 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
154 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
623 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1636 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3983 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3514 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
925) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d
             Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
12960 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w1_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4)
           Double -> Double -> Double
forall a. Num a => a -> a -> a
- ( Double
2132 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_7 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
7915 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_6 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
16821 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
35066 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
87490 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3
             Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
141183 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
95993 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
21640
             ) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d_2
             Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
816480 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w1_3)
           Double -> Double -> Double
forall a. Num a => a -> a -> a
- ( Double
11053 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_8 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
53308 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_7 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
117010 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_6 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
163924 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
116188 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4
             Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
258428 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
677042 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
481940 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
105497
             ) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d_3
             Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
14696640 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w1_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_6)
        e3 :: Double
e3 = -( (Double
3592 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_7 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
8375 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_6 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1323 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
29198 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
89578 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3
                Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
154413 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
116063 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
29632
                ) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)
              )
              Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
816480 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w1_2)
           Double -> Double -> Double
forall a. Num a => a -> a -> a
- ( Double
442043 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_9 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2054169 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_8 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
3803094 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_7 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
3470754 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_6 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2141568 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5
             Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2393568 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
19904934 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
34714674 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
23128299 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
5253353
             ) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d
             Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
146966400 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_6 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w1_3)
           Double -> Double -> Double
forall a. Num a => a -> a -> a
- ( Double
116932 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_10 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
819281 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_9 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2378172 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_8 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
4341330 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_7 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
6806004 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_6
             Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
10622748 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
18739500 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
30651894 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
30869976 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2
             Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
15431867 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2919016
             ) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d_2
             Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
146966400 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w1_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_7)
        eta :: Double
eta = Double -> [Double] -> Double
forall a. Num a => a -> [a] -> a
evaluatePolynomialL (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
a) [Double
eta0, Double
e1, Double
e2, Double
e3]
        -- Now we solve eq 4.2 to recover x using Newton iterations
        u :: Double
u       = Double
eta Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
mu Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
eta Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
mu) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log(Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
mu) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
mu
        cross :: Double
cross   = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
mu);
        lower :: Double
lower   = if Double
eta Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
mu then Double
cross else Double
0
        upper :: Double
upper   = if Double
eta Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
mu then Double
1     else Double
cross
        x_guess :: Double
x_guess = (Double
lower Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
upper) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
2
        func :: Double -> (Double, Double)
func Double
x  = ( Double
u Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
log Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
muDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double -> Double
forall a. Floating a => a -> a
log(Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x)
                  , Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
muDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/(Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
x)
                  )
        Root Double
x0 = NewtonParam
-> (Double, Double, Double)
-> (Double -> (Double, Double))
-> Root Double
newtonRaphson NewtonParam
forall a. Default a => a
def{newtonTol=RelTol 1e-8} (Double
lower, Double
x_guess, Double
upper) Double -> (Double, Double)
func
    in Double
x0
  -- For large a and b approximation from AS109 (Carter
  -- approximation). It's reasonably good in this region
  | Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1 Bool -> Bool -> Bool
&& Double
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1 =
      let r :: Double
r = (Double
yDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
3) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
6
          s :: Double
s = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)
          t :: Double
t = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1)
          h :: Double
h = Double
2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t)
          w :: Double
w = Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt(Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
r) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
h Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
s) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
5Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
6 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
3 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
h))
      in Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
exp(Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
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.
  | Double
chi2 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
ratio Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
1 = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
ratio Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
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.
  | Bool
otherwise = case () of
      ()
_| Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
t Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
w  -> (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
a)
       | Bool
otherwise  -> Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
b)
       where
         lna :: Double
lna = Double -> Double
forall a. Floating a => a -> a
log (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
a Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
b)
         lnb :: Double
lnb = Double -> Double
forall a. Floating a => a -> a
log (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
b Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
b)
         t :: Double
t   = Double -> Double
forall a. Floating a => a -> a
exp( Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
lna ) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a
         u :: Double
u   = Double -> Double
forall a. Floating a => a -> a
exp( Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
lnb ) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
b
         w :: Double
w   = Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
u
  where
    -- Formula [AS64 2]
    ratio :: Double
ratio = (Double
4Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
2) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
chi2
    -- Quantile of chi-squared distribution. Formula [AS64 3].
    chi2 :: Double
chi2 = Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
sqrt Double
t) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** Double
3
      where
        t :: Double
t   = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
9 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b)
    -- `y' is Hasting's approximation of p'th quantile of standard
    -- normal distribution.
    y :: Double
y   = Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
- ( Double
2.30753 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.27061 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r )
              Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ( Double
1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ ( Double
0.99229 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
0.04481 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r ) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r )
      where
        r :: Double
r = Double -> Double
forall a. Floating a => a -> a
sqrt (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ - Double
2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
p



----------------------------------------------------------------
-- Sinc function
----------------------------------------------------------------

-- | Compute sinc function @sin(x)\/x@
sinc :: Double -> Double
sinc :: Double -> Double
sinc Double
x
  | Double
ax Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
eps_0 = Double
1
  | Double
ax Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
eps_2 = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x2Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
6
  | Double
ax Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
eps_4 = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x2Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
6 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x2Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
120
  | Bool
otherwise  = Double -> Double
forall a. Floating a => a -> a
sin Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x
  where
    ax :: Double
ax    = Double -> Double
forall a. Num a => a -> a
abs Double
x
    x2 :: Double
x2    = Double
xDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
x
    -- For explanation of choice see `doc/sinc.hs'
    eps_0 :: Double
eps_0 = Double
1.8250120749944284e-8 -- sqrt (6ε/4)
    eps_2 :: Double
eps_2 = Double
1.4284346431400855e-4 --   (30ε)**(1/4) / 2
    eps_4 :: Double
eps_4 = Double
4.043633626430947e-3  -- (1206ε)**(1/6) / 2


----------------------------------------------------------------
-- Logarithm
----------------------------------------------------------------

-- | Compute log(1+x)-x:
log1pmx :: Double -> Double
log1pmx :: Double -> Double
log1pmx Double
x
  | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<  -Double
1        = String -> Double
forall a. HasCallStack => String -> a
error String
"Domain error"
  | Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== -Double
1        = Double
m_neg_inf
  | Double
ax Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
0.95      = Double -> Double
forall a. Floating a => a -> a
log(Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x
  | Double
ax Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
m_epsilon = -(Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
2
  | Bool
otherwise      = - Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Sequence Double -> Double
sumPowerSeries (-Double
x) (Double -> Double
forall a. Fractional a => a -> a
recip (Double -> Double) -> Sequence Double -> Sequence Double
forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> Double -> Sequence Double
forall a. Num a => a -> Sequence a
enumSequenceFrom Double
2)
  where
   ax :: Double
ax = Double -> Double
forall a. Num a => a -> a
abs Double
x

-- | /O(log n)/ Compute the logarithm in base 2 of the given value.
log2 :: Int -> Int
log2 :: Int -> Int
log2 Int
v0
    | Int
v0 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= Int
0   = String -> Int
forall a. String -> a
modErr (String -> Int) -> String -> Int
forall a b. (a -> b) -> a -> b
$ String
"log2: nonpositive input, got " String -> String -> String
forall a. [a] -> [a] -> [a]
++ Int -> String
forall a. Show a => a -> String
show Int
v0
    | Bool
otherwise = Int -> Int -> Int -> Int
go Int
5 Int
0 Int
v0
  where
    go :: Int -> Int -> Int -> Int
go !Int
i !Int
r !Int
v | Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== -Int
1        = Int
r
                | Int
v Int -> Int -> Int
forall a. Bits a => a -> a -> a
.&. Int -> Int
b Int
i Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
/= Int
0 = let si :: Int
si = Vector Int -> Int -> Int
forall a. Unbox a => Vector a -> Int -> a
U.unsafeIndex Vector Int
sv Int
i
                                   in Int -> Int -> Int -> Int
go (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) (Int
r Int -> Int -> Int
forall a. Bits a => a -> a -> a
.|. Int
si) (Int
v Int -> Int -> Int
forall a. Bits a => a -> Int -> a
`shiftR` Int
si)
                | Bool
otherwise      = Int -> Int -> Int -> Int
go (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Int
r Int
v
    b :: Int -> Int
b = Vector Int -> Int -> Int
forall a. Unbox a => Vector a -> Int -> a
U.unsafeIndex Vector Int
bv
    !bv :: Vector Int
bv = [Int] -> Vector Int
forall a. Unbox a => [a] -> Vector a
U.fromList [ Int
0x02, Int
0x0c, Int
0xf0, Int
0xff00
                     , Word -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Word
0xffff0000 :: Word)
                     , Word -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Word
0xffffffff00000000 :: Word)]
    !sv :: Vector Int
sv = [Int] -> Vector Int
forall a. Unbox a => [a] -> Vector a
U.fromList [Int
1,Int
2,Int
4,Int
8,Int
16,Int
32]


----------------------------------------------------------------
-- Factorial
----------------------------------------------------------------

-- | 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 :: Int -> Double
factorial Int
n
  | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0     = String -> Double
forall a. HasCallStack => String -> a
error String
"Numeric.SpecFunctions.factorial: negative input"
  | Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
170   = Double
m_pos_inf
  | Bool
otherwise = Vector Double -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
U.unsafeIndex Vector Double
factorialTable Int
n

-- | Compute the natural logarithm of the factorial function.  Gives
--   16 decimal digits of precision.
logFactorial :: Integral a => a -> Double
logFactorial :: forall a. Integral a => a -> Double
logFactorial a
n
  | a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<  a
0    = String -> Double
forall a. HasCallStack => String -> a
error String
"Numeric.SpecFunctions.logFactorial: negative input"
  -- For smaller inputs we just look up table
  | a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= a
170  = Double -> Double
forall a. Floating a => a -> a
log (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Vector Double -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
U.unsafeIndex Vector Double
factorialTable (a -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n)
  -- Otherwise we use asymptotic Stirling's series. Number of terms
  -- necessary depends on the argument.
  | a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< a
1500  = Double
stirling Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
rx Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
12) Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
360)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
rxDouble -> Double -> Double
forall a. Num a => a -> a -> a
*Double
rx)
  | Bool
otherwise = Double
stirling Double -> Double -> Double
forall a. Num a => a -> a -> a
+ (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
12)Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
rx
  where
    stirling :: Double
stirling = (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
m_ln_sqrt_2_pi
    x :: Double
x        = a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1
    rx :: Double
rx       = Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x
{-# SPECIALIZE logFactorial :: Int -> Double #-}


-- | Calculate the error term of the Stirling approximation.  This is
-- only defined for non-negative values.
--
-- \[
-- \operatorname{stirlingError}(n) = \log(n!) - \log(\sqrt{2\pi n}\frac{n}{e}^n)
-- \]
stirlingError :: Double -> Double
stirlingError :: Double -> Double
stirlingError Double
n
  | Double
n Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
15.0   = case Double -> (Int, Double)
forall b. Integral b => Double -> (b, Double)
forall a b. (RealFrac a, Integral b) => a -> (b, a)
properFraction (Double
nDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
n) of
                    (Int
i,Double
0) -> Vector Double
sfe Vector Double -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
`U.unsafeIndex` Int
i
                    (Int, Double)
_     -> Double -> Double
logGamma (Double
nDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1.0) Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
nDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
0.5) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
-
                             Double
m_ln_sqrt_2_pi
  | Double
n Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
500     = Double -> [Double] -> Double
forall a. Num a => a -> [a] -> a
evaluateOddPolynomialL (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
n) [Double
s0,-Double
s1]
  | Double
n Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
80      = Double -> [Double] -> Double
forall a. Num a => a -> [a] -> a
evaluateOddPolynomialL (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
n) [Double
s0,-Double
s1,Double
s2]
  | Double
n Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> Double
35      = Double -> [Double] -> Double
forall a. Num a => a -> [a] -> a
evaluateOddPolynomialL (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
n) [Double
s0,-Double
s1,Double
s2,-Double
s3]
  | Bool
otherwise   = Double -> [Double] -> Double
forall a. Num a => a -> [a] -> a
evaluateOddPolynomialL (Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
n) [Double
s0,-Double
s1,Double
s2,-Double
s3,Double
s4]
  where
    s0 :: Double
s0 = Double
0.083333333333333333333        -- 1/12
    s1 :: Double
s1 = Double
0.00277777777777777777778      -- 1/360
    s2 :: Double
s2 = Double
0.00079365079365079365079365   -- 1/1260
    s3 :: Double
s3 = Double
0.000595238095238095238095238  -- 1/1680
    s4 :: Double
s4 = Double
0.0008417508417508417508417508 -- 1/1188
    sfe :: Vector Double
sfe = [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
U.fromList [ Double
0.0,
                Double
0.1534264097200273452913848,   Double
0.0810614667953272582196702,
                Double
0.0548141210519176538961390,   Double
0.0413406959554092940938221,
                Double
0.03316287351993628748511048,  Double
0.02767792568499833914878929,
                Double
0.02374616365629749597132920,  Double
0.02079067210376509311152277,
                Double
0.01848845053267318523077934,  Double
0.01664469118982119216319487,
                Double
0.01513497322191737887351255,  Double
0.01387612882307074799874573,
                Double
0.01281046524292022692424986,  Double
0.01189670994589177009505572,
                Double
0.01110455975820691732662991,  Double
0.010411265261972096497478567,
                Double
0.009799416126158803298389475, Double
0.009255462182712732917728637,
                Double
0.008768700134139385462952823, Double
0.008330563433362871256469318,
                Double
0.007934114564314020547248100, Double
0.007573675487951840794972024,
                Double
0.007244554301320383179543912, Double
0.006942840107209529865664152,
                Double
0.006665247032707682442354394, Double
0.006408994188004207068439631,
                Double
0.006171712263039457647532867, Double
0.005951370112758847735624416,
                Double
0.005746216513010115682023589, Double
0.005554733551962801371038690 ]


----------------------------------------------------------------
-- Combinatorics
----------------------------------------------------------------

-- |
-- Quickly compute the natural logarithm of /n/ @`choose`@ /k/, with
-- no checking.
--
-- Less numerically stable:
--
-- > exp $ lg (n+1) - lg (k+1) - lg (n-k+1)
-- >   where lg = logGamma . fromIntegral
logChooseFast :: Double -> Double -> Double
logChooseFast :: Double -> Double -> Double
logChooseFast Double
n Double
k = -Double -> Double
forall a. Floating a => a -> a
log (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double -> Double
logBeta (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1)

-- | Calculate binomial coefficient using exact formula
chooseExact :: Int -> Int -> Double
Int
n chooseExact :: Int -> Int -> Double
`chooseExact` Int
k
  = (Double -> Int -> Double) -> Double -> Vector Int -> Double
forall b a. Unbox b => (a -> b -> a) -> a -> Vector b -> a
U.foldl' Double -> Int -> Double
forall {p}. Integral p => Double -> p -> Double
go Double
1 (Vector Int -> Double) -> Vector Int -> Double
forall a b. (a -> b) -> a -> b
$ Int -> Int -> Vector Int
forall a. (Unbox a, Enum a) => a -> a -> Vector a
U.enumFromTo Int
1 Int
k
  where
    go :: Double -> p -> Double
go Double
a p
i      = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
nk Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
j) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
j
        where j :: Double
j = p -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral p
i :: Double
    nk :: Double
nk = Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
n Int -> Int -> Int
forall a. Num a => a -> a -> a
- Int
k)

-- | Compute logarithm of the binomial coefficient.
logChoose :: Int -> Int -> Double
Int
n logChoose :: Int -> Int -> Double
`logChoose` Int
k
    | Int
k  Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n    = (-Double
1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
0
      -- For very large N exact algorithm overflows double so we
      -- switch to beta-function based one
    | Int
k' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
50 Bool -> Bool -> Bool
&& (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
20000000) = Double -> Double
forall a. Floating a => a -> a
log (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Int -> Int -> Double
chooseExact Int
n Int
k'
    | Bool
otherwise                 = Double -> Double -> Double
logChooseFast (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k)
  where
    k' :: Int
k' = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
k (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
k)

-- | Compute the binomial coefficient /n/ @\``choose`\`@ /k/. For
-- values of /k/ > 50, this uses an approximation for performance
-- reasons.  The approximation is accurate to 12 decimal places in the
-- worst case
--
-- Example:
--
-- > 7 `choose` 3 == 35
choose :: Int -> Int -> Double
Int
n choose :: Int -> Int -> Double
`choose` Int
k
    | Int
k  Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n         = Double
0
    | Int
k' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
50        = Int -> Int -> Double
chooseExact Int
n Int
k'
    | Double
approx Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
max64 = Int64 -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int64 -> Double) -> (Double -> Int64) -> Double -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Int64
forall {a}. RealFrac a => a -> Int64
round64 (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double
approx
    | Bool
otherwise      = Double
approx
  where
    k' :: Int
k'             = Int -> Int -> Int
forall a. Ord a => a -> a -> a
min Int
k (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-Int
k)
    approx :: Double
approx         = Double -> Double
forall a. Floating a => a -> a
exp (Double -> Double) -> Double -> Double
forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double
logChooseFast (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) (Int -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k')
    max64 :: Double
max64          = Int64 -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int64
forall a. Bounded a => a
maxBound :: Int64)
    round64 :: a -> Int64
round64 a
x      = a -> Int64
forall b. Integral b => a -> b
forall a b. (RealFrac a, Integral b) => a -> b
round a
x :: Int64

-- | Compute ψ(/x/), the first logarithmic derivative of the gamma
--   function.
--
-- \[
-- \psi(x) = \frac{d}{dx} \ln \left(\Gamma(x)\right) = \frac{\Gamma'(x)}{\Gamma(x)}
-- \]
--
-- Uses Algorithm AS 103 by Bernardo, based on Minka's C implementation.
digamma :: Double -> Double
digamma :: Double -> Double
digamma Double
x
    | Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
x                  = Double
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.
    | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
&& Int64 -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Double -> Int64
forall b. Integral b => Double -> b
forall a b. (RealFrac a, Integral b) => a -> b
truncate Double
x :: Int64) Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
x = Double
m_neg_inf
    -- Jeffery's reflection formula
    | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0     = Double -> Double
digamma (Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double -> Double
forall a. Floating a => a -> a
tan (Double -> Double
forall a. Num a => a -> a
negate Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x)
    | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= Double
1e-6 = - Double
γ Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
trigamma1 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x
    | Double
x' Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
c    = Double
r
    -- De Moivre's expansion
    | Bool
otherwise = let s :: Double
s = Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
x'
                  in  Double -> [Double] -> Double
forall a. Num a => a -> [a] -> a
evaluateEvenPolynomialL Double
s
                        [   Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Floating a => a -> a
log Double
x' Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s
                        , - Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
12
                        ,   Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
120
                        , - Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
252
                        ,   Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
240
                        , - Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
132
                        ,  Double
391Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
32760
                        ]
  where
    γ :: Double
γ  = Double
m_eulerMascheroni
    c :: Double
c  = Double
12
    -- Reduce to digamma (x + n) where (x + n) >= c
    (Double
r, Double
x') = Double -> Double -> (Double, Double)
reduce Double
0 Double
x
      where
        reduce :: Double -> Double -> (Double, Double)
reduce !Double
s Double
y
          | Double
y Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
c     = Double -> Double -> (Double, Double)
reduce (Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
y) (Double
y Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
1)
          | Bool
otherwise = (Double
s, Double
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 :: Vector Double
coefW = [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
U.fromList [ Double
0.0055657196642445571, Double
0.012915947284065419, Double
0.020181515297735382
                   , Double
0.027298621498568734,  Double
0.034213810770299537, Double
0.040875750923643261
                   , Double
0.047235083490265582,  Double
0.053244713977759692, Double
0.058860144245324798
                   , Double
0.064039797355015485,  Double
0.068745323835736408, Double
0.072941885005653087
                   , Double
0.076598410645870640,  Double
0.079687828912071670, Double
0.082187266704339706
                   , Double
0.084078218979661945,  Double
0.085346685739338721, Double
0.085983275670394821
                   ]
coefY :: Vector Double
coefY = [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
U.fromList [ Double
0.0021695375159141994, Double
0.011413521097787704, Double
0.027972308950302116
                   , Double
0.051727015600492421,  Double
0.082502225484340941, Double
0.12007019910960293
                   , Double
0.16415283300752470,   Double
0.21442376986779355,  Double
0.27051082840644336
                   , Double
0.33199876341447887,   Double
0.39843234186401943,  Double
0.46931971407375483
                   , Double
0.54413605556657973,   Double
0.62232745288031077,  Double
0.70331500465597174
                   , Double
0.78649910768313447,   Double
0.87126389619061517,  Double
0.95698180152629142
                   ]
{-# NOINLINE coefW #-}
{-# NOINLINE coefY #-}

trigamma1 :: Double
trigamma1 :: Double
trigamma1 = Double
1.6449340668482264365 -- pi**2 / 6

modErr :: String -> a
modErr :: forall a. String -> a
modErr String
msg = String -> a
forall a. HasCallStack => String -> a
error (String -> a) -> String -> a
forall a b. (a -> b) -> a -> b
$ String
"Numeric.SpecFunctions." String -> String -> String
forall a. [a] -> [a] -> [a]
++ String
msg

factorialTable :: U.Vector Double
{-# NOINLINE factorialTable #-}
factorialTable :: Vector Double
factorialTable = Int -> [Double] -> Vector Double
forall a. Unbox a => Int -> [a] -> Vector a
U.fromListN Int
171
  [ Double
1.0
  , Double
1.0
  , Double
2.0
  , Double
6.0
  , Double
24.0
  , Double
120.0
  , Double
720.0
  , Double
5040.0
  , Double
40320.0
  , Double
362880.0
  , Double
3628800.0
  , Double
3.99168e7
  , Double
4.790016e8
  , Double
6.2270208e9
  , Double
8.71782912e10
  , Double
1.307674368e12
  , Double
2.0922789888e13
  , Double
3.55687428096e14
  , Double
6.402373705728e15
  , Double
1.21645100408832e17
  , Double
2.43290200817664e18
  , Double
5.109094217170944e19
  , Double
1.1240007277776077e21
  , Double
2.5852016738884974e22
  , Double
6.204484017332394e23
  , Double
1.5511210043330984e25
  , Double
4.032914611266056e26
  , Double
1.0888869450418352e28
  , Double
3.0488834461171384e29
  , Double
8.841761993739702e30
  , Double
2.6525285981219103e32
  , Double
8.222838654177922e33
  , Double
2.631308369336935e35
  , Double
8.683317618811886e36
  , Double
2.9523279903960412e38
  , Double
1.0333147966386144e40
  , Double
3.719933267899012e41
  , Double
1.3763753091226343e43
  , Double
5.23022617466601e44
  , Double
2.0397882081197442e46
  , Double
8.159152832478977e47
  , Double
3.3452526613163803e49
  , Double
1.4050061177528798e51
  , Double
6.041526306337383e52
  , Double
2.6582715747884485e54
  , Double
1.1962222086548019e56
  , Double
5.5026221598120885e57
  , Double
2.5862324151116818e59
  , Double
1.2413915592536073e61
  , Double
6.082818640342675e62
  , Double
3.0414093201713376e64
  , Double
1.5511187532873822e66
  , Double
8.065817517094388e67
  , Double
4.2748832840600255e69
  , Double
2.308436973392414e71
  , Double
1.2696403353658275e73
  , Double
7.109985878048634e74
  , Double
4.0526919504877214e76
  , Double
2.3505613312828785e78
  , Double
1.386831185456898e80
  , Double
8.32098711274139e81
  , Double
5.075802138772247e83
  , Double
3.146997326038793e85
  , Double
1.9826083154044399e87
  , Double
1.2688693218588415e89
  , Double
8.24765059208247e90
  , Double
5.44344939077443e92
  , Double
3.647111091818868e94
  , Double
2.4800355424368305e96
  , Double
1.711224524281413e98
  , Double
1.197857166996989e100
  , Double
8.504785885678623e101
  , Double
6.1234458376886085e103
  , Double
4.470115461512684e105
  , Double
3.307885441519386e107
  , Double
2.4809140811395396e109
  , Double
1.88549470166605e111
  , Double
1.4518309202828586e113
  , Double
1.1324281178206297e115
  , Double
8.946182130782974e116
  , Double
7.15694570462638e118
  , Double
5.797126020747368e120
  , Double
4.753643337012841e122
  , Double
3.9455239697206583e124
  , Double
3.314240134565353e126
  , Double
2.81710411438055e128
  , Double
2.422709538367273e130
  , Double
2.1077572983795275e132
  , Double
1.8548264225739844e134
  , Double
1.650795516090846e136
  , Double
1.4857159644817613e138
  , Double
1.352001527678403e140
  , Double
1.2438414054641305e142
  , Double
1.1567725070816416e144
  , Double
1.087366156656743e146
  , Double
1.0329978488239058e148
  , Double
9.916779348709496e149
  , Double
9.619275968248211e151
  , Double
9.426890448883246e153
  , Double
9.332621544394413e155
  , Double
9.332621544394415e157
  , Double
9.425947759838358e159
  , Double
9.614466715035125e161
  , Double
9.902900716486179e163
  , Double
1.0299016745145626e166
  , Double
1.0813967582402908e168
  , Double
1.1462805637347082e170
  , Double
1.2265202031961378e172
  , Double
1.3246418194518288e174
  , Double
1.4438595832024934e176
  , Double
1.5882455415227428e178
  , Double
1.7629525510902446e180
  , Double
1.974506857221074e182
  , Double
2.2311927486598134e184
  , Double
2.543559733472187e186
  , Double
2.9250936934930154e188
  , Double
3.393108684451898e190
  , Double
3.9699371608087206e192
  , Double
4.68452584975429e194
  , Double
5.574585761207606e196
  , Double
6.689502913449126e198
  , Double
8.094298525273443e200
  , Double
9.875044200833601e202
  , Double
1.214630436702533e205
  , Double
1.5061417415111406e207
  , Double
1.8826771768889257e209
  , Double
2.372173242880047e211
  , Double
3.0126600184576594e213
  , Double
3.856204823625804e215
  , Double
4.974504222477286e217
  , Double
6.466855489220473e219
  , Double
8.471580690878819e221
  , Double
1.1182486511960041e224
  , Double
1.4872707060906857e226
  , Double
1.9929427461615188e228
  , Double
2.6904727073180504e230
  , Double
3.6590428819525483e232
  , Double
5.012888748274991e234
  , Double
6.917786472619488e236
  , Double
9.615723196941088e238
  , Double
1.3462012475717523e241
  , Double
1.898143759076171e243
  , Double
2.6953641378881624e245
  , Double
3.8543707171800725e247
  , Double
5.5502938327393044e249
  , Double
8.047926057471992e251
  , Double
1.1749972043909107e254
  , Double
1.7272458904546386e256
  , Double
2.5563239178728654e258
  , Double
3.808922637630569e260
  , Double
5.713383956445854e262
  , Double
8.62720977423324e264
  , Double
1.3113358856834524e267
  , Double
2.0063439050956823e269
  , Double
3.0897696138473508e271
  , Double
4.789142901463393e273
  , Double
7.471062926282894e275
  , Double
1.1729568794264143e278
  , Double
1.8532718694937346e280
  , Double
2.946702272495038e282
  , Double
4.714723635992061e284
  , Double
7.590705053947218e286
  , Double
1.2296942187394494e289
  , Double
2.0044015765453023e291
  , Double
3.287218585534296e293
  , Double
5.423910666131589e295
  , Double
9.003691705778436e297
  , Double
1.5036165148649988e300
  , Double
2.526075744973198e302
  , Double
4.269068009004705e304
  , Double
7.257415615307998e306
  ]


-- [NOTE: incompleteGamma.taylorP]
--
-- Incompltete gamma uses several algorithms for different parts of
-- parameter space. Most troublesome is P(a,x) Taylor series
-- [Temme1994,Eq.5.5] which requires to evaluate rather nasty
-- expression:
--
--       x^a             x^a
--  ------------- = -------------
--  exp(x)·Γ(a+1)   exp(x)·a·Γ(a)
--
--  Conditions:
--    | 0.5<x<1.1  = x < 4/3*a
--    | otherwise  = x < a
--
-- For small `a` computation could be performed directly. However for
-- largish values of `a` it's possible some of factor in the
-- expression overflow. Values below take into account ranges for
-- Taylor P approximation:
--
--  · a > 155    - x^a could overflow
--  · a > 1182.5 - exp(x) could overflow
--
-- Usual way to avoid overflow problem is to perform calculations in
-- the log domain. It however doesn't work very well in this case
-- since we encounter catastrophic cancellations and could easily lose
-- up to 6(!) digits for large `a`.
--
-- So we take another approach and use Stirling approximation with
-- correction (logGammaCorrection).
--
--              x^a               / x·e \^a         1
--  ≈ ------------------------- = | --- | · ----------------
--    exp(x)·sqrt(2πa)·(a/e)^a)   \  a  /   exp(x)·sqrt(2πa)
--
-- We're using this approach as soon as logGammaCorrection starts
-- working (a>10) because we don't have implementation for gamma
-- function and exp(logGamma z) results in errors for large a.
--
-- Once we get into region when exp(x) could overflow we rewrite
-- expression above once more:
--
--  / x·e            \^a     1
--  | --- · e^(-x/a) | · ---------
--  \  a             /   sqrt(2πa)
--
-- This approach doesn't work very well but it's still big improvement
-- over calculations in the log domain.