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