enhancements to log/exp series; general zeta_ui; series expansion of log gamma

arb.h

@@ -148,7 +148,10 @@ void _arb_rad_add_ufloat(arb_t y, const ufloat_t err);
 void arb_const_pi_chudnovsky(arb_t x);
 void arb_const_euler_brent_mcmillan(arb_t x);
 void arb_const_zeta3_bsplit(arb_t x);
+
 void arb_zeta_ui_bsplit(arb_t x, ulong s);
+void arb_zeta_ui_mpfr(arb_t x, ulong n);
+void arb_zeta_ui(arb_t x, ulong n);
 
 /* fmpz extras */
 

arb/log_ui.c

@@ -63,4 +63,6 @@ arb_log_ui(arb_t x, ulong n)
     }
 
     arb_set_mpfr(x, t, 1);
+
+    mpfr_clear(t);
 }

arb/test/t-zeta_ui.c

@@ -0,0 +1,81 @@
+/*=============================================================================
+
+    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 "arb.h"
+
+int main()
+{
+    long iter;
+    flint_rand_t state;
+
+    printf("zeta_ui....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000; iter++)
+    {
+        arb_t r;
+        ulong n;
+        mpfr_t s;
+        long effective_prec;
+
+        arb_init(r, 1 + n_randint(state, 1 << n_randint(state, 14)));
+        mpfr_init2(s, arb_prec(r) + 100);
+        arb_randtest(r, state, 10);
+
+        do { n = n_randint(state, 1 << n_randint(state, 10)); } while (n == 1);
+
+        arb_zeta_ui(r, n);
+        mpfr_zeta_ui(s, n, MPFR_RNDN);
+
+        if (!arb_contains_mpfr(r, s))
+        {
+            printf("FAIL: containment\n\n");
+            printf("n = %lu\n\n", n);
+            printf("r = "); arb_debug(r); printf("\n\n");
+            printf("s = "); mpfr_printf("%.275Rf\n", s); printf("\n\n");
+            abort();
+        }
+
+        effective_prec = fmpz_bits(arb_midref(r)) - fmpz_bits(arb_radref(r));
+
+        if (!fmpz_is_zero(arb_radref(r)) && effective_prec < arb_prec(r) - 4)
+        {
+            printf("FAIL: poor accuracy\n\n");
+            printf("r = "); arb_debug(r); printf("\n\n");
+            abort();
+        }
+
+        arb_clear(r);
+        mpfr_clear(s);
+    }
+
+    flint_randclear(state);
+    _fmpz_cleanup();
+    mpfr_free_cache();
+    printf("PASS\n");
+    return EXIT_SUCCESS;
+}

arb/zeta_ui.c

@@ -0,0 +1,104 @@
+/*=============================================================================
+
+    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 "arb.h"
+#include "arith.h"
+
+static void
+mpfr_zeta_ui_bernoulli(mpfr_t x, ulong n)
+{
+    fmpq_t b;
+    mpfr_t t;
+    mpz_t f;
+
+    if (n % 2)
+        abort();
+
+    fmpq_init(b);
+    mpfr_init2(t, mpfr_get_prec(x) + FLINT_BIT_COUNT(n) + 2);
+    mpz_init(f);
+
+    bernoulli_number(b, n);
+
+    fmpq_get_mpfr(x, b, MPFR_RNDD);
+    mpfr_const_pi(t, MPFR_RNDD);
+    mpfr_mul_2exp(t, t, 1, MPFR_RNDD);
+    mpfr_pow_ui(t, t, n, MPFR_RNDD);
+    mpz_fac_ui(f, n);
+    mpfr_div_z(t, t, f, MPFR_RNDD);
+    mpfr_mul(x, x, t, MPFR_RNDD);
+    mpfr_abs(x, x, MPFR_RNDD);
+    mpfr_div_2exp(x, x, 1, MPFR_RNDD);
+
+    mpfr_clear(t);
+    fmpq_clear(b);
+    mpz_clear(f);
+}
+
+void
+arb_zeta_ui(arb_t x, ulong n)
+{
+    long prec = arb_prec(x);
+
+    /* asymptotic case */
+    if (n == 0 || n > 0.7 * prec)
+    {
+        arb_zeta_ui_mpfr(x, n);
+    }
+    else if (n == 3)
+    {
+        arb_const_zeta3_bsplit(x);
+    }
+    /* small even n */
+    else if ((n % 2 == 0) && (n < 40 + 0.11*prec))
+    {
+        mpfr_t t;
+        mpfr_init2(t, prec + 10);
+        mpfr_zeta_ui_bernoulli(t, n);
+        arb_set_mpfr(x, t, 100);  /* XXX */
+        mpfr_clear(t);
+    }
+    /* small odd n, extremely high precision */
+    /* FIXME: arb: n < prec * 0.0006 */
+    else if (n < prec * 0.0006)
+    {
+        arb_zeta_ui_bsplit(x, n);
+    }
+    /* large n */
+    else if (prec > 20 && n > 0.4 * pow(prec, 0.8))
+    {
+        mpfr_t t;
+        mpfr_init2(t, prec + 10);
+        _zeta_inv_euler_product(t, n, 0);
+        mpfr_si_div(t, 1, t, MPFR_RNDD);
+        arb_set_mpfr(x, t, 100);  /* XXX */
+        mpfr_clear(t);
+    }
+    /* fallback */
+    else
+    {
+        arb_zeta_ui_mpfr(x, n);
+    }
+}

