{-# 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

-- |Gamma function.  Minimal definition is ether 'gamma' or 'lnGamma'.
class (Eq a, Floating a, Factorial a) => Gamma a where
    -- |The gamma function:  gamma z == integral from 0 to infinity of
    -- @\t -> t**(z-1) * exp (negate t)@
    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)))


    -- |Natural log of the gamma function
    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)
    
    -- |Natural log of the factorial function
    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 -- limited only by time and space
            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 -- limited only by time and space
            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 -- limited only by time and space
            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 -- limited only by time and space
            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)


-- |Incomplete gamma functions.
class Gamma a => IncGamma a where
    -- |Lower gamma function: lowerGamma s x == integral from 0 to x of 
    -- @\t -> t**(s-1) * exp (negate t)@
    lowerGamma :: a -> a -> a
    -- |Natural log of lower gamma function
    lnLowerGamma :: a -> a -> a 
    -- |Regularized lower incomplete gamma function: lowerGamma s x / gamma s
    p :: a -> a -> a
    
    -- |Upper gamma function: lowerGamma s x == integral from x to infinity of 
    -- @\t -> t**(s-1) * exp (negate t)@
    upperGamma :: a -> a -> a
    -- |Natural log of upper gamma function
    lnUpperGamma :: a -> a -> a
    -- |Regularized upper incomplete gamma function: upperGamma s x / gamma s
    q :: a -> a -> a

-- |This instance uses the Double instance.
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)

-- |I have not yet come up with a good strategy for evaluating these 
-- functions for negative @x@.  They can be rather numerically unstable.
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)