add complex digamma function

fmpcb.h

@@ -531,6 +531,7 @@ void fmpcb_cos(fmpcb_t r, const fmpcb_t z, long prec);
 void fmpcb_sin_cos(fmpcb_t s, fmpcb_t c, const fmpcb_t z, long prec);
 
 void fmpcb_sin_pi(fmpcb_t r, const fmpcb_t z, long prec);
+void fmpcb_sin_cos_pi(fmpcb_t s, fmpcb_t c, const fmpcb_t z, long prec);
 
 void fmpcb_pow(fmpcb_t r, const fmpcb_t x, const fmpcb_t y, long prec);
 
@@ -543,6 +544,8 @@ void fmpcb_gamma(fmpcb_t y, const fmpcb_t x, long prec);
 void fmpcb_rgamma(fmpcb_t y, const fmpcb_t x, long prec);
 void fmpcb_lgamma(fmpcb_t y, const fmpcb_t x, long prec);
 
+void fmpcb_digamma(fmpcb_t y, const fmpcb_t x, long prec);
+
 void fmpcb_rising_ui(fmpcb_t y, const fmpcb_t x, ulong n, long prec);
 
 

fmpcb/digamma.c

@@ -0,0 +1,71 @@
+/*=============================================================================
+
+    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) 2013 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb.h"
+#include "gamma.h"
+
+void
+fmpcb_digamma(fmpcb_t y, const fmpcb_t x, long prec)
+{
+    int reflect;
+    long r, n, wp;
+    fmpcb_t t, u, v;
+
+    wp = prec + FLINT_BIT_COUNT(prec);
+
+    gamma_stirling_choose_param_fmpcb(&reflect, &r, &n, x, 1, 1, wp);
+
+    fmpcb_init(t);
+    fmpcb_init(u);
+    fmpcb_init(v);
+
+    /* psi(x) = psi((1-x)+r) - h(1-x,r) - pi*cot(pi*x) */
+    if (reflect)
+    {
+        fmpcb_sub_ui(t, x, 1, wp);
+        fmpcb_neg(t, t);
+        fmpcb_sin_cos_pi(u, v, x, wp);
+        fmpcb_div(v, v, u, wp);  /* todo: write a complex cotangent function */
+        fmprb_const_pi(fmpcb_realref(u), wp);
+        fmpcb_mul_fmprb(v, v, fmpcb_realref(u), wp);
+        gamma_harmonic_sum_fmpcb_ui_bsplit(u, t, r, wp);
+        fmpcb_add(v, v, u, wp);
+        fmpcb_add_ui(t, t, r, wp);
+        gamma_stirling_eval_series_fmpcb(u, t, n, 1, wp);
+        fmpcb_sub(y, u, v, wp);
+    }
+    else
+    {
+        fmpcb_add_ui(t, x, r, wp);
+        gamma_stirling_eval_series_fmpcb(u, t, n, 1, wp);
+        gamma_harmonic_sum_fmpcb_ui_bsplit(t, x, r, wp);
+        fmpcb_sub(y, u, t, prec);
+    }
+
+    fmpcb_clear(t);
+    fmpcb_clear(u);
+    fmpcb_clear(v);
+}
+

fmpcb/sin_cos_pi.c

