diff --git a/Makefile.in b/Makefile.in
index 25dc6d1a..a8ca68b3 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -102,5 +102,5 @@ build/%.lo: %.c
 build/%.o: %.c
 	$(CC) -fPIC $(CFLAGS) $(INCS) -c $< -o $@
 
-BUILD_DIRS = fmpr fmprb fmprb_poly fmprb_mat fmpz_holonomic fmpcb fmpcb_poly fmpcb_mat bernoulli
+BUILD_DIRS = fmpr fmprb fmprb_poly fmprb_mat fmpz_holonomic fmpcb fmpcb_poly fmpcb_mat bernoulli hypgeom
 
diff --git a/doc/source/bernoulli.rst b/doc/source/bernoulli.rst
index 5b9c8279..e43f4981 100644
--- a/doc/source/bernoulli.rst
+++ b/doc/source/bernoulli.rst
@@ -1,4 +1,4 @@
-**bernoulli.h** -- support code for Bernoulli numbers
+**bernoulli.h** -- support for Bernoulli numbers
 ===============================================================================
 
 Generation of Bernoulli numbers
diff --git a/doc/source/fmpcb.rst b/doc/source/fmpcb.rst
index fd401204..1b926ea8 100644
--- a/doc/source/fmpcb.rst
+++ b/doc/source/fmpcb.rst
@@ -371,11 +371,11 @@ Special functions
 
     .. math ::
 
-        \begin{align}
-        \sum_{k=N}^U f(k) = & \int_N^U f(t) dt + \frac{1}{2} \left(f(N) + f(U)\right) \\
-                           & + \sum_{k=1}^{M} \frac{B_{2k}}{(2k)!} \left( f^{(2k-1)}(U) - f^{(2k-1)}(N) \right) \\
-                          & - \int_N^U \frac{\tilde B_{2M}(t)}{(2M)!} f^{(2M)}(t) dt
-        \end{align}
+        \sum_{k=N}^U f(k) = \int_N^U f(t) dt + \frac{1}{2} \left(f(N) + f(U)\right)
+
+                           + \sum_{k=1}^{M} \frac{B_{2k}}{(2k)!} \left( f^{(2k-1)}(U) - f^{(2k-1)}(N) \right)
+
+                          - \int_N^U \frac{\tilde B_{2M}(t)}{(2M)!} f^{(2M)}(t) dt
 
     where `f` is a sufficiently differentiable function (for example,
     analytic), `B_n` is a Bernoulli number, and
