diff --git a/doc/source/credits.rst b/doc/source/credits.rst
index e2db7923..dc70edae 100644
--- a/doc/source/credits.rst
+++ b/doc/source/credits.rst
@@ -68,6 +68,8 @@ Bibliography
 
 .. [NIST2012] National Institute of Standards and Technology, *Digital Library of Mathematical Functions* (2012), http://dlmf.nist.gov/
 
+.. [Olv1997] \F. Olver, *Asymptotics and special functions*, AKP Classics, AK Peters Ltd., Wellesley, MA, 1997. Reprint of the 1974 original.
+
 .. [PS1973] \M. S. Paterson and L. J. Stockmeyer, "On the number of nonscalar multiplications necessary to evaluate polynomials", SIAM J. Comput (1973)
 
 .. [Smi2001] \D. M. Smith, "Algorithm: Fortran 90 Software for Floating-Point Multiple Precision Arithmetic, Gamma and Related Functions", Transactions on Mathematical Software 27 (2001) 377-387, http://myweb.lmu.edu/dmsmith/toms2001.pdf
diff --git a/doc/source/gamma.rst b/doc/source/gamma.rst
index da704e07..1c6dd3ac 100644
--- a/doc/source/gamma.rst
+++ b/doc/source/gamma.rst
@@ -113,6 +113,50 @@ Evaluation using the Stirling series
     differentiated term and the exponent `2n` changes to `2n+1`
     (section 5.11 in [NIST2012]_).
 
+.. function :: void gamma_stirling_bound_phase(fmpr_t bound, const fmpcb_t z, long prec)
+
+    Sets *bound* to an upper bound for the phase factor
+    `b = 1/\cos(\operatorname{arg}(z)/2)` which appears in the error bound
+    for the Stirling series. By trigonometric identities, assuming
+    that `z = x+yi`, we have `b = \sqrt{1 + t^2}` where
+
+    .. math ::
+
+        t = \frac{y}{\sqrt{x^2 + y^2} + x} = \frac{\sqrt{x^2 + y^2} - x}{y}
+
+    We bound `x` and `y` such that `|\operatorname{arg}(x+yi)|`
+    is maximized, and then evaluate `t` with the choice of square root
+    expression that avoids cancellation, using directional rounding throughout.
+
+.. function :: void gamma_stirling_bound_fmprb(fmpr_struct * err, const fmprb_t z, long k0, long knum, long n)
+
+.. function :: void gamma_stirling_bound_fmpcb(fmpr_struct * err, const fmpcb_t z, long k0, long knum, long n)
+
+    Computes bounds for the truncation error in the Stirling series
+    when summed up to term `n - 1` inclusive. An exact expression for the
+    truncation error is given (see [Olv1997]_ pp. 293-295) by
+
+    .. math ::
+
+        R_n(z) = \int_0^{\infty} \frac{B_{2n} - {\tilde B}_{2n}(x)}{2n(x+z)^{2n}} dt.
+
+    We optionally evaluate the bound for several terms in the
+    Taylor series: considering `R_n(z+t) \in \mathbb{C}[[t]]`, we
+    compute bounds for the coefficient of `t^k` for *knum* consecutive
+    values of *k* starting with *k0*.
+    Using the fact that the numerator of the integrand is bounded in
+    absolute value by `2 |B_{2n}|`, and using the bound for `|x+z|`
+    given by [Olv1997]_, we obtain
+
+    .. math ::
+
+        |[t^k] R_n(z+t)| \le 2 |B_{2n}|
+            \frac{\Gamma(2n+k-1)}{\Gamma(k+1) \Gamma(2n+1)}
+            \; |z| \; (b / |z|)^{2n+k}
+
+    where `b` is the phase factor implemented by
+    :func:`gamma_stirling_bound_phase`.
+
 
 Rising factorials
 --------------------------------------------------------------------------------
