improve lgamma_series; remove log_gamma_series

doc/source/fmpcb.rst

@@ -15,11 +15,11 @@ Types, macros and constants
 
 .. type:: fmpcb_ptr
 
-   Alias for :code:`fmpcb_struct *`, used for vectors of numbers.
+   Alias for ``fmpcb_struct *``, used for vectors of numbers.
 
 .. type:: fmpcb_srcptr
 
-   Alias for :code:`const fmpcb_struct *`, used for vectors of numbers
+   Alias for ``const fmpcb_struct *``, used for vectors of numbers
    when passed as constant input to functions.
 
 .. macro:: fmpcb_realref(x)

doc/source/fmprb.rst

@@ -56,11 +56,11 @@ Types, macros and constants
 
 .. type:: fmprb_ptr
 
-   Alias for :code:`fmprb_struct *`, used for vectors of numbers.
+   Alias for ``fmprb_struct *``, used for vectors of numbers.
 
 .. type:: fmprb_srcptr
 
-   Alias for :code:`const fmprb_struct *`, used for vectors of numbers
+   Alias for ``const fmprb_struct *``, used for vectors of numbers
    when passed as constant input to functions.
 
 .. macro:: FMPRB_RAD_PREC

doc/source/fmprb_poly.rst

@@ -656,10 +656,20 @@ Special functions
     The underscore version does not support aliasing, and requires
     the lengths to be nonzero.
 
-.. function:: void fmprb_poly_log_gamma_series(fmprb_poly_t f, long n, long prec)
+.. function:: void _fmprb_poly_lgamma_series(fmprb_ptr res, fmprb_ptr h, long hlen, long n, long prec)
 
-    Sets `f` to the series expansion of `\log(\Gamma(1-x))`, truncated to
-    length `n`.
+.. function:: void fmprb_poly_lgamma_series(fmprb_poly_t res, const fmprb_poly_t h, long n, long prec)
+
+    Sets *res* to the series expansion of `\log \Gamma(h(x))`, truncated to
+    length *n*.
+
+    This function first computes the expansion with respect to the constant
+    term of *h*, and then calls :func:`_fmprb_poly_compose_series`.
+    It uses an expansion in terms of the Riemann zeta function if the
+    constant term of *h* is a small integer, and Stirling's series otherwise.
+
+    The underscore method supports aliasing of the input and output
+    arrays. It requires that *hlen* and *n* are greater than zero.
 
 .. function:: void _fmprb_poly_rfac_series_ui(fmprb_ptr res, fmprb_srcptr f, long flen, ulong r, long trunc, long prec)
 

fmprb_poly.h

