diff --git a/gamma.cabal b/gamma.cabal --- a/gamma.cabal +++ b/gamma.cabal @@ -1,5 +1,5 @@ name: gamma -version: 0.7.0.1 +version: 0.9.0.1 stability: provisional cabal-version: >= 1.6 @@ -13,11 +13,12 @@ category: Math, Numerical synopsis: Gamma function and related functions. description: Approximations of the gamma function, incomplete gamma - functions, beta function, and factorials. + functions, and factorials. tested-with: GHC == 6.10.4, GHC == 6.12.1, GHC == 6.12.3, - GHC == 7.0.1, GHC == 7.0.2 + GHC == 7.0.1, GHC == 7.0.2, + GHC == 7.2.1 extra-source-files: extras/*.hs @@ -28,7 +29,9 @@ Library hs-source-dirs: src - exposed-modules: Math.Gamma + ghc-options: -Wall + exposed-modules: Math.Factorial + Math.Gamma Math.Gamma.Incomplete Math.Gamma.Stirling Math.Gamma.Lanczos @@ -36,4 +39,4 @@ continued-fractions >= 0.9.1, converge, template-haskell, - vector >= 0.5 && < 0.8 + vector >= 0.5 diff --git a/src/Math/Factorial.hs b/src/Math/Factorial.hs new file mode 100644 --- /dev/null +++ b/src/Math/Factorial.hs @@ -0,0 +1,33 @@ +{-# LANGUAGE FlexibleInstances #-} +module Math.Factorial where + +import Data.Complex +import qualified Data.Vector.Unboxed as V +import GHC.Float (double2Float) + +-- |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 (Complex Float) where + factorial = (:+ 0) . 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 +instance Factorial (Complex Double) where + factorial = (:+ 0) . factorial diff --git a/src/Math/Gamma.hs b/src/Math/Gamma.hs --- a/src/Math/Gamma.hs +++ b/src/Math/Gamma.hs @@ -3,15 +3,14 @@ ( Gamma(..) , Factorial(..) , IncGamma(..) - , beta ) where import Math.Gamma.Lanczos import Math.Gamma.Incomplete +import Math.Factorial import Data.Complex -import Data.List (sortBy, findIndex) -import Data.Ord (comparing) +import Data.List (findIndex) import GHC.Float (float2Double, double2Float) import qualified Data.Vector.Unboxed as V import Language.Haskell.TH (litE, Lit(IntegerL)) @@ -19,7 +18,7 @@ import Math.Sequence.Converge -- |Gamma function. Minimal definition is ether 'gamma' or 'lnGamma'. -class Floating a => Gamma a where +class (Floating a, Factorial 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 @@ -50,7 +49,7 @@ | x >= floatGammaInfCutoff = 1/0 | otherwise = case properFraction x of (n,0) | n < 1 -> 0/0 - | otherwise -> factorial (n-1) + | otherwise -> factorial (n-1 :: Integer) _ | x < (-20) -> let s = pi / sin (pi * x) in signum s * exp (log (abs s) - lnGamma (1-x)) | otherwise -> reflect (gammaLanczos g cs) x @@ -95,7 +94,7 @@ | x >= doubleGammaInfCutoff = 1/0 | otherwise = case properFraction x of (n,0) | n < 1 -> 0/0 - | otherwise -> factorial (n-1) + | otherwise -> factorial (n-1 :: Integer) _ | 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 @@ -147,16 +146,6 @@ 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 @@ -296,29 +285,3 @@ 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)) diff --git a/src/Math/Gamma.hs-boot b/src/Math/Gamma.hs-boot --- a/src/Math/Gamma.hs-boot +++ b/src/Math/Gamma.hs-boot @@ -1,6 +1,8 @@ module Math.Gamma where -class Floating a => Gamma a where +import Math.Factorial + +class (Floating a, Factorial 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 diff --git a/src/Math/Gamma/Incomplete.hs b/src/Math/Gamma/Incomplete.hs --- a/src/Math/Gamma/Incomplete.hs +++ b/src/Math/Gamma/Incomplete.hs @@ -108,9 +108,16 @@ ---- various utility functions ---- +evalSign :: Floating a => (a,a) -> a evalSign (s,x) = log s + x + +signLog :: Floating a => a -> (a,a) signLog x = (signum x, log (abs x)) + +addSignLog :: (Num a, Num b) => (a,b) -> (a,b) -> (a,b) addSignLog (xS,xL) (yS,yL) = (xS*yS, xL+yL) + +negateSignLog :: (Num b) => (a,b) -> (a,b) negateSignLog (s,l) = (s, negate l) -- |Special case of Kummer's confluent hypergeometric function, used @@ -118,25 +125,31 @@ -- -- m_1_sp1 s z = M(1;s+1;z) -- +m_1_sp1 :: Fractional a => a -> a -> a m_1_sp1 s z = converge . scanl (+) 0 . scanl (*) 1 $ [z / x | x <- iterate (1+) (s+1)] +log_m_1_sp1 :: Floating a => a -> a -> (a,a) log_m_1_sp1 s z = converge (concat (log_m_1_sp1_convergents s z)) +log_m_1_sp1_convergents :: Floating a => a -> a -> [[(a,a)]] log_m_1_sp1_convergents s z = modifiedLentzWith signLog addSignLog negateSignLog 1e-30 $ sumPartialProducts (1:[z / x | x <- iterate (1+) (s+1)]) +interleave :: [a] -> [a] -> [a] 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 :: (Gamma a, Ord a) => a -> a -> a 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 :: (Floating a, Ord a) => a -> a -> a pow_x_s_div_exp_x s x | x > 0 = exp (log x * s - x) | otherwise = x ** s / exp x diff --git a/src/Math/Gamma/Lanczos.hs b/src/Math/Gamma/Lanczos.hs --- a/src/Math/Gamma/Lanczos.hs +++ b/src/Math/Gamma/Lanczos.hs @@ -42,18 +42,22 @@ z = zp1 - 1 {-# INLINE a #-} -a [] z = error "Math.Gamma.Lanczos.a: empty coefficient list" +a :: Fractional t => [t] -> t -> t +a [] _ = error "Math.Gamma.Lanczos.a: empty coefficient list" a cs z = head cs + sum [c / (z + k) | c <- tail cs | k <- iterate (1+) 1] +fractionalPart :: RealFloat a => a -> a +fractionalPart x = case properFraction x of + (i,f) -> let _ = i :: Int in 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 gamma z - | z > 0.5 = gamma z - | otherwise = case properFraction z of - (_,0) -> 0/0 - _ -> pi / (sin (pi * z) * gamma (1-z)) + | z > 0.5 = gamma z + | fractionalPart z == 0 = 0 + | otherwise = 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. @@ -62,7 +66,7 @@ reflectC gamma z | realPart z > 0.5 = gamma z | imagPart z == 0 - && snd (properFraction (realPart z)) == 0 + && fractionalPart (realPart z) == 0 = 0/0 | otherwise = pi / (sin (pi * z) * gamma (1-z)) @@ -71,10 +75,9 @@ {-# 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) + | z > 0.5 = lnGamma z + | fractionalPart z == 0 = log (0/0) + | otherwise = 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. @@ -83,7 +86,7 @@ reflectLnC lnGamma z | realPart z > 0.5 = lnGamma z | imagPart z == 0 - && snd (properFraction (realPart z)) == 0 + && fractionalPart (realPart z) == 0 = log (0/0) | otherwise = log pi - log (sin (pi * z)) - lnGamma (1-z) diff --git a/src/Math/Gamma/Stirling.hs b/src/Math/Gamma/Stirling.hs --- a/src/Math/Gamma/Stirling.hs +++ b/src/Math/Gamma/Stirling.hs @@ -11,11 +11,15 @@ -- 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 +lnGammaStirling ks z + = (z - 0.5) * log z + - z + + 0.5 * log (2*pi) + + sum (zipWith (/) ks (risingPowers (z+1))) {-# INLINE risingPowers #-} -risingPowers x = scanl1 (*) (iterate (1+) x) +risingPowers :: Num a => a -> [a] +risingPowers = scanl1 (*) . iterate (1+) -- |The c_n series in the convergent version of Stirling's approximation given -- on wikipedia at @@ -34,17 +38,17 @@ | 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 + | 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 + table = [V.generate (i+1) $ s' i | i <- [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 + s'' n'' k'' = table !! n'' V.! k'' -- |The (unsigned) Stirling numbers of the first kind. abs_s :: Int -> Int -> Integer @@ -52,17 +56,17 @@ | 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 + | 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 + 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 + 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 @@ -70,12 +74,15 @@ -- 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) +terms prec z = stepsNeeded (eps z) (f z) where cs' = cs - f z = scanl1 (+) [c / q | c <- cs' | q <- risingPowers (z+1)] + f x = scanl1 (+) (zipWith (/) cs' (risingPowers (x+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)] + eps x = prec * abs ((x - 0.5) * log x - x + 0.5 * log (2*pi)) + stepsNeeded e xs = go 1 xs + where + go n (x:y:zs) + | abs(x-y)<=e = n + | otherwise = go (n+1) (y:zs) + go n _ = n -- this case should be impossible