arb/zeta_ui_mpfr.c

@@ -0,0 +1,42 @@
+/*=============================================================================
+
+    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 "arb.h"
+
+void
+arb_zeta_ui_mpfr(arb_t x, ulong n)
+{
+    mpfr_t t;
+
+    mpfr_init2(t, FLINT_MAX(arb_prec(x), FLINT_BITS));
+
+    if (n == 1)
+        mpfr_const_euler(t, MPFR_RNDN);
+    else
+        mpfr_zeta_ui(t, n, MPFR_RNDN);
+
+    arb_set_mpfr(x, t, 1);
+    mpfr_clear(t);
+}

arb_poly.h

@@ -113,7 +113,13 @@ void arb_poly_inv_series(arb_poly_t z, const arb_poly_t x, long n);
 void arb_poly_derivative(arb_poly_t z, const arb_poly_t x);
 void arb_poly_integral(arb_poly_t z, const arb_poly_t x);
 
-void arb_poly_log_series(arb_poly_t z, const arb_poly_t x, long n);
-void arb_poly_exp_series(arb_poly_t z, const arb_poly_t x, long n);
+void arb_poly_log_series(arb_poly_t z, const arb_poly_t x, long n, int exact_const);
+void arb_poly_exp_series(arb_poly_t z, const arb_poly_t x, long n, int exact_const);
+
+void _arb_poly_set_coeff_same_exp(arb_poly_t z, long i, const arb_t c);
+
+
+void arb_poly_gamma_log_series(arb_poly_t z, ulong n);
+
 
 #endif

arb_poly/exp_series.c

@@ -26,7 +26,7 @@
 #include "arb_poly.h"
 
 void
-arb_poly_exp_series(arb_poly_t z, const arb_poly_t x, long n)
+arb_poly_exp_series(arb_poly_t z, const arb_poly_t x, long n, int exact_const)
 {
     long a[FLINT_BITS];
     long i;
@@ -52,19 +52,22 @@ arb_poly_exp_series(arb_poly_t z, const arb_poly_t x, long n)
 
     arb_poly_init(c, arb_poly_prec(z));
 
-    /* first coefficient */
-    arb_init(r, arb_poly_prec(z));
-    fmpz_set(arb_midref(r), arb_poly_coeffs(x));
-    fmpz_set(arb_radref(r), arb_poly_radref(x));
-    fmpz_set(arb_expref(r), arb_poly_expref(x));
+    if (!exact_const)
+    {
+        /* first coefficient */
+        arb_init(r, arb_poly_prec(z));
+        fmpz_set(arb_midref(r), arb_poly_coeffs(x));
+        fmpz_set(arb_radref(r), arb_poly_radref(x));
+        fmpz_set(arb_expref(r), arb_poly_expref(x));
 
-    arb_exp(r, r);
-    _arb_poly_fit_length(c, 1);
-    fmpz_set(arb_poly_coeffs(c), arb_midref(r));
-    fmpz_set(arb_poly_radref(c), arb_radref(r));
-    fmpz_set(arb_poly_expref(c), arb_expref(r));
-    c->length = 1;
-    arb_clear(r);
+        arb_exp(r, r);
+        _arb_poly_fit_length(c, 1);
+        fmpz_set(arb_poly_coeffs(c), arb_midref(r));
+        fmpz_set(arb_poly_radref(c), arb_radref(r));
+        fmpz_set(arb_poly_expref(c), arb_expref(r));
+        c->length = 1;
+        arb_clear(r);
+    }
 
     arb_poly_set_si(z, 1);
 
@@ -72,14 +75,15 @@ arb_poly_exp_series(arb_poly_t z, const arb_poly_t x, long n)
     {
         n = a[i];
 
-        arb_poly_log_series(t, z, n);
+        arb_poly_log_series(t, z, n, 1);
         arb_poly_neg(t, t);
         arb_poly_add(t, t, v);
         arb_poly_mullow(t, z, t, n);
         arb_poly_add(z, z, t);
     }
 
-    arb_poly_mul(z, z, c);
+    if (!exact_const)
+        arb_poly_mul(z, z, c);
 
     arb_poly_clear(t);
     arb_poly_clear(v);

arb_poly/gamma_log_series.c