@@ -492,7 +492,9 @@ void _fmprb_poly_tan_series(fmprb_ptr g, fmprb_srcptr h, long hlen, long len, lo
 
 void fmprb_poly_tan_series(fmprb_poly_t g, const fmprb_poly_t h, long n, long prec);
 
-void fmprb_poly_log_gamma_series(fmprb_poly_t z, long n, long prec);
+void _fmprb_poly_lgamma_series(fmprb_ptr res, fmprb_ptr h, long hlen, long len, long prec);
+
+void fmprb_poly_lgamma_series(fmprb_poly_t res, const fmprb_poly_t f, long n, long prec);
 
 void _fmprb_poly_rfac_series_ui(fmprb_ptr res, fmprb_srcptr f, long flen, ulong r, long trunc, long prec);
 

fmprb_poly/lgamma_series.c

@@ -25,39 +25,114 @@
 
 #include "fmprb_poly.h"
 #include "gamma.h"
+#include "zeta.h"
+
+int
+fmpr_cmpabs_ui(const fmpr_t x, ulong y)
+{
+    fmpr_t t;
+    int res;
+    fmpr_init(t);
+    fmpr_set_ui(t, y);
+    res = fmpr_cmpabs(x, t);
+    fmpr_clear(t);
+    return res;
+}
+
+
+static __inline__ int
+fmprb_is_int(const fmprb_t x)
+{
+    return fmprb_is_zero(x) ||
+        (fmprb_is_exact(x) &&
+                 fmpz_sgn(fmpr_expref(fmprb_midref(x))) >= 0);
+}
+
+static __inline__ void
+_fmprb_vec_indeterminate(fmprb_ptr vec, long len)
+{
+    long i;
+    for (i = 0; i < len; i++)
+    {
+        fmpr_nan(fmprb_midref(vec + i));
+        fmpr_pos_inf(fmprb_radref(vec + i));
+    }
+}
+
+static __inline__ void
+_log_rfac_series(fmprb_ptr t, const fmprb_t x, long r, long len, long prec)
+{
+    fmprb_struct f[2];
+    long rflen;
+
+    fmprb_init(f);
+    fmprb_init(f + 1);
+    fmprb_set(f, x);
+    fmprb_one(f + 1);
+
+    rflen = FLINT_MIN(len, r + 1);
+    _fmprb_poly_rfac_series_ui(t, f, FLINT_MIN(2, len), r, rflen, prec);
+    _fmprb_poly_log_series(t, t, rflen, len, prec);
+
+    fmprb_clear(f);
+    fmprb_clear(f + 1);
+}
 
 void
 _fmprb_poly_lgamma_series(fmprb_ptr res, fmprb_ptr h, long hlen, long len, long prec)
 {
     int reflect;
-    long r, n, rflen, wp;
+    long i, r, n, wp;
     fmprb_t zr;
-    fmprb_ptr t, u, f;
+    fmprb_ptr t, u;
 
     hlen = FLINT_MIN(hlen, len);
     wp = prec + FLINT_BIT_COUNT(prec);
-    gamma_stirling_choose_param_fmprb(&reflect, &r, &n, h, 0, 0, wp);
 
     t = _fmprb_vec_init(len);
     u = _fmprb_vec_init(len);
-    f = _fmprb_vec_init(2);
     fmprb_init(zr);
 
-    /* log(gamma(z)) = log(gamma(z+r)) - log(rf(z,r)) */
-    fmprb_add_ui(zr, h, r, wp);
-
-    /* u = log(gamma(z+r)) */
-    gamma_stirling_eval_fmprb_series(u, zr, n, len, wp);
-
-    /* subtract t = log(rf(z,r)) */
-    if (r != 0)
+    /* use zeta values at small integers */
+    if (fmprb_is_int(h) && (fmpr_cmpabs_ui(fmprb_midref(h), prec / 2) < 0))
     {
-        fmprb_set(f, h);
-        fmprb_one(f + 1);
-        rflen = FLINT_MIN(len, r + 1);
-        _fmprb_poly_rfac_series_ui(t, f, FLINT_MIN(2, len), r, rflen, prec);
-        _fmprb_poly_log_series(t, t, rflen, len, wp);
-        _fmprb_vec_sub(u, u, t, len, prec);
+        r = fmpr_get_si(fmprb_midref(h), FMPR_RND_DOWN);
+
+        if (r <= 0)
+        {
+            _fmprb_vec_indeterminate(res, len);
+        }
+        else
+        {
+            fmprb_zero(t);
+
+            if (len > 1) fmprb_const_euler(u + 1, wp);
+            if (len > 2) zeta_ui_vec(u + 2, 2, len - 2, wp);
+            for (i = 2; i < len; i++)
+                fmprb_div_ui(u + i, u + i, i, wp);
+            for (i = 1; i < len; i += 2)
+                fmprb_neg(u + i, u + i);
+
+            if (r != 1)
+            {
+                fmprb_one(zr);
+                _log_rfac_series(t, zr, r - 1, len, wp);
+                _fmprb_vec_add(u, u, t, len, wp);
+            }
+        }
+    }
+    else
+    {
+        /* otherwise use Stirling series */
+        gamma_stirling_choose_param_fmprb(&reflect, &r, &n, h, 0, 0, wp);
+        fmprb_add_ui(zr, h, r, wp);
+        gamma_stirling_eval_fmprb_series(u, zr, n, len, wp);
+
+        if (r != 0)
+        {
+            _log_rfac_series(t, h, r, len, wp);
+            _fmprb_vec_sub(u, u, t, len, wp);
+        }
     }
 
     /* compose with nonconstant part */
@@ -68,7 +143,6 @@ _fmprb_poly_lgamma_series(fmprb_ptr res, fmprb_ptr h, long hlen, long len, long
     fmprb_clear(zr);
     _fmprb_vec_clear(t, len);
     _fmprb_vec_clear(u, len);
-    _fmprb_vec_clear(f, 2);
 }
 
 void

fmprb_poly/log_gamma_series.c

@@ -1,47 +0,0 @@
-/*=============================================================================
-
-    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 "fmprb_poly.h"
-#include "zeta.h"
-
-/* series expansion for log(gamma(1-x)) at x = 0 */
-
-void
-fmprb_poly_log_gamma_series(fmprb_poly_t z, long n, long prec)
-{
-    long i;
-
-    fmprb_poly_fit_length(z, n);
-    _fmprb_poly_set_length(z, n);
-
-    if (n > 0) fmprb_zero(z->coeffs);
-    if (n > 1) fmprb_const_euler(z->coeffs + 1, prec);
-    if (n > 2) zeta_ui_vec(z->coeffs + 2, 2, n - 2, prec);
-
-    for (i = 2; i < n; i++)
-        fmprb_div_ui(z->coeffs + i, z->coeffs + i, i, prec);
-
-    _fmprb_poly_normalise(z);
-}

gamma/taylor_precompute.c

@@ -52,7 +52,9 @@ gamma_taylor_precompute(long num, long prec)
 
         /* TODO: cleanup */
         fmprb_poly_init(A);
-        fmprb_poly_log_gamma_series(A, num, wp);
+        fmprb_poly_set_coeff_si(A, 0, 1);
+        fmprb_poly_set_coeff_si(A, 1, -1);
+        fmprb_poly_lgamma_series(A, A, num, wp);
         fmprb_poly_neg(A, A);
         fmprb_poly_exp_series(A, A, num, wp);