diff --git a/doc/source/hypgeom.rst b/doc/source/hypgeom.rst
new file mode 100644
index 00000000..78dc8e24
--- /dev/null
+++ b/doc/source/hypgeom.rst
@@ -0,0 +1,195 @@
+**hypgeom.h** -- support for hypergeometric series
+===============================================================================
+
+This module provides functions for high-precision evaluation of series
+of the form
+
+.. math ::
+
+    \sum_{k=0}^{n-1} \frac{A(k)}{B(k)} \prod_{j=1}^k \frac{P(k)}{Q(k)} z^k
+
+where `A, B, P, Q` are polynomials. The present version only supports
+`A, B, P, Q \in \mathbb{Z}[k]` (represented using the
+FLINT *fmpz_poly_t* type). This module also provides functions
+for high-precision evaluation of infinite series (`n \to \infty`),
+with automatic, rigorous error bounding.
+
+Note that we can standardize to `A = B = 1` by
+setting `\tilde P(k) = P(k) A(k) B(k-1), \tilde Q(k) = Q(k) A(k-1) B(k)`.
+However, separating out `A` and `B` is convenient and improves
+efficiency during evaluation.
+
+
+Strategy for error bounding
+-------------------------------------------------------------------------------
+
+We wish to evaluate `S(z) = \sum_{k=0}^{\infty} T(k) z^k` where `T(k)`
+satisfies `T(0) = 1` and
+
+.. math ::
+
+    T(k) = R(k) T(k-1) = \left( \frac{P(k)}{Q(k)} \right) T(k-1)
+
+for given polynomials
+
+.. math ::
+
+    P(k) = a_p k^p + a_{p-1} k^{p-1} + \ldots a_0
+
+    Q(k) = b_q k^q + b_{q-1} k^{q-1} + \ldots b_0.
+
+For convergence, we require `p < q`, or `p = q` with `|z| |a_p| < |b_q|`.
+We also assume that `P(k)` and `Q(k)` have no roots among the positive
+integers (if there are positive integer roots, the sum is either finite
+or undefined). With these conditions satisfied, our goal is to find a
+parameter `n \ge 0` such that
+
+.. math ::
+
+    \left\lvert \sum_{k=n}^{\infty} T(k) z^k \right\rvert \le 2^{-d}.
+
+We can rewrite the hypergeometric term ratio as
+
+.. math ::
+
+    z R(k) = z \frac{P(k)}{Q(k)} =
+        z \left( \frac{a_p}{b_q} \right) \frac{1}{k^{q-p}} F(k)
+
+where
+
+.. math ::
+
+    F(k) = \frac{
+    1 + \tilde a_{1} / k + \tilde a_{2} / k^2 + \ldots + \tilde a_q / k^p
+    }{
+    1 + \tilde b_{1} / k + \tilde b_{2} / k^2 + \ldots + \tilde b_q / k^q
+    } = 1 + O(1/k)
+
+and where `\tilde a_i = a_{p-i} / a_p`, `\tilde b_i = b_{q-i} / b_q`.
+Next, we define
+
+.. math ::
+
+    C = \max_{1 \le i \le p} |\tilde a_i|^{(1/i)},
+    \quad D = \max_{1 \le i \le q} |\tilde b_i|^{(1/i)}.
+
+Now, if `k > C`, the magnitude of the numerator of `F(k)` is
+bounded from above by
+
+.. math ::
+
+    1 + \sum_{i=1}^p \left(\frac{C}{k}\right)^i \le 1 + \frac{C}{k-C}
+
+and if `k > 2D`, the magnitude of the denominator of `F(k)` is bounded
+from below by
+
+.. math ::
+
+    1 - \sum_{i=1}^q \left(\frac{D}{k}\right)^i \ge 1 + \frac{D}{D-k}.
+
+Putting the inequalities together gives the following bound,
+valid for `k > K = \max(C, 2D)`:
+
+.. math ::
+
+    |F(k)| \le \frac{k (k-D)}{(k-C)(k-2D)} = \left(1 + \frac{C}{k-C} \right)
+    \left(1 + \frac{D}{k-2D} \right).
+
+Let `r = q-p` and `\tilde z = |z a_p / b_q|`. Assuming
+`k > \max(C, 2D, {\tilde z}^{1/r})`, we have
+
+.. math ::
+
+    |z R(k)| \le G(k) = \frac{\tilde z F(k)}{k^r}
+
+where `G(k)` is monotonically decreasing. Now we just need to find an
+`n` such that `G(n) < 1` and for which `|T(n)| / (1 - G(n)) \le 2^{-d}`.
+This can be done by computing a floating-point guess for `n` then
+trying successively larger values.
+
+This strategy leaves room for some improvement. For example, if
+`\tilde b_1` is positive and large, the bound `B` becomes very pessimistic
+(a larger positive `\tilde b_1` causes faster convergence,
+not slower convergence).
+
+
+Types, macros and constants
+-------------------------------------------------------------------------------
+
+.. type:: hypgeom_struct
+
+.. type:: hypgeom_t
+
+    Stores polynomials *A*, *B*, *P*, *Q* and precomputed bounds,
+    representing a fixed hypergeometric series.
+
+
+Memory management
+-------------------------------------------------------------------------------
+
+.. function:: void hypgeom_init(hypgeom_t hyp)
+
+.. function:: void hypgeom_clear(hypgeom_t hyp)
+
+
+Error bounding
+-------------------------------------------------------------------------------
+
+.. function:: long hypgeom_estimate_terms(const fmpr_t z, int r, long d)
+
+    Computes an approximation of the largest `n` such
+    that `|z|^n/(n!)^r = 2^{-d}`, giving a first-order estimate of the
+    number of terms needed to approximate the sum of a hypergeometric
+    series of weight `r \ge 0` and argument `z` to an absolute
+    precision of `d \ge 0` bits. If `r = 0`, the direct solution of the
+    equation is given by `n = (\log(1-z) - d \log 2) / \log z`.
+    If `r > 0`, using `\log n! \approx n \log n - n` gives an equation
+    that can be solved in terms of the Lambert *W*-function as
+    `n = (d \log 2) / (r\,W\!(t))` where
+    `t = (d \log 2) / (e r z^{1/r})`.
+
+    The evaluation is done using double precision arithmetic.
+    The function aborts if the computed value of `n` is greater
+    than or equal to LONG_MAX / 2.
+
+.. function:: long hypgeom_bound(fmpr_t error, int r, long C, long D, long K, const fmpr_t TK, const fmpr_t z, long prec)
+
+    Computes a truncation parameter sufficient to achieve *prec* bits
+    of absolute accuracy, according to the strategy described above.
+    The input consists of `r`, `C`, `D`, `K`, precomputed bound for `T(K)`,
+    and `\tilde z = z (a_p / b_q)`, such that for `k > K`, the hypergeometric
+    term ratio is bounded by
+
+    .. math ::
+
+        \frac{\tilde z}{k^r} \frac{k(k-D)}{(k-C)(k-2D)}.
+
+    Given this information, we compute a `\varepsilon` and an
+    integer `n` such that
+    `\left| \sum_{k=n}^{\infty} T(k) \right| \le \varepsilon \le 2^{-\mathrm{prec}}`.
+    The output variable *error* is set to the value of `\varepsilon`,
+    and `n` is returned.
+
+.. function:: void hypgeom_precompute(hypgeom_t hyp)
+
+    Precomputes the bounds data `C`, `D`, `K` and an upper bound for `T(K)`.
+
+Summation
+-------------------------------------------------------------------------------
+
+.. function:: void hypgeom_fmprb_sum(fmprb_t P, fmprb_t Q, const hypgeom_t hyp, const long n, long prec)
+
+    Computes `P, Q` such that `P / Q = \sum_{k=0}^{n-1} T(k)` where `T(k)`
+    is defined by *hyp*,
+    using binary splitting and a working precision of *prec* bits.
+
+.. function:: void hypgeom_fmprb_infsum(fmprb_t P, fmprb_t Q, hypgeom_t hyp, long tol, long prec)
+
+    Computes `P, Q` such that `P / Q = \sum_{k=0}^{\infty} T(k)` where `T(k)`
+    is defined by *hyp*, using binary splitting and
+    working precision of *prec* bits.
+    The number of terms is chosen automatically to bound the
+    truncation error by at most `2^{-\mathrm{tol}}`.
+    The bound for the truncation error is included in the output
+    as part of *P*.
+
diff --git a/doc/source/index.rst b/doc/source/index.rst
index 4fe807ad..bf9df5aa 100644
--- a/doc/source/index.rst
+++ b/doc/source/index.rst
@@ -26,7 +26,7 @@ Module documentation
 ::::::::::::::::::::
 
 .. toctree::
