module Math.Gamma.Stirling (lnGammaStirling, cs, s, abs_s, terms) where
import qualified Data.Vector as V
lnGammaStirling :: Floating a => [a] -> a -> a
lnGammaStirling cs z = (z 0.5) * log z z + 0.5 * log (2*pi) + sum [c / q | c <- cs | q <- risingPowers (z+1)]
where
risingPowers x = scanl1 (*) (iterate (1+) x)
cs :: (Fractional a, Ord a) => [a]
cs = map c [1..]
c :: (Fractional a, Ord a) => Int -> a
c n = 0.5 * recip n' * sum [k' * fromInteger (abs_s n k) / ((k' + 1) * (k' + 2)) | k <- [1..n], let k' = fromIntegral k]
where n' = fromIntegral n
s :: Int -> Int -> Integer
s n k
| n < 0 = error "s n k: n < 0"
| k < 0 = error "s n k: k < 0"
| k > n = error "s n k: k > n"
| otherwise = s n k
where
table = [V.generate (n+1) $ \k -> s n k | n <- [0..]]
s 0 0 = 1
s _ 0 = 0
s n k
| n == k = 1
| otherwise = s (n1) (k1) (toInteger n1) * s (n1) k
where
s n k = table !! n V.! k
abs_s :: Int -> Int -> Integer
abs_s n k
| n < 0 = error "abs_s n k: n < 0"
| k < 0 = error "abs_s n k: k < 0"
| k > n = error "abs_s n k: k > n"
| otherwise = abs_s n k
where
table = [V.generate (n+1) $ \k -> abs_s n k | n <- [0..]]
abs_s 0 0 = 1
abs_s _ 0 = 0
abs_s n k
| n == k = 1
| otherwise = abs_s (n1) (k1) + (toInteger n1) * abs_s (n1) k
where
abs_s n k = table !! n V.! k
terms :: (Num t, Floating a, Ord a) => a -> a -> t
terms prec z = converge (eps z) (f z)
where
cs' = cs
f z = scanl1 (+) [c / q | c <- cs' | q <- risingPowers (z+1)]
eps z = prec * abs ((z 0.5) * log z z + 0.5 * log (2*pi))
converge eps xs = go 1 xs where go n (x:y:zs) | abs(xy)<=eps = n | otherwise = go (n+1) (y:zs)
f z = scanl1 (+) [c / q | c <- cs | q <- risingPowers (z+1)]