@@ -0,0 +1,59 @@
+/*=============================================================================
+
+    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 "fmpcb.h"
+
+void
+fmpcb_sin_cos_pi(fmpcb_t s, fmpcb_t c, const fmpcb_t z, long prec)
+{
+#define a fmpcb_realref(z)
+#define b fmpcb_imagref(z)
+
+    fmprb_t sa, ca, sb, cb;
+
+    fmprb_init(sa);
+    fmprb_init(ca);
+    fmprb_init(sb);
+    fmprb_init(cb);
+
+    fmprb_const_pi(sb, prec);
+    fmprb_mul(sb, sb, b, prec);
+
+    fmprb_sin_cos_pi(sa, ca, a, prec);
+    fmprb_sinh_cosh(sb, cb, sb, prec);
+
+    fmprb_mul(fmpcb_realref(s), sa, cb, prec);
+    fmprb_mul(fmpcb_imagref(s), sb, ca, prec);
+
+    fmprb_mul(fmpcb_realref(c), ca, cb, prec);
+    fmprb_mul(fmpcb_imagref(c), sa, sb, prec);
+    fmprb_neg(fmpcb_imagref(c), fmpcb_imagref(c));
+
+    fmprb_clear(sa);
+    fmprb_clear(ca);
+    fmprb_clear(sb);
+    fmprb_clear(cb);
+}
+

fmpcb/test/t-digamma.c

@@ -0,0 +1,89 @@
+/*=============================================================================
+
+    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, 2013 Fredrik Johansson
+
+******************************************************************************/
+
+#include "fmpcb.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("digamma....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 3000; iter++)
+    {
+        fmpcb_t a, b, c;
+        long prec1, prec2;
+
+        prec1 = 2 + n_randint(state, 1000);
+        prec2 = prec1 + 30;
+
+        fmpcb_init(a);
+        fmpcb_init(b);
+        fmpcb_init(c);
+
+        fmprb_randtest_precise(fmpcb_realref(a), state, 1 + n_randint(state, 1000), 3);
+        fmprb_randtest_precise(fmpcb_imagref(a), state, 1 + n_randint(state, 1000), 3);
+
+        fmpcb_digamma(b, a, prec1);
+        fmpcb_digamma(c, a, prec2);
+
+        if (!fmpcb_overlaps(b, c))
+        {
+            printf("FAIL: overlap\n\n");
+            printf("a = "); fmpcb_print(a); printf("\n\n");
+            printf("b = "); fmpcb_print(b); printf("\n\n");
+            printf("c = "); fmpcb_print(c); printf("\n\n");
+            abort();
+        }
+
+        /* check digamma(z+1) = digamma(z) + 1/z */
+        fmpcb_inv(c, a, prec1);
+        fmpcb_add(b, b, c, prec1);
+        fmpcb_add_ui(c, a, 1, prec1);
+        fmpcb_digamma(c, c, prec1);
+
+        if (!fmpcb_overlaps(b, c))
+        {
+            printf("FAIL: functional equation\n\n");
+            printf("a = "); fmpcb_print(a); printf("\n\n");
+            printf("b = "); fmpcb_print(b); printf("\n\n");
+            printf("c = "); fmpcb_print(c); printf("\n\n");
+            abort();
+        }
+
+        fmpcb_clear(a);
+        fmpcb_clear(b);
+        fmpcb_clear(c);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}

gamma.h

