add functions specific to Gram points

acb_dirichlet.h

@@ -180,6 +180,8 @@ void _acb_dirichlet_exact_zeta_nzeros(fmpz_t res, const arf_t t);
 void acb_dirichlet_zeta_nzeros(arb_t res, const arb_t t, slong prec);
 void acb_dirichlet_backlund_s(arb_t res, const arb_t t, slong prec);
 void acb_dirichlet_backlund_s_bound(mag_t res, const arb_t t);
+void acb_dirichlet_zeta_nzeros_gram(fmpz_t res, const fmpz_t n);
+slong acb_dirichlet_backlund_s_gram(const fmpz_t n);
 
 ACB_DIRICHLET_INLINE void
 acb_dirichlet_hardy_z_zero(arb_t res, const fmpz_t n, slong prec)

acb_dirichlet/backlund_s_gram.c

@@ -0,0 +1,33 @@
+/*
+    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"
+
+slong
+acb_dirichlet_backlund_s_gram(const fmpz_t n)
+{
+    slong res = 0;
+    if (fmpz_cmp_si(n, -1) < 0)
+    {
+        flint_printf("n must be >= -1\n");
+        flint_abort();
+    }
+    else
+    {
+        fmpz_t k;
+        fmpz_init(k);
+        acb_dirichlet_zeta_nzeros_gram(k, n);
+        fmpz_sub(k, k, n);
+        res = fmpz_get_si(k) - 1;
+        fmpz_clear(k);
+    }
+    return res;
+}

acb_dirichlet/isolate_hardy_z_zero.c

@@ -1017,16 +1017,18 @@ _separated_gram_list(zz_node_ptr *pu, zz_node_ptr *pv, const fmpz_t n)
 }
 
 /*
- * Isolate up to len zeros, starting from the nth zero.
+ * Get a list of points that fully separate zeros of Z.
- * Return the number of isolated zeros.
+ * Used for both zero isolation and for N(t).
+ * n is the index of a Hardy Z-function zero of interest.
+ * *pu and *pv are the first and last node of an output list.
+ * *pU and *pV are the first and last nodes of a sublist with
+ * fully separated zeros, within the *pu -- *pv list.
  */
-static slong
+static void
-_isolate_hardy_z_zeros(arf_interval_ptr res, const fmpz_t n, slong len)
+_separated_list(zz_node_ptr *pU, zz_node_ptr *pV,
+        zz_node_ptr *pu, zz_node_ptr *pv, const fmpz_t n)
 {
     zz_node_ptr U, V, u, v;
-    slong count;
-
-    /* Get a list of points that fully separate zeros of Z. */
     if (fmpz_cmp_si(n, GRAMS_LAW_MAX) <= 0)
     {
         _separated_gram_list(&u, &v, n);
@@ -1044,6 +1046,39 @@ _isolate_hardy_z_zeros(arf_interval_ptr res, const fmpz_t n, slong len)
         _separated_turing_list(&U, &V, &u, &v, n);
     }
 
+    if (U == NULL || V == NULL)
+    {
+        flint_printf("U and V must not be NULL\n");
+        flint_abort();
+    }
+    if (!zz_node_is_good_gram_node(U) || !zz_node_is_good_gram_node(V))
+    {
+        flint_printf("U and V must be good Gram points\n");
+        flint_abort();
+    }
+    if (U == V)
+    {
+        flint_printf("the list must include at least one interval\n");
+        flint_abort();
+    }
+
+    *pU = U;
+    *pV = V;
+    *pu = u;
+    *pv = v;
+}
+
+/*
+ * Isolate up to len zeros, starting from the nth zero.
+ * Return the number of isolated zeros.
+ */
+static slong
+_isolate_hardy_z_zeros(arf_interval_ptr res, const fmpz_t n, slong len)
+{
+    zz_node_ptr U, V, u, v;
+    slong count;
+
+    _separated_list(&U, &V, &u, &v, n);
     count = count_up_separated_zeros(res, U, V, n, len);
 
     while (u)
@@ -1196,6 +1231,7 @@ gram_index(fmpz_t res, const arf_t t)
             }
             prec *= 2;
         }
+        acb_clear(z);
     }
 }
 
@@ -1228,38 +1264,7 @@ _exact_zeta_multi_nzeros(fmpz *res, arf_srcptr points, slong len)
     fmpz_add_ui(n, n, 2);
 
     /* Get a list of points that fully separate zeros of Z. */
