lgamma series expansion (work in progress; may still have bugs)

fmpr.h

@@ -79,6 +79,8 @@ typedef struct
 fmpr_struct;
 
 typedef fmpr_struct fmpr_t[1];
+typedef fmpr_struct * fmpr_ptr;
+typedef const fmpr_struct * fmpr_srcptr;
 
 #define fmpr_manref(x) (&(x)->man)
 #define fmpr_expref(x) (&(x)->exp)
@@ -692,5 +694,28 @@ void fmpr_pow_sloppy_ui(fmpr_t y, const fmpr_t b, ulong e,
 void fmpr_pow_sloppy_si(fmpr_t y, const fmpr_t b, long e,
     long prec, fmpr_rnd_t rnd);
 
+/* vector functions */
+
+static __inline__ fmpr_ptr
+_fmpr_vec_init(long n)
+{
+    long i;
+    fmpr_ptr v = flint_malloc(sizeof(fmpr_struct) * n);
+
+    for (i = 0; i < n; i++)
+        fmpr_init(v + i);
+
+    return v;
+}
+
+static __inline__ void
+_fmpr_vec_clear(fmpr_ptr v, long n)
+{
+    long i;
+    for (i = 0; i < n; i++)
+        fmpr_clear(v + i);
+    flint_free(v);
+}
+
 #endif
 

fmprb.h

@@ -739,4 +739,12 @@ _fmprb_vec_set_powers(fmprb_ptr xs, const fmprb_t x, long len, long prec)
     }
 }
 
+static __inline__ void
+_fmprb_vec_add_error_fmpr_vec(fmprb_ptr res, fmpr_srcptr err, long len)
+{
+    long i;
+    for (i = 0; i < len; i++)
+        fmprb_add_error_fmpr(res + i, err + i);
+}
+
 #endif

fmprb_poly/lgamma_series.c

@@ -0,0 +1,113 @@
+/*=============================================================================
+
+    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 "fmprb_poly.h"
+#include "gamma.h"
+
+void
+_fmprb_poly_lgamma_cpx_series(fmprb_ptr res, const fmprb_t z, long len, long prec)
+{
+    int reflect;
+    long r, n, rflen, wp;
+    fmprb_t zr;
+    fmprb_ptr t, u, f;
+
+    wp = prec + FLINT_BIT_COUNT(prec);
+
+    gamma_stirling_choose_param_fmprb(&reflect, &r, &n, z, 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, z, r, wp);
+
+    /* u = log(gamma(z+r)) */
+    gamma_stirling_eval_fmprb_series(u, zr, n, len, wp);
+
+    /* t = rising factorial */
+    fmprb_set(f, z);
+    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(res, u, t, len, prec);
+
+    fmprb_clear(zr);
+    _fmprb_vec_clear(t, len);
+    _fmprb_vec_clear(u, len);
+    _fmprb_vec_clear(f, 2);
+}
+
+void
+fmprb_poly_lgamma_cpx_series(fmprb_poly_t res, const fmprb_t z, long num, long prec)
+{
+    fmprb_poly_fit_length(res, num);
+    _fmprb_poly_lgamma_cpx_series(res->coeffs, z, num, prec);
+    _fmprb_poly_set_length(res, num);
+    _fmprb_poly_normalise(res);
+}
+
+void
+_fmprb_poly_lgamma_series(fmprb_ptr res, fmprb_ptr h, long hlen, long len, long prec)
+{
+    fmprb_ptr t, u;
+    fmprb_t z;
+
+    hlen = FLINT_MIN(hlen, len);
+
+    t = _fmprb_vec_init(hlen);
+    u = _fmprb_vec_init(len);
+    fmprb_init(z);
+
+    fmprb_set(z, h);
+    _fmprb_vec_set(t + 1, h + 1, hlen - 1);
+
+    _fmprb_poly_lgamma_cpx_series(u, z, len, prec);
+    _fmprb_poly_compose_series(res, u, len, t, hlen, len, prec);
+
+    fmprb_clear(z);
+    _fmprb_vec_clear(t, hlen);
+    _fmprb_vec_clear(u, len);
+}
+
+void
+fmprb_poly_lgamma_series(fmprb_poly_t res, const fmprb_poly_t f, long n, long prec)
+{
+    if (f->length == 0 || n == 0)
+    {
+        printf("fmprb_poly_lgamma_series: require n > 0 and nonzero input\n");
+        abort();
+    }
+
+    fmprb_poly_fit_length(res, n);
+    _fmprb_poly_lgamma_series(res->coeffs, f->coeffs, f->length, n, prec);
+    _fmprb_poly_set_length(res, n);
+    _fmprb_poly_normalise(res);
+}
+

