diff --git a/acb_poly/gamma_series.c b/acb_poly/gamma_series.c
index 05ac42f7..8852df53 100644
--- a/acb_poly/gamma_series.c
+++ b/acb_poly/gamma_series.c
@@ -45,7 +45,7 @@ bsplit(acb_ptr Q, acb_ptr T, const acb_t z, long a, long b, long num, long prec)
             acb_set(Q, z);
             if (num > 1) acb_one(Q + 1);
             if (num > 2) acb_zero(Q + 2);
-        }             
+        }
         else
         {   /* (z + t)^2 */
             acb_mul(Q, z, z, prec);  /* TODO: precompute */
diff --git a/acb_poly/lgamma_series.c b/acb_poly/lgamma_series.c
index 6860137d..dc250c09 100644
--- a/acb_poly/lgamma_series.c
+++ b/acb_poly/lgamma_series.c
@@ -109,14 +109,57 @@ _acb_poly_lgamma_series(acb_ptr res, acb_srcptr h, long hlen, long len, long pre
     acb_init(zr);
 
     /* use Stirling series */
-    acb_gamma_stirling_choose_param(&reflect, &r, &n, h, 0, 0, wp);
-    acb_add_ui(zr, h, r, wp);
-    _acb_poly_gamma_stirling_eval(u, zr, n, len, wp);
+    acb_gamma_stirling_choose_param(&reflect, &r, &n, h, 1, 0, wp);
 
-    if (r != 0)
+    if (reflect)
     {
-        _log_rising_ui_series(t, h, r, len, wp);
-        _acb_vec_sub(u, u, t, len, wp);
+        /* log gamma(h+x) = log rf(1-(h+x), r) - log gamma(1-(h+x)+r) - log sin(pi (h+x)) + log(pi) */
+        if (r != 0) /* otherwise t = 0 */
+        {
+            acb_sub_ui(u, h, 1, wp);
+            acb_neg(u, u);
+            _log_rising_ui_series(t, u, r, len, wp);
+            for (i = 1; i < len; i += 2)
+                acb_neg(t + i, t + i);
+        }
+
+        acb_sub_ui(u, h, 1, wp);
+        acb_neg(u, u);
+        acb_add_ui(zr, u, r, wp);
+        _acb_poly_gamma_stirling_eval(u, zr, n, len, wp);
+        for (i = 1; i < len; i += 2)
+            acb_neg(u + i, u + i);
+
+        _acb_vec_sub(t, t, u, len, wp);
+
+        /* log(sin) is unstable with large imaginary parts;
+           cot_pi is implemented in a numerically stable way */
+        acb_set(u, h);
+        acb_one(u + 1);
+        _acb_poly_cot_pi_series(u, u, 2, len - 1, wp);
+        _acb_poly_integral(u, u, len, wp);
+        acb_const_pi(u, wp);
+        _acb_vec_scalar_mul(u + 1, u + 1, len - 1, u, wp);
+        acb_log_sin_pi(u, h, wp);
+
+        _acb_vec_sub(u, t, u, len, wp);
+
+        acb_const_pi(t, wp); /* todo: constant for log pi */
+        acb_log(t, t, wp);
+        acb_add(u, u, t, wp);
+    }
+    else
+    {
+        /* log gamma(x) = log gamma(x+r) - log rf(x,r) */
+
+        acb_add_ui(zr, h, r, wp);
+        _acb_poly_gamma_stirling_eval(u, zr, n, len, wp);
+
+        if (r != 0)
+        {
+            _log_rising_ui_series(t, h, r, len, wp);
+            _acb_vec_sub(u, u, t, len, wp);
+        }
     }
 
     /* compose with nonconstant part */
diff --git a/acb_poly/test/t-lgamma_series.c b/acb_poly/test/t-lgamma_series.c
index 114cbc34..47b07200 100644
--- a/acb_poly/test/t-lgamma_series.c
+++ b/acb_poly/test/t-lgamma_series.c
@@ -35,6 +35,35 @@ int main()
 
     flint_randinit(state);
 
+    /* special accuracy test case */
+    {
+        acb_poly_t a;
+        acb_t c;
+
+        acb_init(c);
+        acb_poly_init(a);
+
+        arb_set_str(acb_realref(c), "-20.25", 53);
+        arb_set_str(acb_imagref(c), "1e1000", 53);
+
+        acb_poly_set_coeff_acb(a, 0, c);
+        acb_poly_set_coeff_si(a, 1, 1);
+
+        acb_poly_lgamma_series(a, a, 3, 53);
+
+        if (acb_rel_accuracy_bits(a->coeffs) < 40 ||
+            acb_rel_accuracy_bits(a->coeffs + 1) < 40 ||
+            acb_rel_accuracy_bits(a->coeffs + 2) < 40)
+        {
+            printf("FAIL: accuracy (reflection formula)\n\n");
+            acb_poly_printd(a, 15); printf("\n\n");
+            abort();
+        }
+
+        acb_poly_clear(a);
+        acb_clear(c);
+    }
+
     for (iter = 0; iter < 500; iter++)
     {
         long m, n1, n2, rbits1, rbits2, rbits3;
@@ -53,9 +82,9 @@ int main()
         acb_poly_init(c);
         acb_poly_init(d);
 
-        acb_poly_randtest(a, state, m, rbits1, 3);
-        acb_poly_randtest(b, state, m, rbits1, 3);
-        acb_poly_randtest(c, state, m, rbits1, 3);
+        acb_poly_randtest(a, state, m, rbits1, 10);
+        acb_poly_randtest(b, state, m, rbits1, 10);
+        acb_poly_randtest(c, state, m, rbits1, 10);
 
         acb_poly_lgamma_series(b, a, n1, rbits2);
         acb_poly_lgamma_series(c, a, n2, rbits3);
diff --git a/arb_poly/lgamma_series.c b/arb_poly/lgamma_series.c
index 8ab99d9f..00f80978 100644
--- a/arb_poly/lgamma_series.c
+++ b/arb_poly/lgamma_series.c
@@ -61,6 +61,12 @@ _arb_poly_lgamma_series(arb_ptr res, arb_srcptr h, long hlen, long len, long pre
     arb_t zr;
     arb_ptr t, u;
 
+    if (!arb_is_positive(h))
+    {
+        _arb_vec_indeterminate(res, len);
+        return;
+    }
+
     hlen = FLINT_MIN(hlen, len);
     wp = prec + FLINT_BIT_COUNT(prec);