gamma (empty) → 0.7
raw patch · 8 files changed
+806/−0 lines, 8 filesdep +basedep +continued-fractionsdep +convergesetup-changed
Dependencies added: base, continued-fractions, converge, template-haskell, vector
Files
- Setup.lhs +5/−0
- extras/LanczosConstants.hs +104/−0
- gamma.cabal +38/−0
- src/Math/Gamma.hs +324/−0
- src/Math/Gamma.hs-boot +22/−0
- src/Math/Gamma/Incomplete.hs +143/−0
- src/Math/Gamma/Lanczos.hs +89/−0
- src/Math/Gamma/Stirling.hs +81/−0
+ Setup.lhs view
@@ -0,0 +1,5 @@+#!/usr/bin/env runhaskell++> import Distribution.Simple+> main = defaultMain+
+ extras/LanczosConstants.hs view
@@ -0,0 +1,104 @@+-- This makes use of a not-yet-released matrix library. It could be rewritten+-- to use any of the existing ones on hackage, but I don't know of any of them+-- that support matrices over arbitrary types - they are all focused on+-- efficiently packing the matrices and/or calling foreign libraries+-- (BLAS/GSL/etc.) and do not support any types other than Double, Float, and+-- Complex Double/Float.+-- +-- I am keeping it around anyway and including this file in the source +-- distribution, because with a very small amount of work an end-user +-- could fill in the gaps and use this code to generate their own constants +-- for lanczos gamma function approximations, which one may wish to do if +-- they wanted to implement, say, a gamma function for a very high precision+-- floating point type.+--+-- Note that these really need to be run with significantly higher precision+-- than the target type or truncation error will make the results useless.+-- +-- The algorithm implemented here is by Paul Godfrey, and is described in full+-- at http://www.numericana.com/answer/info/godfrey.htm (as of 21 June 2010).+module LanczosConstants where++import Math.Matrix+import Math.Matrix.Alias++cs g n = vectorToList (applyRat dbc f)+ where+ applyRat :: (Real t, Fractional t) => IMatrix Rational -> IVector t -> IVector t+ applyRat m v = fromRatVec (apply m (toRatVec v))+ fromRatVec :: (Vector v t, Fractional t) => IVector Rational -> v t+ fromRatVec = convertByV fromRational+ toRatVec :: (Vector v t, Real t) => v t -> IVector Rational+ toRatVec = convertByV toRational+ + dbc = dbcMat n+ f = fVec g n++dbcMat n = multRat d (multRat b c)+ where+ multRat :: (Real a, Matrix m1 a, Real b, Matrix m2 b) => m1 a -> m2 b -> IMatrix Rational+ multRat = multiplyWith sum (\d b -> toRational d * toRational b)+ + d = dMat n+ b = bMat n+ c = cMat n++fVec :: (Floating b, Vector v b) => b -> Int -> v b+fVec g n = vector n f+ where+ f a = sqrt (2 / pi)+ * product [fromIntegral i - 0.5 | i <-[1..a]]+ * exp (a' + g + 0.5)+ / (a' + g + 0.5) ** (a' + 0.5)+ where a' = fromIntegral a++cMat :: Int -> IMatrix Rational+cMat n = matrix n n m+ where+ m 0 0 = 1/2+ m i j = fromInteger (c (2*i) (2*j))+ + c 0 0 = 1+ c 1 1 = 1+ c i 0 = negate (c (i-2) 0)+ c i j+ | i == j = 2 * c (i-1) (j-1)+ | i > j = 2 * c (i-1) (j-1) - c (i-2) j+ | otherwise = 0++dMat :: Int -> IAlias Mat Integer+dMat n = AsDiag (IVec (ivector n dFunc)) 0+ where+ dFunc 0 = 1+ dFunc (i+1) = negate (factorial (2*i+2) `div` (2 * factorial i * factorial (i+1)))+ factorial n = product [1..toInteger n]++bMat :: Int -> IMatrix Integer+bMat n = matrixFromList bList+ where+ bList = take n . map (take n) $+ repeat 1 : + [ replicate i 0 ++ bicofs (negate (toInteger i*2))+ | i <- [1..]+ ]+ + + bFunc 0 _ = 1+ bFunc i j+ | i > j = 0+ bFunc i j = bicofs (toInteger (2 * j - 1)) !! i+ + bicofs x = go x 1 1+ where+ go num denom x = x : go (num+signum num) (denom+signum denom) (x * num `div` denom)++-- +-- p g k = sum [c (2*k+1) (2*a+1) * f a | a <- [0..k]]+-- where+-- k' = fromIntegral k+-- f a = +{-# INLINE risingPowers #-}+risingPowers x = scanl1 (*) (iterate (1+) x)++{-# INLINE fallingPowers #-}+fallingPowers x = scanl1 (*) (iterate (subtract 1) x)
+ gamma.cabal view
@@ -0,0 +1,38 @@+name: gamma+version: 0.7+stability: provisional++cabal-version: >= 1.6+build-type: Simple++author: James Cook <mokus@deepbondi.net>+maintainer: James Cook <mokus@deepbondi.net>+license: PublicDomain+homepage: /dev/null++category: Math, Numerical+synopsis: Gamma function and related functions.+description: Approximations of the gamma function, incomplete gamma + functions, beta function, and factorials.++tested-with: GHC == 6.10.4,+ GHC == 6.12.1, GHC == 6.12.3++extra-source-files: extras/*.hs++source-repository head+ type: darcs+ location: http://code.haskell.org/~mokus/gamma+++Library+ hs-source-dirs: src+ exposed-modules: Math.Gamma+ Math.Gamma.Incomplete+ Math.Gamma.Stirling+ Math.Gamma.Lanczos+ build-depends: base >= 3 && <5,+ continued-fractions >= 0.9.1,+ converge,+ template-haskell,+ vector >= 0.5 && < 0.7
+ src/Math/Gamma.hs view
@@ -0,0 +1,324 @@+{-# LANGUAGE FlexibleInstances, TemplateHaskell #-}+module Math.Gamma+ ( Gamma(..)+ , Factorial(..)+ , IncGamma(..)+ , beta+ ) where++import Math.Gamma.Lanczos+import Math.Gamma.Incomplete++import Data.Complex+import Data.List (sortBy, findIndex)+import Data.Ord (comparing)+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 Floating 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 0 = 0/0+ gamma z+ | z == abs z = exp (lnGamma z)+ | otherwise = pi / (sin (pi * z) * exp (lnGamma (1-z)))+++ -- |Natural log of the gamma function+ lnGamma :: a -> a+ lnGamma z = log (gamma z)+ + -- |Natural log of the factorial function+ lnFactorial :: Integral b => b -> a+ lnFactorial n = lnGamma (fromIntegral n+1)++floatGammaInfCutoff :: Double+floatGammaInfCutoff = $( do+ let Just cutoff = findIndex isInfinite (scanl (*) (1::Float) [1..])+ litE (IntegerL (1 + toInteger cutoff))+ )++instance Gamma Float where+ gamma = double2Float . gam . float2Double+ where+ gam x + | x >= floatGammaInfCutoff = 1/0+ | otherwise = case properFraction x of+ (n,0) | n < 1 -> 0/0+ | otherwise -> factorial (n-1)+ _ | x < (-20) -> let s = pi / sin (pi * x)+ in signum s * exp (log (abs s) - lnGamma (1-x))+ | otherwise -> reflect (gammaLanczos g cs) x+ + g = pi+ cs = [ 1.0000000249904433+ , 9.100643759042066+ ,-4.3325519094475+ , 0.12502459858901147+ , 1.1378929685052916e-4+ ,-9.555011214455924e-5+ ]+ + lnGamma = double2Float . reflectLn (lnGammaLanczos g cs) . float2Double+ where+ g = pi+ cs = [ 1.0000000249904433+ , 9.100643759042066+ ,-4.3325519094475+ , 0.12502459858901147+ , 1.1378929685052916e-4+ ,-9.555011214455924e-5+ ]+ + lnFactorial n+ | n' < 0 = error "lnFactorial n: n < 0"+ | n' < toInteger nFacs = facs V.! fromIntegral n+ | otherwise = lnGamma (fromIntegral n+1)+ where+ n' = toInteger n+ nFacs = 2000 -- limited only by time and space+ facs = V.map lnGamma (V.enumFromN 1 nFacs)++doubleGammaInfCutoff :: Double+doubleGammaInfCutoff = $( do+ let Just cutoff = findIndex isInfinite (scanl (*) (1::Double) [1..])+ litE (IntegerL (1 + toInteger cutoff))+ )++instance Gamma Double where+ gamma x + | x >= doubleGammaInfCutoff = 1/0+ | otherwise = case properFraction x of+ (n,0) | n < 1 -> 0/0+ | otherwise -> factorial (n-1)+ _ | x < (-50) -> let s = pi / sin (pi * x)+ in signum s * exp (log (abs s) - lnGammaLanczos g cs (1-x))+ | otherwise -> reflect (gammaLanczos g cs) x+ where+ g = 2*pi+ cs = [ 0.9999999999999858+ , 311.6011750541472+ ,-498.6511904603639+ , 244.08472899976877+ , -38.670364643074194+ , 1.3350900101370549+ , -1.8977221899565682e-3+ , 8.475264614349149e-7+ , 2.59715567376858e-7+ , -2.7166437850607517e-7+ , 6.151114806136299e-8+ ]++ lnGamma = reflectLn (lnGammaLanczos g cs)+ where+ g = exp pi / pi+ cs = [ 1.0000000000000002+ , 1002.5049417114732+ ,-1999.6140446432912+ , 1352.1626218340114+ , -360.6486475548049+ , 33.344988357090685+ , -0.6637188712004668+ , 5.16644552377916e-4+ , 1.684651140163429e-7+ , -1.8148207145896904e-7+ , 6.171532716135051e-8+ , -9.014004881476154e-9+ ]++ lnFactorial n+ | n' < 0 = error "lnFactorial n: n < 0"+ | n' < toInteger nFacs = facs V.! fromIntegral n+ | otherwise = lnGamma (fromIntegral n+1)+ where+ n' = toInteger n+ nFacs = 2000 -- limited only by time and space+ facs = V.map lnGamma (V.enumFromN 1 nFacs)++complexDoubleToFloat :: Complex Double -> Complex Float+complexDoubleToFloat (a :+ b) = double2Float a :+ double2Float b+complexFloatToDouble :: Complex Float -> Complex Double+complexFloatToDouble (a :+ b) = float2Double a :+ float2Double b++instance Gamma (Complex Float) where+ gamma = complexDoubleToFloat . gamma . complexFloatToDouble+ where+ g = pi+ cs = [ 1.0000000249904433+ , 9.100643759042066+ ,-4.3325519094475+ , 0.12502459858901147+ , 1.1378929685052916e-4+ ,-9.555011214455924e-5+ ]+ + lnGamma = complexDoubleToFloat . reflectLnC (lnGammaLanczos g cs) . complexFloatToDouble+ where+ g = pi+ cs = [ 1.0000000249904433+ , 9.100643759042066+ ,-4.3325519094475+ , 0.12502459858901147+ , 1.1378929685052916e-4+ ,-9.555011214455924e-5+ ]++ + lnFactorial n+ | n' < 0 = error "lnFactorial n: n < 0"+ | n' < toInteger nFacs = facs V.! fromIntegral n+ | otherwise = lnGamma (fromIntegral n+1)+ where+ n' = toInteger n+ nFacs = 2000 -- limited only by time and space+ facs = V.map lnGamma (V.enumFromN 1 nFacs)++instance Gamma (Complex Double) where+ gamma = reflectC (gammaLanczos g cs)+ where+ g = 2*pi+ cs = [ 1.0000000000000002+ , 311.60117505414695+ ,-498.65119046033163+ , 244.08472899875767+ , -38.67036462939322+ , 1.3350899103585203+ , -1.8972831806242229e-3+ , -3.935368195357295e-7+ , 2.592464641764731e-6+ , -3.2263565156368265e-6+ , 2.5666169886566876e-6+ , -1.3737776806198937e-6+ , 4.4551204024819644e-7+ , -6.576826592057796e-8+ ]++ lnGamma = reflectLnC (lnGammaLanczos g cs)+ where+ g = exp pi / pi+ cs = [ 1.0000000000000002+ , 1002.5049417114732+ ,-1999.6140446432912+ , 1352.1626218340114+ , -360.6486475548049+ , 33.344988357090685+ , -0.6637188712004668+ , 5.16644552377916e-4+ , 1.684651140163429e-7+ , -1.8148207145896904e-7+ , 6.171532716135051e-8+ , -9.014004881476154e-9+ ]++ lnFactorial n+ | n' < 0 = error "lnFactorial n: n < 0"+ | n' < toInteger nFacs = facs V.! fromIntegral n+ | otherwise = lnGamma (fromIntegral n+1)+ where+ n' = toInteger n+ nFacs = 2000 -- limited only by time and space+ facs = V.map lnGamma (V.enumFromN 1 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 s x = double2Float $ (lowerGamma :: Double -> Double -> Double) (float2Double s) (float2Double x)+ lnLowerGamma s x = double2Float $ (lnLowerGamma :: Double -> Double -> Double) (float2Double s) (float2Double x)+ p s x = double2Float $ (p :: Double -> Double -> Double) (float2Double s) (float2Double x)+ + upperGamma s x = double2Float $ (upperGamma :: Double -> Double -> Double) (float2Double s) (float2Double x)+ lnUpperGamma s x = double2Float $ (lnUpperGamma :: Double -> Double -> Double) (float2Double s) (float2Double x)+ q s x = double2Float $ (q :: Double -> Double -> Double) (float2Double s) (float2Double 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 s x+ | x < 0 = error "lowerGamma: x < 0 is not currently supported."+ | x == 0 = 0+ | x >= s+1 = gamma s - upperGamma s x+ | otherwise = lowerGammaHypGeom s x+ + upperGamma s x+ | x < 0 = error "upperGamma: x < 0 is not currently supported."+ | x == 0 = gamma s+ | x < s+1 = q s x * gamma s+ | otherwise = converge . concat+ $ modifiedLentz 1e-30 (upperGammaCF s x)+ + lnLowerGamma s x+ | x < 0 = error "lnLowerGamma: x < 0 is not currently supported."+ | x == 0 = log 0+ | x >= s+1 = log (p s x) + lnGamma s+ | otherwise = lnLowerGammaHypGeom s x+ + lnUpperGamma s x+ | x < 0 = error "lnUpperGamma: x < 0 is not currently supported."+ | x == 0 = lnGamma s+ | x < s+1 = log (q s x) + lnGamma s+ | otherwise =+ converge (lnUpperGammaConvergents s x)+ + p s x+ | x < 0 = error "p: x < 0 is not currently supported."+ | x == 0 = 0+ | x >= s+1 = 1 - q s x+ | otherwise = pHypGeom s x+ + q s x+ | x < 0 = error "q: x < 0 is not currently supported."+ | x == 0 = 1+ | x < s+1 = 1 - p s x+ | otherwise =+ converge . concat+ $ modifiedLentz 1e-30 (qCF s x)++-- |Factorial function+class Num a => Factorial a where+ factorial :: Integral b => b -> a+ factorial = fromInteger . factorial++instance Factorial Integer where+ factorial n+ | n < 0 = error "factorial: n < 0"+ | otherwise = product [1..toInteger n]++instance Factorial Float where+ factorial = double2Float . factorial+instance Factorial Double where+ factorial n+ | n < 0 = 0/0+ | n < nFacs = facs V.! fromIntegral n+ | otherwise = infinity+ where+ nFacs :: Num a => a+ nFacs = 171 -- any more is pointless, everything beyond here is "Infinity"+ facs = V.scanl (*) 1 (V.enumFromN 1 nFacs)+ infinity = facs V.! nFacs++-- |The beta function: @beta z w@ == @gamma z * gamma w / gamma (z+w)@+beta :: Gamma a => a -> a -> a+beta z w = exp (lnGamma z + lnGamma w - lnGamma (z+w))
+ src/Math/Gamma.hs-boot view
@@ -0,0 +1,22 @@+module Math.Gamma where++class Floating 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 0 = 0/0+ gamma z+ | z == abs z = exp (lnGamma z)+ | otherwise = pi / (sin (pi * z) * exp (lnGamma (1-z)))+++ -- |Natural log of the gamma function+ lnGamma :: a -> a+ lnGamma z = log (gamma z)+ + -- |Natural log of the factorial function+ lnFactorial :: Integral b => b -> a+ lnFactorial n = lnGamma (fromIntegral n+1)++instance Gamma Float where+instance Gamma Double where
+ src/Math/Gamma/Incomplete.hs view
@@ -0,0 +1,143 @@+{-# LANGUAGE ParallelListComp #-}+module Math.Gamma.Incomplete+ ( lowerGammaCF, pCF+ , lowerGammaHypGeom, lnLowerGammaHypGeom, pHypGeom+ , upperGammaCF, lnUpperGammaConvergents, qCF+ ) where++import {-# SOURCE #-} Math.Gamma+import Math.ContinuedFraction+import Math.Sequence.Converge++-- |Continued fraction representation of the lower incomplete gamma function.+lowerGammaCF :: (Floating a, Ord a) => a -> a -> Math.ContinuedFraction.CF a+lowerGammaCF s z = gcf 0+ [ (p,q)+ | p <- pow_x_s_div_exp_x s z+ : interleave+ [negate spn * z | spn <- iterate (1+) s]+ [n * z | n <- iterate (1+) 1]+ | q <- iterate (1+) s+ ]++-- |Lower incomplete gamma function, computed using Kummer's confluent+-- hypergeometric function M(a;b;x). Specifically, this uses the identity:+-- +-- gamma(s,x) = x**s * exp (-x) / s * M(1; 1+s; x)+-- +-- From Abramowitz & Stegun (6.5.12).+--+-- Recommended for use when x < s+1+lowerGammaHypGeom :: Floating b => b -> b -> b+lowerGammaHypGeom 0 0 = 0/0+lowerGammaHypGeom s x = x ** s * exp (negate x) / s * m_1_sp1 s x++-- |Natural logarithm of lower gamma function, based on the same identity as+-- 'lowerGammaHypGeom' and evaluated carefully to avoid overflow and underflow.+-- Recommended for use when x < s+1+lnLowerGammaHypGeom :: Floating a => a -> a -> a+lnLowerGammaHypGeom 0 0 = 0/0+lnLowerGammaHypGeom s x + = log ((signum x)**s * sign_m / signum s)+ + s*log (abs x) - x - log (abs s) + log_m+ where+ (sign_m, log_m) = log_m_1_sp1 s x++-- |Continued fraction representation of the regularized lower incomplete gamma function.+pCF :: (Gamma a, Ord a, Enum a) => a -> a -> CF a+pCF s x = gcf 0+ [ (p,q)+ | p <- pow_x_s_div_gamma_s_div_exp_x s x+ : interleave+ [negate spn * x | spn <- [s..]]+ [n * x | n <- [1..]]+ | q <- [s..]+ ]++-- |Regularized lower incomplete gamma function, computed using Kummer's+-- confluent hypergeometric function. Uses same identity as 'lowerGammaHypGeom'.+-- +-- Recommended for use when x < s+1+pHypGeom :: (Gamma a, Ord a) => a -> a -> a+pHypGeom 0 0 = 0/0+pHypGeom s x+ | s < 0+ = signum x ** s * sin (pi*s) / (-pi)+ * exp (s * log (abs x) - x + lnGamma (-s)) * m_1_sp1 s x++ | s == 0 || x == 0+ = 0++ | otherwise+ = signum x ** s * exp (s * log (abs x) - x - lnGamma (s+1)) * m_1_sp1 s x++-- |Continued fraction representation of the regularized upper incomplete gamma function.+-- Recommended for use when x >= s+1+qCF :: (Gamma a, Ord a, Enum a) => a -> a -> CF a+qCF s x = gcf 0+ [ (p,q)+ | p <- pow_x_s_div_gamma_s_div_exp_x s x+ : zipWith (*) [1..] (iterate (subtract 1) (s-1))+ | q <- [n + x - s | n <- [1,3..]]+ ]++-- |Continued fraction representation of the upper incomplete gamma function.+-- Recommended for use when x >= s+1+upperGammaCF :: (Floating a, Ord a) => a -> a -> CF a+upperGammaCF s z = gcf 0+ [ (p,q)+ | p <- pow_x_s_div_exp_x s z+ : zipWith (*) (iterate (1+) 1) (iterate (subtract 1) (s-1))+ | q <- [n + z - s | n <- iterate (2+) 1]+ ]++-- |Natural logarithms of the convergents of the upper gamma function, +-- evaluated carefully to avoid overflow and underflow.+-- Recommended for use when x >= s+1+lnUpperGammaConvergents :: Floating a => a -> a -> [a]+lnUpperGammaConvergents s x = map (a -) (concat (eval theCF)) + where + eval = map (map evalSign) . modifiedLentzWith signLog addSignLog negateSignLog 1e-30+ + a = s * log x - x+ theCF = gcf (x + 1 - s)+ [ (p,q)+ | p <- zipWith (*) (iterate (1+) 1) (iterate (subtract 1) (s-1))+ | q <- [n + x - s | n <- iterate (2+) 3]+ ]++---- various utility functions ----++evalSign (s,x) = log s + x+signLog x = (signum x, log (abs x))+addSignLog (xS,xL) (yS,yL) = (xS*yS, xL+yL)+negateSignLog (s,l) = (s, negate l)++-- |Special case of Kummer's confluent hypergeometric function, used+-- in lower gamma functions.+-- +-- m_1_sp1 s z = M(1;s+1;z)+-- +m_1_sp1 s z = converge . scanl (+) 0 . scanl (*) 1 $+ [z / x | x <- iterate (1+) (s+1)]++log_m_1_sp1 s z = converge (concat (log_m_1_sp1_convergents s z))++log_m_1_sp1_convergents s z+ = modifiedLentzWith signLog addSignLog negateSignLog 1e-30+ $ sumPartialProducts (1:[z / x | x <- iterate (1+) (s+1)])++interleave [] _ = []+interleave _ [] = []+interleave (x:xs) ys = x:interleave ys xs++-- A common subexpression appearing in both 'pCF' and 'qCF'.+pow_x_s_div_gamma_s_div_exp_x s x + | x > 0 = exp (log x * s - x - lnGamma s)+ | otherwise = x ** s / (exp x * gamma s)++-- The corresponding subexpression from 'lowerGammaCF' and 'upperGammaCF'+pow_x_s_div_exp_x s x + | x > 0 = exp (log x * s - x)+ | otherwise = x ** s / exp x+
+ src/Math/Gamma/Lanczos.hs view
@@ -0,0 +1,89 @@+{-# 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 _ [] _ = error "gammaLanczos: empty coefficient list"+gammaLanczos g cs zp1+ = sqrt (2*pi) * x ** (zp1 - 0.5) * exp (negate x) * a cs z+ where+ x = zp1 + (g - 0.5)+ z = zp1 - 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 _ [] _ = error "lnGammaLanczos: empty coefficient list"+lnGammaLanczos g cs zp1 + = log (sqrt (2*pi)) + log x * (zp1 - 0.5) - x + log (a cs z)+ where + x = zp1 + (g - 0.5)+ z = zp1 - 1++{-# INLINE a #-}+a [] z = error "Math.Gamma.Lanczos.a: empty coefficient list"+a cs z = head cs + sum [c / (z + k) | c <- tail cs | k <- iterate (1+) 1]++-- |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 gamma z+ | z > 0.5 = gamma z+ | otherwise = case properFraction z of+ (_,0) -> 0/0+ _ -> pi / (sin (pi * z) * gamma (1-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 gamma z+ | realPart z > 0.5 = gamma z+ | imagPart z == 0+ && snd (properFraction (realPart z)) == 0+ = 0/0+ | otherwise = pi / (sin (pi * z) * gamma (1-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 lnGamma z+ | z > 0.5 = lnGamma z+ | otherwise = case properFraction z of+ (_,0) -> log (0/0)+ _ -> log pi - log (sin (pi * z)) - lnGamma (1-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 lnGamma z+ | realPart z > 0.5 = lnGamma z+ | imagPart z == 0+ && snd (properFraction (realPart z)) == 0+ = log (0/0)+ | otherwise = log pi - log (sin (pi * z)) - lnGamma (1-z)+
+ src/Math/Gamma/Stirling.hs view
@@ -0,0 +1,81 @@+{-# LANGUAGE ParallelListComp #-}+-- |Stirling's approximation to the gamma function and utility functions for+-- selecting coefficients.+module Math.Gamma.Stirling (lnGammaStirling, cs, s, abs_s, terms) where++import qualified Data.Vector as V++-- |Convergent when Re(z) > 0. The first argument is the c_n series to use +-- ('cs' is an ineffecient but generic definition of the full infinite series.+-- Some precomputed finite prefix of 'cs' should be fed to this function, the +-- length of which will determine the accuracy achieved.)+{-# INLINE lnGammaStirling #-}+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++{-# INLINE risingPowers #-}+risingPowers x = scanl1 (*) (iterate (1+) x)++-- |The c_n series in the convergent version of Stirling's approximation given+-- on wikipedia at+-- http:\/\/en.wikipedia.org\/wiki\/Stirling%27s_approximation#A_convergent_version_of_Stirling.27s_formula+-- as fetched on 11 June 2010.+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++-- |The (signed) Stirling numbers of the first kind.+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 (n-1) (k-1) - (toInteger n-1) * s (n-1) k+ where+ s n k = table !! n V.! k++-- |The (unsigned) Stirling numbers of the first kind.+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 (n-1) (k-1) + (toInteger n-1) * abs_s (n-1) k+ where+ abs_s n k = table !! n V.! k++-- |Compute the number of terms required to achieve a given precision for a+-- given value of z. The mamimum will typically (always?) be around 1, and +-- seems to be more or less independent of the precision desired (though not +-- of the machine epsilon - essentially, near zero I think this method is+-- extremely numerically unstable).+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 is 0 at z=0.86639115674955 and z=2.087930091329227)+ 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(x-y)<=eps = n | otherwise = go (n+1) (y:zs)++f z = scanl1 (+) [c / q | c <- cs | q <- risingPowers (z+1)]