-    if (fmpz_cmp_si(n, GRAMS_LAW_MAX) <= 0)
+    _separated_list(&U, &V, &u, &v, n);
-    {
-        _separated_gram_list(&u, &v, n);
-        U = u;
-        V = v;
-    }
-    else if (fmpz_cmp_si(n, ROSSERS_RULE_MAX) <= 0)
-    {
-        _separated_rosser_list(&u, &v, n);
-        U = u;
-        V = v;
-    }
-    else
-    {
-        _separated_turing_list(&U, &V, &u, &v, n);
-    }
-
-    if (U == NULL || V == NULL)
-    {
-        flint_printf("U and V must not be NULL\n");
-        flint_abort();
-    }
-    if (!zz_node_is_good_gram_node(U) || !zz_node_is_good_gram_node(V))
-    {
-        flint_printf("U and V must be good Gram points\n");
-        flint_abort();
-    }
-    if (U == V)
-    {
-        flint_printf("the list must include at least one interval\n");
-        flint_abort();
-    }
 
     p = U;
     fmpz_add_ui(N, U->gram, 1);
@@ -1322,6 +1327,7 @@ _exact_zeta_multi_nzeros(fmpz *res, arf_srcptr points, slong len)
 
     arb_clear(x);
     fmpz_clear(n);
+    fmpz_clear(N);
 
     return i;
 }
@@ -1445,3 +1451,69 @@ acb_dirichlet_zeta_nzeros(arb_t res, const arb_t t, slong prec)
 
     arb_set_round(res, res, prec);
 }
+
+void
+acb_dirichlet_zeta_nzeros_gram(fmpz_t res, const fmpz_t n)
+{
+    zz_node_ptr U, V, u, v, p;
+    fmpz_t k, N;
+
+    if (fmpz_cmp_si(n, -1) < 0)
+    {
+        flint_printf("n must be >= -1\n");
+        flint_abort();
+    }
+
+    fmpz_init(k);
+    fmpz_init(N);
+
+    /*
+     * Get a list of points that fully separate zeros of Z.
+     * The kth zero is expected to be between the k-2 and the k-1 gram points.
+     */
+    fmpz_add_ui(k, n, 2);
+    _separated_list(&U, &V, &u, &v, k);
+
+    p = U;
+    fmpz_add_ui(N, U->gram, 1);
+    fmpz_set_si(res, -1);
+    while (1)
+    {
+        if (p == NULL)
+        {
+            flint_printf("prematurely reached the end of the list\n");
+            flint_abort();
+        }
+        if (zz_node_is_gram_node(p) && fmpz_equal(n, p->gram))
+        {
+            fmpz_set(res, N);
+            break;
+        }
+        if (zz_node_sgn(p) != zz_node_sgn(p->next))
+        {
+            fmpz_add_ui(N, N, 1);
+        }
+        if (p == V)
+        {
+            break;
+        }
+        p = p->next;
+    }
+
+    if (fmpz_sgn(res) < 0)
+    {
+        flint_printf("failed to find the gram point\n");
+        flint_abort();
+    }
+
+    while (u)
+    {
+        v = u;
+        u = u->next;
+        zz_node_clear(v);
+        flint_free(v);
+    }
+
+    fmpz_clear(k);
+    fmpz_clear(N);
+}

acb_dirichlet/test/t-backlund_s_gram.c

