{-# LANGUAGE FlexibleInstances, TemplateHaskell #-}
module Math.Gamma
( Gamma(..)
, Factorial(..)
, IncGamma(..)
) where
import Math.Gamma.Lanczos
import Math.Gamma.Incomplete
import Math.Factorial
import Data.Complex
import Data.List (findIndex)
import GHC.Float (float2Double, double2Float)
import qualified Data.Vector.Unboxed as V
import Language.Haskell.TH (litE, Lit(IntegerL))
import Math.ContinuedFraction
import Math.Sequence.Converge
class (Eq a, Floating a, Factorial a) => Gamma a where
gamma :: a -> a
gamma a
0 = a
0a -> a -> a
forall a. Fractional a => a -> a -> a
/a
0
gamma a
z
| a
z a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a -> a
forall a. Num a => a -> a
abs a
z = a -> a
forall a. Floating a => a -> a
exp (a -> a
forall a. Gamma a => a -> a
lnGamma a
z)
| Bool
otherwise = a
forall a. Floating a => a
pi a -> a -> a
forall a. Fractional a => a -> a -> a
/ (a -> a
forall a. Floating a => a -> a
sin (a
forall a. Floating a => a
pi a -> a -> a
forall a. Num a => a -> a -> a
* a
z) a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
exp (a -> a
forall a. Gamma a => a -> a
lnGamma (a
1a -> a -> a
forall a. Num a => a -> a -> a
-a
z)))
lnGamma :: a -> a
lnGamma a
z = a -> a
forall a. Floating a => a -> a
log (a -> a
forall a. Gamma a => a -> a
gamma a
z)
lnFactorial :: Integral b => b -> a
lnFactorial b
n = a -> a
forall a. Gamma a => a -> a
lnGamma (b -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral b
na -> a -> a
forall a. Num a => a -> a -> a
+a
1)
floatGammaInfCutoff :: Double
floatGammaInfCutoff :: Double
floatGammaInfCutoff = $( do
let Just cutoff = findIndex isInfinite (scanl (*) (1::Float) [1..])
litE (IntegerL (1 + toInteger cutoff))
)
instance Gamma Float where
gamma :: Float -> Float
gamma = Double -> Float
double2Float (Double -> Float) -> (Float -> Double) -> Float -> Float
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Double -> Double
gam (Double -> Double) -> (Float -> Double) -> Float -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Float -> Double
float2Double
where
gam :: Double -> Double
gam Double
x
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
floatGammaInfCutoff = Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
0
| Bool
otherwise = case Double -> (Integer, Double)
forall a b. (RealFrac a, Integral b) => a -> (b, a)
properFraction Double
x of
(Integer
n,Double
0) | Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
1 -> Double
0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
0
| Bool
otherwise -> Integer -> Double
forall a b. (Factorial a, Integral b) => b -> a
factorial (Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1 :: Integer)
(Integer, Double)
_ | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< (-Double
20) -> let s :: Double
s = 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
sin (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x)
in Double -> Double
forall a. Num a => a -> a
signum Double
s 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 -> Double
forall a. Num a => a -> a
abs Double
s) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double
forall a. Gamma a => a -> a
lnGamma (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
x))
| Bool
otherwise -> (Double -> Double) -> Double -> Double
forall a. (RealFloat a, Ord a) => (a -> a) -> a -> a
reflect (Double -> [Double] -> Double -> Double
forall a. Floating a => a -> [a] -> a -> a
gammaLanczos Double
g [Double]
cs) Double
x
g :: Double
g = Double
forall a. Floating a => a
pi
cs :: [Double]
cs = [ Double
1.0000000249904433
, Double
9.100643759042066
,-Double
4.3325519094475
, Double
0.12502459858901147
, Double
1.1378929685052916e-4
,-Double
9.555011214455924e-5
]
lnGamma :: Float -> Float
lnGamma = Double -> Float
double2Float (Double -> Float) -> (Float -> Double) -> Float -> Float
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Double -> Double) -> Double -> Double
forall a. (RealFloat a, Ord a) => (a -> a) -> a -> a
reflectLn (Double -> [Double] -> Double -> Double
forall a. Floating a => a -> [a] -> a -> a
lnGammaLanczos Double
g [Double]
cs) (Double -> Double) -> (Float -> Double) -> Float -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Float -> Double
float2Double
where
g :: Double
g = Double
forall a. Floating a => a
pi
cs :: [Double]
cs = [ Double
1.0000000249904433
, Double
9.100643759042066
,-Double
4.3325519094475
, Double
0.12502459858901147
, Double
1.1378929685052916e-4
,-Double
9.555011214455924e-5
]
lnFactorial :: b -> Float
lnFactorial b
n
| Integer
n' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 = [Char] -> Float
forall a. HasCallStack => [Char] -> a
error [Char]
"lnFactorial n: n < 0"
| Integer
n' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
nFacs = Vector Float
facs Vector Float -> Int -> Float
forall a. Unbox a => Vector a -> Int -> a
V.! b -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral b
n
| Bool
otherwise = Float -> Float
forall a. Gamma a => a -> a
lnGamma (b -> Float
forall a b. (Integral a, Num b) => a -> b
fromIntegral b
nFloat -> Float -> Float
forall a. Num a => a -> a -> a
+Float
1)
where
n' :: Integer
n' = b -> Integer
forall a. Integral a => a -> Integer
toInteger b
n
nFacs :: Int
nFacs = Int
2000
facs :: Vector Float
facs = (Float -> Float) -> Vector Float -> Vector Float
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map Float -> Float
forall a. Gamma a => a -> a
lnGamma (Float -> Int -> Vector Float
forall a. (Unbox a, Num a) => a -> Int -> Vector a
V.enumFromN Float
1 Int
nFacs)
doubleGammaInfCutoff :: Double
doubleGammaInfCutoff :: Double
doubleGammaInfCutoff = $( do
let Just cutoff = findIndex isInfinite (scanl (*) (1::Double) [1..])
litE (IntegerL (1 + toInteger cutoff))
)
instance Gamma Double where
gamma :: Double -> Double
gamma Double
x
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
doubleGammaInfCutoff = Double
1Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
0
| Bool
otherwise = case Double -> (Integer, Double)
forall a b. (RealFrac a, Integral b) => a -> (b, a)
properFraction Double
x of
(Integer
n,Double
0) | Integer
n Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
1 -> Double
0Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/Double
0
| Bool
otherwise -> Integer -> Double
forall a b. (Factorial a, Integral b) => b -> a
factorial (Integer
nInteger -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1 :: Integer)
(Integer, Double)
_ | Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< (-Double
50) -> let s :: Double
s = 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
sin (Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double
x)
in Double -> Double
forall a. Num a => a -> a
signum Double
s 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 -> Double
forall a. Num a => a -> a
abs Double
s) Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> [Double] -> Double -> Double
forall a. Floating a => a -> [a] -> a -> a
lnGammaLanczos Double
g [Double]
cs (Double
1Double -> Double -> Double
forall a. Num a => a -> a -> a
-Double
x))
| Bool
otherwise -> (Double -> Double) -> Double -> Double
forall a. (RealFloat a, Ord a) => (a -> a) -> a -> a
reflect (Double -> [Double] -> Double -> Double
forall a. Floating a => a -> [a] -> a -> a
gammaLanczos Double
g [Double]
cs) Double
x
where
g :: Double
g = Double
2Double -> Double -> Double
forall a. Num a => a -> a -> a
*Double
forall a. Floating a => a
pi
cs :: [Double]
cs = [ Double
0.9999999999999858
, Double
311.6011750541472
,-Double
498.6511904603639
, Double
244.08472899976877
, -Double
38.670364643074194
, Double
1.3350900101370549
, -Double
1.8977221899565682e-3
, Double
8.475264614349149e-7
, Double
2.59715567376858e-7
, -Double
2.7166437850607517e-7
, Double
6.151114806136299e-8
]
lnGamma :: Double -> Double
lnGamma = (Double -> Double) -> Double -> Double
forall a. (RealFloat a, Ord a) => (a -> a) -> a -> a
reflectLn (Double -> [Double] -> Double -> Double
forall a. Floating a => a -> [a] -> a -> a
lnGammaLanczos Double
g [Double]
cs)
where
g :: Double
g = Double -> Double
forall a. Floating a => a -> a
exp Double
forall a. Floating a => a
pi Double -> Double -> Double
forall a. Fractional a => a -> a -> a
/ Double
forall a. Floating a => a
pi
cs :: [Double]
cs = [ Double
1.0000000000000002
, Double
1002.5049417114732
,-Double
1999.6140446432912
, Double
1352.1626218340114
, -Double
360.6486475548049
, Double
33.344988357090685
, -Double
0.6637188712004668
, Double
5.16644552377916e-4
, Double
1.684651140163429e-7
, -Double
1.8148207145896904e-7
, Double
6.171532716135051e-8
, -Double
9.014004881476154e-9
]
lnFactorial :: b -> Double
lnFactorial b
n
| Integer
n' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 = [Char] -> Double
forall a. HasCallStack => [Char] -> a
error [Char]
"lnFactorial n: n < 0"
| Integer
n' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
nFacs = Vector Double
facs Vector Double -> Int -> Double
forall a. Unbox a => Vector a -> Int -> a
V.! b -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral b
n
| Bool
otherwise = Double -> Double
forall a. Gamma a => a -> a
lnGamma (b -> Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral b
nDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1)
where
n' :: Integer
n' = b -> Integer
forall a. Integral a => a -> Integer
toInteger b
n
nFacs :: Int
nFacs = Int
2000
facs :: Vector Double
facs = (Double -> Double) -> Vector Double -> Vector Double
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map Double -> Double
forall a. Gamma a => a -> a
lnGamma (Double -> Int -> Vector Double
forall a. (Unbox a, Num a) => a -> Int -> Vector a
V.enumFromN Double
1 Int
nFacs)
complexDoubleToFloat :: Complex Double -> Complex Float
complexDoubleToFloat :: Complex Double -> Complex Float
complexDoubleToFloat (Double
a :+ Double
b) = Double -> Float
double2Float Double
a Float -> Float -> Complex Float
forall a. a -> a -> Complex a
:+ Double -> Float
double2Float Double
b
complexFloatToDouble :: Complex Float -> Complex Double
complexFloatToDouble :: Complex Float -> Complex Double
complexFloatToDouble (Float
a :+ Float
b) = Float -> Double
float2Double Float
a Double -> Double -> Complex Double
forall a. a -> a -> Complex a
:+ Float -> Double
float2Double Float
b
instance Gamma (Complex Float) where
gamma :: Complex Float -> Complex Float
gamma = Complex Double -> Complex Float
complexDoubleToFloat (Complex Double -> Complex Float)
-> (Complex Float -> Complex Double)
-> Complex Float
-> Complex Float
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Complex Double -> Complex Double
forall a. Gamma a => a -> a
gamma (Complex Double -> Complex Double)
-> (Complex Float -> Complex Double)
-> Complex Float
-> Complex Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Complex Float -> Complex Double
complexFloatToDouble
lnGamma :: Complex Float -> Complex Float
lnGamma = Complex Double -> Complex Float
complexDoubleToFloat (Complex Double -> Complex Float)
-> (Complex Float -> Complex Double)
-> Complex Float
-> Complex Float
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (Complex Double -> Complex Double)
-> Complex Double -> Complex Double
forall a.
RealFloat a =>
(Complex a -> Complex a) -> Complex a -> Complex a
reflectLnC (Complex Double
-> [Complex Double] -> Complex Double -> Complex Double
forall a. Floating a => a -> [a] -> a -> a
lnGammaLanczos Complex Double
g [Complex Double]
cs) (Complex Double -> Complex Double)
-> (Complex Float -> Complex Double)
-> Complex Float
-> Complex Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. Complex Float -> Complex Double
complexFloatToDouble
where
g :: Complex Double
g = Complex Double
forall a. Floating a => a
pi
cs :: [Complex Double]
cs = [ Complex Double
1.0000000249904433
, Complex Double
9.100643759042066
,-Complex Double
4.3325519094475
, Complex Double
0.12502459858901147
, Complex Double
1.1378929685052916e-4
,-Complex Double
9.555011214455924e-5
]
lnFactorial :: b -> Complex Float
lnFactorial b
n
| Integer
n' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 = [Char] -> Complex Float
forall a. HasCallStack => [Char] -> a
error [Char]
"lnFactorial n: n < 0"
| Integer
n' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
nFacs = Vector (Complex Float)
facs Vector (Complex Float) -> Int -> Complex Float
forall a. Unbox a => Vector a -> Int -> a
V.! b -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral b
n
| Bool
otherwise = Complex Float -> Complex Float
forall a. Gamma a => a -> a
lnGamma (b -> Complex Float
forall a b. (Integral a, Num b) => a -> b
fromIntegral b
nComplex Float -> Complex Float -> Complex Float
forall a. Num a => a -> a -> a
+Complex Float
1)
where
n' :: Integer
n' = b -> Integer
forall a. Integral a => a -> Integer
toInteger b
n
nFacs :: Int
nFacs = Int
2000
facs :: Vector (Complex Float)
facs = (Complex Float -> Complex Float)
-> Vector (Complex Float) -> Vector (Complex Float)
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map Complex Float -> Complex Float
forall a. Gamma a => a -> a
lnGamma (Complex Float -> Int -> Vector (Complex Float)
forall a. (Unbox a, Num a) => a -> Int -> Vector a
V.enumFromN Complex Float
1 Int
nFacs)
instance Gamma (Complex Double) where
gamma :: Complex Double -> Complex Double
gamma = (Complex Double -> Complex Double)
-> Complex Double -> Complex Double
forall a.
RealFloat a =>
(Complex a -> Complex a) -> Complex a -> Complex a
reflectC (Complex Double
-> [Complex Double] -> Complex Double -> Complex Double
forall a. Floating a => a -> [a] -> a -> a
gammaLanczos Complex Double
g [Complex Double]
cs)
where
g :: Complex Double
g = Complex Double
2Complex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
*Complex Double
forall a. Floating a => a
pi
cs :: [Complex Double]
cs = [ Complex Double
1.0000000000000002
, Complex Double
311.60117505414695
,-Complex Double
498.65119046033163
, Complex Double
244.08472899875767
, -Complex Double
38.67036462939322
, Complex Double
1.3350899103585203
, -Complex Double
1.8972831806242229e-3
, -Complex Double
3.935368195357295e-7
, Complex Double
2.592464641764731e-6
, -Complex Double
3.2263565156368265e-6
, Complex Double
2.5666169886566876e-6
, -Complex Double
1.3737776806198937e-6
, Complex Double
4.4551204024819644e-7
, -Complex Double
6.576826592057796e-8
]
lnGamma :: Complex Double -> Complex Double
lnGamma = (Complex Double -> Complex Double)
-> Complex Double -> Complex Double
forall a.
RealFloat a =>
(Complex a -> Complex a) -> Complex a -> Complex a
reflectLnC (Complex Double
-> [Complex Double] -> Complex Double -> Complex Double
forall a. Floating a => a -> [a] -> a -> a
lnGammaLanczos Complex Double
g [Complex Double]
cs)
where
g :: Complex Double
g = Complex Double -> Complex Double
forall a. Floating a => a -> a
exp Complex Double
forall a. Floating a => a
pi Complex Double -> Complex Double -> Complex Double
forall a. Fractional a => a -> a -> a
/ Complex Double
forall a. Floating a => a
pi
cs :: [Complex Double]
cs = [ Complex Double
1.0000000000000002
, Complex Double
1002.5049417114732
,-Complex Double
1999.6140446432912
, Complex Double
1352.1626218340114
, -Complex Double
360.6486475548049
, Complex Double
33.344988357090685
, -Complex Double
0.6637188712004668
, Complex Double
5.16644552377916e-4
, Complex Double
1.684651140163429e-7
, -Complex Double
1.8148207145896904e-7
, Complex Double
6.171532716135051e-8
, -Complex Double
9.014004881476154e-9
]
lnFactorial :: b -> Complex Double
lnFactorial b
n
| Integer
n' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Integer
0 = [Char] -> Complex Double
forall a. HasCallStack => [Char] -> a
error [Char]
"lnFactorial n: n < 0"
| Integer
n' Integer -> Integer -> Bool
forall a. Ord a => a -> a -> Bool
< Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
nFacs = Vector (Complex Double)
facs Vector (Complex Double) -> Int -> Complex Double
forall a. Unbox a => Vector a -> Int -> a
V.! b -> Int
forall a b. (Integral a, Num b) => a -> b
fromIntegral b
n
| Bool
otherwise = Complex Double -> Complex Double
forall a. Gamma a => a -> a
lnGamma (b -> Complex Double
forall a b. (Integral a, Num b) => a -> b
fromIntegral b
nComplex Double -> Complex Double -> Complex Double
forall a. Num a => a -> a -> a
+Complex Double
1)
where
n' :: Integer
n' = b -> Integer
forall a. Integral a => a -> Integer
toInteger b
n
nFacs :: Int
nFacs = Int
2000
facs :: Vector (Complex Double)
facs = (Complex Double -> Complex Double)
-> Vector (Complex Double) -> Vector (Complex Double)
forall a b. (Unbox a, Unbox b) => (a -> b) -> Vector a -> Vector b
V.map Complex Double -> Complex Double
forall a. Gamma a => a -> a
lnGamma (Complex Double -> Int -> Vector (Complex Double)
forall a. (Unbox a, Num a) => a -> Int -> Vector a
V.enumFromN Complex Double
1 Int
nFacs)
class Gamma a => IncGamma a where
lowerGamma :: a -> a -> a
lnLowerGamma :: a -> a -> a
p :: a -> a -> a
upperGamma :: a -> a -> a
lnUpperGamma :: a -> a -> a
q :: a -> a -> a
instance IncGamma Float where
lowerGamma :: Float -> Float -> Float
lowerGamma Float
s Float
x = Double -> Float
double2Float (Double -> Float) -> Double -> Float
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double
forall a. IncGamma a => a -> a -> a
lowerGamma :: Double -> Double -> Double) (Float -> Double
float2Double Float
s) (Float -> Double
float2Double Float
x)
lnLowerGamma :: Float -> Float -> Float
lnLowerGamma Float
s Float
x = Double -> Float
double2Float (Double -> Float) -> Double -> Float
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double
forall a. IncGamma a => a -> a -> a
lnLowerGamma :: Double -> Double -> Double) (Float -> Double
float2Double Float
s) (Float -> Double
float2Double Float
x)
p :: Float -> Float -> Float
p Float
s Float
x = Double -> Float
double2Float (Double -> Float) -> Double -> Float
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double
forall a. IncGamma a => a -> a -> a
p :: Double -> Double -> Double) (Float -> Double
float2Double Float
s) (Float -> Double
float2Double Float
x)
upperGamma :: Float -> Float -> Float
upperGamma Float
s Float
x = Double -> Float
double2Float (Double -> Float) -> Double -> Float
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double
forall a. IncGamma a => a -> a -> a
upperGamma :: Double -> Double -> Double) (Float -> Double
float2Double Float
s) (Float -> Double
float2Double Float
x)
lnUpperGamma :: Float -> Float -> Float
lnUpperGamma Float
s Float
x = Double -> Float
double2Float (Double -> Float) -> Double -> Float
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double
forall a. IncGamma a => a -> a -> a
lnUpperGamma :: Double -> Double -> Double) (Float -> Double
float2Double Float
s) (Float -> Double
float2Double Float
x)
q :: Float -> Float -> Float
q Float
s Float
x = Double -> Float
double2Float (Double -> Float) -> Double -> Float
forall a b. (a -> b) -> a -> b
$ (Double -> Double -> Double
forall a. IncGamma a => a -> a -> a
q :: Double -> Double -> Double) (Float -> Double
float2Double Float
s) (Float -> Double
float2Double Float
x)
instance IncGamma Double where
lowerGamma :: Double -> Double -> Double
lowerGamma Double
s Double
x
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 = [Char] -> Double
forall a. HasCallStack => [Char] -> a
error [Char]
"lowerGamma: x < 0 is not currently supported."
| Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 = Double
0
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
sDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1 = Double -> Double
forall a. Gamma a => a -> a
gamma Double
s Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double -> Double
forall a. IncGamma a => a -> a -> a
upperGamma Double
s Double
x
| Bool
otherwise = Double -> Double -> Double
forall b. (Eq b, Floating b) => b -> b -> b
lowerGammaHypGeom Double
s Double
x
upperGamma :: Double -> Double -> Double
upperGamma Double
s Double
x
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 = [Char] -> Double
forall a. HasCallStack => [Char] -> a
error [Char]
"upperGamma: x < 0 is not currently supported."
| Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 = Double -> Double
forall a. Gamma a => a -> a
gamma Double
s
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
sDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1 = Double -> Double -> Double
forall a. IncGamma a => a -> a -> a
q Double
s Double
x Double -> Double -> Double
forall a. Num a => a -> a -> a
* Double -> Double
forall a. Gamma a => a -> a
gamma Double
s
| Bool
otherwise = [Double] -> Double
forall a. Eq a => [a] -> a
converge ([Double] -> Double)
-> ([[Double]] -> [Double]) -> [[Double]] -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat
([[Double]] -> Double) -> [[Double]] -> Double
forall a b. (a -> b) -> a -> b
$ Double -> CF Double -> [[Double]]
forall a. (Fractional a, Eq a) => a -> CF a -> [[a]]
modifiedLentz Double
1e-30 (Double -> Double -> CF Double
forall a. (Floating a, Ord a) => a -> a -> CF a
upperGammaCF Double
s Double
x)
lnLowerGamma :: Double -> Double -> Double
lnLowerGamma Double
s Double
x
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 = [Char] -> Double
forall a. HasCallStack => [Char] -> a
error [Char]
"lnLowerGamma: x < 0 is not currently supported."
| Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 = Double -> Double
forall a. Floating a => a -> a
log Double
0
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
sDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1 = Double -> Double
forall a. Floating a => a -> a
log (Double -> Double -> Double
forall a. IncGamma a => a -> a -> a
p Double
s Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Gamma a => a -> a
lnGamma Double
s
| Bool
otherwise = Double -> Double -> Double
forall b. (Eq b, Floating b) => b -> b -> b
lnLowerGammaHypGeom Double
s Double
x
lnUpperGamma :: Double -> Double -> Double
lnUpperGamma Double
s Double
x
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 = [Char] -> Double
forall a. HasCallStack => [Char] -> a
error [Char]
"lnUpperGamma: x < 0 is not currently supported."
| Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 = Double -> Double
forall a. Gamma a => a -> a
lnGamma Double
s
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
sDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1 = Double -> Double
forall a. Floating a => a -> a
log (Double -> Double -> Double
forall a. IncGamma a => a -> a -> a
q Double
s Double
x) Double -> Double -> Double
forall a. Num a => a -> a -> a
+ Double -> Double
forall a. Gamma a => a -> a
lnGamma Double
s
| Bool
otherwise =
[Double] -> Double
forall a. Eq a => [a] -> a
converge (Double -> Double -> [Double]
forall a. (Eq a, Floating a) => a -> a -> [a]
lnUpperGammaConvergents Double
s Double
x)
p :: Double -> Double -> Double
p Double
s Double
x
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 = [Char] -> Double
forall a. HasCallStack => [Char] -> a
error [Char]
"p: x < 0 is not currently supported."
| Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 = Double
0
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
>= Double
sDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1 = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double -> Double
forall a. IncGamma a => a -> a -> a
q Double
s Double
x
| Bool
otherwise = Double -> Double -> Double
forall a. (Gamma a, Ord a) => a -> a -> a
pHypGeom Double
s Double
x
q :: Double -> Double -> Double
q Double
s Double
x
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
0 = [Char] -> Double
forall a. HasCallStack => [Char] -> a
error [Char]
"q: x < 0 is not currently supported."
| Double
x Double -> Double -> Bool
forall a. Eq a => a -> a -> Bool
== Double
0 = Double
1
| Double
x Double -> Double -> Bool
forall a. Ord a => a -> a -> Bool
< Double
sDouble -> Double -> Double
forall a. Num a => a -> a -> a
+Double
1 = Double
1 Double -> Double -> Double
forall a. Num a => a -> a -> a
- Double -> Double -> Double
forall a. IncGamma a => a -> a -> a
p Double
s Double
x
| Bool
otherwise =
[Double] -> Double
forall a. Eq a => [a] -> a
converge ([Double] -> Double)
-> ([[Double]] -> [Double]) -> [[Double]] -> Double
forall b c a. (b -> c) -> (a -> b) -> a -> c
. [[Double]] -> [Double]
forall (t :: * -> *) a. Foldable t => t [a] -> [a]
concat
([[Double]] -> Double) -> [[Double]] -> Double
forall a b. (a -> b) -> a -> b
$ Double -> CF Double -> [[Double]]
forall a. (Fractional a, Eq a) => a -> CF a -> [[a]]
modifiedLentz Double
1e-30 (Double -> Double -> CF Double
forall a. (Gamma a, Ord a, Enum a) => a -> a -> CF a
qCF Double
s Double
x)