diff --git a/acb_hypgeom.h b/acb_hypgeom.h
index 24a419eb..2d715dfc 100644
--- a/acb_hypgeom.h
+++ b/acb_hypgeom.h
@@ -64,44 +64,38 @@ void acb_hypgeom_pfq_series_direct(acb_poly_t res,
 void acb_hypgeom_u_asymp(acb_t res, const acb_t a, const acb_t b,
     const acb_t z, long n, long prec);
 
+void acb_hypgeom_u_1f1_series(acb_poly_t res,
+    const acb_poly_t a, const acb_poly_t b, const acb_poly_t z,
+    long len, long prec);
+
+void acb_hypgeom_u_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, long prec);
+void acb_hypgeom_u(acb_t res, const acb_t a, const acb_t b, const acb_t z, long prec);
+
 int acb_hypgeom_u_use_asymp(const acb_t z, long prec);
 
 void acb_hypgeom_erf_1f1a(acb_t res, const acb_t z, long prec);
-
 void acb_hypgeom_erf_1f1b(acb_t res, const acb_t z, long prec);
-
 void acb_hypgeom_erf_asymp(acb_t res, const acb_t z, long prec, long prec2);
-
 void acb_hypgeom_erf(acb_t res, const acb_t z, long prec);
 
 void acb_hypgeom_bessel_j_0f1(acb_t res, const acb_t nu, const acb_t z, long prec);
-
 void acb_hypgeom_bessel_j_asymp(acb_t res, const acb_t nu, const acb_t z, long prec);
-
 void acb_hypgeom_bessel_j(acb_t res, const acb_t nu, const acb_t z, long prec);
 
 void acb_hypgeom_bessel_k_0f1(acb_t res, const acb_t nu, const acb_t z, long prec);
-
 void acb_hypgeom_bessel_k_0f1_series(acb_poly_t res, const acb_poly_t n, const acb_poly_t z, long len, long prec);
-
 void acb_hypgeom_bessel_k_asymp(acb_t res, const acb_t nu, const acb_t z, long prec);
-
 void acb_hypgeom_bessel_k(acb_t res, const acb_t nu, const acb_t z, long prec);
 
 void acb_hypgeom_gamma_upper_asymp(acb_t res, const acb_t s, const acb_t z, int modified, long prec);
-
 void acb_hypgeom_gamma_upper_1f1a(acb_t res, const acb_t s, const acb_t z, int modified, long prec);
-
 void acb_hypgeom_gamma_upper_1f1b(acb_t res, const acb_t s, const acb_t z, int modified, long prec);
-
 void acb_hypgeom_gamma_upper_singular(acb_t res, long s, const acb_t z, int modified, long prec);
-
 void acb_hypgeom_gamma_upper(acb_t res, const acb_t s, const acb_t z, int modified, long prec);
 
 void acb_hypgeom_expint(acb_t res, const acb_t s, const acb_t z, long prec);
 
 void acb_hypgeom_erfc(acb_t res, const acb_t z, long prec);
-
 void acb_hypgeom_erfi(acb_t res, const acb_t z, long prec);
 
 void acb_hypgeom_ei_asymp(acb_t res, const acb_t z, long prec);
diff --git a/acb_hypgeom/si.c b/acb_hypgeom/si.c
index a389bef2..c970691a 100644
--- a/acb_hypgeom/si.c
+++ b/acb_hypgeom/si.c
@@ -67,7 +67,7 @@ acb_hypgeom_si_asymp(acb_t res, const acb_t z, long prec)
 
     if (arb_is_zero(acb_realref(z)))
     {
-        /* the function is real */
+        /* the function is imaginary */
         arb_zero(acb_realref(t));
     }
     else if (arb_is_positive(acb_realref(z)))
