diff --git a/CHANGES.md b/CHANGES.md new file mode 100644 --- /dev/null +++ b/CHANGES.md @@ -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]. diff --git a/README.md b/README.md new file mode 100644 --- /dev/null +++ b/README.md @@ -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). diff --git a/Setup.lhs b/Setup.lhs --- a/Setup.lhs +++ b/Setup.lhs @@ -1,5 +1,6 @@ #!/usr/bin/env runhaskell > import Distribution.Simple + > main = defaultMain diff --git a/gamma.cabal b/gamma.cabal --- a/gamma.cabal +++ b/gamma.cabal @@ -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 -maintainer: James Cook +author: James Cook +maintainer: James Cook , + Alexandre Rodrigues Baldé 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 diff --git a/src/Math/Factorial.hs b/src/Math/Factorial.hs --- a/src/Math/Factorial.hs +++ b/src/Math/Factorial.hs @@ -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 diff --git a/src/Math/Gamma.hs b/src/Math/Gamma.hs --- a/src/Math/Gamma.hs +++ b/src/Math/Gamma.hs @@ -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 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 @@ -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 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 @@ -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 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 @@ -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 diff --git a/test-suite/Math/GammaTests.hs b/test-suite/Math/GammaTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/GammaTests.hs @@ -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 + ] diff --git a/test-suite/Math/IncGammaTests.hs b/test-suite/Math/IncGammaTests.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/IncGammaTests.hs @@ -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 + ] + ] diff --git a/test-suite/Math/Reference.hs b/test-suite/Math/Reference.hs new file mode 100644 --- /dev/null +++ b/test-suite/Math/Reference.hs @@ -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') diff --git a/test-suite/Test.hs b/test-suite/Test.hs new file mode 100644 --- /dev/null +++ b/test-suite/Test.hs @@ -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 + ]