gamma 0.9.0.2 → 0.10.0.0

raw patch · 13 files changed

+547/−34 lines, 13 filesdep +QuickCheckdep +erfdep +gammadep ~basedep ~continued-fractionssetup-changednew-uploaderPVP ok

version bump matches the API change (PVP)

Dependencies added: QuickCheck, erf, gamma, numbers, test-framework, test-framework-quickcheck2

Dependency ranges changed: base, continued-fractions

API changes (from Hackage documentation)

Files

• CHANGES.md +11/−0
• README.md +31/−0
• Setup.lhs +1/−0
• gamma.cabal +38/−18
• src/Math/Factorial.hs +5/−3
• src/Math/Gamma.hs +36/−8
• src/Math/Gamma/Incomplete.hs +4/−3
• src/Math/Gamma/Lanczos.hs +1/−1
• src/Math/Gamma/Stirling.hs +6/−1
• test-suite/Math/GammaTests.hs +248/−0
• test-suite/Math/IncGammaTests.hs +104/−0
• test-suite/Math/Reference.hs +49/−0
• test-suite/Test.hs +13/−0

@@ -0,0 +1,11 @@+0.10.0.0+=====++* Temporarily move project's repository from [mokus0](https://github.com/mokus0/gamma) to+  [rockbmb](https://github.com/rockbmb/gamma2). This was a duplication, not a fork, so+  users can create issues in the latter repository as well.+* Enable Travis CI for the project's new repository.+* Increase lower bound for the Cabal version to be used when compiling this+  project from `>=1.6` to `>=1.10`.+* Change the supported GHC versions from the interval [6.10.4 .. 7.4.1] to+  [8.0.2 .. 8.6.1].

