Merge pull request #282 from p15-git-acc/platt_local_hardy_z_zeros

platt local zeros

acb_dirichlet.h

@@ -252,7 +252,7 @@ void acb_dirichlet_platt_ws_precomp_init(acb_dirichlet_platt_ws_precomp_t pre,
         slong A, const arb_t H, slong sigma, slong prec);
 void acb_dirichlet_platt_ws_precomp_clear(acb_dirichlet_platt_ws_precomp_t pre);
 void acb_dirichlet_platt_ws_interpolation_precomp(arb_t res,
-    acb_dirichlet_platt_ws_precomp_t pre, const arb_t t0,
+    const acb_dirichlet_platt_ws_precomp_t pre, const arb_t t0,
     arb_srcptr p, const fmpz_t T, slong A, slong B, slong Ns_max,
     const arb_t H, slong sigma, slong prec);
 void acb_dirichlet_platt_ws_interpolation(arb_t res, const arb_t t0,
@@ -288,6 +288,12 @@ void acb_dirichlet_platt_lemma_B2(arb_t out, slong K, const arb_t h,
 void acb_dirichlet_platt_multieval(arb_ptr out, const fmpz_t T, slong A,
     slong B, const arb_t h, slong J, slong K, slong sigma, slong prec);
 
+slong _acb_dirichlet_platt_local_hardy_z_zeros(
+    arb_ptr res, const fmpz_t n, slong len,
+    const fmpz_t T, slong A, slong B,
+    const arb_t h, slong J, slong K, slong sigma_grid,
+    slong Ns_max, const arb_t H, slong sigma_interp, slong prec);
+
 /* Discrete Fourier Transform */
 
 void acb_dirichlet_dft_index(acb_ptr w, acb_srcptr v, const dirichlet_group_t G, slong prec);

acb_dirichlet/isolate_hardy_z_zero.c

@@ -1093,8 +1093,8 @@ _isolate_hardy_z_zeros(arf_interval_ptr res, const fmpz_t n, slong len)
 }
 
 /* Isolate len zeros, starting from the nth zero. */
-static void
-isolate_hardy_z_zeros(arf_interval_ptr res, const fmpz_t n, slong len)
+void
+acb_dirichlet_isolate_hardy_z_zeros(arf_interval_ptr res, const fmpz_t n, slong len)
 {
     if (len <= 0)
     {
@@ -1392,7 +1392,7 @@ acb_dirichlet_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, slong prec)
     {
         slong i;
         arf_interval_ptr p = _arf_interval_vec_init(len);
-        isolate_hardy_z_zeros(p, n, len);
+        acb_dirichlet_isolate_hardy_z_zeros(p, n, len);
         for (i = 0; i < len; i++)
         {
             _acb_dirichlet_refine_hardy_z_zero(res + i, &p[i].a, &p[i].b, prec);

acb_dirichlet/platt_ws_interpolation.c

@@ -390,7 +390,7 @@ acb_dirichlet_platt_ws_precomp_clear(acb_dirichlet_platt_ws_precomp_t pre)
 }
 
 void acb_dirichlet_platt_ws_interpolation_precomp(arb_t res,
-    acb_dirichlet_platt_ws_precomp_t pre, const arb_t t0,
+    const acb_dirichlet_platt_ws_precomp_t pre, const arb_t t0,
     arb_srcptr p, const fmpz_t T, slong A, slong B, slong Ns_max,
     const arb_t H, slong sigma, slong prec)
 {

acb_dirichlet/test/t-platt_local_hardy_z_zeros.c

@@ -0,0 +1,87 @@
+/*
+    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()
+{
+    /* Check a specific combination of parameter values that is relatively fast
+     * to evaluate and that has relatively tight bounds. */
+    slong A, B, J, K, sigma_grid, Ns_max, sigma_interp;
+    arb_t h, H;
+    fmpz_t T, n;
+    arb_ptr pa, pb;
+    slong count, i;
+    slong maxcount = 50;
+    slong prec = 128;
+
+    flint_printf("platt_local_hardy_z_zeros....");
+    fflush(stdout);
+
+    arb_init(h);
+    arb_init(H);
+    fmpz_init(T);
+    fmpz_init(n);
+    pa = _arb_vec_init(maxcount);
+    pb = _arb_vec_init(maxcount);
+
+    fmpz_set_si(n, 10142);
+
+    /* parameters related to the location/resolution/width of the grid */
+    fmpz_set_si(T, 10000);
+    A = 8;
+    B = 128;
+
+    /* tuning parameters for the evaluation of grid points */
+    J = 1000;
+    K = 30;
+    sigma_grid = 63;
+    arb_set_d(h, 4.5);
+
+    /* tuning parameters for interpolation on the grid */
+    Ns_max = 200;
+    sigma_interp = 21;
+    arb_one(H);
+
+    count = _acb_dirichlet_platt_local_hardy_z_zeros(pa, n, maxcount,
+            T, A, B, h, J, K, sigma_grid, Ns_max, H, sigma_interp, prec);
+    acb_dirichlet_hardy_z_zeros(pb, n, count, prec);
+    if (count != maxcount)
+    {
+        flint_printf("FAIL: not enough zeros were isolated\n\n");
+        flint_printf("count = %wd  maxcount = %wd\n\n", count, maxcount);
+        flint_abort();
+    }
+
+    for (i = 0; i < count; i++)
+    {
+        if (!arb_overlaps(pa+i, pb+i))
+        {
+            flint_printf("FAIL: overlap\n\n");
+            flint_printf("observed[%wd] = ", i);
+            arb_printd(pa+i, 20); flint_printf("\n\n");
+            flint_printf("expected[%wd] = ", i);
+            arb_printd(pb+i, 20); flint_printf("\n\n");
+            flint_abort();
+        }
+    }
+
+    arb_clear(h);
+    arb_clear(H);
+    fmpz_clear(T);
+    fmpz_clear(n);
+    _arb_vec_clear(pa, maxcount);
+    _arb_vec_clear(pb, maxcount);
+
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

doc/source/acb_dirichlet.rst

@@ -777,3 +777,12 @@ and formulas described by David J. Platt in [Pla2017]_.
     interpolation. *sigma* is an odd positive integer tuning parameter
     `\sigma \in 2\mathbb{Z}_{>0}+1` used in computing error bounds.
 
+.. function:: slong _acb_dirichlet_platt_local_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, const fmpz_t T, slong A, slong B, const arb_t h, slong J, slong K, slong sigma_grid, slong Ns_max, const arb_t H, slong sigma_interp, slong prec)
+
+    Sets the entries of *res* to at most *len* consecutive zeros of the
+    Hardy Z-function, beginning with the *n*-th zero. The number of zeros
+    isolated near *T* is returned. Requires positive *n*.
+    Internally this function uses Platt's grid evaluation of the scaled
+    Lambda function, and the final several parameters have the same meanings
+    as in the functions :func:`acb_dirichlet_platt_multieval`
+    and :func:`acb_dirichlet_platt_ws_interpolation`.