{-# LANGUAGE BangPatterns, ScopedTypeVariables #-}
{-# OPTIONS_HADDOCK hide #-}
module Numeric.SpecFunctions.Internal
( module Numeric.SpecFunctions.Internal
, Compat.log1p
, Compat.expm1
) where
import Control.Applicative
import Data.Bits ((.&.), (.|.), shiftR)
import Data.Int (Int64)
import Data.Word (Word)
import Data.Default.Class
import qualified Data.Vector.Unboxed as U
import Data.Vector.Unboxed ((!))
import Text.Printf
import Numeric.Polynomial.Chebyshev (chebyshevBroucke)
import Numeric.Polynomial (evaluatePolynomial, evaluatePolynomialL, evaluateEvenPolynomialL
,evaluateOddPolynomialL)
import Numeric.RootFinding (Root(..), newtonRaphson, NewtonParam(..), Tolerance(..))
import Numeric.Series
import Numeric.MathFunctions.Constants
import Numeric.SpecFunctions.Compat (log1p)
import qualified Numeric.SpecFunctions.Compat as Compat
erf :: Double -> Double
erf :: Double -> Double
erf = Double -> Double
Compat.erf
{-# INLINE erf #-}
erfc :: Double -> Double
erfc :: Double -> Double
erfc = Double -> Double
Compat.erfc
{-# INLINE erfc #-}
invErf :: Double
-> Double
invErf :: Double -> Double
invErf p :: Double
p
| Double
p Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== 1 = Double
m_pos_inf
| Double
p Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== -1 = Double
m_neg_inf
| Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< 1 Bool -> Bool -> Bool
&& Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> -1 = if Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 0 then Double
r else -Double
r
| Bool
otherwise = [Char] -> Double
forall a. HasCallStack => [Char] -> a
error "invErf: p must in [-1,1] range"
where
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
$ 1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
pp
step :: Double -> Double
step x :: 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
invErfc :: Double
-> Double
invErfc :: Double -> Double
invErfc p :: Double
p
| Double
p Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== 2 = Double
m_neg_inf
| Double
p Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== 0 = Double
m_pos_inf
| Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>0 Bool -> Bool -> Bool
&& Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< 2 = if Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= 1 then Double
r else -Double
r
| Bool
otherwise = [Char] -> Double
forall a. [Char] -> a
modErr ([Char] -> Double) -> [Char] -> Double
forall a b. (a -> b) -> a -> b
$ "invErfc: p must be in [0,2] got " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Double -> [Char]
forall a. Show a => a -> [Char]
show Double
p
where
pp :: Double
pp | Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= 1 = Double
p
| Bool
otherwise = 2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 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
pp
step :: Double -> Double
step x :: 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
guessInvErfc :: Double -> Double
guessInvErfc :: Double -> Double
guessInvErfc p :: Double
p
= -0.70711 Double -> Double -> Double
forall a. Num a => a -> a -> a
* ((2.30753 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* 0.27061) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* (0.99229 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* 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
$ -2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log( 0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
p)
invErfcHalleyStep :: Double -> Double -> Double
invErfcHalleyStep :: Double -> Double -> Double
invErfcHalleyStep err :: Double
err x :: 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
/ (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)
data L = L {-# UNPACK #-} !Double {-# UNPACK #-} !Double
logGamma :: Double -> Double
logGamma :: Double -> Double
logGamma z :: Double
z
| Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= 0 = Double
m_pos_inf
| 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 (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)
| Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< 0.5 = Double -> Double -> Double
lgamma1_15 Double
z (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- 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
< 1 = Double -> Double -> Double
lgamma15_2 Double
z (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- 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
<= 1.5 = Double -> Double -> Double
lgamma1_15 (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- 1) (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- 2)
| Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< 2 = Double -> Double -> Double
lgamma15_2 (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- 1) (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- 2)
| Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< 15 = Double -> Double
lgammaSmall Double
z
| Bool
otherwise = Double -> Double
lanczosApprox Double
z
logGammaL :: Double -> Double
logGammaL :: Double -> Double
logGammaL = Double -> Double
logGamma
{-# DEPRECATED logGammaL "Use logGamma instead" #-}
lgamma1_15 :: Double -> Double -> Double
lgamma1_15 :: Double -> Double -> Double
lgamma1_15 zm1 :: Double
zm1 zm2 :: 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 = 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
[ 0.490622454069039543534e-1
, -0.969117530159521214579e-1
, -0.414983358359495381969e0
, -0.406567124211938417342e0
, -0.158413586390692192217e0
, -0.240149820648571559892e-1
, -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
[ 1
, 0.302349829846463038743e1
, 0.348739585360723852576e1
, 0.191415588274426679201e1
, 0.507137738614363510846e0
, 0.577039722690451849648e-1
, 0.195768102601107189171e-2
]
{-# NOINLINE tableLogGamma_1_15Q #-}
lgamma15_2 :: Double -> Double -> Double
lgamma15_2 :: Double -> Double -> Double
lgamma15_2 zm1 :: Double
zm1 zm2 :: 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 = 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
[ -0.292329721830270012337e-1
, 0.144216267757192309184e0
, -0.142440390738631274135e0
, 0.542809694055053558157e-1
, -0.850535976868336437746e-2
, 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
[ 1
, -0.150169356054485044494e1
, 0.846973248876495016101e0
, -0.220095151814995745555e0
, 0.25582797155975869989e-1
, -0.100666795539143372762e-2
, -0.827193521891290553639e-6
]
{-# NOINLINE tableLogGamma_15_2Q #-}
lgamma2_3 :: Double -> Double
lgamma2_3 :: Double -> Double
lgamma2_3 z :: 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
+ 1)
zm2 :: Double
zm2 = Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- 2
y :: Double
y = 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
[ -0.180355685678449379109e-1
, 0.25126649619989678683e-1
, 0.494103151567532234274e-1
, 0.172491608709613993966e-1
, -0.259453563205438108893e-3
, -0.541009869215204396339e-3
, -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
[ 1
, 0.196202987197795200688e1
, 0.148019669424231326694e1
, 0.541391432071720958364e0
, 0.988504251128010129477e-1
, 0.82130967464889339326e-2
, 0.224936291922115757597e-3
, -0.223352763208617092964e-6
]
{-# NOINLINE tableLogGamma_2_3Q #-}
lgammaSmall :: Double -> Double
lgammaSmall :: Double -> Double
lgammaSmall = Double -> Double -> Double
go 0
where
go :: Double -> Double -> Double
go acc :: Double
acc z :: Double
z | Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< 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
- 1
lanczosApprox :: Double -> Double
lanczosApprox :: Double -> Double
lanczosApprox z :: 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
- 0.5) Double -> Double -> Double
forall a. Num a => a -> a -> a
- 1) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
z Double -> Double -> Double
forall a. Num a => a -> a -> a
- 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 = 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
[ (56906521.91347156388090791033559122686859 , 0)
, (103794043.1163445451906271053616070238554 , 39916800)
, (86363131.28813859145546927288977868422342 , 120543840)
, (43338889.32467613834773723740590533316085 , 150917976)
, (14605578.08768506808414169982791359218571 , 105258076)
, (3481712.15498064590882071018964774556468 , 45995730)
, (601859.6171681098786670226533699352302507 , 13339535)
, (75999.29304014542649875303443598909137092 , 2637558)
, (6955.999602515376140356310115515198987526 , 357423)
, (449.9445569063168119446858607650988409623 , 32670)
, (19.51992788247617482847860966235652136208 , 1925)
, (0.5098416655656676188125178644804694509993 , 66)
, (0.006061842346248906525783753964555936883222 , 1)
]
evalRatio :: U.Vector (Double,Double) -> Double -> Double
evalRatio :: Vector (Double, Double) -> Double -> Double
evalRatio coef :: Vector (Double, Double)
coef x :: Double
x
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 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 0 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 0 0) Vector (Double, Double)
coef
where
fini :: L -> Double
fini (L num :: Double
num den :: Double
den) = Double
num Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
den
stepR :: (Double, Double) -> L -> L
stepR (a :: Double
a,b :: Double
b) (L num :: Double
num den :: 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 num :: Double
num den :: Double
den) (a :: Double
a,b :: 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
logGammaCorrection :: Double -> Double
logGammaCorrection :: Double -> Double
logGammaCorrection x :: Double
x
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< 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
* 2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 1) Vector Double
coeffs Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x
| Bool
otherwise = 1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* 12)
where
big :: Double
big = 94906265.62425156
t :: Double
t = 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 [
0.1666389480451863247205729650822e+0,
-0.1384948176067563840732986059135e-4,
0.9810825646924729426157171547487e-8,
-0.1809129475572494194263306266719e-10,
0.6221098041892605227126015543416e-13,
-0.3399615005417721944303330599666e-15,
0.2683181998482698748957538846666e-17
]
incompleteGamma :: Double
-> Double
-> Double
incompleteGamma :: Double -> Double -> Double
incompleteGamma a :: Double
a x :: Double
x
| Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= 0 Bool -> Bool -> Bool
|| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< 0 = [Char] -> Double
forall a. HasCallStack => [Char] -> a
error
([Char] -> Double) -> [Char] -> Double
forall a b. (a -> b) -> a -> b
$ "incompleteGamma: Domain error z=" [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Double -> [Char]
forall a. Show a => a -> [Char]
show Double
a [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ " x=" [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Double -> [Char]
forall a. Show a => a -> [Char]
show Double
x
| Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== 0 = 0
| Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
m_pos_inf = 1
| 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
> 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
* (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
+ 1))
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< 0.5 = case () of
_| (-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
< 1.1 = case () of
_| 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
> 20 Bool -> Bool -> Bool
&& Bool
useTemme = Double
uniformExpansion
| Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- (1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (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
> 200 Bool -> Bool -> Bool
&& 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
< 0.4)
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
(/) 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
+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
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
+1))
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
(/) 1 (Double -> Sequence Double
forall a. Num a => a -> Sequence a
enumSequenceFrom 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)
contFraction :: Double
contFraction = 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 = (\k :: 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
+ 2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 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 0
uniformExpansion :: Double
uniformExpansion =
let
fm :: U.Vector Double
fm :: Vector Double
fm = [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
U.fromList [ 1.00000000000000000000e+00
,-3.33333333333333370341e-01
, 8.33333333333333287074e-02
,-1.48148148148148153802e-02
, 1.15740740740740734316e-03
, 3.52733686067019369930e-04
,-1.78755144032921825352e-04
, 3.91926317852243766954e-05
,-2.18544851067999240532e-06
,-1.85406221071515996597e-06
, 8.29671134095308545622e-07
,-1.76659527368260808474e-07
, 6.70785354340149841119e-09
, 1.02618097842403069078e-08
,-4.38203601845335376897e-09
, 9.14769958223679020897e-10
,-2.55141939949462514346e-11
,-5.83077213255042560744e-11
, 2.43619480206674150369e-11
,-5.02766928011417632057e-12
, 1.10043920319561347525e-13
, 3.37176326240098513631e-13
]
y :: Double
y = - Double -> Double
log1pmx Double
mu
eta :: Double
eta = Double -> Double
forall a. Floating a => a -> a
sqrt (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
loop :: Double -> Double -> Double -> Int -> Double
loop !Double
_ !Double
_ u :: Double
u 0 = Double
u
loop bm1 :: Double
bm1 bm0 :: Double
bm0 u :: Double
u i :: 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
+ 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
-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
-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
-2)) 0 (Int
nInt -> Int -> Int
forall a. Num a => a -> a -> a
-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 1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/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
/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 (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
invIncompleteGamma :: Double
-> Double
-> Double
invIncompleteGamma :: Double -> Double -> Double
invIncompleteGamma a :: Double
a p :: Double
p
| Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= 0 =
[Char] -> Double
forall a. [Char] -> a
modErr ([Char] -> Double) -> [Char] -> Double
forall a b. (a -> b) -> a -> b
$ [Char] -> Double -> Double -> [Char]
forall r. PrintfType r => [Char] -> r
printf "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
< 0 Bool -> Bool -> Bool
|| Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 1 =
[Char] -> Double
forall a. [Char] -> a
modErr ([Char] -> Double) -> [Char] -> Double
forall a b. (a -> b) -> a -> b
$ [Char] -> Double -> Double -> [Char]
forall r. PrintfType r => [Char] -> r
printf "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
== 0 = 0
| Double
p Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== 1 = 1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ 0
| Bool
otherwise = Int -> Double -> Double
loop 0 Double
guess
where
loop :: Int -> Double -> Double
loop :: Int -> Double -> Double
loop i :: Int
i x :: Double
x
| Int
i Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
>= 12 = Double
x'
| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
f' = 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
+ 1) Double
x'
where
f :: Double
f = Double -> Double -> Double
incompleteGamma Double
a Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p
f' :: Double
f' | Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 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'
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
- 1)
dx :: Double
dx = Double
u Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double -> Double
forall a. Ord a => a -> a -> a
min 1.0 Double
corr)
x' :: Double
x' | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
dx = 0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x
| Bool
otherwise = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
dx
guess :: Double
guess
| Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 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
$ -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
< 0.5 then Double
p else 1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
p)
x1 :: Double
x1 = (2.30753 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* 0.27061) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* (0.99229 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
t Double -> Double -> Double
forall a. Num a => a -> a -> a
* 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
< 0.5 then -Double
x1 else Double
x1
in Double -> Double -> Double
forall a. Ord a => a -> a -> a
max 1e-3 (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* (1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/(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
/ (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
** 3)
| Bool
otherwise =
let t :: Double
t = 1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
* (0.253 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
*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
** (1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
a)
else 1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Floating a => a -> a
log( 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
/ (1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
t))
a1 :: Double
a1 = Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- 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
- 1) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
gln )
gln :: Double
gln = Double -> Double
logGamma Double
a
eps :: Double
eps = 1e-8
logBeta
:: Double
-> Double
-> Double
logBeta :: Double -> Double -> Double
logBeta a :: Double
a b :: Double
b
| Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< 0 = Double
m_NaN
| Double
p Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== 0 = Double
m_pos_inf
| Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= 10 = Double
allStirling
| Double
q Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= 10 = Double
twoStirling
| 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
allStirling :: Double
allStirling
= Double -> Double
forall a. Floating a => a -> a
log Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
* (-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
- 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)
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
- 0.5) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log1p(-Double
ppq)
incompleteBeta :: Double
-> Double
-> Double
-> Double
incompleteBeta :: Double -> Double -> Double -> Double
incompleteBeta p :: Double
p q :: Double
q = Double -> Double -> Double -> Double -> Double
incompleteBeta_ (Double -> Double -> Double
logBeta Double
p Double
q) Double
p Double
q
incompleteBeta_ :: Double
-> Double
-> Double
-> Double
-> Double
incompleteBeta_ :: Double -> Double -> Double -> Double -> Double
incompleteBeta_ beta :: Double
beta p :: Double
p q :: Double
q x :: Double
x
| Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= 0 Bool -> Bool -> Bool
|| Double
q Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= 0 =
[Char] -> Double
forall a. [Char] -> a
modErr ([Char] -> Double) -> [Char] -> Double
forall a b. (a -> b) -> a -> b
$ [Char] -> Double -> Double -> Double -> [Char]
forall r. PrintfType r => [Char] -> r
printf "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
< 0 Bool -> Bool -> Bool
|| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 1 Bool -> Bool -> Bool
|| Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
x =
[Char] -> Double
forall a. [Char] -> a
modErr ([Char] -> Double) -> [Char] -> Double
forall a b. (a -> b) -> a -> b
$ [Char] -> Double -> Double -> Double -> [Char]
forall r. PrintfType r => [Char] -> r
printf "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
== 0 Bool -> Bool -> Bool
|| Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== 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 = 1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double -> Double -> Double -> Double
incompleteBetaWorker Double
beta Double
q Double
p (1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x)
incompleteBetaApprox :: Double -> Double -> Double -> Double -> Double
incompleteBetaApprox :: Double -> Double -> Double -> Double -> Double
incompleteBetaApprox beta :: Double
beta p :: Double
p q :: Double
q x :: Double
x
| Double
ans Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 0 = 1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
ans
| Bool
otherwise = -Double
ans
where
p1 :: Double
p1 = Double
p Double -> Double -> Double
forall a. Num a => a -> a -> a
- 1
q1 :: Double
q1 = Double
q Double -> Double -> Double
forall a. Num a => a -> a -> a
- 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)
xu :: Double
xu = Double -> Double -> Double
forall a. Ord a => a -> a -> a
max 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
- 10Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
t) (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- 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
+ 1) )
go :: Double -> Double -> Double
go y :: Double
y w :: 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(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 )
incompleteBetaWorker :: Double -> Double -> Double -> Double -> Double
incompleteBetaWorker :: Double -> Double -> Double -> Double -> Double
incompleteBetaWorker beta :: Double
beta p :: Double
p q :: Double
q x :: Double
x
| Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 3000 Bool -> Bool -> Bool
&& Double
q Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 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 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)) 1 1 1
where
eps :: Double
eps = 1e-15
cx :: Double
cx = 1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x
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
- 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
- 1)
loop :: Double -> Int -> Double -> Double -> Double -> Double
loop !Double
psq (Int
ns :: Int) ai :: Double
ai term :: Double
term betain :: 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
- 1) (Double
ai Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 1) Double
term' Double
betain'
where
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
> 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
== 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
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
< 0 then Double
psq Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 1 else Double
psq
invIncompleteBeta :: Double
-> Double
-> Double
-> Double
invIncompleteBeta :: Double -> Double -> Double -> Double
invIncompleteBeta p :: Double
p q :: Double
q a :: Double
a
| Double
p Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= 0 Bool -> Bool -> Bool
|| Double
q Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= 0 =
[Char] -> Double
forall a. [Char] -> a
modErr ([Char] -> Double) -> [Char] -> Double
forall a b. (a -> b) -> a -> b
$ [Char] -> Double -> Double -> Double -> [Char]
forall r. PrintfType r => [Char] -> r
printf "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
< 0 Bool -> Bool -> Bool
|| Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 1 =
[Char] -> Double
forall a. [Char] -> a
modErr ([Char] -> Double) -> [Char] -> Double
forall a b. (a -> b) -> a -> b
$ [Char] -> Double -> Double -> Double -> [Char]
forall r. PrintfType r => [Char] -> r
printf "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
== 0 Bool -> Bool -> Bool
|| Double
a Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== 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 beta :: Double
beta a :: Double
a b :: Double
b p :: Double
p = Int -> Double -> Double
forall t. (Ord t, Num t) => t -> Double -> Double
loop (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
- 1
b1 :: Double
b1 = Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- 1
loop :: t -> Double -> Double
loop !t
i !Double
x
| Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== 0 Bool -> Bool -> Bool
|| Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== 1 = Double
x
| Double -> Bool
forall a. RealFloat a => a -> Bool
isInfinite Double
f' = Double
x
| t
i t -> t -> Bool
forall a. Ord a => a -> a -> Bool
>= 10 = Double
x
| Double -> Double
forall a. Num a => a -> a
abs Double
dx Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= 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
+1) Double
x'
where
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'
corr :: Double
corr | Double
d Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 1 = 1
| Double
d Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< -1 = -1
| Double -> Bool
forall a. RealFloat a => a -> Bool
isNaN Double
d = 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
/ (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
/ (1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
corr)
x' :: Double
x' | Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< 0 = Double
x Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ 2
| Double
z Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 1 = (Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ 2
| Bool
otherwise = Double
z
where z :: Double
z = Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
dx
invIncBetaGuess :: Double -> Double -> Double -> Double -> Double
invIncBetaGuess :: Double -> Double -> Double -> Double -> Double
invIncBetaGuess beta :: Double
beta a :: Double
a b :: Double
b p :: Double
p
| Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< 1 Bool -> Bool -> Bool
&& Double
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< 1 =
let x_infl :: Double
x_infl = (1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
a) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (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
** (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
/ (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
* (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
** (1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
b) in 1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
xgDouble -> Double -> Double
forall a. Fractional a => a -> a -> a
/(1Double -> Double -> Double
forall a. Num a => a -> a -> a
+Double
xg)
in Double
x
| Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= 6 Bool -> Bool -> Bool
&& Double
aDouble -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>1 Bool -> Bool -> Bool
&& Double
bDouble -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>1 =
let x_infl :: Double
x_infl = (Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- 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
- 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 = 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((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
| Double
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< 5 Bool -> Bool -> Bool
&& Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= 1 =
let x :: Double
x | Double
pDouble -> Double -> Double
forall a. Floating a => a -> a -> a
**(1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
a) Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< 0.5 = (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 -> Double
forall a. Floating a => a -> a
exp Double
beta) Double -> Double -> Double
forall a. Floating a => a -> a -> a
** (1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
a)
| Bool
otherwise = 1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- (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
**(1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
b)
in Double
x
| Double
aDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 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
> 4 =
let
eta0 :: Double
eta0 = Double -> Double -> Double
invIncompleteGamma Double
b (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
w :: Double
w = Double -> Double
forall a. Floating a => a -> a
sqrt(1 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double
mu)
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
+ 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
e1 :: Double
e1 = (Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 2) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- 1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (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
+ 9 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 21 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 5) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d
Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (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
- 13 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 69 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 167 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 46) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d_2
Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (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
- (7 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 21 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 70 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 26 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 93 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- 31) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d_3
Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (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
- (75 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_6 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 202 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 188 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 888 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 1345 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 118 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 138) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d_4
Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (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 = (28 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 131 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 402 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 581 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 208) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- 1)
Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (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
- (35 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_6 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 154 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 623 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 1636 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 3983 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 3514 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- 925) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d
Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (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
- ( 2132 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_7 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 7915 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_6 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 16821 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 35066 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 87490 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3
Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 141183 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 95993 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 21640
) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d_2
Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (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
- ( 11053 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_8 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 53308 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_7 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 117010 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_6 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 163924 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 116188 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4
Double -> Double -> Double
forall a. Num a => a -> a -> a
- 258428 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 677042 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 481940 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- 105497
) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d_3
Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (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 = -( (3592 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_7 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 8375 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_6 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 1323 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 29198 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 89578 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3
Double -> Double -> Double
forall a. Num a => a -> a -> a
- 154413 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 116063 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- 29632
) Double -> Double -> Double
forall a. Num a => a -> a -> a
* (Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- 1)
)
Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (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
- ( 442043 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_9 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 2054169 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_8 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 3803094 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_7 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 3470754 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_6 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 2141568 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5
Double -> Double -> Double
forall a. Num a => a -> a -> a
- 2393568 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 19904934 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 34714674 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 23128299 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
- 5253353
) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d
Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (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
- ( 116932 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_10 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 819281 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_9 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 2378172 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_8 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 4341330 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_7 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 6806004 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_6
Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 10622748 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_5 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 18739500 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_4 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 30651894 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_3 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 30869976 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w_2
Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 15431867 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 2919016
) Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
d_2
Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (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 (1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
a) [Double
eta0, Double
e1, Double
e2, Double
e3]
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
+ (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(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 = 1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (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 0
upper :: Double
upper = if Double
eta Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
mu then 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
/ 2
func :: Double -> (Double, Double)
func x :: 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(1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x)
, 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
/(1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
x)
)
Root x0 :: Double
x0 = NewtonParam
-> (Double, Double, Double)
-> (Double -> (Double, Double))
-> Root Double
newtonRaphson NewtonParam
forall a. Default a => a
def{newtonTol :: Tolerance
newtonTol=Double -> Tolerance
RelTol 1e-8} (Double
lower, Double
x_guess, Double
upper) Double -> (Double, Double)
func
in Double
x0
| Double
a Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 1 Bool -> Bool -> Bool
&& Double
b Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 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
- 3) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ 6
s :: Double
s = 1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
- 1)
t :: Double
t = 1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- 1)
h :: Double
h = 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
+ 5Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/6 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (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(2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
w))
| Double
chi2 Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 0 Bool -> Bool -> Bool
&& Double
ratio Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 1 = 1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- 2 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (Double
ratio Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 1)
| 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
** (1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
a)
| Bool
otherwise -> 1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* (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
** (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
ratio :: Double
ratio = (4Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
a Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
- 2) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
chi2
chi2 :: Double
chi2 = 2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b Double -> Double -> Double
forall a. Num a => a -> a -> a
* (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
** 3
where
t :: Double
t = 1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ (9 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
b)
y :: Double
y = Double
r Double -> Double -> Double
forall a. Num a => a -> a -> a
- ( 2.30753 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 0.27061 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
r )
Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ ( 1.0 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ ( 0.99229 Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 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
$ - 2 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Floating a => a -> a
log Double
p
sinc :: Double -> Double
sinc :: Double -> Double
sinc x :: Double
x
| Double
ax Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
eps_0 = 1
| Double
ax Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
eps_2 = 1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x2Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/6
| Double
ax Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
eps_4 = 1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double
x2Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/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
/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
eps_0 :: Double
eps_0 = 1.8250120749944284e-8
eps_2 :: Double
eps_2 = 1.4284346431400855e-4
eps_4 :: Double
eps_4 = 4.043633626430947e-3
log1pmx :: Double -> Double
log1pmx :: Double -> Double
log1pmx x :: Double
x
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< -1 = [Char] -> Double
forall a. HasCallStack => [Char] -> a
error "Domain error"
| Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== -1 = Double
m_neg_inf
| Double
ax Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
> 0.95 = Double -> Double
forall a. Floating a => a -> a
log(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
/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 2)
where
ax :: Double
ax = Double -> Double
forall a. Num a => a -> a
abs Double
x
log2 :: Int -> Int
log2 :: Int -> Int
log2 v0 :: Int
v0
| Int
v0 Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
<= 0 = [Char] -> Int
forall a. [Char] -> a
modErr ([Char] -> Int) -> [Char] -> Int
forall a b. (a -> b) -> a -> b
$ "log2: nonpositive input, got " [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ Int -> [Char]
forall a. Show a => a -> [Char]
show Int
v0
| Bool
otherwise = Int -> Int -> Int -> Int
go 5 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
== -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
/= 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
-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
-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 [ 0x02, 0x0c, 0xf0, 0xff00
, Word -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral (0xffff0000 :: Word)
, Word -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral (0xffffffff00000000 :: Word)]
!sv :: Vector Int
sv = [Int] -> Vector Int
forall a. Unbox a => [a] -> Vector a
U.fromList [1,2,4,8,16,32]
factorial :: Int -> Double
factorial :: Int -> Double
factorial n :: Int
n
| Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< 0 = [Char] -> Double
forall a. HasCallStack => [Char] -> a
error "Numeric.SpecFunctions.factorial: negative input"
| Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> 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
logFactorial :: Integral a => a -> Double
logFactorial :: a -> Double
logFactorial n :: a
n
| a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< 0 = [Char] -> Double
forall a. HasCallStack => [Char] -> a
error "Numeric.SpecFunctions.logFactorial: negative input"
| a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<= 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)
| a
n a -> a -> Bool
forall a. Ord a => a -> a -> Bool
< 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
* ((1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/12) Double -> Double -> Double
forall a. Num a => a -> a -> a
- (1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/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
+ (1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/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
- 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
+ 1
rx :: Double
rx = 1 Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
x
{-# SPECIALIZE logFactorial :: Int -> Double #-}
stirlingError :: Double -> Double
stirlingError :: Double -> Double
stirlingError n :: Double
n
| Double
n Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= 15.0 = case Double -> (Int, 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
(i :: Int
i,0) -> Vector Double
sfe Vector Double -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
`U.unsafeIndex` Int
i
_ -> Double -> Double
logGamma (Double
nDouble -> Double -> Double
forall a. Num a => a -> a -> a
+1.0) Double -> Double -> Double
forall a. Num a => a -> a -> a
- (Double
nDouble -> Double -> Double
forall a. Num a => a -> a -> a
+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
> 500 = Double -> [Double] -> Double
forall a. Num a => a -> [a] -> a
evaluateOddPolynomialL (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
> 80 = Double -> [Double] -> Double
forall a. Num a => a -> [a] -> a
evaluateOddPolynomialL (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
> 35 = Double -> [Double] -> Double
forall a. Num a => a -> [a] -> a
evaluateOddPolynomialL (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 (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 = 0.083333333333333333333
s1 :: Double
s1 = 0.00277777777777777777778
s2 :: Double
s2 = 0.00079365079365079365079365
s3 :: Double
s3 = 0.000595238095238095238095238
s4 :: Double
s4 = 0.0008417508417508417508417508
sfe :: Vector Double
sfe = [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
U.fromList [ 0.0,
0.1534264097200273452913848, 0.0810614667953272582196702,
0.0548141210519176538961390, 0.0413406959554092940938221,
0.03316287351993628748511048, 0.02767792568499833914878929,
0.02374616365629749597132920, 0.02079067210376509311152277,
0.01848845053267318523077934, 0.01664469118982119216319487,
0.01513497322191737887351255, 0.01387612882307074799874573,
0.01281046524292022692424986, 0.01189670994589177009505572,
0.01110455975820691732662991, 0.010411265261972096497478567,
0.009799416126158803298389475, 0.009255462182712732917728637,
0.008768700134139385462952823, 0.008330563433362871256469318,
0.007934114564314020547248100, 0.007573675487951840794972024,
0.007244554301320383179543912, 0.006942840107209529865664152,
0.006665247032707682442354394, 0.006408994188004207068439631,
0.006171712263039457647532867, 0.005951370112758847735624416,
0.005746216513010115682023589, 0.005554733551962801371038690 ]
logChooseFast :: Double -> Double -> Double
logChooseFast :: Double -> Double -> Double
logChooseFast n :: Double
n k :: Double
k = -Double -> Double
forall a. Floating a => a -> a
log (Double
n Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 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
+ 1) (Double
k Double -> Double -> Double
forall a. Num a => a -> a -> a
+ 1)
chooseExact :: Int -> Int -> Double
n :: Int
n chooseExact :: Int -> Int -> Double
`chooseExact` k :: 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 a. Integral a => Double -> a -> Double
go 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 1 Int
k
where
go :: Double -> a -> Double
go a :: Double
a i :: a
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 = a -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral a
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)
logChoose :: Int -> Int -> Double
n :: Int
n logChoose :: Int -> Int -> Double
`logChoose` k :: Int
k
| Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n = (-1) Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ 0
| Int
k' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< 50 Bool -> Bool -> Bool
&& (Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< 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)
choose :: Int -> Int -> Double
n :: Int
n choose :: Int -> Int -> Double
`choose` k :: Int
k
| Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n = 0
| Int
k' Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< 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 x :: a
x = a -> Int64
forall a b. (RealFrac a, Integral b) => a -> b
round a
x :: Int64
digamma :: Double -> Double
digamma :: Double -> Double
digamma x :: 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
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
<= 0 Bool -> Bool -> Bool
&& Int64 -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral (Double -> Int64
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
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< 0 = Double -> Double
digamma (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
<= 1e-6 = - Double
γ Double -> Double -> Double
forall a. Num a => a -> a -> a
- 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
| Bool
otherwise = let s :: Double
s = 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
- 0.5 Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
s
, - 1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/12
, 1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/120
, - 1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/252
, 1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/240
, - 1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/132
, 391Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/32760
]
where
γ :: Double
γ = Double
m_eulerMascheroni
c :: Double
c = 12
(r :: Double
r, x' :: Double
x') = Double -> Double -> (Double, Double)
reduce 0 Double
x
where
reduce :: Double -> Double -> (Double, Double)
reduce !Double
s y :: 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
- 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
+ 1)
| Bool
otherwise = (Double
s, Double
y)
coefW,coefY :: U.Vector Double
coefW :: Vector Double
coefW = [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
U.fromList [ 0.0055657196642445571, 0.012915947284065419, 0.020181515297735382
, 0.027298621498568734, 0.034213810770299537, 0.040875750923643261
, 0.047235083490265582, 0.053244713977759692, 0.058860144245324798
, 0.064039797355015485, 0.068745323835736408, 0.072941885005653087
, 0.076598410645870640, 0.079687828912071670, 0.082187266704339706
, 0.084078218979661945, 0.085346685739338721, 0.085983275670394821
]
coefY :: Vector Double
coefY = [Double] -> Vector Double
forall a. Unbox a => [a] -> Vector a
U.fromList [ 0.0021695375159141994, 0.011413521097787704, 0.027972308950302116
, 0.051727015600492421, 0.082502225484340941, 0.12007019910960293
, 0.16415283300752470, 0.21442376986779355, 0.27051082840644336
, 0.33199876341447887, 0.39843234186401943, 0.46931971407375483
, 0.54413605556657973, 0.62232745288031077, 0.70331500465597174
, 0.78649910768313447, 0.87126389619061517, 0.95698180152629142
]
{-# NOINLINE coefW #-}
{-# NOINLINE coefY #-}
trigamma1 :: Double
trigamma1 :: Double
trigamma1 = 1.6449340668482264365
modErr :: String -> a
modErr :: [Char] -> a
modErr msg :: [Char]
msg = [Char] -> a
forall a. HasCallStack => [Char] -> a
error ([Char] -> a) -> [Char] -> a
forall a b. (a -> b) -> a -> b
$ "Numeric.SpecFunctions." [Char] -> [Char] -> [Char]
forall a. [a] -> [a] -> [a]
++ [Char]
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 171
[ 1.0
, 1.0
, 2.0
, 6.0
, 24.0
, 120.0
, 720.0
, 5040.0
, 40320.0
, 362880.0
, 3628800.0
, 3.99168e7
, 4.790016e8
, 6.2270208e9
, 8.71782912e10
, 1.307674368e12
, 2.0922789888e13
, 3.55687428096e14
, 6.402373705728e15
, 1.21645100408832e17
, 2.43290200817664e18
, 5.109094217170944e19
, 1.1240007277776077e21
, 2.5852016738884974e22
, 6.204484017332394e23
, 1.5511210043330984e25
, 4.032914611266056e26
, 1.0888869450418352e28
, 3.0488834461171384e29
, 8.841761993739702e30
, 2.6525285981219103e32
, 8.222838654177922e33
, 2.631308369336935e35
, 8.683317618811886e36
, 2.9523279903960412e38
, 1.0333147966386144e40
, 3.719933267899012e41
, 1.3763753091226343e43
, 5.23022617466601e44
, 2.0397882081197442e46
, 8.159152832478977e47
, 3.3452526613163803e49
, 1.4050061177528798e51
, 6.041526306337383e52
, 2.6582715747884485e54
, 1.1962222086548019e56
, 5.5026221598120885e57
, 2.5862324151116818e59
, 1.2413915592536073e61
, 6.082818640342675e62
, 3.0414093201713376e64
, 1.5511187532873822e66
, 8.065817517094388e67
, 4.2748832840600255e69
, 2.308436973392414e71
, 1.2696403353658275e73
, 7.109985878048634e74
, 4.0526919504877214e76
, 2.3505613312828785e78
, 1.386831185456898e80
, 8.32098711274139e81
, 5.075802138772247e83
, 3.146997326038793e85
, 1.9826083154044399e87
, 1.2688693218588415e89
, 8.24765059208247e90
, 5.44344939077443e92
, 3.647111091818868e94
, 2.4800355424368305e96
, 1.711224524281413e98
, 1.197857166996989e100
, 8.504785885678623e101
, 6.1234458376886085e103
, 4.470115461512684e105
, 3.307885441519386e107
, 2.4809140811395396e109
, 1.88549470166605e111
, 1.4518309202828586e113
, 1.1324281178206297e115
, 8.946182130782974e116
, 7.15694570462638e118
, 5.797126020747368e120
, 4.753643337012841e122
, 3.9455239697206583e124
, 3.314240134565353e126
, 2.81710411438055e128
, 2.422709538367273e130
, 2.1077572983795275e132
, 1.8548264225739844e134
, 1.650795516090846e136
, 1.4857159644817613e138
, 1.352001527678403e140
, 1.2438414054641305e142
, 1.1567725070816416e144
, 1.087366156656743e146
, 1.0329978488239058e148
, 9.916779348709496e149
, 9.619275968248211e151
, 9.426890448883246e153
, 9.332621544394413e155
, 9.332621544394415e157
, 9.425947759838358e159
, 9.614466715035125e161
, 9.902900716486179e163
, 1.0299016745145626e166
, 1.0813967582402908e168
, 1.1462805637347082e170
, 1.2265202031961378e172
, 1.3246418194518288e174
, 1.4438595832024934e176
, 1.5882455415227428e178
, 1.7629525510902446e180
, 1.974506857221074e182
, 2.2311927486598134e184
, 2.543559733472187e186
, 2.9250936934930154e188
, 3.393108684451898e190
, 3.9699371608087206e192
, 4.68452584975429e194
, 5.574585761207606e196
, 6.689502913449126e198
, 8.094298525273443e200
, 9.875044200833601e202
, 1.214630436702533e205
, 1.5061417415111406e207
, 1.8826771768889257e209
, 2.372173242880047e211
, 3.0126600184576594e213
, 3.856204823625804e215
, 4.974504222477286e217
, 6.466855489220473e219
, 8.471580690878819e221
, 1.1182486511960041e224
, 1.4872707060906857e226
, 1.9929427461615188e228
, 2.6904727073180504e230
, 3.6590428819525483e232
, 5.012888748274991e234
, 6.917786472619488e236
, 9.615723196941088e238
, 1.3462012475717523e241
, 1.898143759076171e243
, 2.6953641378881624e245
, 3.8543707171800725e247
, 5.5502938327393044e249
, 8.047926057471992e251
, 1.1749972043909107e254
, 1.7272458904546386e256
, 2.5563239178728654e258
, 3.808922637630569e260
, 5.713383956445854e262
, 8.62720977423324e264
, 1.3113358856834524e267
, 2.0063439050956823e269
, 3.0897696138473508e271
, 4.789142901463393e273
, 7.471062926282894e275
, 1.1729568794264143e278
, 1.8532718694937346e280
, 2.946702272495038e282
, 4.714723635992061e284
, 7.590705053947218e286
, 1.2296942187394494e289
, 2.0044015765453023e291
, 3.287218585534296e293
, 5.423910666131589e295
, 9.003691705778436e297
, 1.5036165148649988e300
, 2.526075744973198e302
, 4.269068009004705e304
, 7.257415615307998e306
]