diff --git a/CHANGELOG.md b/CHANGELOG.md
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,5 +1,9 @@
 # hypergeometric
 
+## 0.1.4.0
+
+  * Add `agm`, `completeElliptic`
+
 ## 0.1.3.0
 
   * Add `bessel1`
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.3.0
+version:         0.1.4.0
 license:         AGPL-3
 license-file:    COPYING
 copyright:       Copyright: (c) 2022 Vanessa McHale
diff --git a/src/Math/Hypergeometric.hs b/src/Math/Hypergeometric.hs
--- a/src/Math/Hypergeometric.hs
+++ b/src/Math/Hypergeometric.hs
@@ -13,20 +13,20 @@
 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 #-}
+{-# SPECIALIZE ncdf :: Float -> Float #-}
 -- | CDF of the standard normal \( N(0,1) \)
 ncdf :: (Eq a, Floating a) => a -> a
 ncdf z = (1/2) * (1 + erf (z / sqrt 2))
 
 {-# SPECIALIZE erf :: Double -> Double #-}
+{-# SPECIALIZE erf :: Float -> Float #-}
 -- | [erf](https://mathworld.wolfram.com/Erf.html)
 erf :: (Eq a, Floating a) => a -> a
 erf z = (2 * z * exp (-z^(2::Int)) / sqrt pi) * hypergeometric [1] [3/2] (z^(2::Int))
 
 {-# SPECIALIZE hypergeometric :: [Double] -> [Double] -> Double -> Double #-}
+{-# SPECIALIZE hypergeometric :: [Float] -> [Float] -> Float -> Float #-}
 -- | \( _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.
diff --git a/src/Math/SpecialFunction.hs b/src/Math/SpecialFunction.hs
--- a/src/Math/SpecialFunction.hs
+++ b/src/Math/SpecialFunction.hs
@@ -3,6 +3,8 @@
                             , bessel1
                             , gamma
                             , gammaln
+                            , agm
+                            , completeElliptic
                             , fcdf
                             , chisqcdf
                             , tcdf
@@ -11,6 +13,7 @@
 import           Math.Hypergeometric
 
 {-# SPECIALIZE tcdf :: Double -> Double -> Double #-}
+{-# SPECIALIZE tcdf :: Float -> Float -> Float #-}
 -- | Converges if and only if \(|x| < \sqrt{\nu} \)
 --
 -- @since 0.1.2.0
@@ -20,6 +23,8 @@
      -> 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)/𝜈)
 
+{-# SPECIALIZE bessel1 :: Double -> Double -> Double #-}
+{-# SPECIALIZE bessel1 :: Float -> Float -> Float #-}
 -- | Bessel functions of the first kind, \( J_\alpha(x)\).
 --
 -- @since 0.1.3.0
@@ -29,8 +34,30 @@
         -> a
 bessel1 𝛼 x = ((x/2)**𝛼/gamma(𝛼+1))*hypergeometric [] [𝛼+1] (-(x^(2::Int))/4)
 
+{-# SPECIALIZE completeElliptic :: Double -> Double #-}
+{-# SPECIALIZE completeElliptic :: Float -> Float #-}
+-- | [Complete elliptic integral of the first kind](https://mathworld.wolfram.com/CompleteEllipticIntegraloftheFirstKind.html)
+--
+-- @since 0.1.4.0
+completeElliptic :: (Ord a, Floating a) => a -> a
+completeElliptic k = pi/(2*agm 1 (sqrt (1-k^(2::Int))))
+
+{-# SPECIALIZE agm :: Double -> Double -> Double #-}
+{-# SPECIALIZE agm :: Float -> Float -> Float #-}
+-- | [Arithmetic-geometric mean](https://mathworld.wolfram.com/Arithmetic-GeometricMean.html)
+--
+-- @since 0.1.4.0
+agm :: (Ord a, Floating a) => a -> a -> a
+agm a b =
+    let a' = (a+b)/2
+        b' = sqrt(a*b)
+    in if (a'-b')/b'<1e-15 && (b'-a')/b'<1e-15
+        then a'
+        else agm a' b'
+
 -- | @since 0.1.2.0
 {-# SPECIALIZE chisqcdf :: Double -> Double -> Double #-}
+{-# SPECIALIZE chisqcdf :: Float -> Float -> Float #-}
 chisqcdf :: (Floating a, Ord a)
          => a -- ^ \(r\) (degrees of freedom)
          -> a -- ^ \(\chi^2\)
@@ -38,6 +65,7 @@
 chisqcdf r x = incgamma (0.5*r) (0.5*x) / gamma (0.5*r)
 
 {-# SPECIALIZE incgamma :: Double -> Double -> Double #-}
+{-# SPECIALIZE incgamma :: Float -> Float -> Float #-}
 -- | \(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)
@@ -48,6 +76,7 @@
 -- (1 2 H. _1.1) 1 hangs indefinitely
 
 {-# SPECIALIZE incbeta :: Double -> Double -> Double -> Double #-}
+{-# SPECIALIZE incbeta :: Float -> Float -> Float -> Float #-}
 -- | Incomplete beta function, \(|z|<1\)
 --
 -- Calculated with \(B(z;a,b)=\displaystyle\frac{z^a}{a}{}_2F_1(a, 1-b; a+1; z)\)
@@ -61,11 +90,13 @@
 incbeta z a b = z**a/a * hypergeometric [a,1-b] [a+1] z
 
 {-# SPECIALIZE regbeta :: Double -> Double -> Double -> Double #-}
+{-# SPECIALIZE regbeta :: Float -> Float -> Float -> Float #-}
 -- | \(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 #-}
+{-# SPECIALIZE fcdf :: Float -> Float -> Float -> Float #-}
 -- | @since 0.1.2.0
 fcdf :: (Floating a, Ord a)
      => a -- ^ \(n\)
@@ -75,6 +106,7 @@
 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 #-}
+{-# SPECIALIZE beta :: Float -> Float -> Float #-}
 -- | \(B(x, y) = \displaystyle\frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}\)
 --
 -- This uses 'gammaln' under the hood to extend its domain somewhat.
@@ -84,6 +116,7 @@
 beta x y = exp (betaln x y)
 
 {-# SPECIALIZE betaln :: Double -> Double -> Double #-}
+{-# SPECIALIZE betaln :: Float -> Float -> Float #-}
 betaln :: (Floating a, Ord a) => a -> a -> a
 betaln x y = gammaln x + gammaln y - gammaln (x+y)
 
@@ -94,6 +127,7 @@
 gamma = exp . gammaln
 
 {-# SPECIALIZE gammaln :: Double -> Double #-}
+{-# SPECIALIZE gammaln :: Float -> Float #-}
 -- | \(\text{log} (\Gamma(z))\)
 --
 -- Lanczos approximation.