add Taylor precomp method for fast parameter multi-evaluation of Hurwitz zeta function

acb_dirichlet.h

@@ -27,6 +27,24 @@
 extern "C" {
 #endif
 
+typedef struct
+{
+    acb_struct s;
+    mag_struct err;
+    acb_ptr coeffs;
+    slong A;
+    slong N;
+    slong K;
+}
+acb_dirichlet_hurwitz_precomp_struct;
+
+typedef acb_dirichlet_hurwitz_precomp_struct acb_dirichlet_hurwitz_precomp_t[1];
+
+void acb_dirichlet_hurwitz_precomp_init(acb_dirichlet_hurwitz_precomp_t pre, const acb_t s, slong A, slong K, slong N, slong prec);
+void acb_dirichlet_hurwitz_precomp_clear(acb_dirichlet_hurwitz_precomp_t pre);
+void acb_dirichlet_hurwitz_precomp_bound(mag_t res, const acb_t s, slong A, slong K, slong N);
+void acb_dirichlet_hurwitz_precomp_eval(acb_t res, const acb_dirichlet_hurwitz_precomp_t pre, ulong p, ulong q, slong prec);
+
 void _acb_dirichlet_euler_product_real_ui(arb_t res, ulong s,
     const signed char * chi, int mod, int reciprocal, slong prec);
 

acb_dirichlet/hurwitz_precomp_bound.c

@@ -0,0 +1,88 @@
+/*
+    Copyright (C) 2016 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_dirichlet.h"
+
+void
+acb_dirichlet_hurwitz_precomp_bound(mag_t res, const acb_t s,
+    slong A, slong K, slong N)
+{
+    acb_t s1;
+    mag_t x, t, TK, C;
+    slong sigmaK;
+    arf_t u;
+
+    if (A < 1 || K < 1 || N < 1)
+    {
+        mag_inf(res);
+        return;
+    }
+
+    /* sigmaK = re(s) + K, floor bound */
+    arf_init(u);
+    arf_set_mag(u, arb_radref(acb_realref(s)));
+    arf_sub(u, arb_midref(acb_realref(s)), u, MAG_BITS, ARF_RND_FLOOR);
+    arf_add_ui(u, u, K, MAG_BITS, ARF_RND_FLOOR);
+    if (arf_cmp_ui(u, 2) < 0 || arf_cmp_2exp_si(u, FLINT_BITS - 2) > 0)
+    {
+        mag_inf(res);
+        arf_clear(u);
+        return;
+    }
+    sigmaK = arf_get_si(u, ARF_RND_FLOOR);
+    arf_clear(u);
+
+    acb_init(s1);
+    mag_init(x);
+    mag_init(t);
+    mag_init(TK);
+    mag_init(C);
+
+    /* With N grid points, we will have |x| <= 1 / (2N). */
+    mag_one(x);
+    mag_div_ui(x, x, 2 * N);
+
+    /* T(K) = |x|^K |(s)_K| / K! * [1/A^(sigma+K) + ...] */
+    mag_pow_ui(TK, x, K);
+    acb_rising_ui_get_mag(t, s, K);
+    mag_mul(TK, TK, t);
+    mag_rfac_ui(t, K);
+    mag_mul(TK, TK, t);
+    /* Note: here we assume that mag_hurwitz_zeta_uiui uses an error bound
+       that is at least as large as the one used in the proof. */
+    mag_hurwitz_zeta_uiui(t, sigmaK, A);
+    mag_mul(TK, TK, t);
+
+    /* C = |x|/A (1 + 1/(K+sigma+A-1)) (1 + |s-1|/(K+1)) */
+
+    mag_div_ui(C, x, A);
+
+    mag_one(t);
+    mag_div_ui(t, t, sigmaK + A - 1);
+    mag_add_ui(t, t, 1);
+    mag_mul(C, C, t);
+
+    acb_sub_ui(s1, s, 1, MAG_BITS);
+    acb_get_mag(t, s1);
+    mag_div_ui(t, t, K + 1);
+    mag_add_ui(t, t, 1);
+    mag_mul(C, C, t);
+
+    mag_geom_series(t, C, 0);
+    mag_mul(res, TK, t);
+
+    acb_clear(s1);
+    mag_clear(x);
+    mag_clear(t);
+    mag_clear(TK);
+    mag_clear(C);
+}
+