gamma.h

@@ -79,6 +79,8 @@ void gamma_stirling_coeff(fmprb_t b, ulong k, int digamma, long prec);
 void gamma_stirling_eval_series_fmprb(fmprb_t s, const fmprb_t z, long nterms, int digamma, long prec);
 void gamma_stirling_eval_series_fmpcb(fmpcb_t s, const fmpcb_t z, long nterms, int digamma, long prec);
 
+void gamma_stirling_eval_fmprb_series(fmprb_ptr res, const fmprb_t z, long n, long num, long prec);
+
 
 void gamma_rising_fmprb_ui_bsplit_simple(fmprb_t y, const fmprb_t x, ulong n, long prec);
 void gamma_rising_fmprb_ui_bsplit_eight(fmprb_t y, const fmprb_t x, ulong n, long prec);

gamma/stirling_eval_fmprb_series.c

@@ -0,0 +1,175 @@
+/*=============================================================================
+
+    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(fmprb_ptr Q, fmprb_ptr T, const fmprb_t z, long a, long b, long num, long prec)
+{
+    if (b - a == 1)
+    {
+        gamma_stirling_coeff(T, a, 0, prec);
+
+        if (a == 1)
+        {   /* (z + t) */
+            fmprb_set(Q, z);
+            if (num > 1) fmprb_one(Q + 1);
+            if (num > 2) fmprb_zero(Q + 2);
+        }             
+        else
+        {   /* (z + t)^2 */
+            fmprb_mul(Q, z, z, prec);  /* TODO: precompute */
+            if (num > 1) fmprb_mul_2exp_si(Q + 1, z, 1);
+            if (num > 2) fmprb_one(Q + 2);
+        }
+    }
+    else
+    {
+        long m, n1, n2, q1len, q2len, t1len, t2len, qlen, tlen, alloc;
+        fmprb_ptr Q1, T1, Q2, T2;
+
+        m = a + (b - a) / 2;
+
+        n1 = m - a;
+        n2 = b - m;
+        q1len = FLINT_MIN(2 * n1 + 1, num);
+        t1len = FLINT_MIN(2 * n1 - 1, num);
+        q2len = FLINT_MIN(2 * n2 + 1, num);
+        t2len = FLINT_MIN(2 * n2 - 1, num);
+        qlen = FLINT_MIN(q1len + q2len - 1, num);
+        tlen = FLINT_MIN(t1len + q2len - 1, num);
+
+        alloc = q1len + q2len + t1len + t2len;
+        Q1 = _fmprb_vec_init(alloc);
+        Q2 = Q1 + q1len;
+        T1 = Q2 + q2len;
+        T2 = T1 + t1len;
+
+        bsplit(Q1, T1, z, a, m, num, prec);
+        bsplit(Q2, T2, z, m, b, num, prec);
+
+        _fmprb_poly_mullow(Q, Q2, q2len, Q1, q1len, qlen, prec);
+        _fmprb_poly_mullow(T, Q2, q2len, T1, t1len, tlen, prec);
+        _fmprb_poly_add(T, T, tlen, T2, t2len, prec);
+
+        _fmprb_vec_clear(Q1, alloc);
+    }
+}
+
+void
+_fmprb_poly_mullow_cpx(fmprb_ptr res, fmprb_srcptr src, long len, const fmprb_t c, long trunc, long prec)
+{
+    long i;
+
+    if (len < trunc)
+        fmprb_set(res + len, src + len - 1);
+
+    for (i = len - 1; i > 0; i--)
+    {
+        fmprb_mul(res + i, src + i, c, prec);
+        fmprb_add(res + i, res + i, src + i - 1, prec);
+    }
+
+    fmprb_mul(res, src, c, prec);
+}
+
+void
+_fmprb_poly_log_cpx_series(fmprb_ptr res, const fmprb_t c, long num, long prec)
+{
+    long i;
+
+    for (i = 0; i < num; i++)
+    {
+        if (i == 0)
+            fmprb_log(res + i, c, prec);
+        else if (i == 1)
+            fmprb_ui_div(res + i, 1, c, prec);
+        else
+            fmprb_mul(res + i, res + i - 1, res + 1, prec);
+    }
+
+    for (i = 2; i < num; i++)
+    {
+        fmprb_div_ui(res + i, res + i, i, prec);
+
+        if (i % 2 == 0)
+            fmprb_neg(res + i, res + i);
+    }
+}
+
+void
+gamma_stirling_eval_fmprb_series(fmprb_ptr res, const fmprb_t z, long n, long num, long prec)
+{
+    long tlen, qlen;
+    fmprb_ptr T, Q;
+    fmpr_ptr err;
+    fmprb_t c;
+
+    T = _fmprb_vec_init(num);
+    Q = _fmprb_vec_init(num);
+    err = _fmpr_vec_init(num);
+    fmprb_init(c);
+
+    gamma_stirling_bound_fmprb(err, z, 0, num, n);
+
+    if (n <= 1)
+    {
+        _fmprb_vec_zero(res, num);
+    }
+    else
+    {
+        tlen = FLINT_MIN(2 * (n - 1) + 1, num);
+        qlen = FLINT_MIN(2 * (n - 1) - 1, num);
+        bsplit(Q, T, z, 1, n, num, prec);
+        _fmprb_poly_div_series(res, T, tlen, Q, qlen, num, prec);
+    }
+
+    /* ((z-1/2) + t) * log(z+t) */
+    _fmprb_poly_log_cpx_series(T, z, num, prec);
+    fmprb_one(c);
+    fmprb_mul_2exp_si(c, c, -1);
+    fmprb_sub(c, z, c, prec);
+    _fmprb_poly_mullow_cpx(T, T, num, c, num, prec);
+
+    /* constant term */
+    fmprb_const_log_sqrt2pi(c, prec);
+    fmprb_add(T, T, c, prec);
+
+    /* subtract (z+t) */
+    fmprb_sub(T, T, z, prec);
+    if (num > 1)
+        fmprb_sub_ui(T + 1, T + 1, 1, prec);
+
+    _fmprb_vec_add(res, res, T, num, prec);
+
+    _fmprb_vec_add_error_fmpr_vec(res, err, num);
+
+    _fmprb_vec_clear(T, num);
+    _fmprb_vec_clear(Q, num);
+    _fmpr_vec_clear(err, num);
+    fmprb_clear(c);
+}
+

