diff --git a/doc/source/acb.rst b/doc/source/acb.rst
index d82fe654..2d12998b 100644
--- a/doc/source/acb.rst
+++ b/doc/source/acb.rst
@@ -812,6 +812,8 @@ Polylogarithms
 Arithmetic-geometric mean
 -------------------------------------------------------------------------------
 
+See :ref:`algorithms_agm` for implementation details.
+
 .. function:: void acb_agm1(acb_t m, const acb_t z, slong prec)
 
     Sets *m* to the arithmetic-geometric mean `M(z) = \operatorname{agm}(1,z)`,
diff --git a/doc/source/agm.rst b/doc/source/agm.rst
new file mode 100644
index 00000000..ef850d70
--- /dev/null
+++ b/doc/source/agm.rst
@@ -0,0 +1,94 @@
+.. _algorithms_agm:
+
+Algorithms for the arithmetic-geometric mean
+===============================================================================
+
+With complex variables, it is convenient to work with the univariate
+function `M(z) = \operatorname{agm}(1,z)`. The general case is given by
+`\operatorname{agm}(a,b) = a M(1,b/a)`.
+
+Functional equation
+------------------------------------------------------------------------------
+
+If the real part of *z* initially is not completely nonnegative, we
+apply the functional equation `M(z) = (z+1) M(u) / 2`
+where `u = \sqrt{z} / (z+1)`.
+
+Note that *u* has nonnegative real part, absent rounding error.
+It is not a problem for correctness if rounding makes the interval
+contain negative points, as this just inflates the final result.
+
+For the derivative, the functional equation becomes
+`M'(z) = [M(u) - (z-1) M'(u) / ((1+z) \sqrt{z})] / 2`.
+
+AGM iteration
+------------------------------------------------------------------------------
+
+Once *z* is in the right half plane, we can apply the AGM iteration
+(`2a_{n+1} = a_n + b_n, b_{n+1}^2 = a_n b_n`) directly.
+The correct square root is given by `\sqrt{a} \sqrt{b}`,
+which is computed as `\sqrt{ab}, i \sqrt{-ab}, -i \sqrt{-ab}, \sqrt{a} \sqrt{b}`
+respectively if both *a* and *b* have positive real part, nonnegative
+imaginary part, nonpositive imaginary part, or otherwise.
+
+It is shown in [Dup2006]_, p. 87 that, for *z* with nonnegative real part,
+`|M(z) - a_n| \le |a_n - b_n|`.
+
+A small improvement would be to switch to a series
+expansion for the last one or two steps.
+
+First derivative
+------------------------------------------------------------------------------
+
+Assuming that *z* is exact and that `|\arg(z)| \le 3 \pi / 4`,
+we compute `(M(z), M'(z))` simultaneously using a finite difference.
+
+The basic inequality we need is `|M(z)| \le \max(1, |z|)`, which is
+an immediate consequence of the AGM iteration.
+
+By Cauchy's integral formula, `|M^{(k)}(z) / k!| \le C D^k` where
+`C = \max(1, |z| + r)` and `D = 1/r`, for any `0 < r < |z|` (we
+choose *r* to be of the order `|z| / 4`). Taylor expansion now gives
+
+.. math ::
+
+    \left|\frac{M(z+h) - M(z)}{h} - M'(z)\right| \le \frac{C D^2 h}{1 - D h}
+
+assuming that *h* is chosen so that it satisfies `h D < 1`.
+
+The forward finite difference requires two function evaluations
+at doubled precision. It would be more efficient to use a central difference
+(this could be implemented in the future).
+
+When *z* is not exact, we evaluate at the midpoint as above
+and bound the propagated error using derivatives.
+Again by Cauchy's integral formula, we have
+
+.. math ::
+
+    |M'(z+\varepsilon)| \le \frac{\max(1, |z|+|\varepsilon|+r)}{r}
+
+    |M''(z+\varepsilon)| \le \frac{2 \max(1, |z|+|\varepsilon|+r)}{r^2}
+
+assuming that the circle centered on *z* with radius `|\varepsilon| + r`
+does not cross the negative half axis. We choose *r* of order `|z| / 2`
+and verify that all assumptions hold.
+
+Higher derivatives
+-------------------------------------------------------------------------------
+
+The function `W(z) = 1 / M(z)` is D-finite. The coefficients of
+`W(z+x) = \sum_{k=0}^{\infty} c_k x^k` satisfy
+
+.. math ::
+
+    -(k+2)(k+3) z (z^2-1) c_{k+3} = (k+2)^2 (3z^2-1) c_{k+2} + (3k(k+3)+7)z c_{k+1} + (k+1)^2 c_{k}
+
+in general, and
+
+.. math ::
+
+    -(k+2)^2 c_{k+2} = (3k(k+3)+7) c_{k+1} + (k+1)^2 c_{k}
+
+when `z = 1`.
+
diff --git a/doc/source/credits.rst b/doc/source/credits.rst
index 3a4c3874..532735eb 100644
--- a/doc/source/credits.rst
+++ b/doc/source/credits.rst
@@ -118,6 +118,8 @@ Bibliography
 
 .. [CP2005] \R. Crandall and C. Pomerance, *Prime Numbers: A Computational Perspective*, second edition, Springer (2005).
 
+.. [Dup2006] \R. Dupont. "Moyenne arithmético-géométrique, suites de Borchardt et applications." These de doctorat, École polytechnique, Palaiseau (2006). http://http://www.lix.polytechnique.fr/Labo/Regis.Dupont/these_soutenance.pdf
+
 .. [EM2004] \O. Espinosa and V. Moll, "A generalized polygamma function", Integral Transforms and Special Functions (2004), 101-115.
 
 .. [Fil1992] \S. Fillebrown, "Faster Computation of Bernoulli Numbers", Journal of Algorithms 13 (1992) 431-445
diff --git a/doc/source/index.rst b/doc/source/index.rst
index a25f377a..b4a23194 100644
--- a/doc/source/index.rst
+++ b/doc/source/index.rst
@@ -161,6 +161,7 @@ lengthy to reproduce in the documentation for each module.
    gamma.rst
    polylogarithms.rst
    hypergeometric.rst
+   agm.rst
 
 Credits and references
 ::::::::::::::::::::::::