diff --git a/acb_hypgeom.h b/acb_hypgeom.h
index b40424d7..57d646b3 100644
--- a/acb_hypgeom.h
+++ b/acb_hypgeom.h
@@ -41,6 +41,8 @@ int acb_hypgeom_gamma_taylor(acb_t res, const acb_t x, int reciprocal, slong pre
 void acb_hypgeom_gamma(acb_t y, const acb_t x, slong prec);
 void acb_hypgeom_rgamma(acb_t y, const acb_t x, slong prec);
 
+void acb_hypgeom_lgamma(acb_t y, const acb_t x, slong prec);
+
 void acb_hypgeom_pfq_bound_factor(mag_t C,
     acb_srcptr a, slong p, acb_srcptr b, slong q, const acb_t z, ulong n);
 
diff --git a/acb_hypgeom/lgamma.c b/acb_hypgeom/lgamma.c
new file mode 100644
index 00000000..392e24e6
--- /dev/null
+++ b/acb_hypgeom/lgamma.c
@@ -0,0 +1,262 @@
+/*
+    Copyright (C) 2014, 2015, 2021 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_hypgeom.h"
+#include "arb_hypgeom.h"
+
+void acb_hypgeom_gamma_stirling_choose_param(int * reflect, slong * r, slong * n,
+    const acb_t z, int use_reflect, int digamma, slong prec);
+
+void acb_hypgeom_gamma_stirling_inner(acb_t s, const acb_t z, slong N, slong prec);
+
+static double
+want_taylor(double x, double y, slong prec)
+{
+    if (y < 0.0) y = -y;
+    if (x < 0.0) x = -2.0 * x;
+
+    if ((prec < 128 && y > 4.0) || (prec < 256 && y > 5.0) ||
+        (prec < 512 && y > 8.0) || (prec < 1024 && y > 9.0) || y > 10.0)
+    {
+        return 0;
+    }
+
+    if (x * (1.0 + 0.75 * y) > 8 + 0.15 * prec)
+    {
+        return 0;
+    }
+
+    return 1;
+}
+
+/* Linear fit on [0.5, 1.5] for
+    lambda x: findroot(lambda y: im(loggamma(x+1j*y)) - (n+0.5)*pi */
+
+static const double Atab[] = {
+    4.5835631239879990091,
+    6.4037921417161376741,
+    7.9938623618272375768,
+    9.4449131928216797873,
+    10.802608819487725856,
+    12.0918817314347272,
+};
+
+static const double Btab[] = {
+    -1.1432582881376479127,
+    -0.86248117216701645437,
+    -0.75778990135448922722,
+    -0.69734688055939976228,
+    -0.65626499937495627271,
+    -0.62578331900739100617,
+};
+
+void
+_arb_const_log_pi(arb_t t, slong prec)
+{
+    arb_const_pi(t, prec + 2);
+    arb_log(t, t, prec);
+}
+
+ARB_DEF_CACHED_CONSTANT(arb_const_log_pi, _arb_const_log_pi)
+
+int
+acb_hypgeom_lgamma_taylor(acb_t res, const acb_t z, slong prec)
+{
+    double x, y, acc;
+    slong k, r, wp;
+    acb_t t, u;
+    int reflect;
+
+    /* Assume xerr, yerr <= 1/16 */
+    if (mag_cmp_2exp_si(arb_radref(acb_realref(z)), -4) > 0)
+        return 0;
+
+    if (mag_cmp_2exp_si(arb_radref(acb_imagref(z)), -4) > 0)
+        return 0;
+
+    acc = acb_rel_accuracy_bits(z);
+    acc = FLINT_MAX(acc, 0);
+    wp = FLINT_MIN(prec, acc + 20);
+    wp = FLINT_MAX(wp, 2);
+
+    /* x, y plus eventual rounding error */
+    x = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_NEAR);
+    y = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_NEAR);
+
+    if (!want_taylor(x, y, wp))
+        return 0;
+
+    acb_init(t);
+    acb_init(u);
+
+    /* Reduce real part to (approximately) [0.5, 1.5]. */
+    r = floor(x - 0.5);
+
+    /* Reflection formula is slower but improves accuracy. */
+    reflect = (x < -3.0);
+
+    if (reflect)
+    {
+        acb_neg(u, z);
+        acb_add_si(u, u, 2 + r, 2 * prec + 10);
+        x = 2.0 + r - x;
+        y = -y;
+    }
+    else
+    {
+        acb_sub_si(u, z, r, 2 * prec + 10);
+        x = x - r;
+    }
+
+    for (k = 0; k < 6; k++)
+    {
+        if (fabs(y) <= Atab[k] + Btab[k] * x)
+        {
+            if (!acb_hypgeom_gamma_taylor(t, u, 1, wp))
+            {
+                acb_clear(t);
+                acb_clear(u);
+                return 0;
+            }
+
+            if (k % 2 == 0)
+            {
+                acb_log(t, t, wp);
+                acb_neg(t, t);
+            }
+            else
+            {
+                acb_neg(t, t);
+                acb_log(t, t, wp);
+                acb_neg(t, t);
+            }
+
+            if (k != 0)
+            {
+                arb_t pi;
+                arb_init(pi);
+                arb_const_pi(pi, wp);
+                arb_addmul_si(acb_imagref(t), pi, (y > 0) ? k : -k, wp);
+                arb_clear(pi);
+            }
+
+            if (reflect)
+            {
+                acb_t v;
+                acb_init(v);
+
+                /* loggamma(x) = log(pi) - lsin(x) - loggamma(2+r-x) - logrf(2+r-x, -r-1) */
+
+                acb_hypgeom_log_rising_ui(v, u, -r-1, wp);
+                acb_log_sin_pi(res, z, wp);
+                acb_add(res, res, v, wp);
+                acb_add(res, res, t, wp);
+                acb_neg(res, res);
+
+                arb_const_log_pi(acb_realref(t), wp);
+                arb_zero(acb_imagref(t));
+                acb_add(res, res, t, prec);
+
+                acb_clear(v);
+            }
+            else if (r == 0)
+            {
+                acb_set_round(res, t, prec);
+            }
+            else if (r > 0)
+            {
+                acb_hypgeom_log_rising_ui(res, u, r, wp);
+                acb_add(res, res, t, prec);
+            }
+            else
+            {
+                acb_hypgeom_log_rising_ui(res, z, -r, wp);
+                acb_sub(res, t, res, prec);
+            }
+
+            acb_clear(t);
+            acb_clear(u);
+            return 1;
+        }
+    }
+
+    acb_clear(t);
+    acb_clear(u);
+    return 0;
+}
+
+void
+acb_hypgeom_lgamma(acb_t y, const acb_t x, slong prec)
+{
+    int reflect;
+    slong r, n, wp;
+    acb_t t, u, v;
+    double acc;
+
+    if (acb_is_real(x) && arb_is_positive(acb_realref(x)))
+    {
+        arb_hypgeom_lgamma(acb_realref(y), acb_realref(x), prec);
+        arb_zero(acb_imagref(y));
+        return;
+    }
+
+    if (acb_hypgeom_lgamma_taylor(y, x, prec))
+        return;
+
+    acc = acb_rel_accuracy_bits(x);
+    acc = FLINT_MAX(acc, 0);
+    wp = FLINT_MIN(prec, acc + 20);
+    wp = FLINT_MAX(wp, 2);
+    wp = wp + FLINT_BIT_COUNT(wp);
+
+    acb_hypgeom_gamma_stirling_choose_param(&reflect, &r, &n, x, 1, 0, wp);
+
+    acb_init(t);
+    acb_init(u);
+    acb_init(v);
+
+    if (reflect)
+    {
+        /* log gamma(x) = log rf(1-x, r) - log gamma(1-x+r) - log sin(pi x) + log(pi) */
+        acb_sub_ui(u, x, 1, wp);
+        acb_neg(u, u);
+
+        acb_hypgeom_log_rising_ui(t, u, r, wp);
+
+        acb_add_ui(u, u, r, wp);
+        acb_hypgeom_gamma_stirling_inner(v, u, n, wp);
+        acb_sub(t, t, v, wp);
+
+        acb_log_sin_pi(u, x, wp);
+        acb_sub(t, t, u, wp);
+
+        arb_const_log_pi(acb_realref(u), wp);
+        arb_zero(acb_imagref(u));
+
+        acb_add(y, t, u, wp);
+    }
+    else
+    {
+        /* log gamma(x) = log gamma(x+r) - log rf(x,r) */
+        acb_add_ui(t, x, r, wp);
+        acb_hypgeom_gamma_stirling_inner(u, t, n, wp);
+        acb_hypgeom_log_rising_ui(t, x, r, wp);
+        acb_sub(y, u, t, prec);
+    }
+
+    if (!acb_is_finite(y))
+        acb_indeterminate(y);
+
+    acb_clear(t);
+    acb_clear(u);
+    acb_clear(v);
+}
+
diff --git a/acb_hypgeom/test/t-lgamma.c b/acb_hypgeom/test/t-lgamma.c
new file mode 100644
index 00000000..e961ee60
--- /dev/null
+++ b/acb_hypgeom/test/t-lgamma.c
@@ -0,0 +1,120 @@
+/*
+    Copyright (C) 2021 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_hypgeom.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("lgamma....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
+    {
+        acb_t z, s1, s2, a, b;
+        slong prec, ebits, prec2;
+
+        prec = 2 + n_randint(state, 200);
+
+        if (n_randint(state, 10) == 0)
+            prec = 2 + n_randint(state, 1000);
+
+        if (n_randint(state, 10) == 0)
+            ebits = 100;
+        else
+            ebits = 10;
+        ebits = 2;
+
+        prec2 = 2 + n_randint(state, 200);
+
+        acb_init(z);
+        acb_init(s1);
+        acb_init(s2);
+        acb_init(a);
+        acb_init(b);
+
+        acb_randtest(z, state, prec, ebits);
+        acb_randtest(s1, state, prec, 10);
+        acb_randtest(s2, state, prec, 10);
+
+        if (n_randint(state, 2))
+        {
+            acb_hypgeom_lgamma(s1, z, prec);
+        }
+        else
+        {
+            acb_set(s1, z);
+            acb_hypgeom_lgamma(s1, s1, prec);
+        }
+
+        acb_add_ui(s2, z, 1, prec2);
+        acb_hypgeom_lgamma(s2, s2, prec2);
+        acb_log(a, z, prec2);
+        acb_sub(s2, s2, a, prec2);
+
+        if (!acb_overlaps(s1, s2))
+        {
+            flint_printf("FAIL\n\n");
+            flint_printf("prec = %wd\n\n", prec);
+            flint_printf("z = "); acb_printn(z, 1000, 0); flint_printf("\n\n");
+            flint_printf("s1 = "); acb_printn(s1, 1000, 0); flint_printf("\n\n");
+            flint_printf("s2 = "); acb_printn(s2, 1000, 0); flint_printf("\n\n");
+            acb_sub(s1, s1, s2, prec2);
+            flint_printf("s1 - s2 = "); acb_printd(s1, 1000); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        acb_get_mid(a, z);
+
+        if (n_randint(state, 2))
+        {
+            arf_set_mag(arb_midref(acb_realref(b)), arb_radref(acb_realref(z)));
+            arf_set_mag(arb_midref(acb_imagref(b)), arb_radref(acb_imagref(z)));
+
+            if (n_randint(state, 2))
+                acb_neg(b, b);
+            if (n_randint(state, 2))
+                acb_conj(b, b);
+
+            acb_add(a, a, b, prec);
+        }
+
+        acb_hypgeom_lgamma(s2, a, prec);
+
+        if (!acb_overlaps(s1, s2))
+        {
+            flint_printf("FAIL (2)\n\n");
+            flint_printf("prec = %wd\n\n", prec);
+            flint_printf("z = "); acb_printn(z, 1000, 0); flint_printf("\n\n");
+            flint_printf("a = "); acb_printn(a, 1000, 0); flint_printf("\n\n");
+            flint_printf("s1 = "); acb_printn(s1, 1000, 0); flint_printf("\n\n");
+            flint_printf("s2 = "); acb_printn(s2, 1000, 0); flint_printf("\n\n");
+            acb_sub(s1, s1, s2, prec2);
+            flint_printf("s1 - s2 = "); acb_printd(s1, 1000); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        acb_clear(z);
+        acb_clear(s1);
+        acb_clear(s2);
+        acb_clear(a);
+        acb_clear(b);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
diff --git a/arb_hypgeom.h b/arb_hypgeom.h
index dbd2031a..92160faa 100644
--- a/arb_hypgeom.h
+++ b/arb_hypgeom.h
@@ -59,6 +59,8 @@ int arb_hypgeom_gamma_taylor(arb_t res, const arb_t x, int reciprocal, slong pre
 void arb_hypgeom_gamma(arb_t y, const arb_t x, slong prec);
 void arb_hypgeom_rgamma(arb_t y, const arb_t x, slong prec);
 
+void arb_hypgeom_lgamma(arb_t y, const arb_t x, slong prec);
+
 void arb_hypgeom_pfq(arb_t res, arb_srcptr a, slong p, arb_srcptr b, slong q,
     const arb_t z, int regularized, slong prec);
 
diff --git a/arb_hypgeom/lgamma.c b/arb_hypgeom/lgamma.c
new file mode 100644
index 00000000..995e6465
--- /dev/null
+++ b/arb_hypgeom/lgamma.c
@@ -0,0 +1,75 @@
+/*
+    Copyright (C) 2021 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 "arb_hypgeom.h"
+
+void arb_hypgeom_gamma_stirling_choose_param(int * reflect, slong * r, slong * n, const arb_t x, int use_reflect, int digamma, slong prec);
+int arb_hypgeom_gamma_exact(arb_t res, const arb_t x, int reciprocal, slong prec);
+void arb_hypgeom_gamma_stirling_inner(arb_t s, const arb_t z, slong N, slong prec);
+
+void
+arb_hypgeom_lgamma_stirling(arb_t y, const arb_t x, slong prec)
+{
+    int reflect;
+    slong r, n, wp;
+    arb_t t, u;
+    double acc;
+
+    /* todo: for large x (if exact or accurate enough), increase precision */
+    acc = arb_rel_accuracy_bits(x);
+    acc = FLINT_MAX(acc, 0);
+    wp = FLINT_MIN(prec, acc + 20);
+    wp = FLINT_MAX(wp, 2);
+    wp = wp + FLINT_BIT_COUNT(wp);
+
+    arb_hypgeom_gamma_stirling_choose_param(&reflect, &r, &n, x, 0, 0, wp);
+
+    arb_init(t);
+    arb_init(u);
+
+    /* log(gamma(x)) = log(gamma(x+r)) - log(rf(x,r)) */
+    arb_init(t);
+    arb_init(u);
+
+    arb_add_ui(t, x, r, wp);
+    arb_hypgeom_gamma_stirling_inner(u, t, n, wp);
+    arb_hypgeom_rising_ui_rec(t, x, r, wp);
+    arb_log(t, t, wp);
+    arb_sub(y, u, t, prec);
+
+    arb_clear(t);
+    arb_clear(u);
+}
+
+void
+arb_hypgeom_lgamma(arb_t res, const arb_t x, slong prec)
+{
+    if (!arb_is_positive(x) || !arb_is_finite(x))
+    {
+        arb_indeterminate(res);
+        return;
+    }
+
+    if (arb_hypgeom_gamma_exact(res, x, 0, prec))
+    {
+        arb_log(res, res, prec);
+        return;
+    }
+
+    if (arb_hypgeom_gamma_taylor(res, x, 0, prec))
+    {
+        arb_log(res, res, prec);
+        return;
+    }
+
+    arb_hypgeom_lgamma_stirling(res, x, prec);
+}
+
diff --git a/arb_hypgeom/test/t-lgamma.c b/arb_hypgeom/test/t-lgamma.c
new file mode 100644
index 00000000..dcbb6561
--- /dev/null
+++ b/arb_hypgeom/test/t-lgamma.c
@@ -0,0 +1,118 @@
+/*
+    Copyright (C) 2021 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 "arb_hypgeom.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("lgamma....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 10000 * arb_test_multiplier(); iter++)
+    {
+        arb_t z, s1, s2, a, b;
+        slong prec, ebits, prec2;
+
+        prec = 2 + n_randint(state, 200);
+
+        if (n_randint(state, 10) == 0)
+            prec = 2 + n_randint(state, 1000);
+
+        if (n_randint(state, 10) == 0)
+            ebits = 100;
+        else
+            ebits = 10;
+        ebits = 2;
+
+        prec2 = 2 + n_randint(state, 200);
+
+        arb_init(z);
+        arb_init(s1);
+        arb_init(s2);
+        arb_init(a);
+        arb_init(b);
+
+        arb_randtest(z, state, prec, ebits);
+        arb_randtest(s1, state, prec, 10);
+        arb_randtest(s2, state, prec, 10);
+
+        if (n_randint(state, 2))
+        {
+            arb_hypgeom_lgamma(s1, z, prec);
+        }
+        else
+        {
+            arb_set(s1, z);
+            arb_hypgeom_lgamma(s1, s1, prec);
+        }
+
+        arb_add_ui(s2, z, 1, prec2);
+        arb_hypgeom_lgamma(s2, s2, prec2);
+        arb_log(a, z, prec2);
+        arb_sub(s2, s2, a, prec2);
+
+        if (!arb_overlaps(s1, s2))
+        {
+            flint_printf("FAIL\n\n");
+            flint_printf("prec = %wd\n\n", prec);
+            flint_printf("z = "); arb_printn(z, 1000, 0); flint_printf("\n\n");
+            flint_printf("s1 = "); arb_printn(s1, 1000, 0); flint_printf("\n\n");
+            flint_printf("s2 = "); arb_printn(s2, 1000, 0); flint_printf("\n\n");
+            arb_sub(s1, s1, s2, prec2);
+            flint_printf("s1 - s2 = "); arb_printd(s1, 1000); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        arb_set(a, z);
+        mag_zero(arb_radref(a));
+
+        if (n_randint(state, 2))
+        {
+            arf_set_mag(arb_midref(b), arb_radref(z));
+
+            if (n_randint(state, 2))
+                arb_neg(b, b);
+
+            arb_add(a, a, b, prec);
+        }
+
+        arb_hypgeom_lgamma(s2, a, prec);
+
+        if (!arb_overlaps(s1, s2))
+        {
+            flint_printf("FAIL (2)\n\n");
+            flint_printf("prec = %wd\n\n", prec);
+            flint_printf("z = "); arb_printn(z, 1000, 0); flint_printf("\n\n");
+            flint_printf("a = "); arb_printn(a, 1000, 0); flint_printf("\n\n");
+            flint_printf("s1 = "); arb_printn(s1, 1000, 0); flint_printf("\n\n");
+            flint_printf("s2 = "); arb_printn(s2, 1000, 0); flint_printf("\n\n");
+            arb_sub(s1, s1, s2, prec2);
+            flint_printf("s1 - s2 = "); arb_printd(s1, 1000); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        arb_clear(z);
+        arb_clear(s1);
+        arb_clear(s2);
+        arb_clear(a);
+        arb_clear(b);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
diff --git a/doc/source/acb_hypgeom.rst b/doc/source/acb_hypgeom.rst
index 17849ad1..be983949 100644
--- a/doc/source/acb_hypgeom.rst
+++ b/doc/source/acb_hypgeom.rst
@@ -99,16 +99,21 @@ Gamma function
     and returns 0. If *reciprocal* is set, the reciprocal gamma function is
     computed instead.
 
-.. function:: void acb_hypgeom_gamma(acb_t y, const acb_t x, slong prec)
+.. function:: void acb_hypgeom_gamma(acb_t res, const acb_t x, slong prec)
 
     Sets *res* to the gamma function of *x* computed using a default
     algorithm choice.
 
-.. function:: void acb_hypgeom_rgamma(acb_t y, const acb_t x, slong prec)
+.. function:: void acb_hypgeom_rgamma(acb_t res, const acb_t x, slong prec)
 
     Sets *res* to the reciprocal gamma function of *x* computed using a default
     algorithm choice.
 
+.. function:: void acb_hypgeom_lgamma(acb_t res, const acb_t x, slong prec)
+
+    Sets *res* to the principal branch of the log-gamma function of *x*
+    computed using a default algorithm choice.
+
 
 Convergent series
 -------------------------------------------------------------------------------
diff --git a/doc/source/arb_hypgeom.rst b/doc/source/arb_hypgeom.rst
index 777cb64e..998497b7 100644
--- a/doc/source/arb_hypgeom.rst
+++ b/doc/source/arb_hypgeom.rst
@@ -103,16 +103,20 @@ Gamma function
     and returns 0. If *reciprocal* is set, the reciprocal gamma function is
     computed instead.
 
-.. function:: void arb_hypgeom_gamma(arb_t y, const arb_t x, slong prec)
+.. function:: void arb_hypgeom_gamma(arb_t res, const arb_t x, slong prec)
 
     Sets *res* to the gamma function of *x* computed using a default
     algorithm choice.
 
-.. function:: void arb_hypgeom_rgamma(arb_t y, const arb_t x, slong prec)
+.. function:: void arb_hypgeom_rgamma(arb_t res, const arb_t x, slong prec)
 
     Sets *res* to the reciprocal gamma function of *x* computed using a default
     algorithm choice.
 
+.. function:: void arb_hypgeom_lgamma(arb_t res, const arb_t x, slong prec)
+
+    Sets *res* to the log-gamma function of *x* computed using a default
+    algorithm choice.
 
 Binomial coefficients
 -------------------------------------------------------------------------------