-   :maxdepth: 2
+   :maxdepth: 1
 
    fmpr.rst
    fmprb.rst
@@ -35,8 +35,9 @@ Module documentation
    fmpcb.rst
    fmpcb_poly.rst
    fmpcb_mat.rst
-   fmpz_holonomic.rst
    bernoulli.rst
+   hypgeom.rst
+   fmpz_holonomic.rst
 
 Credits and references
 ::::::::::::::::::::::::
diff --git a/fmprb.h b/fmprb.h
index 44724af5..f6ce6698 100644
--- a/fmprb.h
+++ b/fmprb.h
@@ -288,40 +288,6 @@ void fmprb_rfac_ui_multipoint(fmprb_t y, const fmprb_t x, ulong n, long prec);
 void fmprb_fib_fmpz(fmprb_t f, const fmpz_t n, long prec);
 void fmprb_fib_ui(fmprb_t f, ulong n, long prec);
 
-
-typedef struct
-{
-    fmpz_poly_t A;
-    fmpz_poly_t B;
-    fmpz_poly_t P;
-    fmpz_poly_t Q;
-
-    /* precomputation data */
-    int have_precomputed;
-    long r;
-    long boundA;
-    long boundB;
-    long boundK;
-    fmpr_t MK;
-}
-hypgeom_struct;
-
-typedef hypgeom_struct hypgeom_t[1];
-
-void hypgeom_init(hypgeom_t hyp);
-
-void hypgeom_clear(hypgeom_t hyp);
-
-void hypgeom_precompute(hypgeom_t hyp);
-
-long hypgeom_bound(fmpr_t error, int r,
-    long A, long B, long K, const fmpr_t TK, const fmpr_t z, long prec);
-
-void fmprb_hypgeom_sum(fmprb_t P, fmprb_t Q, const hypgeom_t hyp, const long n, long prec);
-
-void fmprb_hypgeom_infsum(fmprb_t P, fmprb_t Q, hypgeom_t hyp, long target_prec, long prec);
-
-
 void fmprb_const_pi_chudnovsky(fmprb_t x, long prec);
 void fmprb_const_pi(fmprb_t x, long prec);
 
diff --git a/fmprb/const_pi_chudnovsky.c b/fmprb/const_pi_chudnovsky.c
index f6c3b520..1508542d 100644
--- a/fmprb/const_pi_chudnovsky.c
+++ b/fmprb/const_pi_chudnovsky.c
@@ -24,6 +24,7 @@
 ******************************************************************************/
 
 #include "fmprb.h"
