implement the Riemann-Siegel formula

acb_dirichlet.h

@@ -33,6 +33,12 @@ void acb_dirichlet_powsum_term(acb_ptr res, arb_t log_prev, ulong * prev,
 void acb_dirichlet_powsum_sieved(acb_ptr z, const acb_t s, ulong n, slong len, slong prec);
 void acb_dirichlet_powsum_smooth(acb_ptr z, const acb_t s, ulong n, slong len, slong prec);
 
+void acb_dirichlet_zeta_rs_f_coeffs(acb_ptr c, const arb_t p, slong N, slong prec);
+void acb_dirichlet_zeta_rs_d_coeffs(arb_ptr d, const arb_t sigma, slong k, slong prec);
+void acb_dirichlet_zeta_rs_bound(mag_t err, const acb_t s, slong K);
+void acb_dirichlet_zeta_rs_r(acb_t res, const acb_t s, slong K, slong prec);
+void acb_dirichlet_zeta_rs(acb_t res, const acb_t s, slong K, slong prec);
+
 typedef struct
 {
     acb_struct s;

acb_dirichlet/test/t-zeta_rs.c

@@ -0,0 +1,84 @@
+/*
+    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("zeta_rs....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        acb_t z1, z2, s1, s2;
+        slong prec1, prec2, K1;
+
+        acb_init(z1);
+        acb_init(z2);
+        acb_init(s1);
+        acb_init(s2);
+
+        if (n_randint(state, 2))
+            arb_set_d(acb_realref(s1), 0.5);
+        else
+            arb_randtest(acb_realref(s1), state, 2 + n_randint(state, 100), 2);
+
+        arb_randtest(acb_imagref(s1), state, 2 + n_randint(state, 100), 2);
+        arb_add_ui(acb_imagref(s1), acb_imagref(s1), n_randtest(state) % 1000, 100);
+
+        acb_randtest(z1, state, 2 + n_randint(state, 100), 3);
+        acb_randtest(z2, state, 2 + n_randint(state, 100), 3);
+
+        prec1 = 2 + n_randint(state, 150);
+        prec2 = 2 + n_randint(state, 150);
+
+        K1 = n_randint(state, 10);
+
+        if (n_randint(state, 2))
+        {
+            mag_zero(arb_radref(acb_realref(s1)));
+            mag_zero(arb_radref(acb_imagref(s1)));
+        }
+
+        acb_set(s2, s1);
+
+        acb_dirichlet_zeta_rs(z1, s1, K1, prec1);
+        acb_zeta(z2, s2, prec2);  /* should be Euler-Maclaurin */
+
+        if (!acb_overlaps(z1, z2))
+        {
+            flint_printf("FAIL: overlap\n\n");
+            flint_printf("iter = %wd\n", iter);
+            flint_printf("K1 = %wd\n\n", K1);
+            flint_printf("s1 = "); acb_printn(s1, 50, 0); flint_printf("\n\n");
+            flint_printf("s2 = "); acb_printn(s2, 50, 0); flint_printf("\n\n");
+            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_clear(s1);
+        acb_clear(s2);
+        acb_clear(z1);
+        acb_clear(z2);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

acb_dirichlet/test/t-zeta_rs_r.c

@@ -0,0 +1,94 @@
+/*
+    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("zeta_rs_r....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000 * arb_test_multiplier(); iter++)
+    {
+        acb_t z1, z2, s1, s2;
+        slong prec1, prec2, K1, K2;
+
+        acb_init(z1);
+        acb_init(z2);
+        acb_init(s1);
+        acb_init(s2);
+
+        if (n_randint(state, 2))
+            arb_set_d(acb_realref(s1), 0.5);
+        else
+            arb_randtest(acb_realref(s1), state, 2 + n_randint(state, 100), 2);
+
+        arb_randtest(acb_imagref(s1), state, 2 + n_randint(state, 100), 2);
+        arb_abs(acb_imagref(s1), acb_imagref(s1));
+        arb_add_ui(acb_imagref(s1), acb_imagref(s1), n_randtest(state) % 1000000, 100);
+
+        acb_randtest(z1, state, 2 + n_randint(state, 100), 3);
+        acb_randtest(z2, state, 2 + n_randint(state, 100), 3);
+
+        prec1 = 2 + n_randint(state, 150);
+        prec2 = 2 + n_randint(state, 150);
+
+        K1 = n_randint(state, 10);
+        K2 = n_randint(state, 10);
+
+        if (n_randint(state, 2))
+        {
+            mag_zero(arb_radref(acb_realref(s1)));
+            mag_zero(arb_radref(acb_imagref(s1)));
+        }
+
+        if (n_randint(state, 2))
+        {
+            acb_set(s2, s1);
+        }
+        else
+        {
+            acb_add(s2, s1, z1, prec2);
+            acb_sub(s2, s2, z1, prec2);
+        }
+
+        acb_dirichlet_zeta_rs_r(z1, s1, K1, prec1);
+        acb_dirichlet_zeta_rs_r(z2, s2, K2, prec2);
+
+        if (!acb_overlaps(z1, z2))
+        {
+            flint_printf("FAIL: overlap\n\n");
+            flint_printf("iter = %wd\n", iter);
+            flint_printf("K1 = %wd, K2 = %wd\n\n", K1, K2);
+            flint_printf("s1 = "); acb_printn(s1, 50, 0); flint_printf("\n\n");
+            flint_printf("s2 = "); acb_printn(s2, 50, 0); flint_printf("\n\n");
+            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_clear(s1);
+        acb_clear(s2);
+        acb_clear(z1);
+        acb_clear(z2);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}
+

acb_dirichlet/zeta_rs.c

@@ -0,0 +1,84 @@
+/*
+    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"
+
+/* todo: decompose s into midpoint+radius */
+void
+acb_dirichlet_zeta_rs(acb_t res, const acb_t s, slong K, slong prec)
+{
+    acb_t R1, R2, X, t;
+    slong wp;
+
+    if (arf_sgn(arb_midref(acb_imagref(s))) < 0)
+    {
+        acb_init(t);
+        acb_conj(t, s);
+        acb_dirichlet_zeta_rs(res, t, K, prec);
+        acb_conj(res, res);
+        acb_clear(t);
+        return;
+    }
+
+    acb_init(R1);
+    acb_init(R2);
+    acb_init(X);
+    acb_init(t);
+
+    /* rs_r increases the precision internally */
+    wp = prec;
+
+    acb_dirichlet_zeta_rs_r(R1, s, K, wp);
+
+    if (arb_is_exact(acb_realref(s)) &&
+        (arf_cmp_2exp_si(arb_midref(acb_realref(s)), -1) == 0))
+    {
+        acb_conj(R2, R1);
+    }
+    else
+    {
+        /* conj(R(conj(1-s))) */
+        arb_sub_ui(acb_realref(t), acb_realref(s), 1, 10 * wp);
+        arb_neg(acb_realref(t), acb_realref(t));
+        arb_set(acb_imagref(t), acb_imagref(s));
+        acb_dirichlet_zeta_rs_r(R2, t, K, wp);
+        acb_conj(R2, R2);
+    }
+
+    if (acb_is_finite(R1) && acb_is_finite(R2))
+    {
+        wp += 10 + arf_abs_bound_lt_2exp_si(arb_midref(acb_imagref(s)));
+        wp = FLINT_MAX(wp, 10);
+
+        /* X = pi^(s-1/2) gamma((1-s)/2) rgamma(s/2)
+             = (2 pi)^s rgamma(s) / (2 cos(pi s / 2)) */
+        acb_rgamma(X, s, wp);
+        acb_const_pi(t, wp);
+        acb_mul_2exp_si(t, t, 1);
+        acb_pow(t, t, s, wp);
+        acb_mul(X, X, t, wp);
+        acb_mul_2exp_si(t, s, -1);
+        acb_cos_pi(t, t, wp);
+        acb_mul_2exp_si(t, t, 1);
+        acb_div(X, X, t, wp);
+
+        acb_mul(R2, R2, X, wp);
+    }
+
+    /* R1 + X * R2 */
+    acb_add(res, R1, R2, prec);
+
+    acb_clear(R1);
+    acb_clear(R2);
+    acb_clear(X);
+    acb_clear(t);
+}
+

acb_dirichlet/zeta_rs_bound.c

@@ -0,0 +1,90 @@
+/*
+    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"
+
+/* Arias de Reyna, Theorem 4.2 */
+void
+acb_dirichlet_zeta_rs_bound(mag_t err, const acb_t s, slong K)
+{
+    arb_t a;
+    mag_t c1, c2, c3;
+    slong e;
+
+    if (!arb_is_positive(acb_imagref(s)) || K < 1 || !acb_is_finite(s))
+    {
+        mag_inf(err);
+        return;
+    }
+
+    arb_init(a);
+
+    arb_add_ui(a, acb_realref(s), K, MAG_BITS);
+    arb_sub_ui(a, a, 2, MAG_BITS);
+
+    if (!arb_is_nonnegative(acb_realref(s)) && !arb_is_nonnegative(a))
+    {
+        mag_inf(err);
+        arb_clear(a);
+        return;
+    }
+
+    mag_init(c1);
+    mag_init(c2);
+    mag_init(c3);
+
+    /* c1 = 1/7 2^(3 sigma / 2)  re(sigma) >= 0 */
+    /* c1 < 1/2                  re(sigma) < 0  */
+    arf_set_mag(arb_midref(a), arb_radref(acb_realref(s)));
+    arf_add(arb_midref(a), arb_midref(a), arb_midref(acb_realref(s)), MAG_BITS, ARF_RND_CEIL);
+
+    if (arf_sgn(arb_midref(a)) <= 0)
+    {
+        mag_set_ui_2exp_si(c1, 1, -1);
+    }
+    else if (arf_cmp_2exp_si(arb_midref(a), FLINT_BITS - 4) < 0)
+    {
+        mag_one(c1);
+        mag_div_ui(c1, c1, 7);
+        e = arf_get_si(arb_midref(a), ARF_RND_CEIL);
+        mag_mul_2exp_si(c1, c1, (3 * e + 1) / 2);
+        if (mag_cmp_2exp_si(c1, -1) < 0)
+            mag_set_ui_2exp_si(c1, 1, -1);
+    }
+    else
+    {
+        mag_inf(c1);
+    }
+
+    /* c2 = 1 / ((10/11) sqrt(t/(2pi))) = (11/10) sqrt((2pi)/t) */
+    arb_get_mag_lower(c3, acb_imagref(s));
+    mag_const_pi(c2);
+    mag_mul_2exp_si(c2, c2, 1);
+    mag_div(c2, c2, c3);
+    mag_sqrt(c2, c2);
+    mag_mul_ui(c2, c2, 11);
+    mag_div_ui(c2, c2, 10);
+
+    /* c2 = c2^(K+1) */
+    mag_pow_ui(c2, c2, K + 1);
+
+    /* c3 = gamma((K+1)/2) */
+    mag_fac_ui(c3, K / 2);
+
+    mag_mul(err, c1, c2);
+    mag_mul(err, err, c3);
+
+    mag_clear(c1);
+    mag_clear(c2);
+    mag_clear(c3);
+    arb_clear(a);
+}
+