gamma/test/t-stirling_eval_fmprb_series.c

@@ -0,0 +1,85 @@
+/*=============================================================================
+
+    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"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("stirling_eval_fmprb_series....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    /* check self-consistency */
+    for (iter = 0; iter < 1000; iter++)
+    {
+        fmprb_ptr u, v;
+        fmprb_t x;
+        long i, n, len1, len2, prec1, prec2;
+
+        len1 = 1 + n_randint(state, 40);
+        len2 = 1 + n_randint(state, 40);
+
+        prec1 = 2 + n_randint(state, 1000);
+        prec2 = 2 + n_randint(state, 1000);
+
+        n = 1 + n_randint(state, 40);
+
+        u = _fmprb_vec_init(len1);
+        v = _fmprb_vec_init(len2);
+        fmprb_init(x);
+
+        fmprb_randtest(x, state, 2 + n_randint(state, 1000), 20);
+
+        gamma_stirling_eval_fmprb_series(u, x, n, len1, prec1);
+        gamma_stirling_eval_fmprb_series(v, x, n, len2, prec2);
+
+        for (i = 0; i < FLINT_MIN(len1, len2); i++)
+        {
+            if (!fmprb_overlaps(u + i, v + i))
+            {
+                printf("FAIL: overlap\n\n");
+                printf("n = %ld, len1 = %ld, len2 = %ld, i = %ld\n\n", n, len1, len2, i);
+                printf("x = "); fmprb_printd(x, prec1 / 3.33); printf("\n\n");
+                printf("u = "); fmprb_printd(u + i, prec1 / 3.33); printf("\n\n");
+                printf("v = "); fmprb_printd(v + i, prec2 / 3.33); printf("\n\n");
+                abort();
+            }
+        }
+
+        _fmprb_vec_clear(u, len1);
+        _fmprb_vec_clear(v, len2);
+        fmprb_clear(x);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+