{-# LANGUAGE ParallelListComp #-}
module Math.Gamma.Lanczos
( gammaLanczos, lnGammaLanczos
, reflect, reflectC
, reflectLn, reflectLnC
) where
import Data.Complex
{-# 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
{-# 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
{-# 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))
{-# 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))
{-# 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)
{-# 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)