acb_dirichlet/hurwitz_precomp_clear.c

@@ -0,0 +1,21 @@
+/*
+    Copyright (C) 2016 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_dirichlet.h"
+
+void
+acb_dirichlet_hurwitz_precomp_clear(acb_dirichlet_hurwitz_precomp_t pre)
+{
+    acb_clear(&pre->s);
+    mag_clear(&pre->err);
+    _acb_vec_clear(pre->coeffs, pre->N * pre->K);
+}
+

acb_dirichlet/hurwitz_precomp_eval.c

@@ -0,0 +1,59 @@
+/*
+    Copyright (C) 2016 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_dirichlet.h"
+#include "acb_poly.h"
+
+void
+acb_dirichlet_hurwitz_precomp_eval(acb_t res,
+    const acb_dirichlet_hurwitz_precomp_t pre, ulong p, ulong q, slong prec)
+{
+    slong i;
+    acb_t a, t;
+
+    if (p > q)
+        abort();
+
+    acb_init(a);
+    acb_init(t);
+
+    /* todo: avoid integer overflows below */
+    if (p == q)
+        i = pre->N - 1;
+    else
+        i = (pre->N * p) / q;
+    acb_set_si(a, 2 * pre->N * p - q * (2 * i + 1));
+    acb_div_ui(a, a, 2 * q * pre->N, prec);
+
+    /* compute zeta(s,A+p/q) */
+    _acb_poly_evaluate(res, pre->coeffs + i * pre->K, pre->K, a, prec);
+
+    /* error bound */
+    if (acb_is_real(&pre->s))
+        arb_add_error_mag(acb_realref(res), &pre->err);
+    else
+        acb_add_error_mag(res, &pre->err);
+
+    /* zeta(s,p/q) = (p/q)^-s + ... + (p/q+A-1)^-s zeta(s,A+p/q) */
+    for (i = 0; i < pre->A; i++)
+    {
+        acb_set_ui(a, p);
+        acb_div_ui(a, a, q, prec);
+        acb_add_ui(a, a, i, prec);
+        acb_neg(t, &pre->s);
+        acb_pow(a, a, t, prec);
+        acb_add(res, res, a, prec);
+    }
+
+    acb_clear(a);
+    acb_clear(t);
+}
+

acb_dirichlet/hurwitz_precomp_init.c

@@ -0,0 +1,74 @@
+/*
+    Copyright (C) 2016 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_dirichlet.h"
+
+void
+acb_dirichlet_hurwitz_precomp_init(acb_dirichlet_hurwitz_precomp_t pre,
+        const acb_t s, slong A, slong K, slong N, slong prec)
+{
+    slong i, k;
+
+    if (A < 1 || K < 1 || N < 1)
+        abort();
+
+    pre->A = A;
+    pre->K = K;
+    pre->N = N;
+    pre->coeffs = _acb_vec_init(N * K);
+
+    mag_init(&pre->err);
+    acb_init(&pre->s);
+    acb_set(&pre->s, s);
+
+    acb_dirichlet_hurwitz_precomp_bound(&pre->err, s, A, K, N);
+
+    if (mag_is_finite(&pre->err))
+    {
+        acb_t t, a;
+
+        acb_init(t);
+        acb_init(a);
+
+        /* (-1)^k (s)_k / k! */
+        acb_one(pre->coeffs + 0);
+        for (k = 1; k < K; k++)
+        {
+            acb_add_ui(pre->coeffs + k, s, k - 1, prec);
+            acb_mul(pre->coeffs + k, pre->coeffs + k, pre->coeffs + k - 1, prec);
+            acb_div_ui(pre->coeffs + k, pre->coeffs + k, k, prec);
+            acb_neg(pre->coeffs + k, pre->coeffs + k);
+        }
+
+        for (i = 1; i < N; i++)
+            _acb_vec_set(pre->coeffs + i * K, pre->coeffs, K);
+
+        /* zeta(s+k,a) where a = A + (2*i+1)/(2*N) */
+        for (i = 0; i < N; i++)
+        {
+            acb_set_ui(a, 2 * i + 1);
+            acb_div_ui(a, a, 2 * N, prec);
+            acb_add_ui(a, a, A, prec);
+
+            for (k = 0; k < K; k++)
+            {
+                acb_add_ui(t, s, k, prec);
+                acb_hurwitz_zeta(t, t, a, prec);
+                acb_mul(pre->coeffs + i * K + k,
+                        pre->coeffs + i * K + k, t, prec);
+            }
+        }
+
+        acb_clear(t);
+        acb_clear(a);
+    }
+}
+