@@ -0,0 +1,31 @@+Gamma+=========++[![Build Status](https://travis-ci.com/rockbmb/gamma2.svg?branch=master)](https://travis-ci.com/rockbmb/gamma2)++`gamma` is a number theory package written in Haskell that aims to provide its users+with the following functionality:++* A typeclass+  ```haskell+  class Num a => Factorial a+  ```+  that provides a [`factorial`](https://en.wikipedia.org/wiki/Factorial) function.+* A typeclass+  ```haskell+  class (Eq a, Floating a, Factorial a) => Gamma a+  ```+  that provides the functions+  + [`gamma`](https://en.wikipedia.org/wiki/Gamma_function)+  + `lnGamma` (natural logarithm of the gamma function) and+  + `lnFactorial` (natural logarithm of the factorial function).+* A typeclass+  ```haskell+  class Gamma a => IncGamma a+  ```+  for the [incomplete lower and upper gamma functions](https://en.wikipedia.org/wiki/Incomplete_gamma_function).+* A typeclass+  ```haskell+  class Gamma a => GenGamma a+  ```+  to represent the [multivariate gamma function](https://en.wikipedia.org/wiki/Multivariate_gamma_function).

@@ -1,5 +1,6 @@ #!/usr/bin/env runhaskell  > import Distribution.Simple+ > main = defaultMain 

@@ -1,34 +1,33 @@ name:                   gamma-version:                0.9.0.2+version:                0.10.0.0 stability:              provisional -cabal-version:          >= 1.6+cabal-version:          >= 1.10 build-type:             Simple -author:                 James Cook <mokus@deepbondi.net>-maintainer:             James Cook <mokus@deepbondi.net>+author:                 James Cook <mokus ΑΤ deepbondi dοt net>+maintainer:             James Cook <mokus ΑΤ deepbondi dοt net>,+                        Alexandre Rodrigues Baldé <alexandrer_b ΑΤ outlook dοt com> license:                PublicDomain-homepage:               https://github.com/mokus0/gamma+homepage:               https://github.com/rockbmb/gamma2+bug-reports:            https://github.com/rockbmb/gamma2/issues  category:               Math, Numerical synopsis:               Gamma function and related functions.-description:            Approximations of the gamma function, incomplete gamma -                        functions, and factorials.+description:            Approximations of the gamma function, incomplete gamma+                        functions, generalized gamma 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.2.1,-                        GHC == 7.4.1-rc1+tested-with:            GHC==8.0.2, GHC==8.2.2, GHC==8.4.3, GHC==8.6.1 -extra-source-files:     extras/*.hs+extra-source-files:     extras/*.hs,+                        CHANGES.md,+                        README.md  source-repository head   type: git-  location: git://github.com/mokus0/gamma.git-+  location: git://github.com/rockbmb/gamma2.git -Library+library   hs-source-dirs:       src   ghc-options:          -Wall   exposed-modules:      Math.Factorial@@ -36,8 +35,29 @@                         Math.Gamma.Incomplete                         Math.Gamma.Stirling                         Math.Gamma.Lanczos-  build-depends:        base >= 3 && <5,-                        continued-fractions >= 0.9.1,+  build-depends:        base >= 4.9 && <5,+                        continued-fractions >= 0.10,                         converge,                         template-haskell,                         vector >= 0.5+  default-language: Haskell2010++test-suite gamma-tests+  ghc-options:          -threaded+  hs-source-dirs:       test-suite+  default-language:     Haskell2010+  ghc-options:          -Wall+  type:                 exitcode-stdio-1.0+  main-is:              Test.hs+  build-depends:+    base >= 4.9 && <5,+    erf,+    gamma,+    numbers,+    QuickCheck >= 2.1.1,+    test-framework,+    test-framework-quickcheck2+  other-modules:+    Math.GammaTests+    Math.IncGammaTests+    Math.Reference

@@ -1,9 +1,11 @@ {-# LANGUAGE FlexibleInstances #-}-module Math.Factorial where+module Math.Factorial+    ( Factorial (..)+    ) where -import Data.Complex+import Data.Complex                  (Complex (..)) import qualified Data.Vector.Unboxed as V-import GHC.Float (double2Float)+import GHC.Float                     (double2Float)  -- |Factorial function class Num a => Factorial a where

@@ -3,19 +3,20 @@     ( Gamma(..)     , Factorial(..)     , IncGamma(..)+    , GenGamma(..)     ) where  import Math.Gamma.Lanczos import Math.Gamma.Incomplete import Math.Factorial -import Data.Complex-import Data.List (findIndex)-import GHC.Float (float2Double, double2Float)+import Data.Complex                  (Complex (..))+import Data.List                     (findIndex)+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+import Language.Haskell.TH           (litE, Lit(IntegerL))+import Math.ContinuedFraction        (modifiedLentz)+import Math.Sequence.Converge        (converge)  -- |Gamma function.  Minimal definition is ether 'gamma' or 'lnGamma'. class (Eq a, Floating a, Factorial a) => Gamma a where@@ -27,11 +28,10 @@         | 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)@@ -284,4 +284,32 @@         | otherwise =             converge . concat             $ modifiedLentz 1e-30 (qCF s x)++-- |Generalized gamma function (also known as multivariate gamma function).+-- See http://en.wikipedia.org/wiki/Multivariate_gamma_function+-- (Fetched Feb 29, 2012).+class Gamma a => GenGamma a where+    -- |Generalized gamma function.  generalizedGamma p x = (pi ** ((p - 1) / 2)) * gamma x * generalizedGamma (p - 1) (x - 0.5)+    generalizedGamma :: Int -> a -> a+    -- |Natural log of generalizedGamma+    lnGeneralizedGamma :: Int -> a -> a++instance GenGamma Float where+    generalizedGamma   u x = double2Float $ (generalizedGamma   :: Int -> Double -> Double) u (float2Double x)+    lnGeneralizedGamma u x = double2Float $ (lnGeneralizedGamma :: Int -> Double -> Double) u (float2Double x)++instance GenGamma Double where+    generalizedGamma u x+        | u <= 0    = error "generalizedGamma p x: p must be strictly positive."+        | u == 1    = gamma x+        | otherwise = pi ** (0.25 * u' * (u' - 1))+            * product [ gamma (x - 0.5 * fromIntegral j) | j <- [0 .. u - 1]]+        where u' = fromIntegral u++    lnGeneralizedGamma u x+        | u <= 0    = error "lnGeneralizedGamma p x: p must be strictly positive."+        | u == 1    = lnGamma x+        | otherwise = (0.25 * log pi) * (u' * (u' - 1))+            + sum [ lnGamma (x - 0.5 * fromIntegral j) | j <- [0 .. u - 1]]+        where u' = fromIntegral u 

@@ -5,9 +5,10 @@     , upperGammaCF, lnUpperGammaConvergents, qCF     ) where -import {-# SOURCE #-}  Math.Gamma-import Math.ContinuedFraction-import Math.Sequence.Converge+import {-# SOURCE #-} Math.Gamma++import Math.ContinuedFraction    (CF, gcf, modifiedLentzWith, sumPartialProducts)+import Math.Sequence.Converge    (converge)  -- |Continued fraction representation of the lower incomplete gamma function. lowerGammaCF :: (Floating a, Ord a) => a -> a -> Math.ContinuedFraction.CF a

@@ -14,7 +14,7 @@     , reflectLn, reflectLnC     ) where -import Data.Complex+import Data.Complex (Complex, imagPart, realPart)  -- |Compute Lanczos' approximation to the gamma function, using the specified -- constants.  Valid for Re(x) > 0.5.  Use 'reflect' or 'reflectC' to extend

@@ -1,7 +1,12 @@ {-# 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+module Math.Gamma.Stirling+    ( lnGammaStirling+    , cs+    , s+    , abs_s+    , terms) where  import qualified Data.Vector as V 

@@ -0,0 +1,248 @@+{-# LANGUAGE FlexibleContexts #-}+{-# LANGUAGE ImplicitParams   #-}+{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.GammaTests+    ( eps+    , isSane+    , tests+    , (~=)+    ) where++import qualified Math.Gamma                 as G+import Math.Reference                       (err, lgamma, tgamma)++import Data.Complex                         (Complex (..), conjugate, imagPart,+                                             magnitude, realPart)+import Test.Framework                       (Test, testGroup)+import Test.Framework.Providers.QuickCheck2 (testProperty)+import Test.QuickCheck                      (Arbitrary, NonNegative (..), Positive (..),+                                             (==>))++gammaSign :: (RealFrac a, Fractional t) => a -> t+gammaSign x+    | x > 0     = 1+    | otherwise = case properFraction x of+        (_, 0)              -> 0/0+        (n, _)  | odd n     ->   1+                | otherwise ->  -1++gammaArg :: RealFrac a => a -> Bool+gammaArg = not . isNaN . gammaSign++eps :: RealFloat a => a+eps = eps'+    where+        eps' = encodeFloat 1 (1 - floatDigits eps')++infix 4 ~=+(~=) :: (Ord a, Num a, ?eps::a, ?mag::t -> a, Num t)+     => t -> t -> Bool+x ~= y+    =  absErr <= ?eps+    || absErr <= ?eps * min (?mag x) (?mag y)+    where absErr = ?mag (x-y)++isSane :: RealFloat a => a -> Bool+isSane x = all (\f -> not (f x)) [isNaN, isInfinite, isDenormalized]++tests :: [Test]+tests =+    [ testGroup "Float"+        [ testGroup "Native"  (realTests G.gamma (G.lnGamma :: Float -> Float))+        , testGroup "Complex" (complexTests G.gamma (G.lnGamma :: Complex Float -> Complex Float))+        ]+    , testGroup "Double"+        [ testGroup "Native"  (realTests    G.gamma (G.lnGamma :: Double -> Double))+        , testGroup "FFI"     (realTests    tgamma lgamma)+        , testGroup "Complex" (complexTests G.gamma (G.lnGamma :: Complex Double -> Complex Double))+        ]+    ]++realTests :: (Arbitrary a, G.Gamma a, Show a, RealFloat a)+          => (a -> a) -> (a -> a) -> [Test]+realTests gamma lnGamma = +    let ?mag = abs+        ?eps = eps+        ?complex = False+     in [ testGroup "gamma"       (realGammaTests    gamma)+        , testGroup "lnGamma"     (realLogGammaTests gamma lnGamma)+        , testGroup "lnFactorial" (logFactorialTests lnGamma G.lnFactorial)+        ]++complexTests :: (Arbitrary a, G.Gamma (Complex a), G.Gamma a, Show a, RealFloat a)+             => (Complex a -> Complex a) -> (Complex a -> Complex a) -> [Test]+complexTests gamma lnGamma = +    let ?mag = magnitude+        ?eps = eps+        ?complex = True+     in [ testGroup "gamma"   (complexGammaTests    gamma)+        , testGroup "lnGamma" (complexLogGammaTests gamma lnGamma)+        , testGroup "lnFactorial" (logFactorialTests lnGamma G.lnFactorial)+        ]++realGammaTests+    :: (Arbitrary a, ?mag::a -> t, ?complex::Bool, Show a, RealFloat a, RealFloat t)+    => (a -> a) -> [Test]+realGammaTests gamma = +    gammaTests gamma +++    [ testProperty "between factorials" $ \(Positive x) -> +        let gam w = fromInteger (product [1..w-1])+            gamma_x = gamma x `asTypeOf` eps+         in x > 2 && isSane gamma_x+            ==> gam (floor x) <= gamma_x && gamma_x <= gam (ceiling x)+    , testProperty "agrees with factorial" $ \(Positive x) ->+        let gam w = fromInteger (product [1..w-1])+            gamma_x = gamma (fromInteger x)+         in isSane gamma_x ==> +            let ?eps = 16*eps in gam x ~= gamma_x+    , testProperty "agrees with C tgamma" $ \x ->+        let a = gamma x+            b = realToFrac (tgamma (realToFrac x))+         in isSane a ==> +            let ?eps = 512*eps in a ~= b+    , testProperty "agrees with reference implementation" $ \x ->+        let a = gamma x+         in isSane a ==> snd (err gamma x) <= 256*eps+    , testProperty "monotone when x>2" $ \(Positive x) ->+        let x' = x * (1 + 256*eps)+            a = gamma x+            b = gamma x'+         in (x > 2) && (x <= x') && all isSane [a,b] ==> a <= b+    , testProperty "domain check" $ \x ->+        let a = gamma x+         in not (isInfinite a) ==>+            isNaN a == not (gammaArg x)+    , testProperty "sign check" $ \x ->+        let a = gamma x; s = gammaSign x+            signum' z@0 | isNegativeZero z  = -1+                        | otherwise         = 1+            signum' z = signum z+         in gammaArg x && not (isInfinite a) ==>+            signum' a == s+    ]+++complexGammaTests+    :: (Arbitrary a1, G.Gamma a1, ?mag::Complex a1 -> a, ?eps::a,+        ?complex::Bool, Show a1, RealFloat a1, RealFloat a)+    => (Complex a1 -> Complex a1) -> [Test]+complexGammaTests gamma = +    gammaTests gamma +++    [ testProperty "conjugate" $ \x -> +        let gam = gamma x+         in isSane (magnitude gam) ==> conjugate gam ~= gamma (conjugate x)+    , testProperty "real argument" $ \(Positive x) ->+        let z = x :+ 0+            gam = G.gamma x+         in isSane gam ==> +            let ?mag = abs; ?eps = 512 * eps+             in gam ~= realPart (gamma z)+    ]++gammaTests+    :: (Arbitrary t2, ?mag::t2 -> a, ?complex::Bool, Show t2,+       RealFloat a, Floating t2)+    => (t2 -> t2) -> [Test]+gammaTests gamma =+    [ testProperty "increment arg" $ \x ->+        let a = gamma x+            b = gamma (x + 1)+            margin  | ?complex  = 32+                    | otherwise = 32+         in all (isSane . ?mag) [a,b,recip a, recip b]+            ==> ?mag (a - b/x) <= margin * (max 2 (1 + recip (?mag x))) * eps * ?mag a+             || ?mag (a*x - b) <= margin * (max 2 (1 +        ?mag x))  * eps * ?mag b+    , testProperty "reflect" $ \x ->+        ?mag x > 0 ==>+        let a = gamma x+            b = gamma (1 - x)+            c = sin (pi * x) / pi+            margin  | ?complex  = 16+                    | otherwise = 256+         in all (isSane . ?mag) [a,b,c]+            -- There may be tighter bounds to be found but I haven't+            -- been able to derive them yet.+            ==> ?mag (a*b*c-1) <= margin * eps * (1 + ?mag c * (1 + ?mag (a+b)))+             || ?mag (a*b*c-1) <= margin * eps * (?mag (a*b) + ?mag (a*c) + ?mag (b*c))+    ]++logGammaTests+    :: (Arbitrary t, ?mag::t -> t1, ?complex::Bool, Show t,+        RealFloat t1, Ord a1, Ord a, Num a1, Num a, Floating t)+    => (t -> t) -> (t -> t) -> (t -> a) -> (t -> a1) -> [Test]+logGammaTests gamma lnGamma real imag =+    [ testProperty "increment arg" $ \x ->+        let gam = lnGamma (x+1)+         in real x > 0 && isSane (?mag gam)+            ==> +            let ?eps = 32 * eps+             in gam - log x ~= lnGamma x+                || gam ~= log x + lnGamma x+    , testProperty "reflect" $ \x ->+        ?mag x > 0 ==>+        let a = lnGamma x+            b = lnGamma (1 - x)+            c = log pi - c';    c' = log (sin (pi * x))+         in all (isSane . ?mag) [a,b,c] +            ==> +            let ?eps = 512 * eps+             in     a + b ~= c+                 || a ~= c - b+                 || b ~= c - a+                 || a - b - c' ~= log pi+    , testProperty "agrees with log . gamma" $ \x ->+        let a = log b'    ; a' = exp b+            b = lnGamma x ; b' = gamma x+            margin  | ?complex   = 1024+                    | otherwise  = 16+         in   (real x > 1 || abs (imag x) > 1)+            && all (isSane . ?mag) [a,b,a',b'] ==>+            let ?eps = margin * eps+             in a ~= b || a' ~= b'+    ]++realLogGammaTests+    :: (Arbitrary a, Show a, RealFloat a)+    => (a -> a) -> (a -> a) -> [Test]+realLogGammaTests gamma lnGamma = +    let ?mag = abs+        ?complex = False+     in logGammaTests gamma lnGamma id (const 0) +++        [ testProperty "between factorials" $ \(Positive x) -> +            let gam w = sum $ map (log.fromInteger) [1 .. w-1]+                gamma_x = lnGamma x `asTypeOf` eps+             in x > 2 && isSane gamma_x+                ==> gam (floor x) <= gamma_x && gamma_x <= gam (ceiling x)+        , let ?eps = 2 * eps+           in testProperty "agrees with C lgamma" $ \(NonNegative x) ->+            let a = lnGamma x+                b = realToFrac (lgamma (realToFrac x))+             in let ?eps = 512 * eps+                 in isSane a ==> a ~= b+        ]++complexLogGammaTests+    :: (Arbitrary a, G.Gamma a, Show a, RealFloat a)+    => (Complex a -> Complex a) -> (Complex a -> Complex a) -> [Test]+complexLogGammaTests gamma lnGamma = +    let ?mag = magnitude+        ?complex = True+     in logGammaTests gamma lnGamma realPart imagPart +++        [ testProperty "real argument" $ \(Positive x) ->+            let z = x :+ 0+                gam = G.lnGamma x+             in isSane gam ==> +                let ?eps = 8 * eps+                 in (gam :+ 0) ~= lnGamma z+        ]++logFactorialTests+    :: (?mag::t -> a, ?eps::a, RealFloat a, Num t1, Num t)+    => (t1 -> t) -> (Integer -> t) -> [Test]+logFactorialTests lnGamma lnFactorial =+    [ testProperty "agrees with lnGamma" $ \x ->+        let gam = lnGamma (fromInteger x + 1)+         in isSane (?mag gam)+            ==> let ?eps = 64*eps in gam ~= lnFactorial x+    ]

@@ -0,0 +1,104 @@+{-# LANGUAGE ImplicitParams #-}+module Math.IncGammaTests+    ( tests+    ) where++import Math.GammaTests                      (eps, isSane, (~=))++import Data.Number.Erf                      (Erf (..))+import Math.Gamma                           (Gamma (..), IncGamma (..))+import Test.Framework                       (Test, testGroup)+import Test.Framework.Providers.QuickCheck2 (testProperty)+import Test.QuickCheck                      (Arbitrary, NonNegative (..), Positive (..),+                                             (==>))++tests :: [Test]+tests = +    [ testGroup "incomplete gamma"+        [ testGroup "Float"  (incompleteGammaTests (eps :: Float))+        , testGroup "Double" (incompleteGammaTests (eps :: Double))+        ]+    ]++incompleteGammaTests+    :: (Arbitrary a, Show a, RealFloat a, Erf a, IncGamma a)+    => a -> [Test]+incompleteGammaTests epsi =+    let ?mag = abs+     in [ testProperty "lowerGamma + upperGamma" $ \s (NonNegative x) ->+            let a = lowerGamma s x+                b = upperGamma s x+                c = gamma s+             in all isSane [a,b,c] ==> +                let ?eps = 512 * epsi+                 in  a+b ~= c+                    || a ~= c-b+                    || b ~= c-a+        , testProperty "p + q" $ \s (NonNegative x) ->+            let a = p s x+                b = q s x+             in all isSane [a,b] ==>+                let ?eps = 256 * epsi+                 in  a+b ~= 1+                    || a ~= 1-b+                    || b ~= 1-a+        , testGroup "upperGamma"+            [ testProperty "increment s" $ \s (NonNegative x) ->+                let a = upperGamma (s+1) x+                    b = s * upperGamma s x+                    c = x ** s * exp (-x)+                 in all isSane [a,b,c] ==>+                    let ?eps = 1024*epsi+                     in a ~= b+c || a-b ~= c || a-c ~= b+            , testProperty "x = 0" $ \s ->+                let a = upperGamma s 0+                    b = gamma s+                 in all isSane [a,b] ==>+                    let ?eps = 512*epsi+                     in a ~= b+            , testProperty "s = 1, x > 0" $ \(Positive x) ->+                let ?eps = 256*epsi+                 in x /= 0 ==> upperGamma 1 x ~= exp (-x)+--            , testProperty "s = 1, x < 0" $ \(Positive x) ->+--                let ?eps = 256*eps+--                 in x /= 0 ==> upperGamma 1 (-x) ~= exp x+            , testProperty "s = 0.5" $ \(Positive x) ->+                let ?eps = 128 * epsi+                 in upperGamma 0.5 (x*x) / sqrt pi + erf x ~= 1+            ]+        , testGroup "lnUpperGamma"+            [ testProperty "s = 1, x > 0" $ \(Positive x) ->+                let a = lnUpperGamma 1 x+                 in let ?eps = 1024 * epsi+                     in isSane a ==> a ~= (-x)+--            , testProperty "s = 1, x < 0" $ \(Positive x) ->+--                let a = lnUpperGamma 1 (-x)+--                 in let ?eps = 1024 * eps+--                     in isSane a ==> a ~= x+            , testProperty "agrees with upperGamma" $ \s (NonNegative x) ->+                let a = lnUpperGamma s x;   a' = log b+                    b = upperGamma s x;     b' = exp a+                    +                 in all isSane [a,a',b,b'] ==> +                    let ?eps = 128*epsi+                     in a ~= a' || b ~= b'+            ]+        , testGroup "lowerGamma"+            [ testProperty "increment s" $ \s (NonNegative x) ->+                let a = lowerGamma (s+1) x+                    b = s * lowerGamma s x+                    c = x ** s * exp (-x)+                 in not (s==0 || x==0) && all isSane [a,b,c] ==>+                    let ?eps = 1024*epsi+                     in a ~= b-c || a+c ~= b || b-a ~= c+            , testProperty "s = 0.5" $ \(NonNegative x) ->+                let ?eps = 16 * epsi+                 in lowerGamma 0.5 (x*x) / sqrt pi + erfc x ~= 1+            , testProperty "s = 1, x > 0" $ \(Positive x) ->+                let ?eps = 256*epsi+                 in x /= 0 ==> lowerGamma 1 x + exp (-x) ~= 1+--            , testProperty "s = 1, x < 0" $ \(Positive x) ->+--                let ?eps = 256*eps+--                 in x /= 0 ==> lowerGamma 1 (-x) + exp x ~= 1+            ]+        ]

@@ -0,0 +1,49 @@+{-# LANGUAGE ForeignFunctionInterface #-}+{-# OPTIONS_GHC -fno-warn-type-defaults #-}++module Math.Reference+    ( err+    , tgamma+    , lgamma+    ) where++import Data.Number.BigFloat (BigFloat, Prec50)+import Math.Gamma.Lanczos   (gammaLanczos, reflect)++-- type Ref = BigFloat (PrecPlus20 (PrecPlus20 (PrecPlus20 Prec50)))+type Ref = BigFloat Prec50++refGamma' :: Ref -> Ref+refGamma' x+    | x > 0 = gammaLanczos g cs x+    | otherwise = case properFraction x of+        (n,0) -> if n > 0 then fromInteger (product [1..n-1]) else 0/0+        _     -> reflect (gammaLanczos g cs) x++g::Fractional a => a+cs::Fractional a => [a]++-- g = 4*pi+-- cs = [9.99999999999999999999999996859805857094649110255858638007521499626667910350917912841651754219621023562696911332e-1,2.41021935915031575671085869856051212251086781069965411063137445014791649207976600743220302513537180075741255296e5,-9.39683442214927629464202549615391497731809604507151681446997378587954930487740163463037083316646661765306020767e5,1.50949462278130416707216346338262717752378737207648778592042488727849492084941928909939218772224986523524827271e6,-1.29196480388800562971760352435290649542587719650708873592882554634240911710117984734995323885290793478347646330e6,6.36694128920477129400076563202799705143678072113371399208674087118282528558129436201334367030862320422658246991e5,-1.82471418180570786650926611167028035492537813853558673501139223570434908719610942315092027198546110180213204434e5,2.93229025802547026821102620826974735041884874040373833440904730504747656548457832570174326117298554369079127589e4,-2.42104856309484391962703141393650395582185463486998318684467564163688570420245166352198654798010724012815552809e3,8.70263976705814825513391156852897437425645148258979915072545297700880110517285556763650793819190847901635114815e1,-9.90350651051541084034863017450573967496724797720902585856207732922619895739215627256300090899964209233833558059e-1,1.77111462470160503721207454939884167682509869539871034455375282246563600388989404816104050917245493955854296122e-3,-6.07979338168558073081122926060936198873479110416640952420808116982632249326161497364670524879511535821535057852e-8,2.41177003165992726532912316132624759434987463273639173113847102284810676871579183303359325047268888367893750488e-11,-3.29368781819147402947052101243060083111222984364433304208981956498545487091815683441661638260978104468249346673e-11,2.24453115729169763203131076060244943071934132012991691257619864969664661641236008681074860160604017290650304689e-11,-9.42100302390236927123870056289010276539140677886848718235340783874188089528640301130431598050200664465171227938e-12,2.61294791407497062121812006700436259418110616015207527385647574707939568053359405099396866356159720169068766212e-12,-4.32715100268426644428894578150555672315616608786268711314025106092325719775657429204891211631229932892886729036e-13,2.66596243373922858727920730746965698819332806551704556636516706915515953154082281969632365538824477332141826858e-14,1.77321356734488300746927561707710031581163520443393326611478105950839047661826135657152610805311480453596020540e-15]++-- g = 5*pi+-- cs = [9.99999999999999999999999999993886292491620692316160598686599867606051722036264024489964769505713926399764566018e-1,6.26152503787335403652715374769447739961891146555890980978381175134572365593726917595153339118464456879028200299e6,-3.16335440512758530015177280043484225515804122481774619669315038580009080307294118319420059535265919967069772033e7,6.88835095268754493780533335352036097825929449727064875766761598755949365626134351595870710713568770641108217825e7,-8.46956580333456325307500093500724512731750456929320232801425203134598715275876674318720648707327950233819645139e7,6.47329754881842199772339353831484565977006032033099687741144470787747326331686889977791028752629813425598481151e7,-3.19268433191141469791874286821381470267094549847221319695694408120309390209523512867822254931389981528998092501e7,1.02237086392200672155858881540024824965037923176630968557758516768451829832835513494013895382528680690315693015e7,-2.08890272303961042498315400960754400018306410875319695466462560199939031674854106506443962138469812371142769643e6,2.61321424050856457123433663461415238822132698699274898698444195581028556256292204118246298184152640660180308848e5,-1.86315094230899060297422858598853271269991701222133662269812151190073738494135161781746541881517543444088517450e4,6.72898907328811639584664793087658734023965461622305086223114830302184472501053920325357705077659100329375061164e2,-1.00943157793388998698906942778955472820868316720216759533786330365527456664479196396376802036991852763615117162e1,4.38168530207585179518852579327326813565547855423827923455725302731661798688344684949123570662576923314956409652e-2,-2.56700022740805144017893109688146755394089853324102182680986568703519793504112567874052887016745671988158860509e-5,2.36957641603015216306634255574107168474578275256678299025812704092568833214425267053473328220567647312858972842e-10,1.15694083038124833573384461491621092697227999368580657280434988874851263898593841144273580910524072051004097355e-13,-5.65316256384788292348472745512356223934162359491454541731051455366797335011238891664213587003781560857089287039e-14,8.49474213006453212163254972956476450659399535047675685153282680706022067567488560825213520889992328369290274917e-15,2.28070446151983292921304376186430296779950069002580275483997708509382375273278926321059842672372010645544369814e-15,-1.10200839918603550819509828894654530385656333624024070792512173994075523430080371272123208796268388120390605999e-15,1.38549328695392343672202630014049296547598488298392633928436570519822436737509703285623296516943877681573658856e-16]++-- g = 16+-- cs = [1.00000000000000000000000000000197282505854632603079830840695425875357874385280180091442294078665440768378249098e0,8.46526161663305542149634069187867594722046855412304900854881988813185784117115334278649283694595520258911681345e6,-4.36749138271579683480136186837344589623877333072263981316925882855709622493593700542050638104932540688970449938e7,9.74293930563290544819101625509993544487332726502764966578661696058106791782596806612874880755827781636515874517e7,-1.23204021755002536116556671349951637042719838766566833809791162301292663164161640991555755772951701480576547648e8,9.73252768070985015486610374238589008544120710425762737932789956744946587088886115434931192302762714053106117041e7,-4.99300527228931690175325509502458107530242214409813242610195665492694332478957451089840946485109879221508006660e7,1.67717577881504247993599407438579716255280180251024655266088962746483984030648719508692194681618566416623305775e7,-3.63600580914162233698153144687999173036255278366266794893601829961188765108175483037322789095070584870495716554e6,4.90477280895675364498279606024677805104335312501788362116122862426761656206254084389332939810112070827212068184e5,-3.86134008421465254942234421187656226690536946890041785481282535098425957705374976942206419371390945617575261325e4,1.59799318597802486242335864923532052679307039395883719223156368291075399965515365715900128098360767487863914335e3,-2.92431010362234704149212047766530662168493540855402650510149777639810317548261822822398843227566686564791963281e1,1.74363004063413808172424300916410794128403362491029945418010113517139134404174529869116082145355711405933435336e-1,-1.83889135553342355036509760427653033829857278652806446898154652952072080208456364635472214925432683895177618165e-4,7.33270898586375274415891722892287635561670631889884007144730226480410116396830758987236380519411854625863157502e-9,6.45278000390671865545527420895219728252383858385036218254715732532367519947191407848574538847659484635436688368e-14,-6.94479620536106946656438481553865602897856177541063529072614183743481416853611141761279751208940833601449239078e-14,4.05844190171647299618790214616008990498764985510143852349397107297416322977567534872136697380793297187751912405e-14,-1.31733313666599555356105188113267668954033474011729096786502800127065887481914541186166120448932913511034917752e-14,2.41974657267717455470266193816400294930708183176527991681763060922805701926472008462458376871397921380772169745e-15,-2.02601149216618151544716348082831639958354755920147356248364247315675789908515062799667823982813959353893675328e-16]++g = 16+cs = [1.00000000000000000000000000028760575603001265508442130506014379513956597321465491452116725029663713017318281155302254025482646249861597225426258553986702198310282111660814e0,8.46526161663305542149634069187868312963789113165756291473121826624097304045877402442031648610441688439038240651491328941139813852264197305778456591548908224120094894167246e6,-4.36749138271579683480136186837373556118742078558640362066413127517940638355237243621032816244371668677072024055467961309289136667985789351821077916283681043429310492245048e7,9.74293930563290544819101625512895178907295253257275119766847678373728139529380141034950496249773443820391664483106071950636988772945064231618026141173313776653501642702854e7,-1.23204021755002536116556671362722479624891929606934016862583690470559737742517159607918386688863050322288324639470143518965122467234214549961896872832625644312760192739883e8,9.73252768070985015486610377348101550835843155972616405114822066847313715073349332685175985019283585094322001442327102841230844094173203436533130957487267133091342982903057e7,-4.99300527228931690175325556901821725603229843211023737041708258158778148291264448338140176495149974138608744659105384002613928150272910583070070974922893028621993200088839e7,1.67717577881504247993599895442617443938121576155949723063296052767997809323676003179411440697229440626328003480936049898390416473418451439991990895198019408921657042045389e7,-3.63600580914162233698188824414224911677838005543078337270028990910418595747251549012507094099291095605470928507660671913700304145005691137775901293400972065312587564986433e6,4.90477280895675364500197057167325367117745996005366347188604951009746418243949483785855713959732995794028978290947548233313916565466830346275333882723371675976425249785277e5,-3.86134008421465255019809114269109448642105844599649825606338160546076754535370306095846501745502899782090098798450780912185670412952137486962203827925941120665671522549534e4,1.59799318597802488643745988279123446913929676490559975824444279070965368412618612119375002635982295898117052061965727441889307408959072186532547684742798077234919121659107e3,-2.92431010362235278765051929750170216644701333755426298013832714882892025277018457867985013916936895806752284909319872665280990487229040647950671740315052814101554630827301e1,1.74363004063520607047510940102445356013586745700563922221376454404331547889244236047131697225586656833395848432904781434760496518195945993246362481962863428301094957069334e-1,-1.83889135707484897639922506331631327352704695134176610912010000427967001903568852783866323303740089259923460414793348580157969596224404574122968841894479671002508273424642e-4,7.33288074820458553070627189736492168841920567015366487065477349916549899486012064868307264427889502768168730777757015382371906002557890047275383537043961908550989186406985e-9,-8.14003545649332564829862837358721864946644372854188383148183234349174404580929791883572722401611593795796287442225059328665383671935403293667079588128976652734012111434247e-14,2.30444691867487566543942418193243246136284662851351267620236991760404161776920402164496692134690158656181354762481193169147273885408030319356754409112626454347965677288216e-14,-1.60281647995499836731602960261827820681271025273058699260173363290101149359291286878982917703312074881795295123893833653667997796285463313440572145476822619163620650773297e-15,-1.51961151263436286650314929591024470402675431137010398063941973564941174782047280699515350734586316425940330788036418213975493171611105468630277933898275614760821849251065e-16]++foreign import ccall "math.h lgamma" lgamma :: Double -> Double+foreign import ccall "math.h tgamma" tgamma :: Double -> Double++err :: RealFrac a => (a -> a) -> a -> (a,a)+err approxGamma x = (realToFrac absErr, realToFrac relErr)+    where+        x' = realToFrac x+        +        y = realToFrac (approxGamma x)+        y' = refGamma' x'+        absErr = y - y'+        relErr = absErr / max (abs y) (abs y')

@@ -0,0 +1,13 @@+#!/usr/bin/env runhaskell+module Main where++import qualified Math.GammaTests (tests)+import qualified Math.IncGammaTests (tests)++import Test.Framework (defaultMain, testGroup)++main :: IO ()+main = defaultMain +    [ testGroup "Math.Gamma"            Math.GammaTests.tests+    , testGroup "Math.Gamma.Incomplete" Math.IncGammaTests.tests+    ]