{-# LANGUAGE ParallelListComp #-}
-- |Lanczos' approximation to the gamma function, as described at
-- http:\/\/en.wikipedia.org\/wiki\/Lanczos_approximation
-- (fetched 11 June 2010).
-- 
-- Constants to be supplied by user.  There is a file \"extras/LanczosConstants.hs\"
-- in the source repository that implements a technique by Paul Godfrey for
-- calculating the coefficients.  It is not included in the distribution yet 
-- because it makes use of a linear algebra library I have not yet released 
-- (though I eventually intend to).
module Math.Gamma.Lanczos
    ( gammaLanczos, lnGammaLanczos
    , reflect, reflectC
    , reflectLn, reflectLnC
    ) where

import Data.Complex

-- |Compute Lanczos' approximation to the gamma function, using the specified
-- constants.  Valid for Re(x) > 0.5.  Use 'reflect' or 'reflectC' to extend
-- to the whole real line or complex plane, respectively.
{-# INLINE gammaLanczos #-}
gammaLanczos :: Floating a => a -> [a] -> a -> a
gammaLanczos :: a -> [a] -> a -> a
gammaLanczos a
_ [] a
_ = [Char] -> a
forall a. HasCallStack => [Char] -> a
error [Char]
"gammaLanczos: empty coefficient list"
gammaLanczos a
g [a]
cs a
zp1
    = a -> a
forall a. Floating a => a -> a
sqrt (a
2a -> a -> a
forall a. Num a => a -> a -> a
*a
forall a. Floating a => a
pi) a -> a -> a
forall a. Num a => a -> a -> a
* a
x a -> a -> a
forall a. Floating a => a -> a -> a
** (a
zp1 a -> a -> a
forall a. Num a => a -> a -> a
- a
0.5) a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Floating a => a -> a
exp (a -> a
forall a. Num a => a -> a
negate a
x) a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> a -> a
forall t. Fractional t => [t] -> t -> t
a [a]
cs a
z
    where
        x :: a
x = a
zp1 a -> a -> a
forall a. Num a => a -> a -> a
+ (a
g a -> a -> a
forall a. Num a => a -> a -> a
- a
0.5)
        z :: a
z = a
zp1 a -> a -> a
forall a. Num a => a -> a -> a
- a
1

-- |Compute Lanczos' approximation to the natural logarithm of the gamma
-- function, using the specified constants.  Valid for Re(x) > 0.5.  Use
-- 'reflectLn' or 'reflectLnC' to extend to the whole real line or complex
-- plane, respectively.
{-# INLINE lnGammaLanczos #-}
lnGammaLanczos :: Floating a => a -> [a] -> a -> a
lnGammaLanczos :: a -> [a] -> a -> a
lnGammaLanczos a
_ [] a
_ = [Char] -> a
forall a. HasCallStack => [Char] -> a
error [Char]
"lnGammaLanczos: empty coefficient list"
lnGammaLanczos a
g [a]
cs a
zp1 
    = a -> a
forall a. Floating a => a -> a
log (a -> a
forall a. Floating a => a -> a
sqrt (a
2a -> a -> a
forall a. Num a => a -> a -> a
*a
forall a. Floating a => a
pi)) a -> a -> a
forall a. Num a => a -> a -> a
+ a -> a
forall a. Floating a => a -> a
log a
x a -> a -> a
forall a. Num a => a -> a -> a
* (a
zp1 a -> a -> a
forall a. Num a => a -> a -> a
- a
0.5) a -> a -> a
forall a. Num a => a -> a -> a
- a
x a -> a -> a
forall a. Num a => a -> a -> a
+ a -> a
forall a. Floating a => a -> a
log ([a] -> a -> a
forall t. Fractional t => [t] -> t -> t
a [a]
cs a
z)
    where 
        x :: a
x = a
zp1 a -> a -> a
forall a. Num a => a -> a -> a
+ (a
g a -> a -> a
forall a. Num a => a -> a -> a
- a
0.5)
        z :: a
z = a
zp1 a -> a -> a
forall a. Num a => a -> a -> a
- a
1

{-# INLINE a #-}
a :: Fractional t => [t] -> t -> t
a :: [t] -> t -> t
a [] t
_ = [Char] -> t
forall a. HasCallStack => [Char] -> a
error [Char]
"Math.Gamma.Lanczos.a: empty coefficient list"
a [t]
cs t
z = [t] -> t
forall a. [a] -> a
head [t]
cs t -> t -> t
forall a. Num a => a -> a -> a
+ [t] -> t
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [t
c t -> t -> t
forall a. Fractional a => a -> a -> a
/ (t
z t -> t -> t
forall a. Num a => a -> a -> a
+ t
k) | t
c <- [t] -> [t]
forall a. [a] -> [a]
tail [t]
cs | t
k <- (t -> t) -> t -> [t]
forall a. (a -> a) -> a -> [a]
iterate (t
1t -> t -> t
forall a. Num a => a -> a -> a
+) t
1]

fractionalPart :: RealFloat a => a -> a
fractionalPart :: a -> a
fractionalPart a
x = case a -> (Int, a)
forall a b. (RealFrac a, Integral b) => a -> (b, a)
properFraction a
x of
    (Int
i,a
f) -> let Int
_ = Int
i :: Int in a
f

-- |Extend an approximation of the gamma function from the domain x > 0.5 to
-- the whole real line.
{-# INLINE reflect #-}
reflect :: (RealFloat a, Ord a) => (a -> a) -> a -> a
reflect :: (a -> a) -> a -> a
reflect a -> a
gamma a
z
    | a
z a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0.5               = a -> a
gamma a
z
    | a -> a
forall a. RealFloat a => a -> a
fractionalPart a
z a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = a
0
    | 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
gamma (a
1a -> a -> a
forall a. Num a => a -> a -> a
-a
z))

-- |Extend an approximation of the gamma function from the domain Re(x) > 0.5
-- to the whole complex plane.
{-# INLINE reflectC #-}
reflectC :: RealFloat a => (Complex a -> Complex a) -> Complex a -> Complex a
reflectC :: (Complex a -> Complex a) -> Complex a -> Complex a
reflectC Complex a -> Complex a
gamma Complex a
z
    | Complex a -> a
forall a. Complex a -> a
realPart Complex a
z a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0.5  = Complex a -> Complex a
gamma Complex a
z
    | Complex a -> a
forall a. Complex a -> a
imagPart Complex a
z a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
    Bool -> Bool -> Bool
&& a -> a
forall a. RealFloat a => a -> a
fractionalPart (Complex a -> a
forall a. Complex a -> a
realPart Complex a
z) a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
                        = Complex a
0Complex a -> Complex a -> Complex a
forall a. Fractional a => a -> a -> a
/Complex a
0
    | Bool
otherwise         = Complex a
forall a. Floating a => a
pi Complex a -> Complex a -> Complex a
forall a. Fractional a => a -> a -> a
/ (Complex a -> Complex a
forall a. Floating a => a -> a
sin (Complex a
forall a. Floating a => a
pi Complex a -> Complex a -> Complex a
forall a. Num a => a -> a -> a
* Complex a
z) Complex a -> Complex a -> Complex a
forall a. Num a => a -> a -> a
* Complex a -> Complex a
gamma (Complex a
1Complex a -> Complex a -> Complex a
forall a. Num a => a -> a -> a
-Complex a
z))

-- |Extend an approximation of the natural logarithm of the gamma function 
-- from the domain x > 0.5 to the whole real line.
{-# INLINE reflectLn #-}
reflectLn :: (RealFloat a, Ord a) => (a -> a) -> a -> a
reflectLn :: (a -> a) -> a -> a
reflectLn a -> a
lnGamma a
z
    | a
z a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0.5               = a -> a
lnGamma a
z
    | a -> a
forall a. RealFloat a => a -> a
fractionalPart a
z a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0 = a -> a
forall a. Floating a => a -> a
log (a
0a -> a -> a
forall a. Fractional a => a -> a -> a
/a
0)
    | Bool
otherwise             = a -> a
forall a. Floating a => a -> a
log a
forall a. Floating a => a
pi a -> a -> a
forall a. Num a => a -> a -> a
- a -> a
forall a. Floating a => a -> a
log (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
lnGamma (a
1a -> a -> a
forall a. Num a => a -> a -> a
-a
z)

-- |Extend an approximation of the natural logarithm of the gamma function 
-- from the domain Re(x) > 0.5 to the whole complex plane.
{-# INLINE reflectLnC #-}
reflectLnC :: RealFloat a => (Complex a -> Complex a) -> Complex a -> Complex a
reflectLnC :: (Complex a -> Complex a) -> Complex a -> Complex a
reflectLnC Complex a -> Complex a
lnGamma Complex a
z
    | Complex a -> a
forall a. Complex a -> a
realPart Complex a
z a -> a -> Bool
forall a. Ord a => a -> a -> Bool
> a
0.5  = Complex a -> Complex a
lnGamma Complex a
z
    | Complex a -> a
forall a. Complex a -> a
imagPart Complex a
z a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
    Bool -> Bool -> Bool
&& a -> a
forall a. RealFloat a => a -> a
fractionalPart (Complex a -> a
forall a. Complex a -> a
realPart Complex a
z) a -> a -> Bool
forall a. Eq a => a -> a -> Bool
== a
0
                        = Complex a -> Complex a
forall a. Floating a => a -> a
log (Complex a
0Complex a -> Complex a -> Complex a
forall a. Fractional a => a -> a -> a
/Complex a
0)
    | Bool
otherwise = Complex a -> Complex a
forall a. Floating a => a -> a
log Complex a
forall a. Floating a => a
pi Complex a -> Complex a -> Complex a
forall a. Num a => a -> a -> a
- Complex a -> Complex a
forall a. Floating a => a -> a
log (Complex a -> Complex a
forall a. Floating a => a -> a
sin (Complex a
forall a. Floating a => a
pi Complex a -> Complex a -> Complex a
forall a. Num a => a -> a -> a
* Complex a
z)) Complex a -> Complex a -> Complex a
forall a. Num a => a -> a -> a
- Complex a -> Complex a
lnGamma (Complex a
1Complex a -> Complex a -> Complex a
forall a. Num a => a -> a -> a
-Complex a
z)