diff --git a/acb_hypgeom/test/t-u.c b/acb_hypgeom/test/t-u.c
new file mode 100644
index 00000000..ed6e9e17
--- /dev/null
+++ b/acb_hypgeom/test/t-u.c
@@ -0,0 +1,200 @@
+/*=============================================================================
+
+    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) 2014 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_hypgeom.h"
+
+void
+acb_randtest_maybe_half_int(acb_t x, flint_rand_t state, long prec, long size)
+{
+    if (n_randint(state, 8) == 0)
+    {
+        fmpz_t t;
+        fmpz_init(t);
+        fmpz_randtest(t, state, 1 + n_randint(state, prec));
+        arb_set_fmpz(acb_realref(x), t);
+        arb_zero(acb_imagref(x));
+        acb_mul_2exp_si(x, x, -1);
+        fmpz_clear(t);
+    }
+    else
+    {
+        acb_randtest(x, state, prec, size);
+    }
+}
+
+void
+acb_hypgeom_u_asymp_proper(acb_t res, const acb_t a, const acb_t b, const acb_t z, long prec)
+{
+    acb_t t;
+    acb_init(t);
+    acb_pow(t, z, a, prec);
+    acb_hypgeom_u_asymp(res, a, b, z, -1, prec);
+    acb_div(res, res, t, prec);
+    acb_clear(t);
+}
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("u....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 2000; iter++)
+    {
+        acb_t a0, a1, a2, b, z, w0, w1, w2, t, u;
+        long prec0, prec1, prec2;
+
+        acb_init(a0);
+        acb_init(a1);
+        acb_init(a2);
+        acb_init(b);
+        acb_init(z);
+        acb_init(w0);
+        acb_init(w1);
+        acb_init(w2);
+        acb_init(t);
+        acb_init(u);
+
+        prec0 = 2 + n_randint(state, 700);
+        prec1 = 2 + n_randint(state, 700);
+        prec2 = 2 + n_randint(state, 700);
+
+        acb_randtest_maybe_half_int(a0, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
+        acb_randtest_maybe_half_int(b, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
+        acb_randtest(z, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
+        acb_randtest(w0, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
+        acb_randtest(w1, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
+        acb_randtest(w2, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100));
+
+        acb_add_ui(a1, a0, 1, prec0);
+        acb_add_ui(a2, a0, 2, prec0);
+
+        switch (n_randint(state, 3))
+        {
+            case 0:
+                acb_hypgeom_u_asymp_proper(w0, a0, b, z, prec0);
+                break;
+            case 1:
+                acb_hypgeom_u_1f1(w0, a0, b, z, prec0);
+                break;
+            default:
+                acb_hypgeom_u(w0, a0, b, z, prec0);
+        }
+
+        switch (n_randint(state, 3))
+        {
+            case 0:
+                acb_hypgeom_u_asymp_proper(w1, a0, b, z, prec1);
+                break;
+            case 1:
+                acb_hypgeom_u_1f1(w1, a0, b, z, prec1);
+                break;
+            default:
+                acb_hypgeom_u(w1, a0, b, z, prec1);
+        }
+
+        if (!acb_overlaps(w0, w1))
+        {
+            printf("FAIL: consistency\n\n");
+            printf("a = "); acb_printd(a0, 30); printf("\n\n");
+            printf("b = "); acb_printd(b, 30); printf("\n\n");
+            printf("z = "); acb_printd(z, 30); printf("\n\n");
+            printf("w0 = "); acb_printd(w0, 30); printf("\n\n");
+            printf("w1 = "); acb_printd(w1, 30); printf("\n\n");
+            abort();
+        }
+
+        switch (n_randint(state, 3))
+        {
+            case 0:
+                acb_hypgeom_u_asymp_proper(w1, a1, b, z, prec1);
+                break;
+            case 1:
+                acb_hypgeom_u_1f1(w1, a1, b, z, prec1);
+                break;
+            default:
+                acb_hypgeom_u(w1, a1, b, z, prec1);
+        }
+
+        switch (n_randint(state, 3))
+        {
+            case 0:
+                acb_hypgeom_u_asymp_proper(w2, a2, b, z, prec2);
+                break;
+            case 1:
+                acb_hypgeom_u_1f1(w2, a2, b, z, prec2);
+                break;
+            default:
+                acb_hypgeom_u(w2, a2, b, z, prec2);
+        }
+
+        acb_set(t, w0);
+
+        acb_mul_2exp_si(u, a0, 1);
+        acb_sub(u, u, b, prec0);
+        acb_add(u, u, z, prec0);
+        acb_add_ui(u, u, 2, prec0);
+
+        acb_submul(t, w1, u, prec0);
+
+        acb_sub(u, a2, b, prec0);
+        acb_mul(u, u, a1, prec0);
+        acb_addmul(t, w2, u, prec0);
+
+        if (!acb_contains_zero(t))
+        {
+            printf("FAIL: contiguous relation\n\n");
+            printf("a = "); acb_printd(a0, 30); printf("\n\n");
+            printf("b = "); acb_printd(b, 30); printf("\n\n");
+            printf("z = ");  acb_printd(z, 30); printf("\n\n");
+            printf("w0 = "); acb_printd(w0, 30); printf("\n\n");
+            printf("w1 = "); acb_printd(w1, 30); printf("\n\n");
+            printf("w2 = "); acb_printd(w2, 30); printf("\n\n");
+            printf("t = "); acb_printd(t, 30); printf("\n\n");
+            abort();
+        }
+
+        acb_clear(a0);
+        acb_clear(a1);
+        acb_clear(a2);
+        acb_clear(b);
+        acb_clear(z);
+        acb_clear(w0);
+        acb_clear(w1);
+        acb_clear(w2);
+        acb_clear(t);
+        acb_clear(u);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+
diff --git a/acb_hypgeom/u.c b/acb_hypgeom/u.c
new file mode 100644
index 00000000..6a37bf32
--- /dev/null
+++ b/acb_hypgeom/u.c
@@ -0,0 +1,201 @@
+/*=============================================================================
+
+    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) 2015 Fredrik Johansson
+
+******************************************************************************/
+
+#include "acb_hypgeom.h"
+
+void
+acb_hypgeom_u_1f1_series(acb_poly_t res,
+    const acb_poly_t a, const acb_poly_t b, const acb_poly_t z,
+    long len, long prec)
+{
+    acb_poly_t s, u, A, B;
+    acb_poly_struct aa[3];
+    arb_t c;
+    long wlen;
+    int singular;
+
+    acb_poly_init(s);
+    acb_poly_init(u);
+    acb_poly_init(A);
+    acb_poly_init(B);
+    acb_poly_init(aa + 0);
+    acb_poly_init(aa + 1);
+    acb_poly_init(aa + 2);
+    arb_init(c);
+
+    singular = (b->length == 0) || acb_is_int(b->coeffs);
+    wlen = len + (singular != 0);
+
+    /* A = rgamma(a-b+1) * 1F~1(a,b,z) */
+    acb_poly_sub(u, a, b, prec);
+    acb_poly_add_si(u, u, 1, prec);
+    acb_poly_rgamma_series(A, u, wlen, prec);
+
+    /* todo: handle a = 1 efficiently */
+    acb_poly_set(aa, a);
+    acb_poly_set(aa + 1, b);
+    acb_poly_one(aa + 2);
+    acb_hypgeom_pfq_series_direct(s, aa, 1, aa + 1, 2, z, -1, wlen, prec, 1);
+    acb_poly_mullow(A, A, s, wlen, prec);
+
+    /* B = rgamma(a) * 1F~1(a-b+1,2-b,z) * z^(1-b) */
+    acb_poly_set(aa, u);
+    acb_poly_add_si(aa + 1, b, -2, prec);
+    acb_poly_neg(aa + 1, aa + 1);
+    acb_hypgeom_pfq_series_direct(s, aa, 1, aa + 1, 2, z, -1, wlen, prec, 1);
+    acb_poly_rgamma_series(B, a, wlen, prec);
+    acb_poly_mullow(B, B, s, wlen, prec);
+
+    acb_poly_add_si(u, b, -1, prec);
+    acb_poly_neg(u, u);
+    acb_poly_pow_series(s, z, u, wlen, prec);
+    acb_poly_mullow(B, B, s, wlen, prec);
+
+    acb_poly_sub(A, A, B, prec);
+
+    /* multiply by pi csc(pi b) */
+    acb_poly_sin_pi_series(B, b, wlen, prec);
+
+    if (singular)
+    {
+        acb_poly_shift_right(A, A, 1);
+        acb_poly_shift_right(B, B, 1);
+    }
+
+    acb_poly_div_series(res, A, B, len, prec);
+
+    arb_const_pi(c, prec);
+    _acb_vec_scalar_mul_arb(res->coeffs, res->coeffs, res->length, c, prec);
+
+    acb_poly_clear(s);
+    acb_poly_clear(u);
+    acb_poly_clear(A);
+    acb_poly_clear(B);
+    acb_poly_clear(aa + 0);
+    acb_poly_clear(aa + 1);
+    acb_poly_clear(aa + 2);
+    arb_clear(c);
+}
+
+void
+acb_hypgeom_u_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, long prec)
+{
+    if (acb_is_int(b))
+    {
+        acb_poly_t aa, bb, zz;
+
+        acb_poly_init(aa);
+        acb_poly_init(bb);
+        acb_poly_init(zz);
+
+        acb_poly_set_acb(aa, a);
+        acb_poly_set_coeff_acb(bb, 0, b);
+        acb_poly_set_coeff_si(bb, 1, 1);
+        acb_poly_set_acb(zz, z);
+
+        acb_hypgeom_u_1f1_series(zz, aa, bb, zz, 1, prec);
+        acb_poly_get_coeff_acb(res, zz, 0);
+
+        acb_poly_clear(aa);
+        acb_poly_clear(bb);
+        acb_poly_clear(zz);
+    }
+    else
+    {
+        acb_t t, u, v;
+        acb_struct aa[3];
+
+        acb_init(t);
+        acb_init(u);
+        acb_init(v);
+        acb_init(aa + 0);
+        acb_init(aa + 1);
+        acb_init(aa + 2);
+
+        acb_set(aa, a);
+        acb_set(aa + 1, b);
+        acb_one(aa + 2);
+        acb_hypgeom_pfq_direct(u, aa, 1, aa + 1, 2, z, -1, prec);
+        acb_sub(aa, a, b, prec);
+        acb_add_ui(aa, aa, 1, prec);
+        acb_sub_ui(aa + 1, b, 2, prec);
+        acb_neg(aa + 1, aa + 1);
+        acb_hypgeom_pfq_direct(v, aa, 1, aa + 1, 2, z, -1, prec);
+
+        acb_sub_ui(aa + 1, b, 1, prec);
+
+        /* rgamma(a-b+1) * gamma(1-b) * u */
+        acb_rgamma(t, aa, prec);
+        acb_mul(u, u, t, prec);
+        acb_neg(t, aa + 1);
+        acb_gamma(t, t, prec);
+        acb_mul(u, u, t, prec);
+
+        /* rgamma(a) * gamma(b-1) * z^(1-b) * v */
+        acb_rgamma(t, a, prec);
+        acb_mul(v, v, t, prec);
+        acb_gamma(t, aa + 1, prec);
+        acb_mul(v, v, t, prec);
+        acb_neg(t, aa + 1);
+        acb_pow(t, z, t, prec);
+        acb_mul(v, v, t, prec);
+
+        acb_add(res, u, v, prec);
+
+        acb_clear(t);
+        acb_clear(u);
+        acb_clear(v);
+        acb_clear(aa + 0);
+        acb_clear(aa + 1);
+        acb_clear(aa + 2);
+    }
+}
+
+void
+acb_hypgeom_u(acb_t res, const acb_t a, const acb_t b, const acb_t z, long prec)
+{
+    acb_t t;
+    acb_init(t);
+
+    acb_sub(t, a, b, prec);
+    acb_add_ui(t, t, 1, prec);
+
+    if ((acb_is_int(a) && arf_sgn(arb_midref(acb_realref(a))) <= 0) ||
+        (acb_is_int(t) && arf_sgn(arb_midref(acb_realref(t))) <= 0) ||
+        acb_hypgeom_u_use_asymp(z, prec))
+    {
+        acb_neg(t, a);
+        acb_pow(t, z, t, prec);
+        acb_hypgeom_u_asymp(res, a, b, z, -1, prec);
+        acb_mul(res, res, t, prec);
+    }
+    else
+    {
+        acb_hypgeom_u_1f1(res, a, b, z, prec);
+    }
+
+    acb_clear(t);
+}
+
diff --git a/doc/source/acb_hypgeom.rst b/doc/source/acb_hypgeom.rst
index 03a68a14..aeef7a7c 100644
--- a/doc/source/acb_hypgeom.rst
+++ b/doc/source/acb_hypgeom.rst
@@ -223,6 +223,31 @@ in terms of the following auxiliary quantities
     to evaluate `U(a,b,z)` to *prec* accurate bits (assuming that *a* and *b*
     are small).
 