@@ -0,0 +1,68 @@
+/*=============================================================================
+
+    This file is part of FLINT.
+
+    FLINT 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.
+
+    FLINT 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 FLINT; 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 "arb_poly.h"
+
+/* series expansion for log(gamma(1-x)) at x = 0 */
+
+void
+arb_poly_gamma_log_series(arb_poly_t z, ulong n)
+{
+    arb_t c;
+    long i, prec;
+
+    prec = arb_poly_prec(z) + FLINT_BIT_COUNT(n);
+
+    arb_init(c, prec);
+
+    arb_poly_zero(z);
+    _arb_poly_fit_length(z, n);
+    z->length = n;
+
+    fmpz_set_si(arb_poly_expref(z), -prec);
+
+    for (i = 0; i < n; i++)
+    {
+        if (i == 0)
+        {
+            arb_zero(c);
+        }
+        else if (i == 1)
+        {
+            arb_const_euler_brent_mcmillan(c);
+        }
+        else
+        {
+            /* printf("%ld of %ld\n", i, n); */
+            arb_zeta_ui(c, i);
+            fmpz_tdiv_q_ui(arb_midref(c), arb_midref(c), i);
+            fmpz_add_ui(arb_radref(c), arb_radref(c), 1);
+        }
+
+        _arb_poly_set_coeff_same_exp(z, i, c);
+    }
+
+    arb_clear(c);
+}

arb_poly/log_series.c

@@ -26,43 +26,7 @@
 #include "arb_poly.h"
 
 void
-_arb_poly_set_coeff_same_exp(arb_poly_t z, const arb_t c)
-{
-    long e1, e2;
-    fmpz_t rad;
-
-    e1 = z->exp;
-    e2 = c->exp;
-
-    if (z->length < 1)
-    {
-        _arb_poly_fit_length(z, 1);
-        fmpz_set(arb_expref(z), arb_expref(c));
-        z->length = 1;
-    }
-
-    fmpz_init(rad);
-
-    if (e1 > e2)
-    {
-        fmpz_tdiv_q_2exp(z->coeffs, arb_midref(c), e1 - e2);
-        fmpz_cdiv_q_2exp(rad, arb_radref(c), e1 - e2);
-        fmpz_add_ui(rad, rad, 1);
-    }
-    else
-    {
-        fmpz_mul_2exp(z->coeffs, arb_midref(c), e2 - e1);
-        fmpz_mul_2exp(rad, arb_radref(c), e2 - e1);
-    }
-
-    if (fmpz_cmp(arb_poly_radref(z), rad) < 0)
-        fmpz_set(arb_poly_radref(z), rad);
-
-    fmpz_clear(rad);
-}
-
-void
-arb_poly_log_series(arb_poly_t z, const arb_poly_t x, long n)
+arb_poly_log_series(arb_poly_t z, const arb_poly_t x, long n, int exact_const)
 {
     arb_poly_t t, u;
     arb_t c;
@@ -72,18 +36,24 @@ arb_poly_log_series(arb_poly_t z, const arb_poly_t x, long n)
     arb_init(c, arb_poly_prec(z));
 
     /* first term */
-    fmpz_set(arb_midref(c), x->coeffs);
-    fmpz_set(arb_radref(c), arb_poly_radref(x));
-    fmpz_set(arb_expref(c), arb_poly_expref(x));
-
-    arb_log(c, c);
+    if (!exact_const)
+    {
+        fmpz_set(arb_midref(c), x->coeffs);
+        fmpz_set(arb_radref(c), arb_poly_radref(x));
+        fmpz_set(arb_expref(c), arb_poly_expref(x));
+        arb_log(c, c);
+    }
 
     arb_poly_inv_series(t, x, n);
+
     arb_poly_derivative(u, x);
     arb_poly_mullow(t, t, u, n - 1);
     arb_poly_integral(z, t);
 
-    _arb_poly_set_coeff_same_exp(z, c);
+    if (!exact_const)
+    {
+        _arb_poly_set_coeff_same_exp(z, 0, c);
+    }
 
     arb_clear(c);
     arb_poly_clear(t);

arb_poly/set_coeff_same_exp.c

@@ -0,0 +1,63 @@
+/*=============================================================================
+
+    This file is part of FLINT.
+
+    FLINT 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.
+
+    FLINT 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 FLINT; 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 "arb_poly.h"
+
+void
+_arb_poly_set_coeff_same_exp(arb_poly_t z, long i, const arb_t c)
+{
+    long e1, e2, j;
+    fmpz_t rad;
+
+    e1 = z->exp;
+    e2 = c->exp;
+
+    if (z->length < i + 1)
+    {
+        _arb_poly_fit_length(z, i + 1);
+        for (j = z->length; j < i; j++)
+            fmpz_zero(z->coeffs + j);
+        z->length = i + 1;
+    }
+
+    fmpz_init(rad);
+
+    if (e1 > e2)
+    {
+        fmpz_tdiv_q_2exp(z->coeffs + i, arb_midref(c), e1 - e2);
+        fmpz_cdiv_q_2exp(rad, arb_radref(c), e1 - e2);
+        fmpz_add_ui(rad, rad, 1);
+    }
+    else
+    {
+        fmpz_mul_2exp(z->coeffs + i, arb_midref(c), e2 - e1);
+        fmpz_mul_2exp(rad, arb_radref(c), e2 - e1);
+    }
+
+    if (fmpz_cmp(arb_poly_radref(z), rad) < 0)
+        fmpz_set(arb_poly_radref(z), rad);
+
+    fmpz_clear(rad);
+}