+#include "hypgeom.h"
 
 void
 fmprb_const_pi_chudnovsky(fmprb_t s, long prec)
diff --git a/fmprb/const_zeta3_bsplit.c b/fmprb/const_zeta3_bsplit.c
index 9f8ed7d2..b6d1eead 100644
--- a/fmprb/const_zeta3_bsplit.c
+++ b/fmprb/const_zeta3_bsplit.c
@@ -24,7 +24,7 @@
 ******************************************************************************/
 
 #include "fmprb.h"
-
+#include "hypgeom.h"
 
 void
 fmprb_const_zeta3_bsplit(fmprb_t s, long prec)
diff --git a/hypgeom.h b/hypgeom.h
new file mode 100644
index 00000000..cccca5cc
--- /dev/null
+++ b/hypgeom.h
@@ -0,0 +1,67 @@
+/*=============================================================================
+
+    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
+
+******************************************************************************/
+
+#ifndef HYPGEOM_H
+#define HYPGEOM_H
+
+#include "fmprb.h"
+#include "fmpz_poly.h"
+
+typedef struct
+{
+    fmpz_poly_t A;
+    fmpz_poly_t B;
+    fmpz_poly_t P;
+    fmpz_poly_t Q;
+
+    /* precomputation data */
+    int have_precomputed;
+    long r;
+    long boundC;
+    long boundD;
+    long boundK;
+    fmpr_t MK;
+}
+hypgeom_struct;
+
+typedef hypgeom_struct hypgeom_t[1];
+
+void hypgeom_init(hypgeom_t hyp);
+
+void hypgeom_clear(hypgeom_t hyp);
+
+void hypgeom_precompute(hypgeom_t hyp);
+
+long hypgeom_estimate_terms(const fmpr_t z, int r, long prec);
+
+long hypgeom_bound(fmpr_t error, int r,
+    long C, long D, long K, const fmpr_t TK, const fmpr_t z, long prec);
+
+void fmprb_hypgeom_sum(fmprb_t P, fmprb_t Q, const hypgeom_t hyp, const long n, long prec);
+
+void fmprb_hypgeom_infsum(fmprb_t P, fmprb_t Q, hypgeom_t hyp, long target_prec, long prec);
+
+#endif
+
diff --git a/hypgeom/Makefile b/hypgeom/Makefile
new file mode 100644
index 00000000..c7e54132
--- /dev/null
+++ b/hypgeom/Makefile
@@ -0,0 +1,43 @@
+SOURCES = $(wildcard *.c)
+
+OBJS = $(patsubst %.c, $(BUILD_DIR)/%.o, $(SOURCES))
+
+LIB_OBJS = $(patsubst %.c, $(BUILD_DIR)/%.lo, $(SOURCES))
+
+TEST_SOURCES = $(wildcard test/*.c)
+
+PROF_SOURCES = $(wildcard profile/*.c)
+
+TUNE_SOURCES = $(wildcard tune/*.c)
+
+TESTS = $(patsubst %.c, %, $(TEST_SOURCES))
+
+PROFS = $(patsubst %.c, %, $(PROF_SOURCES))
+
+TUNE = $(patsubst %.c, %, $(TUNE_SOURCES))
+
+all: $(OBJS)
+
+library: $(LIB_OBJS)
+
+profile:
+	$(foreach prog, $(PROFS), $(CC) -O2 -std=c99 $(INCS) $(prog).c ../profiler.o -o $(BUILD_DIR)/$(prog) $(LIBS) -lflint;)
+        
+tune: $(TUNE_SOURCES)
+	$(foreach prog, $(TUNE), $(CC) -O2 -std=c99 $(INCS) $(prog).c -o $(BUILD_DIR)/$(prog) $(LIBS) -lflint;)
+
+$(BUILD_DIR)/%.o: %.c
+	$(CC) $(CFLAGS) -c $(INCS) $< -o $@
+
+$(BUILD_DIR)/%.lo: %.c
+	$(CC) -fPIC $(CFLAGS) $(INCS) -c $< -o $@
+
+clean:
+	rm -rf $(BUILD_DIR)	
+
+check: library
+	$(foreach prog, $(TESTS), $(CC) $(CFLAGS) $(INCS) $(prog).c -o $(BUILD_DIR)/$(prog) $(LIBS) -lflint;)
+	$(foreach prog, $(TESTS), $(BUILD_DIR)/$(prog);)
+
+.PHONY: profile clean check all
+
diff --git a/fmprb/hypgeom_bound.c b/hypgeom/bound.c
similarity index 50%
rename from fmprb/hypgeom_bound.c
rename to hypgeom/bound.c
index 28eba7d4..2801d3be 100644
--- a/fmprb/hypgeom_bound.c
+++ b/hypgeom/bound.c
@@ -24,15 +24,12 @@
 ******************************************************************************/
 
 #include <math.h>
-#include "fmprb.h"
 #include "double_extras.h"
+#include "hypgeom.h"
 
 double fmpr_get_d(const fmpr_t x);
 
-#define LOG2 0.69314718055994530942
-#define EXP1 2.7182818284590452354
-
-double d_root(double x, int r)
+static __inline__ double d_root(double x, int r)
 {
     if (r == 1)
         return x;
@@ -41,46 +38,6 @@ double d_root(double x, int r)
     return pow(x, 1. / r);
 }
 
-/*
-Estimate the truncation point to obtain accuracy 2^(-prec) with the
-hypergeometric series |z|^k / (k!)^r.
-*/
-long
-estimate_nterms(double z, int r, long prec)
-{
-    double y;
-
-    z = fabs(z);
-
-    if (z == 0)
-        return 1;
-
-    if (r == 0)
-    {
-        if (z >= 1)
-        {
-            printf("z must be smaller than 1\n");
-            abort();
-        }
-
-        y = (log(1-z) - prec * LOG2) / log(z) + 1;
-    }
-    else
-    {
-        /* Solve k*log(z) - r*(k*log(k)-k) = -prec*log(2)  */
-        y = prec * LOG2 / (d_root(z, r) * EXP1 * r);
-        y =  prec * LOG2 / (r * d_lambertw(y)) + 1;
-    }
-
-    if (y >= LONG_MAX / 2)
-    {
-        printf("error: series will converge too slowly\n");
-        abort();
-    }
-
-    return y;
-}
-
 void
 fmpr_gamma_ui_lbound(fmpr_t x, ulong n, long prec)
 {
@@ -157,151 +114,153 @@ fmpr_gamma_ui_ubound(fmpr_t x, ulong n, long prec)
     }
 }
 