diff --git a/fmpcb.h b/fmpcb.h
index da27cf6b..0a9bcd64 100644
--- a/fmpcb.h
+++ b/fmpcb.h
@@ -241,35 +241,57 @@ fmpcb_set_round_fmprb(fmpcb_t z, const fmprb_t x, long prec)
 static __inline__ void
 fmpcb_get_abs_ubound_fmpr(fmpr_t u, const fmpcb_t z, long prec)
 {
-    fmpr_t v;
-    fmpr_init(v);
+    if (fmprb_is_zero(fmpcb_imagref(z)))
+    {
+        fmprb_get_abs_ubound_fmpr(u, fmpcb_realref(z), prec);
+    }
+    else if (fmprb_is_zero(fmpcb_realref(z)))
+    {
+        fmprb_get_abs_ubound_fmpr(u, fmpcb_imagref(z), prec);
+    }
+    else
+    {
+        fmpr_t v;
+        fmpr_init(v);
 
-    fmprb_get_abs_ubound_fmpr(u, fmpcb_realref(z), prec);
-    fmprb_get_abs_ubound_fmpr(v, fmpcb_imagref(z), prec);
+        fmprb_get_abs_ubound_fmpr(u, fmpcb_realref(z), prec);
+        fmprb_get_abs_ubound_fmpr(v, fmpcb_imagref(z), prec);
 
-    fmpr_mul(u, u, u, prec, FMPR_RND_UP);
-    fmpr_mul(v, v, v, prec, FMPR_RND_UP);
-    fmpr_add(u, u, v, prec, FMPR_RND_UP);
-    fmpr_sqrt(u, u, prec, FMPR_RND_UP);
+        fmpr_mul(u, u, u, prec, FMPR_RND_UP);
+        fmpr_mul(v, v, v, prec, FMPR_RND_UP);
+        fmpr_add(u, u, v, prec, FMPR_RND_UP);
+        fmpr_sqrt(u, u, prec, FMPR_RND_UP);
 
-    fmpr_clear(v);
+        fmpr_clear(v);
+    }
 }
 
 static __inline__ void
 fmpcb_get_abs_lbound_fmpr(fmpr_t u, const fmpcb_t z, long prec)
 {
-    fmpr_t v;
-    fmpr_init(v);
+    if (fmprb_is_zero(fmpcb_imagref(z)))
+    {
+        fmprb_get_abs_lbound_fmpr(u, fmpcb_realref(z), prec);
+    }
+    else if (fmprb_is_zero(fmpcb_realref(z)))
+    {
+        fmprb_get_abs_lbound_fmpr(u, fmpcb_imagref(z), prec);
+    }
+    else
+    {
+        fmpr_t v;
+        fmpr_init(v);
 
-    fmprb_get_abs_lbound_fmpr(u, fmpcb_realref(z), prec);
-    fmprb_get_abs_lbound_fmpr(v, fmpcb_imagref(z), prec);
+        fmprb_get_abs_lbound_fmpr(u, fmpcb_realref(z), prec);
+        fmprb_get_abs_lbound_fmpr(v, fmpcb_imagref(z), prec);
 
-    fmpr_mul(u, u, u, prec, FMPR_RND_DOWN);
-    fmpr_mul(v, v, v, prec, FMPR_RND_DOWN);
-    fmpr_add(u, u, v, prec, FMPR_RND_DOWN);
-    fmpr_sqrt(u, u, prec, FMPR_RND_DOWN);
+        fmpr_mul(u, u, u, prec, FMPR_RND_DOWN);
+        fmpr_mul(v, v, v, prec, FMPR_RND_DOWN);
+        fmpr_add(u, u, v, prec, FMPR_RND_DOWN);
+        fmpr_sqrt(u, u, prec, FMPR_RND_DOWN);
 
-    fmpr_clear(v);
+        fmpr_clear(v);
+    }
 }
 
 static __inline__ void
diff --git a/gamma.h b/gamma.h
index 20a34811..6fac89fe 100644
--- a/gamma.h
+++ b/gamma.h
@@ -71,6 +71,10 @@ void gamma_taylor_fmprb(fmprb_t y, const fmprb_t x, long prec);
 void gamma_stirling_choose_param_fmprb(int * reflect, long * r, long * n, const fmprb_t z, int use_reflect, int digamma, long prec);
 void gamma_stirling_choose_param_fmpcb(int * reflect, long * r, long * n, const fmpcb_t z, int use_reflect, int digamma, long prec);
 