@@ -92,6 +92,10 @@ void gamma_harmonic_sum_fmprb_ui_bsplit_rectangular(fmprb_t y, const fmprb_t x,
 void gamma_harmonic_sum_fmprb_ui_bsplit_simple(fmprb_t y, const fmprb_t x, ulong n, long prec);
 void gamma_harmonic_sum_fmprb_ui_bsplit(fmprb_t y, const fmprb_t x, ulong n, long prec);
 
+void gamma_harmonic_sum_fmpcb_ui_bsplit_rectangular(fmpcb_t y, const fmpcb_t x, ulong n, ulong step, long prec);
+void gamma_harmonic_sum_fmpcb_ui_bsplit_simple(fmpcb_t y, const fmpcb_t x, ulong n, long prec);
+void gamma_harmonic_sum_fmpcb_ui_bsplit(fmpcb_t y, const fmpcb_t x, ulong n, long prec);
+
 void gamma_series_fmpq_hypgeom(fmprb_struct * res, const fmpq_t a, long len, long prec);
 void gamma_small_frac(fmprb_t y, unsigned int p, unsigned int q, long prec);
 

gamma/harmonic_sum_fmpcb_ui_bsplit.c

@@ -0,0 +1,36 @@
+/*=============================================================================
+
+    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) 2013 Fredrik Johansson
+
+******************************************************************************/
+
+#include "gamma.h"
+
+void
+gamma_harmonic_sum_fmpcb_ui_bsplit(fmpcb_t y, const fmpcb_t x, ulong n, long prec)
+{
+    if (n < 8 || fmpcb_bits(x) < prec / 8)
+        gamma_harmonic_sum_fmpcb_ui_bsplit_simple(y, x, n, prec);
+    else
+        gamma_harmonic_sum_fmpcb_ui_bsplit_rectangular(y, x, n, 0, prec);
+}
+

gamma/harmonic_sum_fmpcb_ui_bsplit_rectangular.c

@@ -0,0 +1,106 @@
+/*=============================================================================
+
+    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) 2013 Fredrik Johansson
+
+******************************************************************************/
+
+#include "gamma.h"
+#include "fmprb_poly.h"
+
+static void
+bsplit_step(fmpcb_t p, fmpcb_t q, const fmprb_struct * poly, ulong step,
+        const fmpcb_t x, ulong a, ulong b, long prec)
+{
+    if (b - a == step)
+    {
+        fmpcb_add_ui(p, x, a, prec);
+        _fmprb_poly_evaluate2_fmpcb_rectangular(q, p, poly, step + 1, p, prec);
+    }
+    else
+    {
+        fmpcb_t r, s;
+        ulong m;
+
+        fmpcb_init(r);
+        fmpcb_init(s);
+
+        m = a + ((b - a) / (2 * step)) * step;
+        bsplit_step(p, q, poly, step, x, a, m, prec);
+        bsplit_step(r, s, poly, step, x, m, b, prec);
+
+        fmpcb_mul(p, p, s, prec);
+        fmpcb_mul(r, r, q, prec);
+        fmpcb_add(p, p, r, prec);
+        fmpcb_mul(q, q, s, prec);
+
+        fmpcb_clear(r);
+        fmpcb_clear(s);
+    }
+}
+
+void
+gamma_harmonic_sum_fmpcb_ui_bsplit_rectangular(fmpcb_t y, const fmpcb_t x, ulong n, ulong step, long prec)
+{
+    fmpcb_t t, u;
+    ulong b, s1, s2;
+    long wp;
+
+    wp = FMPR_PREC_ADD(prec, FLINT_BIT_COUNT(n));
+
+    fmpcb_init(t);
+    fmpcb_init(u);
+
+    if (step == 0)
+    {
+        s1 = 2 * sqrt(n);
+        s2 = 10 * pow(prec, 0.25);
+        step = FLINT_MIN(s1, s2);
+    }
+    step = FLINT_MAX(step, 1);
+    b = (n / step) * step;
+
+    if (b != 0)
+    {
+        fmprb_struct * poly;
+        fmprb_struct xpoly[2];
+        poly = _fmprb_vec_init(step + 1);
+        fmprb_init(xpoly + 0);
+        fmprb_init(xpoly + 1);
+        fmprb_one(xpoly + 1);
+        _fmprb_poly_rfac_series_ui(poly, xpoly, 2, step, step + 1, wp);
+        bsplit_step(t, u, poly, step, x, 0, b, wp);
+        fmpcb_div(t, t, u, wp);
+        _fmprb_vec_clear(poly, step + 1);
+    }
+    else
+    {
+        fmpcb_zero(t);
+    }
+
+    fmpcb_add_ui(u, x, b, wp);
+    gamma_harmonic_sum_fmpcb_ui_bsplit(u, u, n - b, wp);
+    fmpcb_add(y, t, u, prec);
+
+    fmpcb_clear(t);
+    fmpcb_clear(u);
+}
+

gamma/harmonic_sum_fmpcb_ui_bsplit_simple.c

@@ -0,0 +1,101 @@
+/*=============================================================================
+
+    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) 2013 Fredrik Johansson
+
+******************************************************************************/
+
+#include "gamma.h"
+
+static void
+bsplit(fmpcb_t p, fmpcb_t q, const fmpcb_t x, ulong a, ulong b, long prec)
+{
+    if (b - a < 8)
+    {
+        ulong k;
+        fmpcb_t t;
+
+        fmpcb_one(p);
+        fmpcb_add_ui(q, x, a, prec);
+
+        fmpcb_init(t);
+
+        for (k = a + 1; k < b; k++)
+        {
+            fmpcb_add_ui(t, x, k, prec);
+            fmpcb_mul(p, p, t, prec);
+            fmpcb_add(p, p, q, prec);
+            fmpcb_mul(q, q, t, prec);
+        }
+
+        fmpcb_clear(t);
+    }
+    else
+    {
+        fmpcb_t r, s;
+        ulong m;
+
+        fmpcb_init(r);
+        fmpcb_init(s);
+
+        m = a + (b - a) / 2;
+        bsplit(p, q, x, a, m, prec);
+        bsplit(r, s, x, m, b, prec);
+
+        fmpcb_mul(p, p, s, prec);
+        fmpcb_mul(r, r, q, prec);
+        fmpcb_add(p, p, r, prec);
+        fmpcb_mul(q, q, s, prec);
+
+        fmpcb_clear(r);
+        fmpcb_clear(s);
+    }
+}
+
+void
+gamma_harmonic_sum_fmpcb_ui_bsplit_simple(fmpcb_t y, const fmpcb_t x, ulong n, long prec)
+{
+    if (n == 0)
+    {
+        fmpcb_zero(y);
+    }
+    else if (n == 1)
+    {
+        fmpcb_inv(y, x, prec);
+    }
+    else
+    {
+        fmpcb_t p, q;
+        long wp;
+
+        wp = FMPR_PREC_ADD(prec, FLINT_BIT_COUNT(n));
+
+        fmpcb_init(p);
+        fmpcb_init(q);
+
+        bsplit(p, q, x, 0, n, wp);
+        fmpcb_div(y, p, q, prec);
+
+        fmpcb_clear(p);
+        fmpcb_clear(q);
+    }
+}
+

gamma/stirling_eval_series_fmpcb.c

@@ -31,7 +31,7 @@
 B_(2n) / (2n (2n-1)) / |z|^(2n-1) * (1/cos(0.5*arg(z))^(2n))
 */
 void
-gamma_stirling_bound_remainder_fmpcb(fmpr_t err, const fmpcb_t z, long n)
+gamma_stirling_bound_remainder_fmpcb(fmpr_t err, const fmpcb_t z, int digamma, long n)
 {
     fmpr_t t;
     fmprb_t b;
@@ -52,13 +52,17 @@ gamma_stirling_bound_remainder_fmpcb(fmpr_t err, const fmpcb_t z, long n)
         return;
     }
     fmpr_ui_div(err, 1, err, FMPRB_RAD_PREC, FMPR_RND_UP);
-    fmpr_pow_sloppy_ui(err, err, 2 * n - 1, FMPRB_RAD_PREC, FMPR_RND_UP);
+
+    if (digamma)
+        fmpr_pow_sloppy_ui(err, err, 2 * n, FMPRB_RAD_PREC, FMPR_RND_UP);
+    else
+        fmpr_pow_sloppy_ui(err, err, 2 * n - 1, FMPRB_RAD_PREC, FMPR_RND_UP);
 
     /* bound coefficient */
     fmprb_init(b);
     fmpr_init(t);
 
-    gamma_stirling_coeff(b, n, 0, FMPRB_RAD_PREC);
+    gamma_stirling_coeff(b, n, digamma, FMPRB_RAD_PREC);
     fmprb_get_abs_ubound_fmpr(t, b, FMPRB_RAD_PREC);
     fmpr_mul(err, err, t, FMPRB_RAD_PREC, FMPR_RND_UP);
 
@@ -68,7 +72,11 @@ gamma_stirling_bound_remainder_fmpcb(fmpr_t err, const fmpcb_t z, long n)
     fmprb_cos(b, b, FMPRB_RAD_PREC);
     fmprb_get_abs_lbound_fmpr(t, b, FMPRB_RAD_PREC);
     fmpr_ui_div(t, 1, t, FMPRB_RAD_PREC, FMPR_RND_UP);
-    fmpr_pow_sloppy_ui(t, t, 2 * n, FMPRB_RAD_PREC, FMPR_RND_UP);
+
+    if (digamma)
+        fmpr_pow_sloppy_ui(t, t, 2 * n + 1, FMPRB_RAD_PREC, FMPR_RND_UP);
+    else
+        fmpr_pow_sloppy_ui(t, t, 2 * n, FMPRB_RAD_PREC, FMPR_RND_UP);
 
     fmpr_mul(err, err, t, FMPRB_RAD_PREC, FMPR_RND_UP);
 
@@ -79,7 +87,7 @@ gamma_stirling_bound_remainder_fmpcb(fmpr_t err, const fmpcb_t z, long n)
 void
 gamma_stirling_eval_series_fmpcb(fmpcb_t s, const fmpcb_t z, long nterms, int digamma, long prec)
 {
-    fmpcb_t t, u, w, v;
+    fmpcb_t t, logz, zinv, zinv2;
     fmprb_t b;
     fmpr_t err;
 
@@ -87,21 +95,22 @@ gamma_stirling_eval_series_fmpcb(fmpcb_t s, const fmpcb_t z, long nterms, int di
     double z_mag, term_mag;
 
     fmpcb_init(t);
-    fmpcb_init(u);
-    fmpcb_init(w);
-    fmpcb_init(v);
+    fmpcb_init(logz);
+    fmpcb_init(zinv);
+    fmpcb_init(zinv2);
     fmprb_init(b);
 
-    fmpcb_log(w, z, prec);
+    fmpcb_log(logz, z, prec);
+    fmpcb_inv(zinv, z, prec);
 
     nterms = FLINT_MAX(nterms, 1);
 
     fmpcb_zero(s);
     if (nterms > 1)
     {
-        fmpcb_inv(t, z, prec);
-        fmpcb_mul(u, t, t, prec);
-        z_mag = fmpr_get_d(fmprb_midref(fmpcb_realref(w)), FMPR_RND_UP) * 1.44269504088896;
+        fmpcb_mul(zinv2, zinv, zinv, prec);
+
+        z_mag = fmpr_get_d(fmprb_midref(fmpcb_realref(logz)), FMPR_RND_UP) * 1.44269504088896;
 
         for (k = nterms - 1; k >= 1; k--)
         {
@@ -111,49 +120,57 @@ gamma_stirling_eval_series_fmpcb(fmpcb_t s, const fmpcb_t z, long nterms, int di
             term_prec = FLINT_MIN(term_prec, prec);
             term_prec = FLINT_MAX(term_prec, 10);
 
-            gamma_stirling_coeff(b, k, 0, term_prec);
+            gamma_stirling_coeff(b, k, digamma, term_prec);
 
             if (prec > 2000)
             {
-                fmpcb_set_round(v, u, term_prec);
-                fmpcb_mul(s, s, v, term_prec);
+                fmpcb_set_round(t, zinv2, term_prec);
+                fmpcb_mul(s, s, t, term_prec);
             }
             else
-            {
-                fmpcb_mul(s, s, u, term_prec);
-            }
+                fmpcb_mul(s, s, zinv2, term_prec);
 
             fmprb_add(fmpcb_realref(s), fmpcb_realref(s), b, term_prec);
         }
 
-        fmpcb_mul(s, s, t, prec);
+        if (digamma)
+            fmpcb_mul(s, s, zinv2, prec);
+        else
+            fmpcb_mul(s, s, zinv, prec);
     }
 
     /* remainder bound */
     fmpr_init(err);
-    gamma_stirling_bound_remainder_fmpcb(err, z, nterms);
+    gamma_stirling_bound_remainder_fmpcb(err, z, digamma, nterms);
     fmprb_add_error_fmpr(fmpcb_realref(s), err);
     fmprb_add_error_fmpr(fmpcb_imagref(s), err);
     fmpr_clear(err);
 
-    /* (z-0.5)*log(z) - z + log(2*pi)/2 */
-    fmprb_one(b);
-    fmprb_mul_2exp_si(b, b, -1);
-    fmprb_set(fmpcb_imagref(t), fmpcb_imagref(z));
-    fmprb_sub(fmpcb_realref(t), fmpcb_realref(z), b, prec);
-
-    fmpcb_mul(t, w, t, prec);
-
-    fmpcb_add(s, s, t, prec);
-    fmpcb_sub(s, s, z, prec);
-
-    fmprb_const_log_sqrt2pi(b, prec);
-    fmprb_add(fmpcb_realref(s), fmpcb_realref(s), b, prec);
+    if (digamma)
+    {
+        fmpcb_neg(s, s);
+        fmpcb_mul_2exp_si(zinv, zinv, -1);
+        fmpcb_sub(s, s, zinv, prec);
+        fmpcb_add(s, s, logz, prec);
+    }
+    else
+    {
+        /* (z-0.5)*log(z) - z + log(2*pi)/2 */
+        fmprb_one(b);
+        fmprb_mul_2exp_si(b, b, -1);
+        fmprb_set(fmpcb_imagref(t), fmpcb_imagref(z));
+        fmprb_sub(fmpcb_realref(t), fmpcb_realref(z), b, prec);
+        fmpcb_mul(t, logz, t, prec);
+        fmpcb_add(s, s, t, prec);
+        fmpcb_sub(s, s, z, prec);
+        fmprb_const_log_sqrt2pi(b, prec);
+        fmprb_add(fmpcb_realref(s), fmpcb_realref(s), b, prec);
+    }
 
     fmpcb_clear(t);
-    fmpcb_clear(u);
-    fmpcb_clear(w);
-    fmpcb_clear(v);
+    fmpcb_clear(logz);
+    fmpcb_clear(zinv);
+    fmpcb_clear(zinv2);
     fmprb_clear(b);
 }
 

gamma/stirling_eval_series_fmprb.c

@@ -138,7 +138,7 @@ gamma_stirling_eval_series_fmprb(fmprb_t s, const fmprb_t z, long nterms, int di
     else
     {
         /* (z-0.5)*log(z) - z + log(2*pi)/2 */
-        fmprb_set_ui(t, 1);
+        fmprb_one(t);
         fmprb_mul_2exp_si(t, t, -1);
         fmprb_sub(t, z, t, prec);
         fmprb_mul(t, logz, t, prec);