-/* FIXME: assumes no overflow when computing n + r */
-static void
-fmpr_rfac_uiui_ubound(fmpr_t x, ulong n, ulong r, long prec)
+/* FIXME: use more than doubles for this */
+long hypgeom_root_bound(const fmpr_t z, int r)
 {
     if (r == 0)
     {
-        fmpr_one(x);
-    }
-    else if (r == 1)
-    {
-        fmpr_set_ui(x, n);
+        return 0;
     }
     else
     {
-        fmpr_t t;
-        fmpr_init(t);
-        fmpr_gamma_ui_ubound(x, n + r, prec);
-        fmpr_gamma_ui_lbound(t, n, prec);
-        fmpr_div(x, x, t, prec, FMPR_RND_UP);
-        fmpr_clear(t);
+        double zd;
+        zd = fmpr_get_d(z);
+        return d_root(zd, r) + 2;
     }
 }
 
-static void
-fmpr_rfac_uiui_lbound(fmpr_t x, ulong n, ulong r, long prec)
-{
-    if (r == 0)
-    {
-        fmpr_one(x);
-    }
-    else if (r == 1)
-    {
-        fmpr_set_ui(x, n);
-    }
-    else
-    {
-        fmpr_t t;
-        fmpr_init(t);
-        fmpr_gamma_ui_lbound(x, n + r, prec);
-        fmpr_gamma_ui_ubound(t, n, prec);
-        fmpr_div(x, x, t, prec, FMPR_RND_DOWN);
-        fmpr_clear(t);
-    }
-}
-
-
 /*
+Given T(K), compute bound for T(n) z^n.
 
-The general term T(k) is z^k / (k!)^r * R(k) where R(k) = 1 + O(1/k).
+We need to multiply by
 
-We have precomputed integers A, B, K such that for all k > K,
-|R(k)| <= F(k) = k(k-B)/((k-A)(k-2B)) = (1+A/(k-A))(1+B/(k-2B)).
+z^n * 1/rf(K+1,m)^r * (rf(K+1,m)/rf(K+1-A,m)) * (rf(K+1-B,m)/rf(K+1-2B,m))
 
+where m = n - K. This is equal to
+
+z^n * 
+
+(K+A)! (K-2B)! (K-B+m)!
+-----------------------    * ((K+m)! / K!)^(1-r)
+(K-B)! (K-A+m)! (K-2B+m)!
 */