+void gamma_stirling_bound_phase(fmpr_t bound, const fmpcb_t z, long prec);
+void gamma_stirling_bound_fmprb(fmpr_struct * err, const fmprb_t x, long k0, long knum, long n);
+void gamma_stirling_bound_fmpcb(fmpr_struct * err, const fmpcb_t z, long k0, long knum, long n);
+
 void gamma_stirling_coeff(fmprb_t b, ulong k, int digamma, long prec);
 void gamma_stirling_eval_series_fmprb(fmprb_t s, const fmprb_t z, long nterms, int digamma, long prec);
 void gamma_stirling_eval_series_fmpcb(fmpcb_t s, const fmpcb_t z, long nterms, int digamma, long prec);
diff --git a/gamma/stirling_bound.c b/gamma/stirling_bound.c
new file mode 100644
index 00000000..a3a240a7
--- /dev/null
+++ b/gamma/stirling_bound.c
@@ -0,0 +1,108 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2012 Fredrik Johansson
+
+******************************************************************************/
+
+#include <math.h>
+#include "gamma.h"
+#include "bernoulli.h"
+
+void fmpr_gamma_ui_lbound(fmpr_t x, ulong n, long prec);
+
+void fmpr_gamma_ui_ubound(fmpr_t x, ulong n, long prec);
+
+/*
+  2 |B_{2n}| G(2n+k-1) / (G(k+1) G(2n+1)) |z| (T |z|^{-1})^(2n+k)
+*/
+void
+gamma_stirling_bound_fmpcb(fmpr_struct * err, const fmpcb_t z, long k0, long knum, long n)
+{
+    fmpr_t c, t, u, v;
+    long i, k, prec = FMPRB_RAD_PREC;
+
+    if (fmprb_contains_zero(fmpcb_imagref(z)) &&
+        fmprb_contains_nonpositive(fmpcb_realref(z)))
+    {
+        for (i = 0; i < knum; i++)
+            fmpr_pos_inf(err + i);
+        return;
+    }
+
+    fmpr_init(c);
+    fmpr_init(t);
+    fmpr_init(u);
+    fmpr_init(v);
+
+    /* t = lower bound for |z| */
+    fmpcb_get_abs_lbound_fmpr(t, z, prec);
+    fmpcb_get_abs_ubound_fmpr(v, z, prec);
+
+    /* c = upper bound for 1/(cos(arg(z)/2) |z|) */
+    gamma_stirling_bound_phase(c, z, prec);
+    fmpr_div(c, c, t, prec, FMPR_RND_UP);
+
+    /* numerator: 2 B_{2n} gamma(2n+k-1) |z| */
+    BERNOULLI_ENSURE_CACHED(2 * n);
+    fmpr_set_round_fmpz(err, fmpq_numref(bernoulli_cache + 2 * n), prec, FMPR_RND_UP);
+    fmpr_abs(err, err);
+    fmpr_div_fmpz(err, err, fmpq_denref(bernoulli_cache + 2 * n), prec, FMPR_RND_UP);
+    fmpr_mul_2exp_si(err, err, 1);
+    fmpr_gamma_ui_ubound(u, 2 * n + k0 - 1, prec);
+    fmpr_mul(err, err, u, prec, FMPR_RND_UP);
+    fmpr_mul(err, err, v, prec, FMPR_RND_UP);
+
+    /* denominator gamma(k+1) gamma(2n+1) */
+    fmpr_gamma_ui_lbound(t, 2 * n + 1, prec);
+    fmpr_gamma_ui_lbound(u, k0 + 1, prec);
+    fmpr_mul(t, t, u, prec, FMPR_RND_DOWN);
+    fmpr_div(err, err, t, prec, FMPR_RND_UP);
+
+    /* multiply by c^(2n+k) */
+    fmpr_pow_sloppy_ui(t, c, 2 * n + k0, prec, FMPR_RND_UP);
+    fmpr_mul(err, err, t, prec, FMPR_RND_UP);
+
+    for (i = 1; i < knum; i++)
+    {
+        /* recurrence factor: c * (2n+k-2) / k */
+        k = k0 + i;
+        fmpr_mul(err + i, err + i - 1, c, prec, FMPR_RND_UP);
+        fmpr_mul_ui(err + i, err + i, 2 * n + k - 2, prec, FMPR_RND_UP);
+        fmpr_div_ui(err + i, err + i, k, prec, FMPR_RND_UP);
+    }
+
+    fmpr_clear(c);
+    fmpr_clear(t);
+    fmpr_clear(u);
+    fmpr_clear(v);
+}
+
+void
+gamma_stirling_bound_fmprb(fmpr_struct * err, const fmprb_t x, long k0, long knum, long n)
+{
+    fmpcb_t z;
+    fmpcb_init(z);
+    fmpcb_set_fmprb(z, x);
+    gamma_stirling_bound_fmpcb(err, z, k0, knum, n);
+    fmpcb_clear(z);
+}
+
diff --git a/gamma/stirling_bound_phase.c b/gamma/stirling_bound_phase.c
new file mode 100644
index 00000000..b5f27e26
--- /dev/null
+++ b/gamma/stirling_bound_phase.c
@@ -0,0 +1,93 @@
+/*=============================================================================
+
+    This file is part of ARB.
+
+    ARB is free software; you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation; either version 2 of the License, or
+    (at your option) any later version.
+
+    ARB is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with ARB; if not, write to the Free Software
+    Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+=============================================================================*/
+/******************************************************************************
+
+    Copyright (C) 2013 Fredrik Johansson
+
+******************************************************************************/
+
+#include "gamma.h"
+
+void
+gamma_stirling_bound_phase(fmpr_t bound, const fmpcb_t z, long prec)
+{
+    fmpr_t x, y, t, u;
+    int xsign;
+
+    fmpr_init(x);
+    fmpr_init(y);
+    fmpr_init(t);
+    fmpr_init(u);
+
+    /* first compute x, y such that |arg(z)| <= arg(x+yi) */
+
+    /* argument increases with smaller real parts */
+    fmpr_sub(x, fmprb_midref(fmpcb_realref(z)),
+        fmprb_radref(fmpcb_realref(z)), prec, FMPR_RND_FLOOR);
+
+    xsign = fmpr_sgn(x);
+
+    if (xsign >= 0)  /* argument increases away from the real axis */
+        fmprb_get_abs_ubound_fmpr(y, fmpcb_imagref(z), prec);
+    else  /* argument increases closer to the real axis */
+        fmprb_get_abs_lbound_fmpr(y, fmpcb_imagref(z), prec);
+
+    if (fmpr_is_zero(y))
+    {
+        if (xsign > 0)
+            fmpr_one(bound);
+        else
+            fmpr_pos_inf(bound);
+    }
+    else
+    {
+        if (xsign >= 0)
+        {
+            /* compute upper bound for t = y / (sqrt(x^2 + y^2) + x) */
+            fmpr_mul(t, x, x, prec, FMPR_RND_DOWN);
+            fmpr_mul(u, y, y, prec, FMPR_RND_DOWN);
+            fmpr_add(t, t, u, prec, FMPR_RND_DOWN);
+            fmpr_sqrt(t, t, prec, FMPR_RND_DOWN);
+            fmpr_add(t, t, x, prec, FMPR_RND_DOWN);
+            fmpr_div(t, y, t, prec, FMPR_RND_UP);
+        }
+        else
+        {
+            /* compute upper bound for t = (sqrt(x^2 + y^2) - x) / y */
+            fmpr_mul(t, x, x, prec, FMPR_RND_UP);
+            fmpr_mul(u, y, y, prec, FMPR_RND_UP);
+            fmpr_add(t, t, u, prec, FMPR_RND_UP);
+            fmpr_sqrt(t, t, prec, FMPR_RND_UP);
+            fmpr_sub(t, t, x, prec, FMPR_RND_UP);
+            fmpr_div(t, t, y, prec, FMPR_RND_UP);
+        }
+
+        /* compute upper bound for sqrt(1 + t^2) */
+        fmpr_mul(t, t, t, prec, FMPR_RND_UP);
+        fmpr_add_ui(t, t, 1, prec, FMPR_RND_UP);
+        fmpr_sqrt(bound, t, prec, FMPR_RND_UP);
+    }
+
+    fmpr_clear(x);
+    fmpr_clear(y);
+    fmpr_clear(t);
+    fmpr_clear(u);
+}
+