+.. function:: void acb_hypgeom_u_1f1_series(acb_poly_t res, const acb_poly_t a, const acb_poly_t b, const acb_poly_t z, long len, long prec)
+
+    Computes `U(a,b,z)` as a power series truncated to length *len*,
+    given `a, b, z \in \mathbb{C}[[x]]`.
+    If `b \in \mathbb{Z}`, it computes one extra derivative and removes
+    the singularity. 
+    As currently implemented, the output is indeterminate if `b` is nonexact
+    and contains an integer.
+
+.. function:: void acb_hypgeom_u_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, long prec)
+
+    Computes `U(a,b,z)` as a sum of two convergent hypergeometric series.
+    If `b \in \mathbb{Z}`, it computes
+    the limit value via :func:`acb_hypgeom_u_1f1_series`.
+    As currently implemented, the output is indeterminate if `b` is nonexact
+    and contains an integer.
+
+.. function:: void acb_hypgeom_u(acb_t res, const acb_t a, const acb_t b, const acb_t z, long prec)
+
+    Computes `U(a,b,z)` using an automatic algorithm choice. The
+    function :func:`acb_hypgeom_u_asymp` is used
+    if `a` or `a-b+1` is a nonpositive integer (in which
+    case the asymptotic series terminates), or if *z* is sufficiently large.
+    Otherwise :func:`acb_hypgeom_u_1f1` is used.
+
 The error function
 -------------------------------------------------------------------------------