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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    -- Adapted from Numerical Recipes §6.2.1

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

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

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

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

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

    -- Approximation of incomplete beta by quadrature.
    --
    -- Note that x =< p/(p+q)
    incompleteBetaApprox :: Double -> Double -> Double -> Double -> Double
    incompleteBetaApprox :: Double -> Double -> Double -> Double -> Double
    incompleteBetaApprox Double
    beta Double
    p Double
    q Double
    x
      | Double
    ans Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    > Double
    0   = Double
    1 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    ans
      | Bool
    otherwise = -Double
    ans
      where
        -- Constants
        p1 :: Double
    p1    = Double
    p Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    1
        q1 :: Double
    q1    = Double
    q Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    1
        mu :: Double
    mu    = Double
    p Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / (Double
    p Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    q)
        lnmu :: Double
    lnmu  = Double -> Double
    forall a. Floating a => a -> a
    log     Double
    mu
        lnmuc :: Double
    lnmuc = Double -> Double
    forall a. Floating a => a -> a
    log1p (-Double
    mu)
        -- Upper limit for integration
        xu :: Double
    xu = Double -> Double -> Double
    forall a. Ord a => a -> a -> a
    max Double
    0 (Double -> Double) -> Double -> Double
    forall a b. (a -> b) -> a -> b
    $ Double -> Double -> Double
    forall a. Ord a => a -> a -> a
    min (Double
    mu Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    10Double -> Double -> Double
    forall a. Num a => a -> a -> a
    *Double
    t) (Double
    x Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    5Double -> Double -> Double
    forall a. Num a => a -> a -> a
    *Double
    t)
           where
             t :: Double
    t = Double -> Double
    forall a. Floating a => a -> a
    sqrt (Double -> Double) -> Double -> Double
    forall a b. (a -> b) -> a -> b
    $ Double
    pDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    *Double
    q Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / ( (Double
    pDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    +Double
    q) Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * (Double
    pDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    +Double
    q) Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * (Double
    p Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    q Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    1) )
        -- Calculate incomplete beta by quadrature
        go :: Double -> Double -> Double
    go Double
    y Double
    w = let t :: Double
    t = Double
    x Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + (Double
    xu Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    x) Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    y
                 in  Double
    w Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double -> Double
    forall a. Floating a => a -> a
    exp( Double
    p1 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * (Double -> Double
    forall a. Floating a => a -> a
    log Double
    t Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    lnmu) Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    q1 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * (Double -> Double
    forall a. Floating a => a -> a
    log(Double
    1Double -> Double -> Double
    forall a. Num a => a -> a -> a
    -Double
    t) Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    lnmuc) )
        s :: Double
    s   = Vector Double -> Double
    forall a. (Unbox a, Num a) => Vector a -> a
    U.sum (Vector Double -> Double) -> Vector Double -> Double
    forall a b. (a -> b) -> a -> b
    $ (Double -> Double -> Double)
    -> Vector Double -> Vector Double -> Vector Double
    forall a b c.
    (Unbox a, Unbox b, Unbox c) =>
    (a -> b -> c) -> Vector a -> Vector b -> Vector c
    U.zipWith Double -> Double -> Double
    go Vector Double
    coefY Vector Double
    coefW
        ans :: Double
    ans = Double
    s Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * (Double
    xu Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    x) Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double -> Double
    forall a. Floating a => a -> a
    exp( Double
    p1 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    lnmu Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    q1 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    lnmuc Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    beta )

    -- Worker for incomplete beta function. It is separate function to
    -- avoid confusion with parameter during parameter swapping
    incompleteBetaWorker :: Double -> Double -> Double -> Double -> Double
    incompleteBetaWorker :: Double -> Double -> Double -> Double -> Double
    incompleteBetaWorker Double
    beta Double
    p Double
    q Double
    x
      -- For very large p and q this method becomes very slow so another
      -- method is used.
      | Double
    p Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    > Double
    3000 Bool -> Bool -> Bool
    && Double
    q Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    > Double
    3000 = Double -> Double -> Double -> Double -> Double
    incompleteBetaApprox Double
    beta Double
    p Double
    q Double
    x
      | Bool
    otherwise            = Double -> Int -> Double -> Double -> Double -> Double
    loop (Double
    pDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    +Double
    q) (Double -> Int
    forall b. Integral b => Double -> b
    forall a b. (RealFrac a, Integral b) => a -> b
    truncate (Double -> Int) -> Double -> Int
    forall a b. (a -> b) -> a -> b
    $ Double
    q Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    cx Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * (Double
    pDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    +Double
    q)) Double
    1 Double
    1 Double
    1
      where
        -- Constants
        eps :: Double
    eps = Double
    1e-15
        cx :: Double
    cx  = Double
    1 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    x
        -- Common multiplies for expansion. Accurate calculation is a bit
        -- tricky. Performing calculation in log-domain leads to slight
        -- loss of precision for small x, while using ** prone to
        -- underflows.
        --
        -- If either beta function of x**p·(1-x)**(q-1) underflows we
        -- switch to log domain. It could waste work but there's no easy
        -- switch criterion.
        factor :: Double
    factor
          | Double
    beta Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    < Double
    m_min_log Bool -> Bool -> Bool
    || Double
    prod Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    < Double
    m_tiny = Double -> Double
    forall a. Floating a => a -> a
    exp( Double
    p Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double -> Double
    forall a. Floating a => a -> a
    log Double
    x Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + (Double
    q Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    1) Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double -> Double
    forall a. Floating a => a -> a
    log Double
    cx Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    beta)
          | Bool
    otherwise                         = Double
    prod Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / Double -> Double
    forall a. Floating a => a -> a
    exp Double
    beta
          where
            prod :: Double
    prod =  Double
    xDouble -> Double -> Double
    forall a. Floating a => a -> a -> a
    **Double
    p Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    cxDouble -> Double -> Double
    forall a. Floating a => a -> a -> a
    **(Double
    q Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    1)
        -- Soper's expansion of incomplete beta function
        loop :: Double -> Int -> Double -> Double -> Double -> Double
    loop !Double
    psq (Int
    ns :: Int) Double
    ai Double
    term Double
    betain
          | Bool
    done      = Double
    betain' Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    factor Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / Double
    p
          | Bool
    otherwise = Double -> Int -> Double -> Double -> Double -> Double
    loop Double
    psq' (Int
    ns Int -> Int -> Int
    forall a. Num a => a -> a -> a
    - Int
    1) (Double
    ai Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    1) Double
    term' Double
    betain'
          where
            -- New values
            term' :: Double
    term'   = Double
    term Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    fact Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / (Double
    p Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    ai)
            betain' :: Double
    betain' = Double
    betain Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    term'
            fact :: Double
    fact | Int
    ns Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    >  Int
    0   = (Double
    q Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    ai) Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    xDouble -> Double -> Double
    forall a. Fractional a => a -> a -> a
    /Double
    cx
                 | Int
    ns Int -> Int -> Bool
    forall a. Eq a => a -> a -> Bool
    == Int
    0   = (Double
    q Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    ai) Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    x
                 | Bool
    otherwise = Double
    psq Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    x
            -- Iterations are complete
            done :: Bool
    done = Double
    db Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    <= Double
    eps Bool -> Bool -> Bool
    && Double
    db Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    <= Double
    epsDouble -> Double -> Double
    forall a. Num a => a -> a -> a
    *Double
    betain' where db :: Double
    db = Double -> Double
    forall a. Num a => a -> a
    abs Double
    term'
            psq' :: Double
    psq' = if Int
    ns Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    < Int
    0 then Double
    psq Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    1 else Double
    psq

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

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

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

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

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

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

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

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

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

    -- | Compute the factorial function /n/!.  Returns +∞ if the input is
    --   above 170 (above which the result cannot be represented by a
    --   64-bit 'Double').
    factorial :: Int -> Double
    factorial :: Int -> Double
    factorial Int
    n
      | Int
    n Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    < Int
    0     = String -> Double
    forall a. HasCallStack => String -> a
    error String
    "Numeric.SpecFunctions.factorial: negative input"
      | Int
    n Int -> Int -> Bool
    forall a. Ord a => a -> a -> Bool
    > Int
    170   = Double
    m_pos_inf
      | Bool
    otherwise = Vector Double -> Int -> Double
    forall a. Unbox a => Vector a -> Int -> a
    U.unsafeIndex Vector Double
    factorialTable Int
    n

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

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

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

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

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

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

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

    -- | Compute ψ(/x/), the first logarithmic derivative of the gamma
    --   function.
    --
    -- \[
    -- \psi(x) = \frac{d}{dx} \ln \left(\Gamma(x)\right) = \frac{\Gamma'(x)}{\Gamma(x)}
    -- \]
    --
    -- Uses Algorithm AS 103 by Bernardo, based on Minka's C implementation.
    digamma :: Double -> Double
    digamma :: Double -> Double
    digamma Double
    x
        | Double -> Bool
    forall a. RealFloat a => a -> Bool
    isNaN Double
    x Bool -> Bool -> Bool
    || Double -> Bool
    forall a. RealFloat a => a -> Bool
    isInfinite Double
    x                  = Double
    m_NaN
        -- FIXME:
        --   This is ugly. We are testing here that number is in fact
        --   integer. It's somewhat tricky question to answer. When ε for
        --   given number becomes 1 or greater every number is represents
        --   an integer. We also must make sure that excess precision
        --   won't bite us.
        | Double
    x Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    <= Double
    0 Bool -> Bool -> Bool
    && Int64 -> Double
    forall a b. (Integral a, Num b) => a -> b
    fromIntegral (Double -> Int64
    forall b. Integral b => Double -> b
    forall a b. (RealFrac a, Integral b) => a -> b
    truncate Double
    x :: Int64) Double -> Double -> Bool
    forall a. Eq a => a -> a -> Bool
    == Double
    x = Double
    m_neg_inf
        -- Jeffery's reflection formula
        | Double
    x Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    < Double
    0     = Double -> Double
    digamma (Double
    1 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    x) Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    forall a. Floating a => a
    pi Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / Double -> Double
    forall a. Floating a => a -> a
    tan (Double -> Double
    forall a. Num a => a -> a
    negate Double
    forall a. Floating a => a
    pi Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    x)
        | Double
    x Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    <= Double
    1e-6 = - Double
    γ Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    1Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    /Double
    x Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    trigamma1 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    x
        | Double
    x' Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    < Double
    c    = Double
    r
        -- De Moivre's expansion
        | Bool
    otherwise = let s :: Double
    s = Double
    1Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    /Double
    x'
                      in  Double -> [Double] -> Double
    forall a. Num a => a -> [a] -> a
    evaluateEvenPolynomialL Double
    s
                            [   Double
    r Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double -> Double
    forall a. Floating a => a -> a
    log Double
    x' Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    0.5 Double -> Double -> Double
    forall a. Num a => a -> a -> a
    * Double
    s
                            , - Double
    1Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    /Double
    12
                            ,   Double
    1Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    /Double
    120
                            , - Double
    1Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    /Double
    252
                            ,   Double
    1Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    /Double
    240
                            , - Double
    1Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    /Double
    132
                            ,  Double
    391Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    /Double
    32760
                            ]
      where
        γ :: Double
    γ  = Double
    m_eulerMascheroni
        c :: Double
    c  = Double
    12
        -- Reduce to digamma (x + n) where (x + n) >= c
        (Double
    r, Double
    x') = Double -> Double -> (Double, Double)
    reduce Double
    0 Double
    x
          where
            reduce :: Double -> Double -> (Double, Double)
    reduce !Double
    s Double
    y
              | Double
    y Double -> Double -> Bool
    forall a. Ord a => a -> a -> Bool
    < Double
    c     = Double -> Double -> (Double, Double)
    reduce (Double
    s Double -> Double -> Double
    forall a. Num a => a -> a -> a
    - Double
    1 Double -> Double -> Double
    forall a. Fractional a => a -> a -> a
    / Double
    y) (Double
    y Double -> Double -> Double
    forall a. Num a => a -> a -> a
    + Double
    1)
              | Bool
    otherwise = (Double
    s, Double
    y)

    ----------------------------------------------------------------
    -- Constants
    ----------------------------------------------------------------

    -- Coefficients for 18-point Gauss-Legendre integration. They are
    -- used in implementation of incomplete gamma and beta functions.
    coefW,coefY :: U.Vector Double
    coefW :: Vector Double
    coefW = [Double] -> Vector Double
    forall a. Unbox a => [a] -> Vector a
    U.fromList [ Double
    0.0055657196642445571, Double
    0.012915947284065419, Double
    0.020181515297735382
                       , Double
    0.027298621498568734,  Double
    0.034213810770299537, Double
    0.040875750923643261
                       , Double
    0.047235083490265582,  Double
    0.053244713977759692, Double
    0.058860144245324798
                       , Double
    0.064039797355015485,  Double
    0.068745323835736408, Double
    0.072941885005653087
                       , Double
    0.076598410645870640,  Double
    0.079687828912071670, Double
    0.082187266704339706
                       , Double
    0.084078218979661945,  Double
    0.085346685739338721, Double
    0.085983275670394821
                       ]
    coefY :: Vector Double
    coefY = [Double] -> Vector Double
    forall a. Unbox a => [a] -> Vector a
    U.fromList [ Double
    0.0021695375159141994, Double
    0.011413521097787704, Double
    0.027972308950302116
                       , Double
    0.051727015600492421,  Double
    0.082502225484340941, Double
    0.12007019910960293
                       , Double
    0.16415283300752470,   Double
    0.21442376986779355,  Double
    0.27051082840644336
                       , Double
    0.33199876341447887,   Double
    0.39843234186401943,  Double
    0.46931971407375483
                       , Double
    0.54413605556657973,   Double
    0.62232745288031077,  Double
    0.70331500465597174
                       , Double
    0.78649910768313447,   Double
    0.87126389619061517,  Double
    0.95698180152629142
                       ]
    {-# NOINLINE coefW #-}
    {-# NOINLINE coefY #-}

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

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

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

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