{-# LANGUAGE ParallelListComp #-}
module Math.Gamma.Stirling (lnGammaStirling, cs, s, abs_s, terms) where
import qualified Data.Vector as V
{-# INLINE lnGammaStirling #-}
lnGammaStirling :: Floating a => [a] -> a -> a
lnGammaStirling :: [a] -> a -> a
lnGammaStirling [a]
ks a
z
= (a
z 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
log a
z
a -> a -> a
forall a. Num a => a -> a -> a
- a
z
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
log (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 (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Fractional a => a -> a -> a
(/) [a]
ks (a -> [a]
forall a. Num a => a -> [a]
risingPowers (a
za -> a -> a
forall a. Num a => a -> a -> a
+a
1)))
{-# INLINE risingPowers #-}
risingPowers :: Num a => a -> [a]
risingPowers :: a -> [a]
risingPowers = (a -> a -> a) -> [a] -> [a]
forall a. (a -> a -> a) -> [a] -> [a]
scanl1 a -> a -> a
forall a. Num a => a -> a -> a
(*) ([a] -> [a]) -> (a -> [a]) -> a -> [a]
forall b c a. (b -> c) -> (a -> b) -> a -> c
. (a -> a) -> a -> [a]
forall a. (a -> a) -> a -> [a]
iterate (a
1a -> a -> a
forall a. Num a => a -> a -> a
+)
cs :: (Fractional a, Ord a) => [a]
cs :: [a]
cs = (Int -> a) -> [Int] -> [a]
forall a b. (a -> b) -> [a] -> [b]
map Int -> a
forall a. (Fractional a, Ord a) => Int -> a
c [Int
1..]
c :: (Fractional a, Ord a) => Int -> a
c :: Int -> a
c Int
n = a
0.5 a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Fractional a => a -> a
recip a
n' a -> a -> a
forall a. Num a => a -> a -> a
* [a] -> a
forall (t :: * -> *) a. (Foldable t, Num a) => t a -> a
sum [a
k' a -> a -> a
forall a. Num a => a -> a -> a
* Integer -> a
forall a. Num a => Integer -> a
fromInteger (Int -> Int -> Integer
abs_s Int
n Int
k) a -> a -> a
forall a. Fractional a => a -> a -> a
/ ((a
k' a -> a -> a
forall a. Num a => a -> a -> a
+ a
1) a -> a -> a
forall a. Num a => a -> a -> a
* (a
k' a -> a -> a
forall a. Num a => a -> a -> a
+ a
2)) | Int
k <- [Int
1..Int
n], let k' :: a
k' = Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
k]
where n' :: a
n' = Int -> a
forall a b. (Integral a, Num b) => a -> b
fromIntegral Int
n
s :: Int -> Int -> Integer
s :: Int -> Int -> Integer
s Int
n Int
k
| Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = [Char] -> Integer
forall a. HasCallStack => [Char] -> a
error [Char]
"s n k: n < 0"
| Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = [Char] -> Integer
forall a. HasCallStack => [Char] -> a
error [Char]
"s n k: k < 0"
| Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n = [Char] -> Integer
forall a. HasCallStack => [Char] -> a
error [Char]
"s n k: k > n"
| Bool
otherwise = Int -> Int -> Integer
s' Int
n Int
k
where
table :: [Vector Integer]
table = [Int -> (Int -> Integer) -> Vector Integer
forall a. Int -> (Int -> a) -> Vector a
V.generate (Int
iInt -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) ((Int -> Integer) -> Vector Integer)
-> (Int -> Integer) -> Vector Integer
forall a b. (a -> b) -> a -> b
$ Int -> Int -> Integer
s' Int
i | Int
i <- [Int
0..]]
s' :: Int -> Int -> Integer
s' Int
0 Int
0 = Integer
1
s' Int
_ Int
0 = Integer
0
s' Int
n' Int
k'
| Int
n' Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
k' = Integer
1
| Bool
otherwise = Int -> Int -> Integer
s'' (Int
n'Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) (Int
k'Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
- (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n'Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Int -> Integer
s'' (Int
n'Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Int
k'
where
s'' :: Int -> Int -> Integer
s'' Int
n'' Int
k'' = [Vector Integer]
table [Vector Integer] -> Int -> Vector Integer
forall a. [a] -> Int -> a
!! Int
n'' Vector Integer -> Int -> Integer
forall a. Vector a -> Int -> a
V.! Int
k''
abs_s :: Int -> Int -> Integer
abs_s :: Int -> Int -> Integer
abs_s Int
n Int
k
| Int
n Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = [Char] -> Integer
forall a. HasCallStack => [Char] -> a
error [Char]
"abs_s n k: n < 0"
| Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
< Int
0 = [Char] -> Integer
forall a. HasCallStack => [Char] -> a
error [Char]
"abs_s n k: k < 0"
| Int
k Int -> Int -> Bool
forall a. Ord a => a -> a -> Bool
> Int
n = [Char] -> Integer
forall a. HasCallStack => [Char] -> a
error [Char]
"abs_s n k: k > n"
| Bool
otherwise = Int -> Int -> Integer
abs_s' Int
n Int
k
where
table :: [Vector Integer]
table = [Int -> (Int -> Integer) -> Vector Integer
forall a. Int -> (Int -> a) -> Vector a
V.generate (Int
n''Int -> Int -> Int
forall a. Num a => a -> a -> a
+Int
1) ((Int -> Integer) -> Vector Integer)
-> (Int -> Integer) -> Vector Integer
forall a b. (a -> b) -> a -> b
$ \Int
k'' -> Int -> Int -> Integer
abs_s' Int
n'' Int
k'' | Int
n'' <- [Int
0..]]
abs_s' :: Int -> Int -> Integer
abs_s' Int
0 Int
0 = Integer
1
abs_s' Int
_ Int
0 = Integer
0
abs_s' Int
n' Int
k'
| Int
n' Int -> Int -> Bool
forall a. Eq a => a -> a -> Bool
== Int
k' = Integer
1
| Bool
otherwise = Int -> Int -> Integer
abs_s'' (Int
n'Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) (Int
k'Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
+ (Int -> Integer
forall a. Integral a => a -> Integer
toInteger Int
n'Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
-Integer
1) Integer -> Integer -> Integer
forall a. Num a => a -> a -> a
* Int -> Int -> Integer
abs_s'' (Int
n'Int -> Int -> Int
forall a. Num a => a -> a -> a
-Int
1) Int
k'
where
abs_s'' :: Int -> Int -> Integer
abs_s'' Int
n'' Int
k'' = [Vector Integer]
table [Vector Integer] -> Int -> Vector Integer
forall a. [a] -> Int -> a
!! Int
n'' Vector Integer -> Int -> Integer
forall a. Vector a -> Int -> a
V.! Int
k''
terms :: (Num t, Floating a, Ord a) => a -> a -> t
terms :: a -> a -> t
terms a
prec a
z = a -> [a] -> t
forall t a. (Num t, Ord a, Num a) => a -> [a] -> t
stepsNeeded (a -> a
eps a
z) (a -> [a]
f a
z)
where
cs' :: [a]
cs' = [a]
forall a. (Fractional a, Ord a) => [a]
cs
f :: a -> [a]
f a
x = (a -> a -> a) -> [a] -> [a]
forall a. (a -> a -> a) -> [a] -> [a]
scanl1 a -> a -> a
forall a. Num a => a -> a -> a
(+) ((a -> a -> a) -> [a] -> [a] -> [a]
forall a b c. (a -> b -> c) -> [a] -> [b] -> [c]
zipWith a -> a -> a
forall a. Fractional a => a -> a -> a
(/) [a]
cs' (a -> [a]
forall a. Num a => a -> [a]
risingPowers (a
xa -> a -> a
forall a. Num a => a -> a -> a
+a
1)))
eps :: a -> a
eps a
x = a
prec a -> a -> a
forall a. Num a => a -> a -> a
* a -> a
forall a. Num a => a -> a
abs ((a
x 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
log a
x a -> a -> a
forall a. Num a => a -> a -> a
- a
x 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
log (a
2a -> a -> a
forall a. Num a => a -> a -> a
*a
forall a. Floating a => a
pi))
stepsNeeded :: a -> [a] -> t
stepsNeeded a
e [a]
xs = t -> [a] -> t
forall t. Num t => t -> [a] -> t
go t
1 [a]
xs
where
go :: t -> [a] -> t
go t
n (a
x:a
y:[a]
zs)
| a -> a
forall a. Num a => a -> a
abs(a
xa -> a -> a
forall a. Num a => a -> a -> a
-a
y)a -> a -> Bool
forall a. Ord a => a -> a -> Bool
<=a
e = t
n
| Bool
otherwise = t -> [a] -> t
go (t
nt -> t -> t
forall a. Num a => a -> a -> a
+t
1) (a
ya -> [a] -> [a]
forall a. a -> [a] -> [a]
:[a]
zs)
go t
n [a]
_ = t
n