{-# 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.Word          (Word)
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 forall a. Eq a => a -> a -> Bool
==  Double
1         = Double
m_pos_inf
  | Double
p forall a. Eq a => a -> a -> Bool
== -Double
1         = Double
m_neg_inf
  | Double
p forall a. Ord a => a -> a -> Bool
< Double
1 Bool -> Bool -> Bool
&& Double
p forall a. Ord a => a -> a -> Bool
> -Double
1 = if Double
p forall a. Ord a => a -> a -> Bool
> Double
0 then Double
r else -Double
r
  | Bool
otherwise       = 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     = forall a. Num a => a -> a
abs Double
p
    r :: Double
r      = Double -> Double
step forall a b. (a -> b) -> a -> b
$ Double -> Double
step forall a b. (a -> b) -> a -> b
$ Double -> Double
guessInvErfc forall a b. (a -> b) -> a -> b
$ Double
1 forall a. Num a => a -> a -> a
- Double
pp
    step :: Double -> Double
step Double
x = Double -> Double -> Double
invErfcHalleyStep (Double
pp 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 forall a. Eq a => a -> a -> Bool
== Double
2        = Double
m_neg_inf
  | Double
p forall a. Eq a => a -> a -> Bool
== Double
0        = Double
m_pos_inf
  | Double
p forall a. Ord a => a -> a -> Bool
>Double
0 Bool -> Bool -> Bool
&& Double
p forall a. Ord a => a -> a -> Bool
< Double
2 = if Double
p forall a. Ord a => a -> a -> Bool
<= Double
1 then Double
r else -Double
r
  | Bool
otherwise     = forall a. String -> a
modErr forall a b. (a -> b) -> a -> b
$ String
"invErfc: p must be in [0,2] got " forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show Double
p
  where
    pp :: Double
pp | Double
p forall a. Ord a => a -> a -> Bool
<= Double
1    = Double
p
       | Bool
otherwise = Double
2 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 forall a b. (a -> b) -> a -> b
$ Double -> Double
step 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 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 forall a. Num a => a -> a -> a
* ((Double
2.30753 forall a. Num a => a -> a -> a
+ Double
t forall a. Num a => a -> a -> a
* Double
0.27061) forall a. Fractional a => a -> a -> a
/ (Double
1 forall a. Num a => a -> a -> a
+ Double
t forall a. Num a => a -> a -> a
* (Double
0.99229 forall a. Num a => a -> a -> a
+ Double
t forall a. Num a => a -> a -> a
* Double
0.04481)) forall a. Num a => a -> a -> a
- Double
t)
  where
    t :: Double
t = forall a. Floating a => a -> a
sqrt forall a b. (a -> b) -> a -> b
$ -Double
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log( Double
0.5 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 forall a. Num a => a -> a -> a
+ Double
err forall a. Fractional a => a -> a -> a
/ (Double
1.12837916709551257 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
exp(-Double
x forall a. Num a => a -> a -> a
* Double
x) forall a. Num a => a -> a -> a
- Double
x 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 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 forall a. Ord a => a -> a -> Bool
< Double
m_sqrt_eps = forall a. Floating a => a -> a
log (Double
1forall a. Fractional a => a -> a -> a
/Double
z 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 forall a. Ord a => a -> a -> Bool
< Double
0.5   = Double -> Double -> Double
lgamma1_15 Double
z (Double
z forall a. Num a => a -> a -> a
- Double
1) forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
log Double
z
  | Double
z forall a. Ord a => a -> a -> Bool
< Double
1     = Double -> Double -> Double
lgamma15_2 Double
z (Double
z forall a. Num a => a -> a -> a
- Double
1) forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
log Double
z
  -- Piecewise polynomial approximations
  | Double
z forall a. Ord a => a -> a -> Bool
<= Double
1.5  = Double -> Double -> Double
lgamma1_15 (Double
z forall a. Num a => a -> a -> a
- Double
1) (Double
z forall a. Num a => a -> a -> a
- Double
2)
  | Double
z forall a. Ord a => a -> a -> Bool
< Double
2     = Double -> Double -> Double
lgamma15_2 (Double
z forall a. Num a => a -> a -> a
- Double
1) (Double
z forall a. Num a => a -> a -> a
- Double
2)
  | Double
z 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 forall a. Num a => a -> a -> a
* Double
y forall a. Num a => a -> a -> a
+ Double
r forall a. Num a => a -> a -> a
* ( forall (v :: * -> *) a. (Vector v a, Num a) => a -> v a -> a
evaluatePolynomial Double
zm1 Vector Double
tableLogGamma_1_15P
                 forall a. Fractional a => a -> a -> a
/ 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 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 = 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 = 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 forall a. Num a => a -> a -> a
* Double
y forall a. Num a => a -> a -> a
+ Double
r forall a. Num a => a -> a -> a
* ( forall (v :: * -> *) a. (Vector v a, Num a) => a -> v a -> a
evaluatePolynomial (-Double
zm2) Vector Double
tableLogGamma_15_2P
                 forall a. Fractional a => a -> a -> a
/ 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 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 = 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 = 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 forall a. Num a => a -> a -> a
* Double
y forall a. Num a => a -> a -> a
+ Double
r forall a. Num a => a -> a -> a
* ( forall (v :: * -> *) a. (Vector v a, Num a) => a -> v a -> a
evaluatePolynomial Double
zm2 Vector Double
tableLogGamma_2_3P
                forall a. Fractional a => a -> a -> a
/ 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 forall a. Num a => a -> a -> a
* (Double
z forall a. Num a => a -> a -> a
+ Double
1)
    zm2 :: Double
zm2 = Double
z 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 = 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 = 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 forall a. Ord a => a -> a -> Bool
< Double
3     = Double
acc forall a. Num a => a -> a -> a
+ Double -> Double
lgamma2_3 Double
z
             | Bool
otherwise = Double -> Double -> Double
go (Double
acc forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
log Double
zm1) Double
zm1
             where
               zm1 :: Double
zm1 = Double
z 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
  = (forall a. Floating a => a -> a
log (Double
z forall a. Num a => a -> a -> a
+ Double
g forall a. Num a => a -> a -> a
- Double
0.5) forall a. Num a => a -> a -> a
- Double
1) forall a. Num a => a -> a -> a
* (Double
z forall a. Num a => a -> a -> a
- Double
0.5)
  forall a. Num a => a -> a -> a
+ 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 = 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 forall a. Ord a => a -> a -> Bool
> Double
1     = L -> Double
fini forall a b. (a -> b) -> a -> b
$ 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 forall a b. (a -> b) -> a -> b
$ 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 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 forall a. Num a => a -> a -> a
* Double
x  forall a. Num a => a -> a -> a
+ Double
a) (Double
den forall a. Num a => a -> a -> a
* Double
x  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 forall a. Num a => a -> a -> a
* Double
rx forall a. Num a => a -> a -> a
+ Double
a) (Double
den forall a. Num a => a -> a -> a
* Double
rx forall a. Num a => a -> a -> a
+ Double
b)
    rx :: Double
rx = 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 forall a. Ord a => a -> a -> Bool
< Double
10    = Double
m_NaN
    | Double
x forall a. Ord a => a -> a -> Bool
< Double
big   = forall (v :: * -> *).
Vector v Double =>
Double -> v Double -> Double
chebyshevBroucke (Double
t forall a. Num a => a -> a -> a
* Double
t forall a. Num a => a -> a -> a
* Double
2 forall a. Num a => a -> a -> a
- Double
1) Vector Double
coeffs forall a. Fractional a => a -> a -> a
/ Double
x
    | Bool
otherwise = Double
1 forall a. Fractional a => a -> a -> a
/ (Double
x forall a. Num a => a -> a -> a
* Double
12)
  where
    big :: Double
big    = Double
94906265.62425156
    t :: Double
t      = Double
10 forall a. Fractional a => a -> a -> a
/ Double
x
    coeffs :: Vector Double
coeffs = 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 forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
|| Double
x forall a. Ord a => a -> a -> Bool
< Double
0 = forall a. HasCallStack => String -> a
error
     forall a b. (a -> b) -> a -> b
$ String
"incompleteGamma: Domain error z=" forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show Double
a forall a. [a] -> [a] -> [a]
++ String
" x=" forall a. [a] -> [a] -> [a]
++ forall a. Show a => a -> String
show Double
x
  | Double
x forall a. Eq a => a -> a -> Bool
== Double
0          = Double
0
  | Double
x 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 forall a. Ord a => a -> a -> Bool
< forall a. Floating a => a -> a
sqrt Double
m_epsilon Bool -> Bool -> Bool
&& Double
a forall a. Ord a => a -> a -> Bool
> Double
1
    = Double
xforall a. Floating a => a -> a -> a
**Double
a forall a. Fractional a => a -> a -> a
/ Double
a forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
exp (Double -> Double
logGamma Double
a) forall a. Num a => a -> a -> a
* (Double
1 forall a. Num a => a -> a -> a
- Double
aforall a. Num a => a -> a -> a
*Double
x forall a. Fractional a => a -> a -> a
/ (Double
a forall a. Num a => a -> a -> a
+ Double
1))
  | Double
x forall a. Ord a => a -> a -> Bool
< Double
0.5 = case () of
    ()
_| (-Double
0.4)forall a. Fractional a => a -> a -> a
/forall a. Floating a => a -> a
log Double
x forall a. Ord a => a -> a -> Bool
< Double
a  -> Double
taylorSeriesP
     | Bool
otherwise         -> Double
taylorSeriesComplQ
  | Double
x forall a. Ord a => a -> a -> Bool
< Double
1.1 = case () of
    ()
_| Double
0.75forall a. Num a => a -> a -> a
*Double
x forall a. Ord a => a -> a -> Bool
< Double
a        -> Double
taylorSeriesP
     | Bool
otherwise         -> Double
taylorSeriesComplQ
  | Double
a forall a. Ord a => a -> a -> Bool
> Double
20 Bool -> Bool -> Bool
&& Bool
useTemme    = Double
uniformExpansion
  | Double
x forall a. Num a => a -> a -> a
- (Double
1 forall a. Fractional a => a -> a -> a
/ (Double
3 forall a. Num a => a -> a -> a
* Double
x)) forall a. Ord a => a -> a -> Bool
< Double
a = Double
taylorSeriesP
  | Bool
otherwise             = Double
contFraction
  where
    mu :: Double
mu = (Double
x forall a. Num a => a -> a -> a
- Double
a) forall a. Fractional a => a -> a -> a
/ Double
a
    useTemme :: Bool
useTemme = (Double
a forall a. Ord a => a -> a -> Bool
> Double
200 Bool -> Bool -> Bool
&& Double
20forall a. Fractional a => a -> a -> a
/Double
a forall a. Ord a => a -> a -> Bool
> Double
muforall a. Num a => a -> a -> a
*Double
mu)
            Bool -> Bool -> Bool
|| (forall a. Num a => a -> a
abs Double
mu 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 forall a. Ord a => a -> a -> Bool
< Double
10     = Double
x forall a. Floating a => a -> a -> a
** Double
a
                   forall a. Fractional a => a -> a -> a
/ (forall a. Floating a => a -> a
exp Double
x forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
exp (Double -> Double
logGamma (Double
a forall a. Num a => a -> a -> a
+ Double
1)))
      | Double
a forall a. Ord a => a -> a -> Bool
< Double
1182.5 = (Double
x forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
exp Double
1 forall a. Fractional a => a -> a -> a
/ Double
a) forall a. Floating a => a -> a -> a
** Double
a
                   forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
exp Double
x
                   forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt (Double
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
a)
                   forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
exp (Double -> Double
logGammaCorrection Double
a)
      | Bool
otherwise  = (Double
x forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
exp Double
1 forall a. Fractional a => a -> a -> a
/ Double
a forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
exp (-Double
xforall a. Fractional a => a -> a -> a
/Double
a)) forall a. Floating a => a -> a -> a
** Double
a
                   forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt (Double
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
a)
                   forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
exp (Double -> Double
logGammaCorrection Double
a)
    taylorSeriesP :: Double
taylorSeriesP
      = Double -> Sequence Double -> Double
sumPowerSeries Double
x (forall b a. (b -> a -> b) -> b -> Sequence a -> Sequence b
scanSequence forall a. Fractional a => a -> a -> a
(/) Double
1 forall a b. (a -> b) -> a -> b
$ forall a. Num a => a -> Sequence a
enumSequenceFrom (Double
aforall a. Num a => a -> a -> a
+Double
1))
      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) (forall b a. (b -> a -> b) -> b -> Sequence a -> Sequence b
scanSequence forall a. Fractional a => a -> a -> a
(/) Double
1 (forall a. Num a => a -> Sequence a
enumSequenceFrom Double
1) forall a. Fractional a => a -> a -> a
/ forall a. Num a => a -> Sequence a
enumSequenceFrom Double
a)
      forall a. Num a => a -> a -> a
* Double
xforall a. Floating a => a -> a -> a
**Double
a forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
exp(Double -> Double
logGamma Double
a)
    -- Legendre continued fractions
    contFraction :: Double
contFraction = Double
1 forall a. Num a => a -> a -> a
- ( forall a. Floating a => a -> a
exp ( forall a. Floating a => a -> a
log Double
x forall a. Num a => a -> a -> a
* Double
a forall a. Num a => a -> a -> a
- Double
x forall a. Num a => a -> a -> a
- Double -> Double
logGamma Double
a )
                       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
kforall a. Num a => a -> a -> a
*(Double
aforall a. Num a => a -> a -> a
-Double
k), Double
x forall a. Num a => a -> a -> a
- Double
a forall a. Num a => a -> a -> a
+ Double
2forall a. Num a => a -> a -> a
*Double
k forall a. Num a => a -> a -> a
+ Double
1)) forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> 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 = 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 = forall a. Floating a => a -> a
sqrt (Double
2 forall a. Num a => a -> a -> a
* Double
y) forall a. Num a => a -> a -> a
* 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 forall a. Unbox a => Vector a -> Int -> a
! Int
i) forall a. Num a => a -> a -> a
+ (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
i forall a. Num a => a -> a -> a
+ Double
1)forall a. Num a => a -> a -> a
*Double
bm1 forall a. Fractional a => a -> a -> a
/ Double
a
                                 u' :: Double
u' = Double
eta forall a. Num a => a -> a -> a
* Double
u forall a. Num a => a -> a -> a
+ Double
t
                             in  Double -> Double -> Double -> Int -> Double
loop Double
bm0 Double
t Double
u' (Int
iforall a. Num a => a -> a -> a
-Int
1)
          s_a :: Double
s_a = let n :: Int
n = forall a. Unbox a => Vector a -> Int
U.length Vector Double
fm
                in  Double -> Double -> Double -> Int -> Double
loop (Vector Double
fm forall a. Unbox a => Vector a -> Int -> a
! (Int
nforall a. Num a => a -> a -> a
-Int
1)) (Vector Double
fm forall a. Unbox a => Vector a -> Int -> a
! (Int
nforall a. Num a => a -> a -> a
-Int
2)) Double
0 (Int
nforall a. Num a => a -> a -> a
-Int
3)
                  forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
exp (Double -> Double
logGammaCorrection Double
a)
      in Double
1forall a. Fractional a => a -> a -> a
/Double
2 forall a. Num a => a -> a -> a
* Double -> Double
erfc(-Double
etaforall a. Num a => a -> a -> a
*forall a. Floating a => a -> a
sqrt(Double
aforall a. Fractional a => a -> a -> a
/Double
2)) forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
exp(-(Double
aforall a. Num a => a -> a -> a
*Double
y)) forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
sqrt (Double
2forall a. Num a => a -> a -> a
*forall a. Floating a => a
piforall a. Num a => a -> a -> a
*Double
a) 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 forall a. Ord a => a -> a -> Bool
<= Double
0         =
      forall a. String -> a
modErr forall a b. (a -> b) -> a -> b
$ forall r. PrintfType r => String -> r
printf String
"invIncompleteGamma: a must be positive. a=%g p=%g" Double
a Double
p
  | Double
p forall a. Ord a => a -> a -> Bool
< Double
0 Bool -> Bool -> Bool
|| Double
p forall a. Ord a => a -> a -> Bool
> Double
1 =
      forall a. String -> a
modErr forall a b. (a -> b) -> a -> b
$ 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 forall a. Eq a => a -> a -> Bool
== Double
0         = Double
0
  | Double
p forall a. Eq a => a -> a -> Bool
== Double
1         = Double
1 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 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.
      | forall a. RealFloat a => a -> Bool
isInfinite Double
f'     = Double
0
      | forall a. Num a => a -> a
abs Double
dx forall a. Ord a => a -> a -> Bool
< Double
eps forall a. Num a => a -> a -> a
* Double
x' = Double
x'
      | Bool
otherwise         = Int -> Double -> Double
loop (Int
i 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 forall a. Num a => a -> a -> a
- Double
p
        -- dγ(a,x)/dx
        f' :: Double
f'   | Double
a forall a. Ord a => a -> a -> Bool
> Double
1     = Double
afac forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
exp( -(Double
x forall a. Num a => a -> a -> a
- Double
a1) forall a. Num a => a -> a -> a
+ Double
a1 forall a. Num a => a -> a -> a
* (forall a. Floating a => a -> a
log Double
x forall a. Num a => a -> a -> a
- Double
lna1))
             | Bool
otherwise = forall a. Floating a => a -> a
exp( -Double
x forall a. Num a => a -> a -> a
+ Double
a1 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log Double
x forall a. Num a => a -> a -> a
- Double
gln)
        u :: Double
u    = Double
f forall a. Fractional a => a -> a -> a
/ Double
f'
        -- Halley correction to Newton-Rapson step
        corr :: Double
corr = Double
u forall a. Num a => a -> a -> a
* (Double
a1 forall a. Fractional a => a -> a -> a
/ Double
x forall a. Num a => a -> a -> a
- Double
1)
        dx :: Double
dx   = Double
u forall a. Fractional a => a -> a -> a
/ (Double
1 forall a. Num a => a -> a -> a
- Double
0.5 forall a. Num a => a -> a -> a
* forall a. Ord a => a -> a -> a
min Double
1.0 Double
corr)
        -- New approximation to x
        x' :: Double
x'   | Double
x forall a. Ord a => a -> a -> Bool
< Double
dx    = Double
0.5 forall a. Num a => a -> a -> a
* Double
x -- Do not go below 0
             | Bool
otherwise = Double
x forall a. Num a => a -> a -> a
- Double
dx
    -- Calculate initial guess for root
    guess :: Double
guess
      --
      | Double
a forall a. Ord a => a -> a -> Bool
> Double
1   =
         let t :: Double
t  = forall a. Floating a => a -> a
sqrt forall a b. (a -> b) -> a -> b
$ -Double
2 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log(if Double
p forall a. Ord a => a -> a -> Bool
< Double
0.5 then Double
p else Double
1 forall a. Num a => a -> a -> a
- Double
p)
             x1 :: Double
x1 = (Double
2.30753 forall a. Num a => a -> a -> a
+ Double
t forall a. Num a => a -> a -> a
* Double
0.27061) forall a. Fractional a => a -> a -> a
/ (Double
1 forall a. Num a => a -> a -> a
+ Double
t forall a. Num a => a -> a -> a
* (Double
0.99229 forall a. Num a => a -> a -> a
+ Double
t forall a. Num a => a -> a -> a
* Double
0.04481)) forall a. Num a => a -> a -> a
- Double
t
             x2 :: Double
x2 = if Double
p forall a. Ord a => a -> a -> Bool
< Double
0.5 then -Double
x1 else Double
x1
         in forall a. Ord a => a -> a -> a
max Double
1e-3 (Double
a forall a. Num a => a -> a -> a
* (Double
1 forall a. Num a => a -> a -> a
- Double
1forall a. Fractional a => a -> a -> a
/(Double
9forall a. Num a => a -> a -> a
*Double
a) forall a. Num a => a -> a -> a
- Double
x2 forall a. Fractional a => a -> a -> a
/ (Double
3 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt Double
a)) 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 forall a. Num a => a -> a -> a
- Double
a forall a. Num a => a -> a -> a
* (Double
0.253 forall a. Num a => a -> a -> a
+ Double
aforall a. Num a => a -> a -> a
*Double
0.12)
         in if Double
p forall a. Ord a => a -> a -> Bool
< Double
t
            then (Double
p forall a. Fractional a => a -> a -> a
/ Double
t) forall a. Floating a => a -> a -> a
** (Double
1 forall a. Fractional a => a -> a -> a
/ Double
a)
            else Double
1 forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
log( Double
1 forall a. Num a => a -> a -> a
- (Double
pforall a. Num a => a -> a -> a
-Double
t) forall a. Fractional a => a -> a -> a
/ (Double
1forall a. Num a => a -> a -> a
-Double
t))
    -- Constants
    a1 :: Double
a1   = Double
a forall a. Num a => a -> a -> a
- Double
1
    lna1 :: Double
lna1 = forall a. Floating a => a -> a
log Double
a1
    afac :: Double
afac = forall a. Floating a => a -> a
exp( Double
a1 forall a. Num a => a -> a -> a
* (Double
lna1 forall a. Num a => a -> a -> a
- Double
1) 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 forall a. Ord a => a -> a -> Bool
< Double
0     = Double
m_NaN
  | Double
p forall a. Eq a => a -> a -> Bool
== Double
0    = Double
m_pos_inf
  | Double
p forall a. Ord a => a -> a -> Bool
>= Double
10   = Double
allStirling
  | Double
q 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 forall a. Num a => a -> a -> a
+ (Double -> Double
logGamma Double
q forall a. Num a => a -> a -> a
- Double -> Double
logGamma Double
pq)
  where
    p :: Double
p   = forall a. Ord a => a -> a -> a
min Double
a Double
b
    q :: Double
q   = forall a. Ord a => a -> a -> a
max Double
a Double
b
    ppq :: Double
ppq = Double
p forall a. Fractional a => a -> a -> a
/ Double
pq
    pq :: Double
pq  = Double
p 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
      = forall a. Floating a => a -> a
log Double
q forall a. Num a => a -> a -> a
* (-Double
0.5)
      forall a. Num a => a -> a -> a
+ Double
m_ln_sqrt_2_pi
      forall a. Num a => a -> a -> a
+ Double -> Double
logGammaCorrection Double
p
      forall a. Num a => a -> a -> a
+ (Double -> Double
logGammaCorrection Double
q forall a. Num a => a -> a -> a
- Double -> Double
logGammaCorrection Double
pq)
      forall a. Num a => a -> a -> a
+ (Double
p forall a. Num a => a -> a -> a
- Double
0.5) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log Double
ppq
      forall a. Num a => a -> a -> a
+ Double
q forall a. Num a => a -> a -> a
* 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
      forall a. Num a => a -> a -> a
+ (Double -> Double
logGammaCorrection Double
q forall a. Num a => a -> a -> a
- Double -> Double
logGammaCorrection Double
pq)
      forall a. Num a => a -> a -> a
+ Double
p
      forall a. Num a => a -> a -> a
- Double
p forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log Double
pq
      forall a. Num a => a -> a -> a
+ (Double
q forall a. Num a => a -> a -> a
- Double
0.5) forall a. Num a => a -> a -> a
* 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 forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
|| Double
q forall a. Ord a => a -> a -> Bool
<= Double
0            =
      forall a. String -> a
modErr forall a b. (a -> b) -> a -> b
$ 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 forall a. Ord a => a -> a -> Bool
<  Double
0 Bool -> Bool -> Bool
|| Double
x forall a. Ord a => a -> a -> Bool
>  Double
1 Bool -> Bool -> Bool
|| forall a. RealFloat a => a -> Bool
isNaN Double
x =
      forall a. String -> a
modErr forall a b. (a -> b) -> a -> b
$ 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 forall a. Eq a => a -> a -> Bool
== Double
0 Bool -> Bool -> Bool
|| Double
x forall a. Eq a => a -> a -> Bool
== Double
1            = Double
x
  | Double
p forall a. Ord a => a -> a -> Bool
>= (Double
pforall a. Num a => a -> a -> a
+Double
q) 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 forall a. Num a => a -> a -> a
- Double -> Double -> Double -> Double -> Double
incompleteBetaWorker Double
beta Double
q Double
p (Double
1 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 forall a. Ord a => a -> a -> Bool
> Double
0   = Double
1 forall a. Num a => a -> a -> a
- Double
ans
  | Bool
otherwise = -Double
ans
  where
    -- Constants
    p1 :: Double
p1    = Double
p forall a. Num a => a -> a -> a
- Double
1
    q1 :: Double
q1    = Double
q forall a. Num a => a -> a -> a
- Double
1
    mu :: Double
mu    = Double
p forall a. Fractional a => a -> a -> a
/ (Double
p forall a. Num a => a -> a -> a
+ Double
q)
    lnmu :: Double
lnmu  = forall a. Floating a => a -> a
log     Double
mu
    lnmuc :: Double
lnmuc = forall a. Floating a => a -> a
log1p (-Double
mu)
    -- Upper limit for integration
    xu :: Double
xu = forall a. Ord a => a -> a -> a
max Double
0 forall a b. (a -> b) -> a -> b
$ forall a. Ord a => a -> a -> a
min (Double
mu forall a. Num a => a -> a -> a
- Double
10forall a. Num a => a -> a -> a
*Double
t) (Double
x forall a. Num a => a -> a -> a
- Double
5forall a. Num a => a -> a -> a
*Double
t)
       where
         t :: Double
t = forall a. Floating a => a -> a
sqrt forall a b. (a -> b) -> a -> b
$ Double
pforall a. Num a => a -> a -> a
*Double
q forall a. Fractional a => a -> a -> a
/ ( (Double
pforall a. Num a => a -> a -> a
+Double
q) forall a. Num a => a -> a -> a
* (Double
pforall a. Num a => a -> a -> a
+Double
q) forall a. Num a => a -> a -> a
* (Double
p forall a. Num a => a -> a -> a
+ Double
q 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 forall a. Num a => a -> a -> a
+ (Double
xu forall a. Num a => a -> a -> a
- Double
x) forall a. Num a => a -> a -> a
* Double
y
             in  Double
w forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
exp( Double
p1 forall a. Num a => a -> a -> a
* (forall a. Floating a => a -> a
log Double
t forall a. Num a => a -> a -> a
- Double
lnmu) forall a. Num a => a -> a -> a
+ Double
q1 forall a. Num a => a -> a -> a
* (forall a. Floating a => a -> a
log(Double
1forall a. Num a => a -> a -> a
-Double
t) forall a. Num a => a -> a -> a
- Double
lnmuc) )
    s :: Double
s   = forall a. (Unbox a, Num a) => Vector a -> a
U.sum forall a b. (a -> b) -> a -> b
$ 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 forall a. Num a => a -> a -> a
* (Double
xu forall a. Num a => a -> a -> a
- Double
x) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
exp( Double
p1 forall a. Num a => a -> a -> a
* Double
lnmu forall a. Num a => a -> a -> a
+ Double
q1 forall a. Num a => a -> a -> a
* Double
lnmuc 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 forall a. Ord a => a -> a -> Bool
> Double
3000 Bool -> Bool -> Bool
&& Double
q 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
pforall a. Num a => a -> a -> a
+Double
q) (forall a b. (RealFrac a, Integral b) => a -> b
truncate forall a b. (a -> b) -> a -> b
$ Double
q forall a. Num a => a -> a -> a
+ Double
cx forall a. Num a => a -> a -> a
* (Double
pforall 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 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 forall a. Ord a => a -> a -> Bool
< Double
m_min_log Bool -> Bool -> Bool
|| Double
prod forall a. Ord a => a -> a -> Bool
< Double
m_tiny = forall a. Floating a => a -> a
exp( Double
p forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log Double
x forall a. Num a => a -> a -> a
+ (Double
q forall a. Num a => a -> a -> a
- Double
1) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log Double
cx forall a. Num a => a -> a -> a
- Double
beta)
      | Bool
otherwise                         = Double
prod forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
exp Double
beta
      where
        prod :: Double
prod =  Double
xforall a. Floating a => a -> a -> a
**Double
p forall a. Num a => a -> a -> a
* Double
cxforall a. Floating a => a -> a -> a
**(Double
q 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' forall a. Num a => a -> a -> a
* Double
factor forall a. Fractional a => a -> a -> a
/ Double
p
      | Bool
otherwise = Double -> Int -> Double -> Double -> Double -> Double
loop Double
psq' (Int
ns forall a. Num a => a -> a -> a
- Int
1) (Double
ai forall a. Num a => a -> a -> a
+ Double
1) Double
term' Double
betain'
      where
        -- New values
        term' :: Double
term'   = Double
term forall a. Num a => a -> a -> a
* Double
fact forall a. Fractional a => a -> a -> a
/ (Double
p forall a. Num a => a -> a -> a
+ Double
ai)
        betain' :: Double
betain' = Double
betain forall a. Num a => a -> a -> a
+ Double
term'
        fact :: Double
fact | Int
ns forall a. Ord a => a -> a -> Bool
>  Int
0   = (Double
q forall a. Num a => a -> a -> a
- Double
ai) forall a. Num a => a -> a -> a
* Double
xforall a. Fractional a => a -> a -> a
/Double
cx
             | Int
ns forall a. Eq a => a -> a -> Bool
== Int
0   = (Double
q forall a. Num a => a -> a -> a
- Double
ai) forall a. Num a => a -> a -> a
* Double
x
             | Bool
otherwise = Double
psq forall a. Num a => a -> a -> a
* Double
x
        -- Iterations are complete
        done :: Bool
done = Double
db forall a. Ord a => a -> a -> Bool
<= Double
eps Bool -> Bool -> Bool
&& Double
db forall a. Ord a => a -> a -> Bool
<= Double
epsforall a. Num a => a -> a -> a
*Double
betain' where db :: Double
db = forall a. Num a => a -> a
abs Double
term'
        psq' :: Double
psq' = if Int
ns forall a. Ord a => a -> a -> Bool
< Int
0 then Double
psq 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 forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
|| Double
q forall a. Ord a => a -> a -> Bool
<= Double
0 =
      forall a. String -> a
modErr forall a b. (a -> b) -> a -> b
$ 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 forall a. Ord a => a -> a -> Bool
<  Double
0 Bool -> Bool -> Bool
|| Double
a forall a. Ord a => a -> a -> Bool
>  Double
1 =
      forall a. String -> a
modErr forall a b. (a -> b) -> a -> b
$ 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 forall a. Eq a => a -> a -> Bool
== Double
0 Bool -> Bool -> Bool
|| Double
a 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 = 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 forall a. Num a => a -> a -> a
- Double
1
    b1 :: Double
b1 = Double
b 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 forall a. Eq a => a -> a -> Bool
== Double
0 Bool -> Bool -> Bool
|| Double
x 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.
      | 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 forall a. Ord a => a -> a -> Bool
>= t
10                      = Double
x
      -- Solution converges
      | forall a. Num a => a -> a
abs Double
dx forall a. Ord a => a -> a -> Bool
<= Double
16 forall a. Num a => a -> a -> a
* Double
m_epsilon forall a. Num a => a -> a -> a
* Double
x = Double
x'
      | Bool
otherwise                    = t -> Double -> Double
loop (t
iforall 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 forall a. Num a => a -> a -> a
- Double
p
        f' :: Double
f'  = forall a. Floating a => a -> a
exp forall a b. (a -> b) -> a -> b
$ Double
a1 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log Double
x forall a. Num a => a -> a -> a
+ Double
b1 forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log1p (-Double
x) forall a. Num a => a -> a -> a
- Double
beta
        u :: Double
u   = Double
f 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 forall a. Ord a => a -> a -> Bool
> Double
1     = Double
1
             | Double
d forall a. Ord a => a -> a -> Bool
< -Double
1    = -Double
1
             | forall a. RealFloat a => a -> Bool
isNaN Double
d   = Double
0
             | Bool
otherwise = Double
d
          where
            d :: Double
d = Double
u forall a. Num a => a -> a -> a
* (Double
a1 forall a. Fractional a => a -> a -> a
/ Double
x forall a. Num a => a -> a -> a
- Double
b1 forall a. Fractional a => a -> a -> a
/ (Double
1 forall a. Num a => a -> a -> a
- Double
x))
        dx :: Double
dx  = Double
u forall a. Fractional a => a -> a -> a
/ (Double
1 forall a. Num a => a -> a -> a
- Double
0.5 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 forall a. Ord a => a -> a -> Bool
< Double
0     = Double
x forall a. Fractional a => a -> a -> a
/ Double
2
            | Double
z forall a. Ord a => a -> a -> Bool
> Double
1     = (Double
x forall a. Num a => a -> a -> a
+ Double
1) forall a. Fractional a => a -> a -> a
/ Double
2
            | Bool
otherwise = Double
z
            where z :: Double
z = Double
x 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 forall a. Ord a => a -> a -> Bool
< Double
1 Bool -> Bool -> Bool
&& Double
b forall a. Ord a => a -> a -> Bool
< Double
1 =
    let x_infl :: Double
x_infl = (Double
1 forall a. Num a => a -> a -> a
- Double
a) forall a. Fractional a => a -> a -> a
/ (Double
2 forall a. Num a => a -> a -> a
- Double
a 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 forall a. Ord a => a -> a -> Bool
< Double
p_infl = let xg :: Double
xg = (Double
a forall a. Num a => a -> a -> a
* Double
p     forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
exp Double
beta) forall a. Floating a => a -> a -> a
** (Double
1forall a. Fractional a => a -> a -> a
/Double
a) in Double
xg forall a. Fractional a => a -> a -> a
/ (Double
1forall a. Num a => a -> a -> a
+Double
xg)
          | Bool
otherwise  = let xg :: Double
xg = (Double
b forall a. Num a => a -> a -> a
* (Double
1forall a. Num a => a -> a -> a
-Double
p) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
exp Double
beta) forall a. Floating a => a -> a -> a
** (Double
1forall a. Fractional a => a -> a -> a
/Double
b) in Double
1 forall a. Num a => a -> a -> a
- Double
xgforall a. Fractional a => a -> a -> a
/(Double
1forall 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
aforall a. Num a => a -> a -> a
+Double
b forall a. Ord a => a -> a -> Bool
<= Double
6 Bool -> Bool -> Bool
&& Double
aforall a. Ord a => a -> a -> Bool
>Double
1 Bool -> Bool -> Bool
&& Double
bforall a. Ord a => a -> a -> Bool
>Double
1 =
    let x_infl :: Double
x_infl = (Double
a forall a. Num a => a -> a -> a
- Double
1) forall a. Fractional a => a -> a -> a
/ (Double
a forall a. Num a => a -> a -> a
+ Double
b 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 forall a. Ord a => a -> a -> Bool
< Double
p_infl = forall a. Floating a => a -> a
exp ((forall a. Floating a => a -> a
log(Double
p forall a. Num a => a -> a -> a
* Double
a) forall a. Num a => a -> a -> a
+ Double
beta) forall a. Fractional a => a -> a -> a
/ Double
a)
          | Bool
otherwise  = Double
1 forall a. Num a => a -> a -> a
- forall a. Floating a => a -> a
exp((forall a. Floating a => a -> a
log((Double
1forall a. Num a => a -> a -> a
-Double
p) forall a. Num a => a -> a -> a
* Double
b) forall a. Num a => a -> a -> a
+ Double
beta) 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 forall a. Ord a => a -> a -> Bool
< Double
5 Bool -> Bool -> Bool
&& Double
a forall a. Ord a => a -> a -> Bool
<= Double
1 =
    let x :: Double
x | Double
pforall a. Floating a => a -> a -> a
**(Double
1forall a. Fractional a => a -> a -> a
/Double
a) forall a. Ord a => a -> a -> Bool
< Double
0.5 = (Double
p forall a. Floating a => a -> a -> a
** (Double
a forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
exp Double
beta)) forall a. Floating a => a -> a -> a
** (Double
1forall a. Fractional a => a -> a -> a
/Double
a)
          | Bool
otherwise      = Double
1 forall a. Num a => a -> a -> a
- (Double
1 forall a. Num a => a -> a -> a
- Double
p forall a. Floating a => a -> a -> a
** (Double
b forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
exp Double
beta))forall a. Floating a => a -> a -> a
**(Double
1forall 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
aforall a. Num a => a -> a -> a
+Double
b forall a. Ord a => a -> a -> Bool
> Double
5 Bool -> Bool -> Bool
&&  Double
aforall a. Fractional a => a -> a -> a
/Double
b 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
1forall a. Num a => a -> a -> a
-Double
p) forall a. Fractional a => a -> a -> a
/ Double
a
        mu :: Double
mu   = Double
b forall a. Fractional a => a -> a -> a
/ Double
a            -- Eq. 4.3
        -- A lot of helpers for calculation of
        w :: Double
w    = forall a. Floating a => a -> a
sqrt(Double
1 forall a. Num a => a -> a -> a
+ Double
mu)     -- Eq. 4.9
        w_2 :: Double
w_2  = Double
w forall a. Num a => a -> a -> a
* Double
w
        w_3 :: Double
w_3  = Double
w_2 forall a. Num a => a -> a -> a
* Double
w
        w_4 :: Double
w_4  = Double
w_2 forall a. Num a => a -> a -> a
* Double
w_2
        w_5 :: Double
w_5  = Double
w_3 forall a. Num a => a -> a -> a
* Double
w_2
        w_6 :: Double
w_6  = Double
w_3 forall a. Num a => a -> a -> a
* Double
w_3
        w_7 :: Double
w_7  = Double
w_4 forall a. Num a => a -> a -> a
* Double
w_3
        w_8 :: Double
w_8  = Double
w_4 forall a. Num a => a -> a -> a
* Double
w_4
        w_9 :: Double
w_9  = Double
w_5 forall a. Num a => a -> a -> a
* Double
w_4
        w_10 :: Double
w_10 = Double
w_5 forall a. Num a => a -> a -> a
* Double
w_5
        d :: Double
d    = Double
eta0 forall a. Num a => a -> a -> a
- Double
mu
        d_2 :: Double
d_2  = Double
d forall a. Num a => a -> a -> a
* Double
d
        d_3 :: Double
d_3  = Double
d_2 forall a. Num a => a -> a -> a
* Double
d
        d_4 :: Double
d_4  = Double
d_2 forall a. Num a => a -> a -> a
* Double
d_2
        w1 :: Double
w1   = Double
w forall a. Num a => a -> a -> a
+ Double
1
        w1_2 :: Double
w1_2 = Double
w1 forall a. Num a => a -> a -> a
* Double
w1
        w1_3 :: Double
w1_3 = Double
w1 forall a. Num a => a -> a -> a
* Double
w1_2
        w1_4 :: Double
w1_4 = Double
w1_2 forall a. Num a => a -> a -> a
* Double
w1_2
        -- Evaluation of eq 4.10
        e1 :: Double
e1 = (Double
w forall a. Num a => a -> a -> a
+ Double
2) forall a. Num a => a -> a -> a
* (Double
w forall a. Num a => a -> a -> a
- Double
1) forall a. Fractional a => a -> a -> a
/ (Double
3 forall a. Num a => a -> a -> a
* Double
w)
           forall a. Num a => a -> a -> a
+ (Double
w_3 forall a. Num a => a -> a -> a
+ Double
9 forall a. Num a => a -> a -> a
* Double
w_2 forall a. Num a => a -> a -> a
+ Double
21 forall a. Num a => a -> a -> a
* Double
w forall a. Num a => a -> a -> a
+ Double
5) forall a. Num a => a -> a -> a
* Double
d
             forall a. Fractional a => a -> a -> a
/ (Double
36 forall a. Num a => a -> a -> a
* Double
w_2 forall a. Num a => a -> a -> a
* Double
w1)
           forall a. Num a => a -> a -> a
- (Double
w_4 forall a. Num a => a -> a -> a
- Double
13 forall a. Num a => a -> a -> a
* Double
w_3 forall a. Num a => a -> a -> a
+ Double
69 forall a. Num a => a -> a -> a
* Double
w_2 forall a. Num a => a -> a -> a
+ Double
167 forall a. Num a => a -> a -> a
* Double
w forall a. Num a => a -> a -> a
+ Double
46) forall a. Num a => a -> a -> a
* Double
d_2
             forall a. Fractional a => a -> a -> a
/ (Double
1620 forall a. Num a => a -> a -> a
* Double
w1_2 forall a. Num a => a -> a -> a
* Double
w_3)
           forall a. Num a => a -> a -> a
- (Double
7 forall a. Num a => a -> a -> a
* Double
w_5 forall a. Num a => a -> a -> a
+ Double
21 forall a. Num a => a -> a -> a
* Double
w_4 forall a. Num a => a -> a -> a
+ Double
70 forall a. Num a => a -> a -> a
* Double
w_3 forall a. Num a => a -> a -> a
+ Double
26 forall a. Num a => a -> a -> a
* Double
w_2 forall a. Num a => a -> a -> a
- Double
93 forall a. Num a => a -> a -> a
* Double
w forall a. Num a => a -> a -> a
- Double
31) forall a. Num a => a -> a -> a
* Double
d_3
             forall a. Fractional a => a -> a -> a
/ (Double
6480 forall a. Num a => a -> a -> a
* Double
w1_3 forall a. Num a => a -> a -> a
* Double
w_4)
           forall a. Num a => a -> a -> a
- (Double
75 forall a. Num a => a -> a -> a
* Double
w_6 forall a. Num a => a -> a -> a
+ Double
202 forall a. Num a => a -> a -> a
* Double
w_5 forall a. Num a => a -> a -> a
+ Double
188 forall a. Num a => a -> a -> a
* Double
w_4 forall a. Num a => a -> a -> a
- Double
888 forall a. Num a => a -> a -> a
* Double
w_3 forall a. Num a => a -> a -> a
- Double
1345 forall a. Num a => a -> a -> a
* Double
w_2 forall a. Num a => a -> a -> a
+ Double
118 forall a. Num a => a -> a -> a
* Double
w forall a. Num a => a -> a -> a
+ Double
138) forall a. Num a => a -> a -> a
* Double
d_4
             forall a. Fractional a => a -> a -> a
/ (Double
272160 forall a. Num a => a -> a -> a
* Double
w1_4 forall a. Num a => a -> a -> a
* Double
w_5)
        e2 :: Double
e2 = (Double
28 forall a. Num a => a -> a -> a
* Double
w_4 forall a. Num a => a -> a -> a
+ Double
131 forall a. Num a => a -> a -> a
* Double
w_3 forall a. Num a => a -> a -> a
+ Double
402 forall a. Num a => a -> a -> a
* Double
w_2 forall a. Num a => a -> a -> a
+ Double
581 forall a. Num a => a -> a -> a
* Double
w forall a. Num a => a -> a -> a
+ Double
208) forall a. Num a => a -> a -> a
* (Double
w forall a. Num a => a -> a -> a
- Double
1)
             forall a. Fractional a => a -> a -> a
/ (Double
1620 forall a. Num a => a -> a -> a
* Double
w1 forall a. Num a => a -> a -> a
* Double
w_3)
           forall a. Num a => a -> a -> a
- (Double
35 forall a. Num a => a -> a -> a
* Double
w_6 forall a. Num a => a -> a -> a
- Double
154 forall a. Num a => a -> a -> a
* Double
w_5 forall a. Num a => a -> a -> a
- Double
623 forall a. Num a => a -> a -> a
* Double
w_4 forall a. Num a => a -> a -> a
- Double
1636 forall a. Num a => a -> a -> a
* Double
w_3 forall a. Num a => a -> a -> a
- Double
3983 forall a. Num a => a -> a -> a
* Double
w_2 forall a. Num a => a -> a -> a
- Double
3514 forall a. Num a => a -> a -> a
* Double
w forall a. Num a => a -> a -> a
- Double
925) forall a. Num a => a -> a -> a
* Double
d
             forall a. Fractional a => a -> a -> a
/ (Double
12960 forall a. Num a => a -> a -> a
* Double
w1_2 forall a. Num a => a -> a -> a
* Double
w_4)
           forall a. Num a => a -> a -> a
- ( Double
2132 forall a. Num a => a -> a -> a
* Double
w_7 forall a. Num a => a -> a -> a
+ Double
7915 forall a. Num a => a -> a -> a
* Double
w_6 forall a. Num a => a -> a -> a
+ Double
16821 forall a. Num a => a -> a -> a
* Double
w_5 forall a. Num a => a -> a -> a
+ Double
35066 forall a. Num a => a -> a -> a
* Double
w_4 forall a. Num a => a -> a -> a
+ Double
87490 forall a. Num a => a -> a -> a
* Double
w_3
             forall a. Num a => a -> a -> a
+ Double
141183 forall a. Num a => a -> a -> a
* Double
w_2 forall a. Num a => a -> a -> a
+ Double
95993 forall a. Num a => a -> a -> a
* Double
w forall a. Num a => a -> a -> a
+ Double
21640
             ) forall a. Num a => a -> a -> a
* Double
d_2
             forall a. Fractional a => a -> a -> a
/ (Double
816480 forall a. Num a => a -> a -> a
* Double
w_5 forall a. Num a => a -> a -> a
* Double
w1_3)
           forall a. Num a => a -> a -> a
- ( Double
11053 forall a. Num a => a -> a -> a
* Double
w_8 forall a. Num a => a -> a -> a
+ Double
53308 forall a. Num a => a -> a -> a
* Double
w_7 forall a. Num a => a -> a -> a
+ Double
117010 forall a. Num a => a -> a -> a
* Double
w_6 forall a. Num a => a -> a -> a
+ Double
163924 forall a. Num a => a -> a -> a
* Double
w_5 forall a. Num a => a -> a -> a
+ Double
116188 forall a. Num a => a -> a -> a
* Double
w_4
             forall a. Num a => a -> a -> a
- Double
258428 forall a. Num a => a -> a -> a
* Double
w_3 forall a. Num a => a -> a -> a
- Double
677042 forall a. Num a => a -> a -> a
* Double
w_2 forall a. Num a => a -> a -> a
- Double
481940 forall a. Num a => a -> a -> a
* Double
w forall a. Num a => a -> a -> a
- Double
105497
             ) forall a. Num a => a -> a -> a
* Double
d_3
             forall a. Fractional a => a -> a -> a
/ (Double
14696640 forall a. Num a => a -> a -> a
* Double
w1_4 forall a. Num a => a -> a -> a
* Double
w_6)
        e3 :: Double
e3 = -( (Double
3592 forall a. Num a => a -> a -> a
* Double
w_7 forall a. Num a => a -> a -> a
+ Double
8375 forall a. Num a => a -> a -> a
* Double
w_6 forall a. Num a => a -> a -> a
- Double
1323 forall a. Num a => a -> a -> a
* Double
w_5 forall a. Num a => a -> a -> a
- Double
29198 forall a. Num a => a -> a -> a
* Double
w_4 forall a. Num a => a -> a -> a
- Double
89578 forall a. Num a => a -> a -> a
* Double
w_3
                forall a. Num a => a -> a -> a
- Double
154413 forall a. Num a => a -> a -> a
* Double
w_2 forall a. Num a => a -> a -> a
- Double
116063 forall a. Num a => a -> a -> a
* Double
w forall a. Num a => a -> a -> a
- Double
29632
                ) forall a. Num a => a -> a -> a
* (Double
w forall a. Num a => a -> a -> a
- Double
1)
              )
              forall a. Fractional a => a -> a -> a
/ (Double
816480 forall a. Num a => a -> a -> a
* Double
w_5 forall a. Num a => a -> a -> a
* Double
w1_2)
           forall a. Num a => a -> a -> a
- ( Double
442043 forall a. Num a => a -> a -> a
* Double
w_9 forall a. Num a => a -> a -> a
+ Double
2054169 forall a. Num a => a -> a -> a
* Double
w_8 forall a. Num a => a -> a -> a
+ Double
3803094 forall a. Num a => a -> a -> a
* Double
w_7 forall a. Num a => a -> a -> a
+ Double
3470754 forall a. Num a => a -> a -> a
* Double
w_6 forall a. Num a => a -> a -> a
+ Double
2141568 forall a. Num a => a -> a -> a
* Double
w_5
             forall a. Num a => a -> a -> a
- Double
2393568 forall a. Num a => a -> a -> a
* Double
w_4 forall a. Num a => a -> a -> a
- Double
19904934 forall a. Num a => a -> a -> a
* Double
w_3 forall a. Num a => a -> a -> a
- Double
34714674 forall a. Num a => a -> a -> a
* Double
w_2 forall a. Num a => a -> a -> a
- Double
23128299 forall a. Num a => a -> a -> a
* Double
w forall a. Num a => a -> a -> a
- Double
5253353
             ) forall a. Num a => a -> a -> a
* Double
d
             forall a. Fractional a => a -> a -> a
/ (Double
146966400 forall a. Num a => a -> a -> a
* Double
w_6 forall a. Num a => a -> a -> a
* Double
w1_3)
           forall a. Num a => a -> a -> a
- ( Double
116932 forall a. Num a => a -> a -> a
* Double
w_10 forall a. Num a => a -> a -> a
+ Double
819281 forall a. Num a => a -> a -> a
* Double
w_9 forall a. Num a => a -> a -> a
+ Double
2378172 forall a. Num a => a -> a -> a
* Double
w_8 forall a. Num a => a -> a -> a
+ Double
4341330 forall a. Num a => a -> a -> a
* Double
w_7 forall a. Num a => a -> a -> a
+ Double
6806004 forall a. Num a => a -> a -> a
* Double
w_6
             forall a. Num a => a -> a -> a
+ Double
10622748 forall a. Num a => a -> a -> a
* Double
w_5 forall a. Num a => a -> a -> a
+ Double
18739500 forall a. Num a => a -> a -> a
* Double
w_4 forall a. Num a => a -> a -> a
+ Double
30651894 forall a. Num a => a -> a -> a
* Double
w_3 forall a. Num a => a -> a -> a
+ Double
30869976 forall a. Num a => a -> a -> a
* Double
w_2
             forall a. Num a => a -> a -> a
+ Double
15431867 forall a. Num a => a -> a -> a
* Double
w forall a. Num a => a -> a -> a
+ Double
2919016
             ) forall a. Num a => a -> a -> a
* Double
d_2
             forall a. Fractional a => a -> a -> a
/ (Double
146966400 forall a. Num a => a -> a -> a
* Double
w1_4 forall a. Num a => a -> a -> a
* Double
w_7)
        eta :: Double
eta = forall a. Num a => a -> [a] -> a
evaluatePolynomialL (Double
1forall 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 forall a. Num a => a -> a -> a
- Double
mu forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log Double
eta forall a. Num a => a -> a -> a
+ (Double
1 forall a. Num a => a -> a -> a
+ Double
mu) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log(Double
1 forall a. Num a => a -> a -> a
+ Double
mu) forall a. Num a => a -> a -> a
- Double
mu
        cross :: Double
cross   = Double
1 forall a. Fractional a => a -> a -> a
/ (Double
1 forall a. Num a => a -> a -> a
+ Double
mu);
        lower :: Double
lower   = if Double
eta forall a. Ord a => a -> a -> Bool
< Double
mu then Double
cross else Double
0
        upper :: Double
upper   = if Double
eta forall a. Ord a => a -> a -> Bool
< Double
mu then Double
1     else Double
cross
        x_guess :: Double
x_guess = (Double
lower forall a. Num a => a -> a -> a
+ Double
upper) forall a. Fractional a => a -> a -> a
/ Double
2
        func :: Double -> (Double, Double)
func Double
x  = ( Double
u forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
log Double
x forall a. Num a => a -> a -> a
+ Double
muforall a. Num a => a -> a -> a
*forall a. Floating a => a -> a
log(Double
1 forall a. Num a => a -> a -> a
- Double
x)
                  , Double
1forall a. Fractional a => a -> a -> a
/Double
x forall a. Num a => a -> a -> a
- Double
muforall a. Fractional a => a -> a -> a
/(Double
1forall a. Num a => a -> a -> a
-Double
x)
                  )
        Root Double
x0 = NewtonParam
-> (Double, Double, Double)
-> (Double -> (Double, Double))
-> Root Double
newtonRaphson forall a. Default a => a
def{newtonTol :: Tolerance
newtonTol=Double -> Tolerance
RelTol Double
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 forall a. Ord a => a -> a -> Bool
> Double
1 Bool -> Bool -> Bool
&& Double
b forall a. Ord a => a -> a -> Bool
> Double
1 =
      let r :: Double
r = (Double
yforall a. Num a => a -> a -> a
*Double
y forall a. Num a => a -> a -> a
- Double
3) forall a. Fractional a => a -> a -> a
/ Double
6
          s :: Double
s = Double
1 forall a. Fractional a => a -> a -> a
/ (Double
2forall a. Num a => a -> a -> a
*Double
a forall a. Num a => a -> a -> a
- Double
1)
          t :: Double
t = Double
1 forall a. Fractional a => a -> a -> a
/ (Double
2forall a. Num a => a -> a -> a
*Double
b forall a. Num a => a -> a -> a
- Double
1)
          h :: Double
h = Double
2 forall a. Fractional a => a -> a -> a
/ (Double
s forall a. Num a => a -> a -> a
+ Double
t)
          w :: Double
w = Double
y forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt(Double
h forall a. Num a => a -> a -> a
+ Double
r) forall a. Fractional a => a -> a -> a
/ Double
h forall a. Num a => a -> a -> a
- (Double
t forall a. Num a => a -> a -> a
- Double
s) forall a. Num a => a -> a -> a
* (Double
r forall a. Num a => a -> a -> a
+ Double
5forall a. Fractional a => a -> a -> a
/Double
6 forall a. Num a => a -> a -> a
- Double
2 forall a. Fractional a => a -> a -> a
/ (Double
3 forall a. Num a => a -> a -> a
* Double
h))
      in Double
a forall a. Fractional a => a -> a -> a
/ (Double
a forall a. Num a => a -> a -> a
+ Double
b forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
exp(Double
2 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 forall a. Ord a => a -> a -> Bool
> Double
0 Bool -> Bool -> Bool
&& Double
ratio forall a. Ord a => a -> a -> Bool
> Double
1 = Double
1 forall a. Num a => a -> a -> a
- Double
2 forall a. Fractional a => a -> a -> a
/ (Double
ratio 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 forall a. Ord a => a -> a -> Bool
< Double
t forall a. Fractional a => a -> a -> a
/ Double
w  -> (Double
a forall a. Num a => a -> a -> a
* Double
p forall a. Num a => a -> a -> a
* Double
w) forall a. Floating a => a -> a -> a
** (Double
1forall a. Fractional a => a -> a -> a
/Double
a)
       | Bool
otherwise  -> Double
1 forall a. Num a => a -> a -> a
- (Double
b forall a. Num a => a -> a -> a
* (Double
1 forall a. Num a => a -> a -> a
- Double
p) forall a. Num a => a -> a -> a
* Double
w) forall a. Floating a => a -> a -> a
** (Double
1forall a. Fractional a => a -> a -> a
/Double
b)
       where
         lna :: Double
lna = forall a. Floating a => a -> a
log forall a b. (a -> b) -> a -> b
$ Double
a forall a. Fractional a => a -> a -> a
/ (Double
aforall a. Num a => a -> a -> a
+Double
b)
         lnb :: Double
lnb = forall a. Floating a => a -> a
log forall a b. (a -> b) -> a -> b
$ Double
b forall a. Fractional a => a -> a -> a
/ (Double
aforall a. Num a => a -> a -> a
+Double
b)
         t :: Double
t   = forall a. Floating a => a -> a
exp( Double
a forall a. Num a => a -> a -> a
* Double
lna ) forall a. Fractional a => a -> a -> a
/ Double
a
         u :: Double
u   = forall a. Floating a => a -> a
exp( Double
b forall a. Num a => a -> a -> a
* Double
lnb ) forall a. Fractional a => a -> a -> a
/ Double
b
         w :: Double
w   = Double
t forall a. Num a => a -> a -> a
+ Double
u
  where
    -- Formula [AS64 2]
    ratio :: Double
ratio = (Double
4forall a. Num a => a -> a -> a
*Double
a forall a. Num a => a -> a -> a
+ Double
2forall a. Num a => a -> a -> a
*Double
b forall a. Num a => a -> a -> a
- Double
2) forall a. Fractional a => a -> a -> a
/ Double
chi2
    -- Quantile of chi-squared distribution. Formula [AS64 3].
    chi2 :: Double
chi2 = Double
2 forall a. Num a => a -> a -> a
* Double
b forall a. Num a => a -> a -> a
* (Double
1 forall a. Num a => a -> a -> a
- Double
t forall a. Num a => a -> a -> a
+ Double
y forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
sqrt Double
t) forall a. Floating a => a -> a -> a
** Double
3
      where
        t :: Double
t   = Double
1 forall a. Fractional a => a -> a -> a
/ (Double
9 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 forall a. Num a => a -> a -> a
- ( Double
2.30753 forall a. Num a => a -> a -> a
+ Double
0.27061 forall a. Num a => a -> a -> a
* Double
r )
              forall a. Fractional a => a -> a -> a
/ ( Double
1.0 forall a. Num a => a -> a -> a
+ ( Double
0.99229 forall a. Num a => a -> a -> a
+ Double
0.04481 forall a. Num a => a -> a -> a
* Double
r ) forall a. Num a => a -> a -> a
* Double
r )
      where
        r :: Double
r = forall a. Floating a => a -> a
sqrt forall a b. (a -> b) -> a -> b
$ - Double
2 forall a. Num a => a -> a -> a
* 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 forall a. Ord a => a -> a -> Bool
< Double
eps_0 = Double
1
  | Double
ax forall a. Ord a => a -> a -> Bool
< Double
eps_2 = Double
1 forall a. Num a => a -> a -> a
- Double
x2forall a. Fractional a => a -> a -> a
/Double
6
  | Double
ax forall a. Ord a => a -> a -> Bool
< Double
eps_4 = Double
1 forall a. Num a => a -> a -> a
- Double
x2forall a. Fractional a => a -> a -> a
/Double
6 forall a. Num a => a -> a -> a
+ Double
x2forall a. Num a => a -> a -> a
*Double
x2forall a. Fractional a => a -> a -> a
/Double
120
  | Bool
otherwise  = forall a. Floating a => a -> a
sin Double
x forall a. Fractional a => a -> a -> a
/ Double
x
  where
    ax :: Double
ax    = forall a. Num a => a -> a
abs Double
x
    x2 :: Double
x2    = Double
xforall 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 forall a. Ord a => a -> a -> Bool
<  -Double
1        = forall a. HasCallStack => String -> a
error String
"Domain error"
  | Double
x forall a. Eq a => a -> a -> Bool
== -Double
1        = Double
m_neg_inf
  | Double
ax forall a. Ord a => a -> a -> Bool
> Double
0.95      = forall a. Floating a => a -> a
log(Double
1 forall a. Num a => a -> a -> a
+ Double
x) forall a. Num a => a -> a -> a
- Double
x
  | Double
ax forall a. Ord a => a -> a -> Bool
< Double
m_epsilon = -(Double
x forall a. Num a => a -> a -> a
* Double
x) forall a. Fractional a => a -> a -> a
/Double
2
  | Bool
otherwise      = - Double
x forall a. Num a => a -> a -> a
* Double
x forall a. Num a => a -> a -> a
* Double -> Sequence Double -> Double
sumPowerSeries (-Double
x) (forall a. Fractional a => a -> a
recip forall (f :: * -> *) a b. Functor f => (a -> b) -> f a -> f b
<$> forall a. Num a => a -> Sequence a
enumSequenceFrom Double
2)
  where
   ax :: Double
ax = 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 forall a. Ord a => a -> a -> Bool
<= Int
0   = forall a. String -> a
modErr forall a b. (a -> b) -> a -> b
$ String
"log2: nonpositive input, got " forall a. [a] -> [a] -> [a]
++ 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 forall a. Eq a => a -> a -> Bool
== -Int
1        = Int
r
                | Int
v forall a. Bits a => a -> a -> a
.&. Int -> Int
b Int
i forall a. Eq a => a -> a -> Bool
/= Int
0 = let si :: Int
si = forall a. Unbox a => Vector a -> Int -> a
U.unsafeIndex Vector Int
sv Int
i
                                   in Int -> Int -> Int -> Int
go (Int
iforall a. Num a => a -> a -> a
-Int
1) (Int
r forall a. Bits a => a -> a -> a
.|. Int
si) (Int
v forall a. Bits a => a -> Int -> a
`shiftR` Int
si)
                | Bool
otherwise      = Int -> Int -> Int -> Int
go (Int
iforall a. Num a => a -> a -> a
-Int
1) Int
r Int
v
    b :: Int -> Int
b = forall a. Unbox a => Vector a -> Int -> a
U.unsafeIndex Vector Int
bv
    !bv :: Vector Int
bv = forall a. Unbox a => [a] -> Vector a
U.fromList [ Int
0x02, Int
0x0c, Int
0xf0, Int
0xff00
                     , forall a b. (Integral a, Num b) => a -> b
fromIntegral (Word
0xffff0000 :: Word)
                     , forall a b. (Integral a, Num b) => a -> b
fromIntegral (Word
0xffffffff00000000 :: Word)]
    !sv :: Vector Int
sv = 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 forall a. Ord a => a -> a -> Bool
< Int
0     = forall a. HasCallStack => String -> a
error String
"Numeric.SpecFunctions.factorial: negative input"
  | Int
n forall a. Ord a => a -> a -> Bool
> Int
170   = Double
m_pos_inf
  | Bool
otherwise = 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 forall a. Ord a => a -> a -> Bool
<  a
0    = forall a. HasCallStack => String -> a
error String
"Numeric.SpecFunctions.logFactorial: negative input"
  -- For smaller inputs we just look up table
  | a
n forall a. Ord a => a -> a -> Bool
<= a
170  = forall a. Floating a => a -> a
log forall a b. (a -> b) -> a -> b
$ forall a. Unbox a => Vector a -> Int -> a
U.unsafeIndex Vector Double
factorialTable (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 forall a. Ord a => a -> a -> Bool
< a
1500  = Double
stirling forall a. Num a => a -> a -> a
+ Double
rx forall a. Num a => a -> a -> a
* ((Double
1forall a. Fractional a => a -> a -> a
/Double
12) forall a. Num a => a -> a -> a
- (Double
1forall a. Fractional a => a -> a -> a
/Double
360)forall a. Num a => a -> a -> a
*Double
rxforall a. Num a => a -> a -> a
*Double
rx)
  | Bool
otherwise = Double
stirling forall a. Num a => a -> a -> a
+ (Double
1forall a. Fractional a => a -> a -> a
/Double
12)forall a. Num a => a -> a -> a
*Double
rx
  where
    stirling :: Double
stirling = (Double
x forall a. Num a => a -> a -> a
- Double
0.5) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log Double
x forall a. Num a => a -> a -> a
- Double
x forall a. Num a => a -> a -> a
+ Double
m_ln_sqrt_2_pi
    x :: Double
x        = forall a b. (Integral a, Num b) => a -> b
fromIntegral a
n forall a. Num a => a -> a -> a
+ Double
1
    rx :: Double
rx       = Double
1 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 forall a. Ord a => a -> a -> Bool
<= Double
15.0   = case forall a b. (RealFrac a, Integral b) => a -> (b, a)
properFraction (Double
nforall a. Num a => a -> a -> a
+Double
n) of
                    (Int
i,Double
0) -> Vector Double
sfe forall a. Unbox a => Vector a -> Int -> a
`U.unsafeIndex` Int
i
                    (Int, Double)
_     -> Double -> Double
logGamma (Double
nforall a. Num a => a -> a -> a
+Double
1.0) forall a. Num a => a -> a -> a
- (Double
nforall a. Num a => a -> a -> a
+Double
0.5) forall a. Num a => a -> a -> a
* forall a. Floating a => a -> a
log Double
n forall a. Num a => a -> a -> a
+ Double
n forall a. Num a => a -> a -> a
-
                             Double
m_ln_sqrt_2_pi
  | Double
n forall a. Ord a => a -> a -> Bool
> Double
500     = forall a. Num a => a -> [a] -> a
evaluateOddPolynomialL (Double
1forall a. Fractional a => a -> a -> a
/Double
n) [Double
s0,-Double
s1]
  | Double
n forall a. Ord a => a -> a -> Bool
> Double
80      = forall a. Num a => a -> [a] -> a
evaluateOddPolynomialL (Double
1forall a. Fractional a => a -> a -> a
/Double
n) [Double
s0,-Double
s1,Double
s2]
  | Double
n forall a. Ord a => a -> a -> Bool
> Double
35      = forall a. Num a => a -> [a] -> a
evaluateOddPolynomialL (Double
1forall a. Fractional a => a -> a -> a
/Double
n) [Double
s0,-Double
s1,Double
s2,-Double
s3]
  | Bool
otherwise   = forall a. Num a => a -> [a] -> a
evaluateOddPolynomialL (Double
1forall 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 = 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 = -forall a. Floating a => a -> a
log (Double
n forall a. Num a => a -> a -> a
+ Double
1) forall a. Num a => a -> a -> a
- Double -> Double -> Double
logBeta (Double
n forall a. Num a => a -> a -> a
- Double
k forall a. Num a => a -> a -> a
+ Double
1) (Double
k 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
  = forall b a. Unbox b => (a -> b -> a) -> a -> Vector b -> a
U.foldl' forall {p}. Integral p => Double -> p -> Double
go Double
1 forall a b. (a -> b) -> a -> b
$ 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 forall a. Num a => a -> a -> a
* (Double
nk forall a. Num a => a -> a -> a
+ Double
j) forall a. Fractional a => a -> a -> a
/ Double
j
        where j :: Double
j = forall a b. (Integral a, Num b) => a -> b
fromIntegral p
i :: Double
    nk :: Double
nk = forall a b. (Integral a, Num b) => a -> b
fromIntegral (Int
n 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  forall a. Ord a => a -> a -> Bool
> Int
n    = (-Double
1) 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' forall a. Ord a => a -> a -> Bool
< Int
50 Bool -> Bool -> Bool
&& (Int
n forall a. Ord a => a -> a -> Bool
< Int
20000000) = forall a. Floating a => a -> a
log forall a b. (a -> b) -> a -> b
$ Int -> Int -> Double
chooseExact Int
n Int
k'
    | Bool
otherwise                 = Double -> Double -> Double
logChooseFast (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k)
  where
    k' :: Int
k' = forall a. Ord a => a -> a -> a
min Int
k (Int
nforall 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  forall a. Ord a => a -> a -> Bool
> Int
n         = Double
0
    | Int
k' forall a. Ord a => a -> a -> Bool
< Int
50        = Int -> Int -> Double
chooseExact Int
n Int
k'
    | Double
approx forall a. Ord a => a -> a -> Bool
< Double
max64 = forall a b. (Integral a, Num b) => a -> b
fromIntegral forall b c a. (b -> c) -> (a -> b) -> a -> c
. forall {a}. RealFrac a => a -> Int64
round64 forall a b. (a -> b) -> a -> b
$ Double
approx
    | Bool
otherwise      = Double
approx
  where
    k' :: Int
k'             = forall a. Ord a => a -> a -> a
min Int
k (Int
nforall a. Num a => a -> a -> a
-Int
k)
    approx :: Double
approx         = forall a. Floating a => a -> a
exp forall a b. (a -> b) -> a -> b
$ Double -> Double -> Double
logChooseFast (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n) (forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k')
    max64 :: Double
max64          = forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall a. Bounded a => a
maxBound :: Int64)
    round64 :: a -> Int64
round64 a
x      = 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
    | forall a. RealFloat a => a -> Bool
isNaN Double
x Bool -> Bool -> 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 forall a. Ord a => a -> a -> Bool
<= Double
0 Bool -> Bool -> Bool
&& forall a b. (Integral a, Num b) => a -> b
fromIntegral (forall a b. (RealFrac a, Integral b) => a -> b
truncate Double
x :: Int64) forall a. Eq a => a -> a -> Bool
== Double
x = Double
m_neg_inf
    -- Jeffery's reflection formula
    | Double
x forall a. Ord a => a -> a -> Bool
< Double
0     = Double -> Double
digamma (Double
1 forall a. Num a => a -> a -> a
- Double
x) forall a. Num a => a -> a -> a
+ forall a. Floating a => a
pi forall a. Fractional a => a -> a -> a
/ forall a. Floating a => a -> a
tan (forall a. Num a => a -> a
negate forall a. Floating a => a
pi forall a. Num a => a -> a -> a
* Double
x)
    | Double
x forall a. Ord a => a -> a -> Bool
<= Double
1e-6 = - Double
γ forall a. Num a => a -> a -> a
- Double
1forall a. Fractional a => a -> a -> a
/Double
x forall a. Num a => a -> a -> a
+ Double
trigamma1 forall a. Num a => a -> a -> a
* Double
x
    | Double
x' forall a. Ord a => a -> a -> Bool
< Double
c    = Double
r
    -- De Moivre's expansion
    | Bool
otherwise = let s :: Double
s = Double
1forall a. Fractional a => a -> a -> a
/Double
x'
                  in  forall a. Num a => a -> [a] -> a
evaluateEvenPolynomialL Double
s
                        [   Double
r forall a. Num a => a -> a -> a
+ forall a. Floating a => a -> a
log Double
x' forall a. Num a => a -> a -> a
- Double
0.5 forall a. Num a => a -> a -> a
* Double
s
                        , - Double
1forall a. Fractional a => a -> a -> a
/Double
12
                        ,   Double
1forall a. Fractional a => a -> a -> a
/Double
120
                        , - Double
1forall a. Fractional a => a -> a -> a
/Double
252
                        ,   Double
1forall a. Fractional a => a -> a -> a
/Double
240
                        , - Double
1forall a. Fractional a => a -> a -> a
/Double
132
                        ,  Double
391forall 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 forall a. Ord a => a -> a -> Bool
< Double
c     = Double -> Double -> (Double, Double)
reduce (Double
s forall a. Num a => a -> a -> a
- Double
1 forall a. Fractional a => a -> a -> a
/ Double
y) (Double
y 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 = 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 = 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 = forall a. HasCallStack => String -> a
error forall a b. (a -> b) -> a -> b
$ String
"Numeric.SpecFunctions." forall a. [a] -> [a] -> [a]
++ String
msg

factorialTable :: U.Vector Double
{-# NOINLINE factorialTable #-}
factorialTable :: Vector Double
factorialTable = 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.