acb_dirichlet/test/t-hurwitz_precomp.c

@@ -0,0 +1,81 @@
+/*
+    Copyright (C) 2016 Fredrik Johansson
+
+    This file is part of Arb.
+
+    Arb is free software: you can redistribute it and/or modify it under
+    the terms of the GNU Lesser General Public License (LGPL) as published
+    by the Free Software Foundation; either version 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "acb_dirichlet.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("hurwitz_precomp....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 500 * arb_test_multiplier(); iter++)
+    {
+        acb_t s, a, z1, z2;
+        ulong p, q;
+        slong prec1, prec2, A, K, N, i;
+        acb_dirichlet_hurwitz_precomp_t pre;
+
+        prec1 = 2 + n_randint(state, 100);
+        prec2 = 2 + n_randint(state, 100);
+        A = 1 + n_randint(state, 10);
+        K = 1 + n_randint(state, 10);
+        N = 1 + n_randint(state, 10);
+
+        acb_init(s);
+        acb_init(a);
+        acb_init(z1);
+        acb_init(z2);
+
+        acb_randtest(s, state, 1 + n_randint(state, 200), 2);
+        acb_dirichlet_hurwitz_precomp_init(pre, s, A, K, N, prec1);
+
+        for (i = 0; i < 10; i++)
+        {
+            q = 1 + n_randint(state, 1000);
+            p = 1 + n_randint(state, q);
+
+            acb_dirichlet_hurwitz_precomp_eval(z1, pre, p, q, prec1);
+
+            acb_set_ui(a, p);
+            acb_div_ui(a, a, q, prec2);
+            acb_hurwitz_zeta(z2, s, a, prec2);
+
+            if (!acb_overlaps(z1, z2))
+            {
+                flint_printf("FAIL! (overlap)");
+                flint_printf("s = "); acb_printn(s, 50, 0); flint_printf("\n\n");
+                flint_printf("A = %wd  K = %wd  N = %wd\n\n", A, K, N);
+                flint_printf("p = %wu  q = %wu\n\n", p, q);
+                flint_printf("z1 = "); acb_printn(z1, 50, 0); flint_printf("\n\n");
+                flint_printf("z2 = "); acb_printn(z2, 50, 0); flint_printf("\n\n");
+                abort();
+            }
+        }
+
+        acb_dirichlet_hurwitz_precomp_clear(pre);
+
+        acb_clear(s);
+        acb_clear(a);
+        acb_clear(z1);
+        acb_clear(z2);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

doc/source/acb_dirichlet.rst

@@ -22,6 +22,44 @@ The code in other modules for computing the Riemann zeta function,
 Hurwitz zeta function and polylogarithm will possibly be migrated to this
 module in the future.
 
+Hurwitz zeta function
+-------------------------------------------------------------------------------
+
+.. type:: acb_dirichlet_hurwitz_precomp_struct
+
+.. type:: acb_dirichlet_hurwitz_precomp_t
+
+.. function:: void acb_dirichlet_hurwitz_precomp_init(acb_dirichlet_hurwitz_precomp_t pre, const acb_t s, ulong A, ulong K, ulong N, slong prec)
+
+    Precomputes a grid of Taylor polynomials for fast evaluation of
+    `\zeta(s,a)` on `a \in (0,1]` with fixed *s*.
+    *A* is the initial shift to apply to *a*, *K* is the number of Taylor terms,
+    *N* is the number of grid points.  The precomputation requires *NK*
+    evaluations of the Hurwitz zeta function, and each subsequent evaluation
+    requires *2K* simple arithmetic operations (polynomial evaluation) plus
+    *A* powers. As *K* grows, the error is at most `O(1/(2AN)^K)`.
+
+    We require that *A*, *K* and *N* are all positive. Moreover, for a finite
+    error bound, we require `K+\operatorname{re}(s) > 1`.
+    To avoid an initial "bump" that steals precision
+    and slows convergence, *AN* should be at least roughly as large as `|s|`,
+    e.g. it is a good idea to have at least `AN > 0.5 |s|`.
+
+.. function:: void acb_dirichlet_hurwitz_precomp_clear(acb_dirichlet_hurwitz_precomp_t pre)
+
+    Clears the precomputed data.
+
+.. function:: void acb_dirichlet_hurwitz_precomp_bound(mag_t res, const acb_t s, ulong A, ulong K, ulong N)
+
+    Computes an upper bound for the truncation error (not accounting for
+    roundoff error) when evaluating `\zeta(s,a)` with precomputation parameters
+    *A*, *K*, *N*, assuming that `0 < a \le 1`.
+    For details, see :ref:`algorithms_hurwitz`.
+
+.. function:: void acb_dirichlet_hurwitz_precomp_eval(acb_t res, const acb_dirichlet_hurwitz_precomp_t pre, ulong p, ulong q, slong prec)
+
+    Evaluates `\zeta(s,p/q)` using precomputed data, assuming that `0 < p/q \le 1`.
+
 Character evaluation
 -------------------------------------------------------------------------------
 

doc/source/gamma.rst

@@ -1,6 +1,6 @@
 .. _algorithms_gamma:
 
-Algorithms for gamma functions
+Algorithms for the gamma function
 ===============================================================================
 
 The Stirling series

doc/source/hurwitz.rst

@@ -0,0 +1,52 @@
+.. _algorithms_hurwitz:
+
+Algorithms for the Hurwitz zeta function
+===============================================================================
+
+Euler-Maclaurin summation
+-------------------------------------------------------------------------------
+
+The Euler-Maclaurin formula allows evaluating the Hurwitz zeta function and
+its derivatives for general complex input. The algorithm is described
+in [Joh2013]_.
+
+Parameter Taylor series
+-------------------------------------------------------------------------------
+
+To evaluate `\zeta(s,a)` for several nearby parameter values, the following
+Taylor expansion is useful:
+
+.. math ::
+
+    \zeta(s,a+x) = \sum_{k=0}^{\infty} (-x)^k \frac{(s)_k}{k!} \zeta(s+k,a)
+
+We assume that `a \ge 1` is real and that `\sigma = \operatorname{re}(s)`
+with `K + \sigma > 1`. The tail is bounded by
+
+.. math ::
+
+    \sum_{k=K}^{\infty} |x|^k \frac{|(s)_k|}{k!} \zeta(\sigma+k,a) \le
+    \sum_{k=K}^{\infty}
+        |x|^k \frac{|(s)_k|}{k!} \left[
+            \frac{1}{a^{\sigma+k}} + \frac{1}{(\sigma+k-1) a^{\sigma+k-1}} \right].
+
+Denote the term on the right by `T(k)`. Then
+
+.. math ::
+
+    \left|\frac{T(k+1)}{T(k)}\right| =
+            \frac{|x|}{a}
+            \frac{(k+\sigma-1)}{(k+\sigma)}
+            \frac{(k+\sigma+a)}{(k+\sigma+a-1)}
+            \frac{|k+s|}{(k+1)}
+        \le
+            \frac{|x|}{a}
+            \left(1 + \frac{1}{K+\sigma+a-1}\right)
+            \left(1 + \frac{|s-1|}{K+1}\right) = C
+
+and if `C < 1`,
+
+.. math ::
+
+    \sum_{k=K}^{\infty} T(k) \le \frac{T(K)}{1-C}.
+

doc/source/index.rst

@@ -191,6 +191,7 @@ lengthy to reproduce in the documentation for each module.
    formulas.rst
    constants.rst
    gamma.rst
+   hurwitz.rst
    polylogarithms.rst
    hypergeometric.rst
    agm.rst