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