diff --git a/acb_dirichlet.h b/acb_dirichlet.h
index df353a91..63d3edb6 100644
--- a/acb_dirichlet.h
+++ b/acb_dirichlet.h
@@ -165,6 +165,8 @@ void acb_dirichlet_hardy_theta_series(acb_poly_t res, const acb_poly_t s, const
 void _acb_dirichlet_hardy_z_series(acb_ptr res, acb_srcptr s, slong slen, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec);
 void acb_dirichlet_hardy_z_series(acb_poly_t res, const acb_poly_t s, const dirichlet_group_t G, const dirichlet_char_t chi, slong len, slong prec);
 
+void acb_dirichlet_gram_point(arb_t res, const fmpz_t n, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
+
 /* Discrete Fourier Transform */
 
 void acb_dirichlet_dft_index(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec);
diff --git a/acb_dirichlet/gram_point.c b/acb_dirichlet/gram_point.c
new file mode 100644
index 00000000..950dc93d
--- /dev/null
+++ b/acb_dirichlet/gram_point.c
@@ -0,0 +1,174 @@
+/*
+    Copyright (C) 2019 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_dirichlet.h"
+
+/*
+Claim: the error is bounded by 1/256 if n <= 1 and (1/64) (log(n)/n) if n >= 2.
+
+A crude lower bound for g_n is 2 pi exp(W(n)), or 8*n/log(n) for n >= 8.
+
+We want to solve pi n = -t/2 log(2 pi/t) - t/2 - pi/8 + epsilon for t (= g_n).
+Using (47) in Brent [https://arxiv.org/abs/1609.03682], |epsilon| <= 1/(8 t) for for t >= 2.
+Also, for x >= 3, |f'(x)| < 0.5 where f(x) = exp(W(x)).
+
+Assume n >= 9, so that (n+1/8)/e >= 3.35. Then inverting gives
+
+t = 2 pi exp[W( [pi n - epsilon + pi/8] / (pi e) ) + 1]
+  = 2 pi e exp[W((n+1/8)/e  - epsilon / (pi e))]
+  = 2 pi e exp[W((n+1/8)/e)] + epsilon2, |epsilon2| <= 1/(8 t) <= (1/64) (log(n)/n)
+
+One can check 0 <= n <= 8 separately.
+*/
+static void
+gram_point_initial(arb_t x, const fmpz_t n, slong prec)
+{
+    arb_t pi, e;
+    mag_t b;
+
+    arb_init(pi);
+    arb_init(e);
+    mag_init(b);
+
+    arb_const_pi(pi, prec);
+    arb_const_e(e, prec);
+
+    /* x = 2*pi*exp(1 + W((n+1/8)/e)) */
+    arb_one(x);
+    arb_mul_2exp_si(x, x, -3);
+    arb_add_fmpz(x, x, n, prec);
+    arb_div(x, x, e, prec);
+    arb_lambertw(x, x, 0, prec);
+    arb_add_ui(x, x, 1, prec);
+    arb_exp(x, x, prec);
+    arb_mul(x, x, pi, prec);
+    arb_mul_2exp_si(x, x, 1);
+
+    if (fmpz_cmp_ui(n, 1) <= 0)
+    {
+        mag_set_ui_2exp_si(b, 1, -8);
+    }
+    else
+    {
+        mag_set_fmpz(b, n);
+        mag_log(b, b);
+        mag_div_fmpz(b, b, n);
+        mag_mul_2exp_si(b, b, -6);
+    }
+
+    arb_add_error_mag(x, b);
+
+    arb_clear(pi);
+    arb_clear(e);
+    mag_clear(b);
+}
+
+void
+acb_dirichlet_gram_point(arb_t res, const fmpz_t n, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
+{
+    slong asymp_accuracy;
+
+    /* Only implemented for n >= 0 and Riemann zeta. */
+    if (fmpz_sgn(n) < 0 || G != NULL || chi != NULL)
+    {
+        arb_indeterminate(res);
+        return;
+    }
+
+    asymp_accuracy = 2 * fmpz_bits(n);
+    asymp_accuracy = FLINT_MIN(asymp_accuracy, prec);
+
+    gram_point_initial(res, n, asymp_accuracy + 20);
+    asymp_accuracy = arb_rel_accuracy_bits(res);
+
+    if (asymp_accuracy < prec)
+    {
+        acb_struct tmp[2];
+        arb_t f, fprime, root;
+        mag_t C, r;
+        slong * steps;
+        slong wp, step;
+
+        acb_init(tmp);
+        acb_init(tmp + 1);
+        arb_init(f);
+        arb_init(fprime);
+        arb_init(root);
+        mag_init(C);
+        mag_init(r);
+        steps = flint_malloc(sizeof(slong) * FLINT_BITS);
+
+        step = 0;
+        steps[step] = prec * 1.05 + 10;
+
+        while (steps[step] / 2 > asymp_accuracy)
+        {
+            steps[step + 1] = steps[step] / 2;
+            step++;
+        }
+
+        arb_set(root, res);
+
+        /* theta''(x) <= C = 1/x, x >= 1 */
+        arb_get_mag_lower(C, root);
+        if (mag_cmp_2exp_si(C, 0) >= 0)
+            mag_inv(C, C);
+        else
+            mag_inf(C);
+
+        arb_set(root, res);
+
+        for ( ; step >= 0; step--)
+        {
+            wp = steps[step] + 10;
+            wp = FLINT_MAX(wp, arb_rel_accuracy_bits(root) + 10);
+
+            /* store radius, set root to the midpoint */
+            mag_set(r, arb_radref(root));
+            mag_zero(arb_radref(root));
+
+            acb_set_arb(tmp, root);
+            acb_dirichlet_hardy_theta(tmp, tmp, NULL, NULL, 2, wp);
+            arb_set(f, acb_realref(tmp));
+            arb_const_pi(acb_imagref(tmp), wp);
+            arb_submul_fmpz(f, acb_imagref(tmp), n, wp);
+
+            arb_set(fprime, acb_realref(tmp + 1));
+
+            /* f'([m+/-r]) = f'(m) +/- f''([m +/- r]) * r */
+            mag_mul(r, C, r);
+            arb_add_error_mag(fprime, r);
+            arb_div(f, f, fprime, wp);
+            arb_sub(root, root, f, wp);
+
+            /* Verify inclusion so that C is still valid. */
+            if (!arb_contains(res, root))
+            {
+                flint_printf("unexpected: no containment computing Gram point\n");
+                arb_set(root, res);
+                break;
+            }
+        }
+
+        arb_set(res, root);
+
+        acb_clear(tmp);
+        acb_clear(tmp + 1);
+        arb_clear(f);
+        arb_clear(fprime);
+        arb_clear(root);
+        mag_clear(C);
+        mag_clear(r);
+        flint_free(steps);
+    }
+
+    arb_set_round(res, res, prec);
+}
diff --git a/acb_dirichlet/test/t-gram_point.c b/acb_dirichlet/test/t-gram_point.c
new file mode 100644
index 00000000..a82a6d9f
--- /dev/null
+++ b/acb_dirichlet/test/t-gram_point.c
@@ -0,0 +1,94 @@
+/*
+    Copyright (C) 2019 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_dirichlet.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("gram_point....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        arb_t x1, x2, v1, v2, pin;
+        acb_t t;
+        fmpz_t n;
+        slong prec1, prec2;
+
+        arb_init(x1);
+        arb_init(x2);
+        arb_init(v1);
+        arb_init(v2);
+        arb_init(pin);
+        acb_init(t);
+        fmpz_init(n);
+
+        fmpz_randtest(n, state, 500);
+        fmpz_abs(n, n);
+        prec1 = 2 + n_randtest(state) % 500;
+        prec2 = 2 + n_randtest(state) % 2000;
+
+        acb_dirichlet_gram_point(x1, n, NULL, NULL, prec1);
+        acb_dirichlet_gram_point(x2, n, NULL, NULL, prec2);
+
+        arb_const_pi(pin, FLINT_MAX(prec1, prec2) + 20);
+        arb_mul_fmpz(pin, pin, n, FLINT_MAX(prec1, prec2) + 20);
+
+        acb_set_arb(t, x1);
+        acb_dirichlet_hardy_theta(t, t, NULL, NULL, 1, prec1 + 20);
+        arb_set(v1, acb_realref(t));
+
+        acb_set_arb(t, x2);
+        acb_dirichlet_hardy_theta(t, t, NULL, NULL, 1, prec2 + 20);
+        arb_set(v2, acb_realref(t));
+
+        if (!arb_overlaps(x1, x2) || !arb_contains(v1, pin) || !arb_contains(v2, pin))
+        {
+            flint_printf("FAIL: overlap\n\n");
+            flint_printf("n = "); fmpz_print(n);
+            flint_printf("   prec1 = %wd  prec2 = %wd\n\n", prec1, prec2);
+            flint_printf("x1 = "); arb_printn(x1, 100, 0); flint_printf("\n\n");
+            flint_printf("x2 = "); arb_printn(x2, 100, 0); flint_printf("\n\n");
+            flint_printf("v1 = "); arb_printn(v1, 100, 0); flint_printf("\n\n");
+            flint_printf("v2 = "); arb_printn(v2, 100, 0); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        if (arb_rel_accuracy_bits(x1) < prec1 - 3 || arb_rel_accuracy_bits(x2) < prec2 - 3)
+        {
+            flint_printf("FAIL: accuracy\n\n");
+            flint_printf("n = "); fmpz_print(n);
+            flint_printf("   prec1 = %wd  prec2 = %wd\n\n", prec1, prec2);
+            flint_printf("acc(x1) = %wd, acc(x2) = %wd\n\n", arb_rel_accuracy_bits(x1), arb_rel_accuracy_bits(x2));
+            flint_printf("x1 = "); arb_printn(x1, 100, 0); flint_printf("\n\n");
+            flint_printf("x2 = "); arb_printn(x2, 100, 0); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        arb_clear(x1);
+        arb_clear(x2);
+        arb_clear(v1);
+        arb_clear(v2);
+        arb_clear(pin);
+        acb_clear(t);
+        fmpz_clear(n);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
diff --git a/doc/source/acb_dirichlet.rst b/doc/source/acb_dirichlet.rst
index 7c8f0705..173a6f03 100644
--- a/doc/source/acb_dirichlet.rst
+++ b/doc/source/acb_dirichlet.rst
@@ -613,3 +613,10 @@ Currently, these methods require *chi* to be a primitive character.
 
     Sets *res* to the power series `Z(t)` where *t* is a given power series, truncating the result to length *len*.
 
+.. function:: void acb_dirichlet_gram_point(arb_t res, const fmpz_t n, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
+
+    Sets *res* to the *n*-th Gram point `g_n`, defined as the solution of
+    `\theta(g_n) = \pi n`. Currently only the Gram points corresponding to the
+    Riemann zeta function are supported and *G* and *chi* must both be set to
+    *NULL*. Requires `n \ge 0`.
+