@@ -0,0 +1,72 @@
+/*
+    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("backlund_s_gram....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 130 + 20 * arb_test_multiplier(); iter++)
+    {
+        arb_t t, x;
+        fmpz_t n;
+        slong S, prec1, prec2;
+
+        arb_init(t);
+        arb_init(x);
+        fmpz_init(n);
+
+        if (iter < 130)
+        {
+            fmpz_set_si(n, iter - 1);
+        }
+        else
+        {
+            fmpz_randtest_unsigned(n, state, 20);
+            fmpz_add_ui(n, n, 129);
+        }
+
+        prec1 = 2 + n_randtest(state) % 100;
+        prec2 = 2 + n_randtest(state) % 100;
+
+        S = acb_dirichlet_backlund_s_gram(n);
+
+        acb_dirichlet_gram_point(t, n, NULL, NULL, prec1);
+        acb_dirichlet_backlund_s(x, t, prec2);
+
+        if (!arb_contains_si(x, S))
+        {
+            flint_printf("FAIL: containment\n\n");
+            flint_printf("n = "); fmpz_print(n);
+            flint_printf("   prec1 = %wd  prec2 = %wd\n\n", prec1, prec2);
+            flint_printf("S = %ld\n\n", S);
+            flint_printf("t = "); arb_printn(t, 100, 0); flint_printf("\n\n");
+            flint_printf("x = "); arb_printn(x, 100, 0); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        arb_clear(t);
+        arb_clear(x);
+        fmpz_clear(n);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

acb_dirichlet/test/t-zeta_nzeros.c

@@ -23,7 +23,7 @@ int main()
 
     for (iter = 0; iter < 130 + 20 * arb_test_multiplier(); iter++)
     {
-        arb_t x, y, t, u, v, a;
+        arb_t x, y, t;
         arf_t f1, f2;
         fmpz_t n, m, k1, k2;
         acb_t z;
@@ -32,9 +32,6 @@ int main()
         arb_init(x);
         arb_init(y);
         arb_init(t);
-        arb_init(u);
-        arb_init(v);
-        arb_init(a);
         acb_init(z);
         arf_init(f1);
         arf_init(f2);
@@ -50,7 +47,7 @@ int main()
         else
         {
             fmpz_randtest_unsigned(n, state, 20);
-            fmpz_add_ui(n, n, 1);
+            fmpz_add_ui(n, n, 131);
         }
 
         prec1 = 2 + n_randtest(state) % 50;
@@ -64,7 +61,7 @@ int main()
         if (!arb_contains_fmpz(x, n) || !arb_contains_fmpz(x, m))
         {
             flint_printf("FAIL: zero containment\n\n");
-            flint_printf("n = "); fmpz_print(n); flint_printf("\n\n");
+            flint_printf("n = "); fmpz_print(n);
             flint_printf("   prec1 = %wd  prec2 = %wd\n\n", prec1, prec2);
             flint_printf("t = "); arb_printn(t, 100, 0); flint_printf("\n\n");
             flint_printf("x = "); arb_printn(x, 100, 0); flint_printf("\n\n");
@@ -116,9 +113,6 @@ int main()
         arb_clear(x);
         arb_clear(y);
         arb_clear(t);
-        arb_clear(u);
-        arb_clear(v);
-        arb_clear(a);
         acb_clear(z);
         arf_clear(f1);
         arf_clear(f2);

acb_dirichlet/test/t-zeta_nzeros_gram.c

@@ -0,0 +1,74 @@
+/*
+    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("zeta_nzeros_gram....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 130 + 20 * arb_test_multiplier(); iter++)
+    {
+        arb_t t, x;
+        fmpz_t N, n;
+        slong prec1, prec2;
+
+        arb_init(t);
+        arb_init(x);
+        fmpz_init(n);
+        fmpz_init(N);
+
+        if (iter < 130)
+        {
+            fmpz_set_si(n, iter - 1);
+        }
+        else
+        {
+            fmpz_randtest_unsigned(n, state, 20);
+            fmpz_add_ui(n, n, 129);
+        }
+
+        prec1 = 2 + n_randtest(state) % 100;
+        prec2 = 2 + n_randtest(state) % 100;
+
+        acb_dirichlet_zeta_nzeros_gram(N, n);
+
+        acb_dirichlet_gram_point(t, n, NULL, NULL, prec1);
+        acb_dirichlet_zeta_nzeros(x, t, prec2);
+
+        if (!arb_contains_fmpz(x, N))
+        {
+            flint_printf("FAIL: containment\n\n");
+            flint_printf("n = "); fmpz_print(n);
+            flint_printf("   prec1 = %wd  prec2 = %wd\n\n", prec1, prec2);
+            flint_printf("N = "); fmpz_print(N); flint_printf("\n\n");
+            flint_printf("t = "); arb_printn(t, 100, 0); flint_printf("\n\n");
+            flint_printf("x = "); arb_printn(x, 100, 0); flint_printf("\n\n");
+            flint_abort();
+        }
+
+        arb_clear(t);
+        arb_clear(x);
+        fmpz_clear(n);
+        fmpz_clear(N);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return EXIT_SUCCESS;
+}

doc/source/acb_dirichlet.rst

@@ -722,3 +722,14 @@ mpmath [Joh2018b]_ by Juan Arias de Reyna, described in [Ari2012]_.
     Compute an upper bound for `|S(t)|` quickly. Theorem 1
     and the bounds in (1.2) in [Tru2014]_ are used.
 
+.. function:: void acb_dirichlet_zeta_nzeros_gram(fmpz_t res, const fmpz_t n)
+
+    Compute `N(g_n)`. That is, compute the number of zeros (counted according
+    to their multiplicities) of `\zeta(s)` in the region
+    `0 < \operatorname{Im}(s) \le g_n` where `g_n` is the *n*-th Gram point.
+    Requires `n \ge -1`.
+
+.. function:: slong acb_dirichlet_backlund_s_gram(const fmpz_t n)
+
+    Compute `S(g_n)` where `g_n` is the *n*-th Gram point. Requires `n \ge -1`.
+