+void
+hypgeom_term_bound(fmpr_t Tn, const fmpr_t TK, long K, long A, long B, int r, const fmpr_t z, long n, long wp)
+{
+    fmpr_t t, u, num, den;
+    long m;
+
+    fmpr_init(t);
+    fmpr_init(u);
+    fmpr_init(num);
+    fmpr_init(den);
+
+    m = n - K;
+    if (m < 0)
+        abort();
+
+    /* TK * z^n */
+    fmpr_pow_sloppy_ui(t, z, n, wp, FMPR_RND_UP);
+    fmpr_mul(num, TK, t, wp, FMPR_RND_UP);
+
+    /* (K+A)! (K-2B)! (K-B+m)!, upper bounding */
+    fmpr_gamma_ui_ubound(t, K+A+1, wp);
+    fmpr_mul(num, num, t, wp, FMPR_RND_UP);
+    fmpr_gamma_ui_ubound(t, K-2*B+1, wp);
+    fmpr_mul(num, num, t, wp, FMPR_RND_UP);
+    fmpr_gamma_ui_ubound(t, K-B+m, wp);
+    fmpr_mul(num, num, t, wp, FMPR_RND_UP);
+
+    /* (K-B)! (K-A+m)! (K-2B+m)!, lower bounding */
+    fmpr_gamma_ui_lbound(den, K-B+1, wp);
+    fmpr_gamma_ui_lbound(t, K-A+m+1, wp);
+    fmpr_mul(den, den, t, wp, FMPR_RND_DOWN);
+    fmpr_gamma_ui_lbound(t, K-2*B+m+1, wp);
+    fmpr_mul(den, den, t, wp, FMPR_RND_DOWN);
+
+    /* ((K+m)! / K!)^(1-r) */
+    if (r == 0)
+    {
+        fmpr_gamma_ui_ubound(t, K+m+1, wp);
+        fmpr_mul(num, num, t, wp, FMPR_RND_UP);
+        fmpr_gamma_ui_lbound(t, K+1, wp);
+        fmpr_mul(den, den, t, wp, FMPR_RND_DOWN);
+    }
+    else if (r != 1)
+    {
+        fmpr_gamma_ui_ubound(t, K+1, wp);
+        fmpr_gamma_ui_lbound(u, K+m+1, wp);
+        fmpr_div(t, t, u, wp, FMPR_RND_UP);
+        fmpr_pow_sloppy_ui(t, t, r-1, wp, FMPR_RND_UP);
+        fmpr_mul(num, num, t, wp, FMPR_RND_UP);
+    }
+
+    fmpr_div(Tn, num, den, wp, FMPR_RND_UP);
+
+    fmpr_clear(t);
+    fmpr_clear(u);
+    fmpr_clear(num);
+    fmpr_clear(den);
+}
+
 long
 hypgeom_bound(fmpr_t error, int r,
     long A, long B, long K, const fmpr_t TK, const fmpr_t z, long prec)
 {
-    fmpr_t Tn, t, u, one, tol;
+    fmpr_t Tn, t, u, one, tol, num, den;
     long wp = FMPRB_RAD_PREC;
-    long n;
-    double zd;
+    long n, m;
 
     fmpr_init(Tn);
     fmpr_init(t);
     fmpr_init(u);
     fmpr_init(one);
     fmpr_init(tol);
+    fmpr_init(num);
+    fmpr_init(den);
 
     fmpr_one(one);
     fmpr_set_ui_2exp_si(tol, 1UL, -prec);
 
-    zd = fmpr_get_d(z);
-
-    n = estimate_nterms(zd, r, prec);
+    /* approximate number of needed terms */
+    n = hypgeom_estimate_terms(z, r, prec);
 
     /* required for 1 + O(1/k) part to be decreasing */
     n = FLINT_MAX(n, K + 1);
 
-    /* required for z^k / (k!)^r to be decreasing
-       (TODO: don't use doubles for this) */
-    if (r > 0)
-    {
-        long nbd = d_root(zd, r) + 2;
-        n = FLINT_MAX(n, nbd);
-    }
+    /* required for z^k / (k!)^r to be decreasing */
+    m = hypgeom_root_bound(z, r);
+    n = FLINT_MAX(n, m);
 
-    /* We are now sure that |R(k)| is either decreasing or strictly
-       smaller than 1 for k >= n, which means that we can bound the tail
-       using a geometric series as soon as soon as |R(k)| < 1  */
+    /*  We now have |R(k)| <= G(k) where G(k) is monotonically decreasing,
+        and can bound the tail using a geometric series as soon
+        as soon as G(k) < 1.  */
 
-    /* Compute an upper bound for T(n) */
-    /* (z^n) / (n!)^r * TK * [(K+1)(K+2)...(n)] * [(K-B+1)(K-B+2)...(n-B)]
-                             ---------------------------------------------
-                             [(K-A+1)(K-A+2)...(n)] * [(K-2B+1)(K-2B+2)...(n-2B)]
-    */
-
-    /* z^n * TK */
-    fmpr_pow_sloppy_ui(Tn, z, n, wp, FMPR_RND_UP);
-    fmpr_mul(Tn, Tn, TK, wp, FMPR_RND_UP);
-
-    /* divide by (n!)^r */
-    if (r != 0)
-    {
-        fmpr_gamma_ui_lbound(t, n + 1, wp);
-        fmpr_ui_div(t, 1UL, t, wp, FMPR_RND_UP);
-        fmpr_pow_sloppy_ui(t, t, r, wp, FMPR_RND_UP);
-        fmpr_mul(Tn, Tn, t, wp, FMPR_RND_UP);
-    }
-
-    fmpr_rfac_uiui_ubound(t, K+1, n-K, wp);
-    fmpr_mul(Tn, Tn, t, wp, FMPR_RND_UP);
-
-    fmpr_rfac_uiui_ubound(t, K-B+1, n-K, wp);
-    fmpr_mul(Tn, Tn, t, wp, FMPR_RND_UP);
-
-    fmpr_rfac_uiui_lbound(t, K-A+1, n-K, wp);
-    fmpr_div(Tn, Tn, t, wp, FMPR_RND_UP);
-
-    fmpr_rfac_uiui_lbound(t, K-2*B+1, n-K, wp);
-    fmpr_div(Tn, Tn, t, wp, FMPR_RND_UP);
+    /* bound T(n-1) */
+    hypgeom_term_bound(Tn, TK, K, A, B, r, z, n-1, wp);
 
     while (1)
     {
-        /* bound for term ratio: z * F(n) / n^r */
-
-        /* F(n) <= n (n-B) / ((n-A) (n-2B)) */
-        fmpr_set_ui(t, n);
-        fmpr_mul_ui(t, t, n - B, wp, FMPR_RND_UP);
-        fmpr_div_ui(t, t, n - A, wp, FMPR_RND_UP);
-        fmpr_div_ui(t, t, n - 2*B, wp, FMPR_RND_UP);
-
-        fmpr_mul(t, t, z, wp, FMPR_RND_UP);
+        /* bound R(n) */
+        fmpr_mul_ui(num, z, n, wp, FMPR_RND_UP);
+        fmpr_mul_ui(num, num, n - B, wp, FMPR_RND_UP);
+        fmpr_set_ui(den, n - A);
+        fmpr_mul_ui(den, den, n - 2*B, wp, FMPR_RND_DOWN);
 
         if (r != 0)
         {
-            fmpr_div_ui(u, one, n, wp, FMPR_RND_UP);
-            fmpr_pow_sloppy_ui(u, u, r, wp, FMPR_RND_UP);
-            fmpr_mul(t, t, u, wp, FMPR_RND_UP);
+            fmpr_set_ui(u, n);
+            fmpr_pow_sloppy_ui(u, u, r, wp, FMPR_RND_DOWN);
+            fmpr_mul(den, den, u, wp, FMPR_RND_DOWN);
         }
 
-        /* bound by geometric series: Tn / (1 - t) */
-        /* where the term ratio must be < 1 */
-        fmpr_sub(u, one, t, wp, FMPR_RND_DOWN);
+        fmpr_div(t, num, den, wp, FMPR_RND_UP);
 
+        /* multiply bound for T(n-1) by bound for R(n) to bound T(n) */
+        fmpr_mul(Tn, Tn, t, wp, FMPR_RND_UP);
+
+        /* geometric series termination check */
+        fmpr_sub(u, one, t, wp, FMPR_RND_DOWN);
         if (fmpr_sgn(u) > 0)
         {
             fmpr_div(u, Tn, u, wp, FMPR_RND_UP);
@@ -314,7 +273,6 @@ hypgeom_bound(fmpr_t error, int r,
         }
 
         /* move on to next term */
-        fmpr_mul(Tn, Tn, t, wp, FMPR_RND_UP);
         n++;
     }
 
@@ -323,6 +281,8 @@ hypgeom_bound(fmpr_t error, int r,
     fmpr_clear(u);
     fmpr_clear(one);
     fmpr_clear(tol);
+    fmpr_clear(num);
+    fmpr_clear(den);
 
     return n;
 }
diff --git a/hypgeom/estimate_terms_d.c b/hypgeom/estimate_terms_d.c
new file mode 100644
index 00000000..a8c4dcda
--- /dev/null
+++ b/hypgeom/estimate_terms_d.c
@@ -0,0 +1,79 @@
+/*=============================================================================
+
+    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 "double_extras.h"
+#include "hypgeom.h"
+
+#define LOG2 0.69314718055994530942
+#define EXP1 2.7182818284590452354
+
+static __inline__ double d_root(double x, int r)
+{
+    if (r == 1)
+        return x;
+    if (r == 2)
+        return sqrt(x);
+    return pow(x, 1. / r);
+}
+
+double fmpr_get_d(const fmpr_t x);
+
+long
+hypgeom_estimate_terms(const fmpr_t z, int r, long prec)
+{
+    double y, t;
+
+    t = fmpr_get_d(z);
+    t = fabs(t);
+
+    if (t == 0)
+        return 1;
+
+    if (r == 0)
+    {
+        if (t >= 1)
+        {
+            printf("z must be smaller than 1\n");
+            abort();
+        }
+
+        y = (log(1-t) - prec * LOG2) / log(t) + 1;
+    }
+    else
+    {
+        y = (prec * LOG2) / (d_root(t, r) * EXP1 * r);
+        y = (prec * LOG2) / (r * d_lambertw(y)) + 1;
+    }
+
+    if (y >= LONG_MAX / 2)
+    {
+        printf("error: series will converge too slowly\n");
+        abort();
+    }
+
+    return y;
+}
+
diff --git a/fmprb/hypgeom_init.c b/hypgeom/init.c
similarity index 98%
rename from fmprb/hypgeom_init.c
rename to hypgeom/init.c
index 15749ed3..9d5fbd01 100644
--- a/fmprb/hypgeom_init.c
+++ b/hypgeom/init.c
@@ -23,7 +23,7 @@
 
 ******************************************************************************/
 
-#include "fmprb.h"
+#include "hypgeom.h"
 
 void
 hypgeom_init(hypgeom_t hyp)
diff --git a/fmprb/hypgeom_precompute.c b/hypgeom/precompute.c
similarity index 96%
rename from fmprb/hypgeom_precompute.c
rename to hypgeom/precompute.c
index c624cf6f..d98de9f3 100644
--- a/fmprb/hypgeom_precompute.c
+++ b/hypgeom/precompute.c
@@ -24,7 +24,7 @@
 ******************************************************************************/
 
 #include <math.h>
-#include "fmprb.h"
+#include "hypgeom.h"
 
 /* Compute a pure ratio P2(k)/Q2(k) for the term
    A(k)/B(k) * [P(1)P(2)...P(k)] / [Q(1)Q(2)...Q(k)] */
@@ -122,9 +122,9 @@ _hypgeom_precompute(hypgeom_t hyp, const fmpz_poly_t P, const fmpz_poly_t Q)
     fmpz_init(t);
 
     hyp->r = fmpz_poly_degree(Q) - fmpz_poly_degree(P);
-    hyp->boundA = hypgeom_root_norm(P);
-    hyp->boundB = hypgeom_root_norm(Q);
-    hyp->boundK = 1 + FLINT_MAX(hyp->boundA, 2 * hyp->boundB);
+    hyp->boundC = hypgeom_root_norm(P);
+    hyp->boundD = hypgeom_root_norm(Q);
+    hyp->boundK = 1 + FLINT_MAX(hyp->boundC, 2 * hyp->boundD);
 
     fmpr_one(hyp->MK);
 
diff --git a/fmprb/hypgeom_sum.c b/hypgeom/sum.c
similarity index 98%
rename from fmprb/hypgeom_sum.c
rename to hypgeom/sum.c
index bc4dac05..c7cf6bf1 100644
--- a/fmprb/hypgeom_sum.c
+++ b/hypgeom/sum.c
@@ -23,7 +23,7 @@
 
 ******************************************************************************/
 
-#include "fmprb.h"
+#include "hypgeom.h"
 
 static __inline__ void
 fmpz_poly_evaluate_si(fmpz_t y, const fmpz_poly_t poly, long x)
@@ -201,7 +201,7 @@ fmprb_hypgeom_infsum(fmprb_t P, fmprb_t Q, hypgeom_t hyp, long target_prec, long
         hyp->have_precomputed = 1;
     }
 
-    n = hypgeom_bound(err, hyp->r, hyp->boundA, hyp->boundB,
+    n = hypgeom_bound(err, hyp->r, hyp->boundC, hyp->boundD,
         hyp->boundK, hyp->MK, z, target_prec);
 
     fmprb_hypgeom_sum(P, Q, hyp, n, prec);