diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,9 @@
 # hypergeometric
 
+## 0.1.1.0
+
+  * Add `Math.SpecialFunction` with `gamma`, `beta` etc.
+
 ## 0.1.0.0
 
 Initial release
diff --git a/hypergeometric.cabal b/hypergeometric.cabal
--- a/hypergeometric.cabal
+++ b/hypergeometric.cabal
@@ -1,11 +1,12 @@
 cabal-version:   1.18
 name:            hypergeometric
-version:         0.1.0.0
+version:         0.1.1.0
 license:         AGPL-3
 license-file:    COPYING
 copyright:       Copyright: (c) 2022 Vanessa McHale
 maintainer:      vamchale@gmail.com
 author:          Vanessa McHale
+bug-reports:     https://github.com/vmchale/hypergeometric/issues
 synopsis:        Hypergeometric functions
 description:
     Haskell implementation of hypergeometric functions and associated statistical functions, viz. erf, normal cdf
@@ -21,7 +22,10 @@
     location: https://github.com/vmchale/hypergeometric
 
 library
-    exposed-modules:  Math.Hypergeometric
+    exposed-modules:
+        Math.Hypergeometric
+        Math.SpecialFunction
+
     hs-source-dirs:   src
     default-language: Haskell2010
     ghc-options:      -Wall
diff --git a/src/Math/Hypergeometric.hs b/src/Math/Hypergeometric.hs
--- a/src/Math/Hypergeometric.hs
+++ b/src/Math/Hypergeometric.hs
@@ -6,6 +6,9 @@
 
 import           Data.Functor ((<$>))
 
+-- choose :: Integral a => a -> a -> a
+-- choose n k = product [(n-k+1) .. n] `quot` factorial (fromIntegral k)
+
 risingFactorial :: Num a => a -> Int -> a
 risingFactorial _ 0 = 1
 risingFactorial a n = (a + fromIntegral n - 1) * risingFactorial a (n-1)
@@ -25,6 +28,9 @@
 
 {-# SPECIALIZE hypergeometric :: [Double] -> [Double] -> Double -> Double #-}
 -- | \( _pF_q(a_1,\ldots,a_p;b_1,\ldots,b_q;z) = \displaystyle\sum_{n=0}^\infty\frac{(a_1)_n\cdots(a_p)_n}{(b_1)_b\cdots(b_q)_n}\frac{z^n}{n!} \)
+--
+-- This iterates until the result stabilizes, so don't use it on
+-- arbitrary-precision types!
 hypergeometric :: (Eq a, Fractional a)
                => [a] -- ^ \( a_1,\ldots,a_p \)
                -> [a] -- ^ \( b_1,\ldots,b_q \)
@@ -32,6 +38,9 @@
                -> a
 hypergeometric as bs z = sumUntilEq
     [ (product (fmap (`risingFactorial` n) as) / product (fmap (`risingFactorial` n) bs)) * (z ^ n) / factorial n | n <- [0..] ]
+    -- [ exp (nth n) | n <- [0..] ]
+    -- where nth n = sum (fmap (log . (`risingFactorial` n)) as) - sum (fmap (log . (`risingFactorial` n)) bs) + fromIntegral n * log z - log (factorial n)
+    -- TODO: Revisit the exponential approach using complex numbers?
 
 sumUntilEq :: (Eq a, Num a) => [a] -> a
 sumUntilEq = sumUntilEqLoop 0
diff --git a/src/Math/SpecialFunction.hs b/src/Math/SpecialFunction.hs
new file mode 100644
--- /dev/null
+++ b/src/Math/SpecialFunction.hs
@@ -0,0 +1,81 @@
+module Math.SpecialFunction ( incbeta
+                            , beta
+                            , gamma
+                            , gammaln
+                            ) where
+
+import           Math.Hypergeometric
+
+-- prop_betamatch :: Double -> Double -> Bool
+-- prop_betamatch x y = x <= 0 || y <= 0 || abs (beta x y - incbeta 1 x y) < 1e-15
+
+{-# SPECIALIZE incbeta :: Double -> Double -> Double -> Double #-}
+-- | Incomplete beta function.
+--
+-- Calculated with \(B(z;a,b)=\displaystyle\frac{z^a}{a}{}_2F_1(a, 1-b; a+1; z)\)
+--
+-- @since 0.1.1.0
+incbeta :: (Floating a, Eq a)
+        => a -- ^ \(z\)
+        -> a -- ^ \(a\)
+        -> a -- ^ \(b\)
+        -> a
+incbeta z a b = z**a/a * hypergeometric [a,1-b] [a+1] z
+
+{-# SPECIALIZE beta :: Double -> Double -> Double #-}
+-- | \(B(x, y) = \displaystyle\frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}\)
+--
+-- This uses 'gammaln' under the hood to extend its domain somewhat.
+--
+-- @since 0.1.1.0
+beta :: (Floating a, Ord a) => a -> a -> a
+beta x y = exp (betaln x y)
+
+{-# SPECIALIZE betaln :: Double -> Double -> Double #-}
+betaln :: (Floating a, Ord a) => a -> a -> a
+betaln x y = gammaln x + gammaln y - gammaln (x+y)
+
+-- | \(\Gamma(z)\)
+--
+-- @since 0.1.1.0
+gamma :: (Floating a, Ord a) => a -> a
+gamma = exp . gammaln
+
+-- gamma from beta:
+-- Γ(z)Γ(1-z) = 𝜋/sin(𝜋z)
+--
+-- THENCE, B(z,1-z)=Γ(z)Γ(1-z)/Γ(1)=...
+
+{-# SPECIALIZE gammaln :: Double -> Double #-}
+-- | \(\text{log} (\Gamma(z))\)
+--
+-- Lanczos approximation.
+-- This is exactly the approach described in Press, William H. et al. /Numerical Recipes/, 3rd ed., extended to work on negative real numbers.
+--
+-- @since 0.1.1.0
+gammaln :: (Floating a, Ord a)
+        => a -- ^ \( z \)
+        -> a
+gammaln 0 = -log 0
+gammaln z | z >= 0.5 = (z' + 1/2) * log (z' + 𝛾 + 1/2) - (z' + 1/2 + 𝛾) + log (sqrt (2*pi) * (c0 + sum series))
+    where series = zipWith (\c x -> c / (z' + fromIntegral x)) coeff [(1::Int)..]
+          c0 = 0.999999999999997092
+          -- constants from Numerical Recipes
+          coeff = [ 57.1562356658629235
+                  , -59.5979603554754912
+                  , 14.1360979747417471
+                  , -0.491913816097620199
+                  , 0.339946499848118887e-4
+                  , 0.465236289270485756e-4
+                  , -0.983744753048795646e-4
+                  , 0.158088703224912494e-3
+                  , -0.210264441724104883e-3
+                  , 0.217439618115212643e-3
+                  , -0.164318106536763890e-3
+                  , 0.844182239838527433e-4
+                  , -0.261908384015814087e-4
+                  , 0.368991826595316234e-5
+                  ]
+          𝛾 = 607/128
+          z' = z-1
+gammaln z = log pi - log (sin (pi * z)) - gammaln (1 - z)