Merge pull request #259 from p15-git-acc/turing_method_bound

add turing_method_bound

acb_dirichlet.h

@@ -167,6 +167,7 @@ void acb_dirichlet_hardy_z_series(acb_poly_t res, const acb_poly_t s, const diri
 
 void acb_dirichlet_gram_point(arb_t res, const fmpz_t n, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec);
 void acb_dirichlet_backlund_s_bound(mag_t res, const arb_t t);
+ulong acb_dirichlet_turing_method_bound(const fmpz_t p);
 
 /* Discrete Fourier Transform */
 

acb_dirichlet/gram_point.c

@@ -12,7 +12,7 @@
 #include "acb_dirichlet.h"
 
 /*
-Claim: the error is bounded by 1/256 if n <= 1 and (1/64) (log(n)/n) if n >= 2.
+Claim: the error is bounded by 1/64 if n <= 1 and (1/64) (log(n)/n) if n >= 2.
 
 A crude lower bound for g_n is 2 pi exp(W(n)), or 8*n/log(n) for n >= 8.
 
@@ -54,7 +54,7 @@ gram_point_initial(arb_t x, const fmpz_t n, slong prec)
 
     if (fmpz_cmp_ui(n, 1) <= 0)
     {
-        mag_set_ui_2exp_si(b, 1, -8);
+        mag_set_ui_2exp_si(b, 1, -6);
     }
     else
     {
@@ -76,8 +76,8 @@ acb_dirichlet_gram_point(arb_t res, const fmpz_t n, const dirichlet_group_t G, c
 {
     slong asymp_accuracy;
 
-    /* Only implemented for n >= 0 and Riemann zeta. */
-    if (fmpz_sgn(n) < 0 || G != NULL || chi != NULL)
+    /* Only implemented for n >= -1 and Riemann zeta. */
+    if (fmpz_cmp_si(n, -1) < 0 || G != NULL || chi != NULL)
     {
         arb_indeterminate(res);
         return;

acb_dirichlet/test/t-gram_point.c

@@ -36,8 +36,8 @@ int main()
         acb_init(t);
         fmpz_init(n);
 
-        fmpz_randtest(n, state, 500);
-        fmpz_abs(n, n);
+        fmpz_randtest_unsigned(n, state, 500);
+        fmpz_sub_ui(n, n, 1);
         prec1 = 2 + n_randtest(state) % 500;
         prec2 = 2 + n_randtest(state) % 2000;
 

acb_dirichlet/test/t-turing_method_bound.c

@@ -0,0 +1,59 @@
+/*
+    Copyright (C) 2019 D.H.J Polymath
+
+    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("turing_method_bound....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 500 * arb_test_multiplier(); iter++)
+    {
+        ulong u, v;
+        fmpz_t a, b;
+
+        fmpz_init(a);
+        fmpz_init(b);
+
+        fmpz_randtest_unsigned(a, state, 500);
+        fmpz_randtest_unsigned(b, state, 500);
+        fmpz_sub_ui(a, a, 1);
+        fmpz_sub_ui(b, b, 1);
+
+        u = acb_dirichlet_turing_method_bound(a);
+        v = acb_dirichlet_turing_method_bound(b);
+
+        if ((fmpz_cmp(a, b) < 0 && u > v) ||
+            (fmpz_cmp(a, b) > 0 && u < v))
+        {
+            flint_printf("FAIL: increasing on p >= -1\n\n");
+            flint_printf("a = "); fmpz_print(a); flint_printf("\n\n");
+            flint_printf("b = "); fmpz_print(b); flint_printf("\n\n");
+            flint_printf("u = %lu\n", u);
+            flint_printf("v = %lu\n", v);
+            flint_abort();
+        }
+
+        fmpz_clear(a);
+        fmpz_clear(b);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

acb_dirichlet/turing_method_bound.c

@@ -0,0 +1,67 @@
+/*
+    Copyright (C) 2019 D.H.J Polymath
+
+    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"
+
+static void
+_mag_div_ui_ui(mag_t res, ulong a, ulong b)
+{
+    mag_set_ui(res, a);
+    mag_div_ui(res, res, b);
+}
+
+ulong
+acb_dirichlet_turing_method_bound(const fmpz_t p)
+{
+    ulong result;
+    arb_t t;
+    fmpz_t k;
+    mag_t m, b1, b2, c;
+
+    arb_init(t);
+    fmpz_init(k);
+    mag_init(m);
+    mag_init(b1);
+    mag_init(b2);
+    mag_init(c);
+
+    acb_dirichlet_gram_point(t, p, NULL, NULL, FLINT_MAX(8, fmpz_bits(p)));
+    arb_get_mag(m, t);
+    mag_log(m, m);
+
+    /* b1 = 0.0061*log(gram(p))^2 + 0.08*log(gram(p)) */
+    _mag_div_ui_ui(c, 61, 10000);
+    mag_mul(b1, c, m);
+    _mag_div_ui_ui(c, 8, 100);
+    mag_add(b1, b1, c);
+    mag_mul(b1, b1, m);
+
+    /* b2 = 0.0031*log(gram(p))^2 + 0.11*log(gram(p)) */
+    _mag_div_ui_ui(c, 31, 10000);
+    mag_mul(b2, c, m);
+    _mag_div_ui_ui(c, 11, 100);
+    mag_add(b2, b2, c);
+    mag_mul(b2, b2, m);
+
+    /* result = ceil(min(b1, b2)) */
+    mag_min(m, b1, b2);
+    mag_get_fmpz(k, m);
+    result = fmpz_get_ui(k);
+
+    arb_clear(t);
+    fmpz_clear(k);
+    mag_clear(m);
+    mag_clear(b1);
+    mag_clear(b2);
+    mag_clear(c);
+
+    return result;
+}

doc/source/acb_dirichlet.rst

@@ -615,13 +615,24 @@ Currently, these methods require *chi* to be a primitive character.
 
 .. function:: void acb_dirichlet_gram_point(arb_t res, const fmpz_t n, const dirichlet_group_t G, const dirichlet_char_t chi, slong prec)
 
-    Sets *res* to the *n*-th Gram point `g_n`, defined as the solution of
-    `\theta(g_n) = \pi n`. Currently only the Gram points corresponding to the
-    Riemann zeta function are supported and *G* and *chi* must both be set to
-    *NULL*. Requires `n \ge 0`.
+    Sets *res* to the *n*-th Gram point `g_n`, defined as the unique solution
+    in `[7, \infty)` of `\theta(g_n) = \pi n`. Currently only the Gram points
+    corresponding to the Riemann zeta function are supported and *G* and *chi*
+    must both be set to *NULL*. Requires `n \ge -1`.
 
 .. function:: void acb_dirichlet_backlund_s_bound(mag_t res, const arb_t t)
 
     Computes an upper bound for `|S(t)|` quickly. Theorem 1
     and the bounds in (1.2) in [Tru2014]_ are used.
 
+.. function:: ulong acb_dirichlet_turing_method_bound(const fmpz_t p)
+
+    Computes an upper bound *B* for the minimum number of consecutive good
+    Gram blocks sufficient to count nontrivial zeros of the Riemann zeta
+    function using Turing's method [Tur1953]_ as updated by [Leh1970]_,
+    [Bre1979]_, and [Tru2011]_.
+
+    Let `N(T)` denote the number of zeros (counted according to their
+    multiplicities) of `\zeta(s)` in the region `0 < \mathcal{I}(s) \le T`.
+    If at least *B* consecutive Gram blocks with union `[g_n, g_p)`
+    satisfy Rosser's rule, then `N(g_n) \le n + 1` and `N(g_p) \ge p + 1`.

doc/source/credits.rst

@@ -161,6 +161,8 @@ Bibliography
 
 .. [Bre1978] \R. P. Brent, "A Fortran multiple-precision arithmetic package", ACM Transactions on Mathematical Software, 4(1):57–70, March 1978.
 
+.. [Bre1979] \R. P. Brent, "On the Zeros of the Riemann Zeta Function in the Critical Strip", Mathematics of Computation 33 (1979), 1361-1372, https://doi.org/10.1090/S0025-5718-1979-0537983-2
+
 .. [Bre2010] \R. P. Brent, "Ramanujan and Euler's Constant", http://wwwmaths.anu.edu.au/~brent/pd/Euler_CARMA_10.pdf
 
 .. [BJ2013] \R. P. Brent and F. Johansson, "A bound for the error term in the Brent-McMillan algorithm", preprint (2013), http://arxiv.org/abs/1312.0039
@@ -227,6 +229,8 @@ Bibliography
 
 .. [Kri2013] \A. Krishnamoorthy and D. Menon, "Matrix Inversion Using Cholesky Decomposition" Proc. of the International Conference on Signal Processing Algorithms, Architectures, Arrangements, and Applications (SPA-2013), pp. 70-72, 2013.
 
+.. [Leh1970] \R. S. Lehman, "On the Distribution of Zeros of the Riemann Zeta-Function", Proc. of the London Mathematical Society 20(3) (1970), 303-320, https://doi.org/10.1112/plms/s3-20.2.303
+
 .. [Miy2010] \S. Miyajima, "Fast enclosure for all eigenvalues in generalized eigenvalue problems", Journal of Computational and Applied Mathematics, 233 (2010), 2994-3004, https://dx.doi.org/10.1016/j.cam.2009.11.048
 
 .. [MPFR2012] The MPFR team, "MPFR Algorithms" (2012), http://www.mpfr.org/algo.html
@@ -251,4 +255,8 @@ Bibliography
 
 .. [Tre2008] \L. N. Trefethen, "Is Gauss Quadrature Better than Clenshaw-Curtis?", SIAM Review, 50:1 (2008), 67-87, https://doi.org/10.1137/060659831
 
+.. [Tru2011] \T. S. Trudgian, "Improvements to Turing's method", Mathematics of Computation 80 (2011), 2259-2279, https://doi.org/10.1090/S0025-5718-2011-02470-1 
+
 .. [Tru2014] \T. S. Trudgian, "An improved upper bound for the argument of the Riemann zeta-function on the critical line II", Journal of Number Theory 134 (2014), 280-292, https://doi.org/10.1016/j.jnt.2013.07.017
+
+.. [Tur1953] \A. M. Turing, "Some Calculations of the Riemann Zeta-Function", Proc. of the London Mathematical Society 3(3) (1953), 99-117, https://doi.org/10.1112/plms/s3-3.1.99