gamma 0.7.0.1 → 0.9.0.1

raw patch · 7 files changed

+108/−84 lines, 7 filesdep ~vectorPVP ok

version bump matches the API change (PVP)

Dependency ranges changed: vector

API changes (from Hackage documentation)

Files

• gamma.cabal +8/−5
• src/Math/Factorial.hs +33/−0
• src/Math/Gamma.hs +5/−42
• src/Math/Gamma.hs-boot +3/−1
• src/Math/Gamma/Incomplete.hs +13/−0
• src/Math/Gamma/Lanczos.hs +14/−11
• src/Math/Gamma/Stirling.hs +32/−25

@@ -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

@@ -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

@@ -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))

@@ -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

@@ -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

@@ -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) 

@@ -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