diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,9 @@
 # hypergeometric
 
+## 0.1.2.0
+
+  * Add `tcdf`, `chisqcdf`, `fcdf`
+
 ## 0.1.1.0
 
   * Add `Math.SpecialFunction` with `gamma`, `beta` etc.
diff --git a/hypergeometric.cabal b/hypergeometric.cabal
--- a/hypergeometric.cabal
+++ b/hypergeometric.cabal
@@ -1,6 +1,6 @@
 cabal-version:   1.18
 name:            hypergeometric
-version:         0.1.1.0
+version:         0.1.2.0
 license:         AGPL-3
 license-file:    COPYING
 copyright:       Copyright: (c) 2022 Vanessa McHale
@@ -9,7 +9,8 @@
 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
+    Haskell implementation of hypergeometric functions and associated statistical and special functions, viz. erf, normal cdf, incomplete beta, F-distribution cdf, \(\chi^2\)-distribution cdf, t-distrubtion cdf.
+    Also includes Lanczos' approximation of the gamma function.
 
 category:        Math, Statistics
 build-type:      Simple
diff --git a/src/Math/Hypergeometric.hs b/src/Math/Hypergeometric.hs
--- a/src/Math/Hypergeometric.hs
+++ b/src/Math/Hypergeometric.hs
@@ -6,9 +6,6 @@
 
 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)
@@ -16,6 +13,9 @@
 factorial :: Num a => Int -> a
 factorial n = product (fromIntegral <$> [1..n])
 
+-- prop_cdf :: (Double -> Double) -> Double -> Bool
+-- prop_cdf f x = f x <= 1
+
 {-# SPECIALIZE ncdf :: Double -> Double #-}
 -- | CDF of the standard normal \( N(0,1) \)
 ncdf :: (Eq a, Floating a) => a -> a
@@ -29,8 +29,7 @@
 {-# 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!
+-- This iterates until the result stabilizes.
 hypergeometric :: (Eq a, Fractional a)
                => [a] -- ^ \( a_1,\ldots,a_p \)
                -> [a] -- ^ \( b_1,\ldots,b_q \)
@@ -38,9 +37,6 @@
                -> 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
--- a/src/Math/SpecialFunction.hs
+++ b/src/Math/SpecialFunction.hs
@@ -2,26 +2,68 @@
                             , beta
                             , gamma
                             , gammaln
+                            , fcdf
+                            , chisqcdf
+                            , tcdf
                             ) 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 tcdf :: Double -> Double -> Double #-}
+-- | Converges if and only if \(|x| < \sqrt{\nu} \)
+--
+-- @since 0.1.2.0
+tcdf :: (Floating a, Ord a)
+     => a -- ^ \(\nu\) (degrees of freedom)
+     -> a -- ^ \(x\)
+     -> a
+tcdf 𝜈 x = 0.5 + x * gamma (0.5*(𝜈+1)) / (sqrt(pi*𝜈) * gamma(𝜈/2)) * hypergeometric [0.5, 0.5*(𝜈+1)] [1.5] (-x^(2::Int)/𝜈)
 
+-- | @since 0.1.2.0
+{-# SPECIALIZE chisqcdf :: Double -> Double -> Double #-}
+chisqcdf :: (Floating a, Ord a)
+         => a -- ^ \(r\) (degrees of freedom)
+         -> a -- ^ \(\chi^2\)
+         -> a
+chisqcdf r x = incgamma (0.5*r) (0.5*x) / gamma (0.5*r)
+
+{-# SPECIALIZE incgamma :: Double -> Double -> Double #-}
+-- | \(a^{-1}x^a{}_1F_1(a;1+a;-x) \)
+incgamma :: (Floating a, Ord a) => a -> a -> a
+incgamma a x = (1/a) * x ** a * hypergeometric [a] [1+a] (-x)
+-- TODO: writeup?
+--
+-- chisqcdf 10 28 works better this way than w/ e^-x ... x
+--
+-- (1 2 H. _1.1) 1 hangs indefinitely
+
 {-# SPECIALIZE incbeta :: Double -> Double -> Double -> Double #-}
--- | Incomplete beta function.
+-- | Incomplete beta function, \(|z|<1\)
 --
 -- 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)
+incbeta :: (Floating a, Ord a)
         => a -- ^ \(z\)
         -> a -- ^ \(a\)
         -> a -- ^ \(b\)
         -> a
 incbeta z a b = z**a/a * hypergeometric [a,1-b] [a+1] z
 
+{-# SPECIALIZE regbeta :: Double -> Double -> Double -> Double #-}
+-- | \(I(z;a,b) = \displaystyle\frac{B(z;a,b)}{B(a,b)}\)
+regbeta :: (Floating a, Ord a) => a -> a -> a -> a
+regbeta z a b = incbeta z a b / beta a b
+
+{-# SPECIALIZE fcdf :: Double -> Double -> Double -> Double #-}
+-- | @since 0.1.2.0
+fcdf :: (Floating a, Ord a)
+     => a -- ^ \(n\)
+     -> a -- ^ \(m\)
+     -> a -- ^ \(x\)
+     -> a
+fcdf n m x = regbeta (n * x / (m + n * x)) (0.5 * n) (0.5 * m) -- we can use hypergeo because nx/(m+nx) < 1
+
 {-# SPECIALIZE beta :: Double -> Double -> Double #-}
 -- | \(B(x, y) = \displaystyle\frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}\)
 --
@@ -40,11 +82,6 @@
 -- @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))\)