<commit_before>/*
    Copyright (C) 2010 <NAME>
    Copyright (C) 2019 <NAME>

    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 "arb_calc.h"

/*
 * For a detailed explanation of the algorithm implemented in this file, see:
 *
 * <NAME>, "Programs for Riemann's zeta function", (<NAME>
 * Vonderen, Ed.) Leven met getallen : liber amicorum ter gelegenheid van de
 * pensionering van Herman te Riele, CWI (2012) 102-112,
 * https://ir.cwi.nl/pub/19724
 */

/*
 * These constants define the largest n such that for every k <= n
 * the kth zero of the Hardy Z function is governed by Gram's law
 * or by Rosser's rule respectively.
 */
static const slong GRAMS_LAW_MAX = 126;
static const slong ROSSERS_RULE_MAX = 13999526;

/*
 * This structure describes a node of a doubly linked list.
 * Each node represents a height t at which the Hardy Z-function
 * Z(t) has been evaluated.
 *
 * As it relates to this data structure, the big picture of the algorithm
 * to isolate the nth zero is roughly:
 *
 * 1) Initialize a two-node list consisting of the two points
 *    predicted by Gram's law to surround the nth zero.
 * 2) Append and prepend as many additional Gram points as are indicated
 *    by the theory behind Turing's method, and occasionally insert points
 *    between existing points for the purpose of certifying blocks of intervals
 *    as 'good'.
 * 3) Repeatedly insert new points into the list, interleaving existing points,
 *    until the number of sign changes indicated by theory is observed.
 *    This is where the algorithm would enter an infinite loop if it encountered
 *    a counterexample to the Riemann hypothesis.
 * 4) Traverse the list until we find the pair of adjacent nodes with opposite
 *    signs of Z(t) that corresponds to the nth zero; the t heights of the
 *    two nodes define the endpoints of an isolating interval.
 * 5) Delete the list.
 */
typedef struct _zz_node_struct
{
    arf_struct t; /* height t where v = Z(t) is evaluated */
    arb_struct v; /* Z(t) */
    fmpz *gram; /* Gram point index or NULL if not a Gram point */
    slong prec; /* precision of evaluation of v */
    struct _zz_node_struct *prev;
    struct _zz_node_struct *next;
}
zz_node_struct;

typedef zz_node_struct zz_node_t[1];
typedef zz_node_struct * zz_node_ptr;
typedef const zz_node_struct * zz_node_srcptr;

static int
zz_node_is_gram_node(const zz_node_t p)
{
    return p->gram != NULL;
}

static int
zz_node_sgn(const zz_node_t p)
{
    int s = arb_sgn_nonzero(&p->v);
    if (!s)
    {
        flint_printf("unexpectedly imprecise evaluation of Z(t)\n");
        flint_abort();
    }
    return s;
}

/* Good Gram points are Gram points where sgn(Z(g(n)))*(-1)^n > 0. */
static int
zz_node_is_good_gram_node(const zz_node_t p)
{
    if (zz_node_is_gram_node(p))
    {
        int s = zz_node_sgn(p);
        if ((s > 0 && fmpz_is_even(p->gram)) ||
            (s < 0 && fmpz_is_odd(p->gram)))
        {
            return 1;
        }
    }
    return 0;
}

static void
zz_node_init(zz_node_t p)
{
    arf_init(&p->t);
    arb_init(&p->v);
    arb_indeterminate(&p->v);
    p->prec = 0;
    p->gram = NULL;
    p->prev = NULL;
    p->next = NULL;
}

static void
zz_node_clear(zz_node_t p)
{
    arf_clear(&p->t);
    arb_clear(&p->v);
    if (p->gram)
    {
        fmpz_clear(p->gram);
        flint_free(p->gram);
    }
    p->prec = 0;
    p->gram = NULL;
    p->prev = NULL;
    p->next = NULL;
}

/*
 * The node p represents the evaluation of the Hardy Z-function
 * at a point t with some precision. This function re-evaluates
 * the Hardy Z-function at t with greater precision, and it also
 * guarantees that the precision is high enough to determine the
 * sign of Z(t).
 */
static int
zz_node_refine(zz_node_t p)
{
    slong default_prec = arf_bits(&p->t) + 8;
    if (p->prec < default_prec)
    {
        p->prec = default_prec;
    }
    else
    {
        p->prec *= 2;
    }
    return _acb_dirichlet_definite_hardy_z(&p->v, &p->t, &p->prec);
}

/*
 * Create a node representing an evaluation of the Hardy Z-function
 * at arbitrary point t. Upon creation the sign of Z(t) will
 * be known with certainty.
 */
static zz_node_ptr
create_non_gram_node(const arf_t t)
{
    zz_node_ptr p = flint_malloc(sizeof(zz_node_struct));
    zz_node_init(p);
    arf_set(&p->t, t);
    zz_node_refine(p);
    return p;
}

/*
 * Create a node representing an evaluation of the Hardy Z-function
 * at the nth Gram point g(n). Upon creation a floating point number t will be
 * assigned to this node, with the property that there are no zeros of the
 * Hardy Z-function between t and the actual value of g(n).
 * The sign of Z(t) will also be known with certainty.
 */
static zz_node_ptr
create_gram_node(const fmpz_t n)
{
    zz_node_ptr p;
    arb_t t, v;
    acb_t z;
    slong prec = fmpz_bits(n) + 8;

    arb_init(t);
    arb_init(v);
    acb_init(z);

    while (1)
    {
        acb_dirichlet_gram_point(t, n, NULL, NULL, prec + fmpz_bits(n));
        acb_set_arb(z, t);
        acb_dirichlet_hardy_z(z, z, NULL, NULL, 1, prec);
        acb_get_real(v, z);
        if (!arb_contains_zero(v))
        {
            break;
        }
        prec *= 2;
    }

    p = flint_malloc(sizeof(zz_node_struct));
    zz_node_init(p);
    p->gram = flint_malloc(sizeof(fmpz));
    fmpz_init(p->gram);

    /* t contains g(n) and does not contain a zero of the Z function */
    fmpz_set(p->gram, n);
    arf_set(&p->t, arb_midref(t));
    arb_set(&p->v, v);
    p->prec = prec;

    arb_clear(t);
    arb_clear(v);
    acb_clear(z);

    return p;
}
<commit_msg>Add a comment motivating the extra precision for Gram points<commit_after>/*
    Copyright (C) 2010 <NAME>
    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"
#include "arb_calc.h"

/*
 * For a detailed explanation of the algorithm implemented in this file, see:
 *
 * <NAME>, "Programs for Riemann's zeta function", (<NAME>
 * Vonderen, Ed.) Leven met getallen : liber amicorum ter gelegenheid van de
 * pensionering van Herman te Riele, CWI (2012) 102-112,
 * https://ir.cwi.nl/pub/19724
 */

/*
 * These constants define the largest n such that for every k <= n
 * the kth zero of the Hardy Z function is governed by Gram's law
 * or by Rosser's rule respectively.
 */
static const slong GRAMS_LAW_MAX = 126;
static const slong ROSSERS_RULE_MAX = 13999526;

/*
 * This structure describes a node of a doubly linked list.
 * Each node represents a height t at which the Hardy Z-function
 * Z(t) has been evaluated.
 *
 * As it relates to this data structure, the big picture of the algorithm
 * to isolate the nth zero is roughly:
 *
 * 1) Initialize a two-node list consisting of the two points
 *    predicted by Gram's law to surround the nth zero.
 * 2) Append and prepend as many additional Gram points as are indicated
 *    by the theory behind Turing's method, and occasionally insert points
 *    between existing points for the purpose of certifying blocks of intervals
 *    as 'good'.
 * 3) Repeatedly insert new points into the list, interleaving existing points,
 *    until the number of sign changes indicated by theory is observed.
 *    This is where the algorithm would enter an infinite loop if it encountered
 *    a counterexample to the Riemann hypothesis.
 * 4) Traverse the list until we find the pair of adjacent nodes with opposite
 *    signs of Z(t) that corresponds to the nth zero; the t heights of the
 *    two nodes define the endpoints of an isolating interval.
 * 5) Delete the list.
 */
typedef struct _zz_node_struct
{
    arf_struct t; /* height t where v = Z(t) is evaluated */
    arb_struct v; /* Z(t) */
    fmpz *gram; /* Gram point index or NULL if not a Gram point */
    slong prec; /* precision of evaluation of v */
    struct _zz_node_struct *prev;
    struct _zz_node_struct *next;
}
zz_node_struct;

typedef zz_node_struct zz_node_t[1];
typedef zz_node_struct * zz_node_ptr;
typedef const zz_node_struct * zz_node_srcptr;

static int
zz_node_is_gram_node(const zz_node_t p)
{
    return p->gram != NULL;
}

static int
zz_node_sgn(const zz_node_t p)
{
    int s = arb_sgn_nonzero(&p->v);
    if (!s)
    {
        flint_printf("unexpectedly imprecise evaluation of Z(t)\n");
        flint_abort();
    }
    return s;
}

/* Good Gram points are Gram points where sgn(Z(g(n)))*(-1)^n > 0. */
static int
zz_node_is_good_gram_node(const zz_node_t p)
{
    if (zz_node_is_gram_node(p))
    {
        int s = zz_node_sgn(p);
        if ((s > 0 && fmpz_is_even(p->gram)) ||
            (s < 0 && fmpz_is_odd(p->gram)))
        {
            return 1;
        }
    }
    return 0;
}

static void
zz_node_init(zz_node_t p)
{
    arf_init(&p->t);
    arb_init(&p->v);
    arb_indeterminate(&p->v);
    p->prec = 0;
    p->gram = NULL;
    p->prev = NULL;
    p->next = NULL;
}

static void
zz_node_clear(zz_node_t p)
{
    arf_clear(&p->t);
    arb_clear(&p->v);
    if (p->gram)
    {
        fmpz_clear(p->gram);
        flint_free(p->gram);
    }
    p->prec = 0;
    p->gram = NULL;
    p->prev = NULL;
    p->next = NULL;
}

/*
 * The node p represents the evaluation of the Hardy Z-function
 * at a point t with some precision. This function re-evaluates
 * the Hardy Z-function at t with greater precision, and it also
 * guarantees that the precision is high enough to determine the
 * sign of Z(t).
 */
static int
zz_node_refine(zz_node_t p)
{
    slong default_prec = arf_bits(&p->t) + 8;
    if (p->prec < default_prec)
    {
        p->prec = default_prec;
    }
    else
    {
        p->prec *= 2;
    }
    return _acb_dirichlet_definite_hardy_z(&p->v, &p->t, &p->prec);
}

/*
 * Create a node representing an evaluation of the Hardy Z-function
 * at arbitrary point t. Upon creation the sign of Z(t) will
 * be known with certainty.
 */
static zz_node_ptr
create_non_gram_node(const arf_t t)
{
    zz_node_ptr p = flint_malloc(sizeof(zz_node_struct));
    zz_node_init(p);
    arf_set(&p->t, t);
    zz_node_refine(p);
    return p;
}

/*
 * Create a node representing an evaluation of the Hardy Z-function
 * at the nth Gram point g(n). Upon creation a floating point number t will be
 * assigned to this node, with the property that there are no zeros of the
 * Hardy Z-function between t and the actual value of g(n).
 * The sign of Z(t) will also be known with certainty.
 */
static zz_node_ptr
create_gram_node(const fmpz_t n)
{
    zz_node_ptr p;
    arb_t t, v;
    acb_t z;
    slong prec = fmpz_bits(n) + 8;

    arb_init(t);
    arb_init(v);
    acb_init(z);

    while (1)
    {
        /* Computing the Gram point to higher precision improves performance
           significantly. The likely explanation (not verified) is that
           when evaluating the Z-function at an inexact ball using the
           Riemann-Siegel formula, error propagation uses a bound for Z'
           that is far from tight. The extra precision compensates
           for this lack of tightness. */
        acb_dirichlet_gram_point(t, n, NULL, NULL, prec + fmpz_bits(n));
        acb_set_arb(z, t);
        acb_dirichlet_hardy_z(z, z, NULL, NULL, 1, prec);
        acb_get_real(v, z);
        if (!arb_contains_zero(v))
        {
            break;
        }
        prec *= 2;
    }

    p = flint_malloc(sizeof(zz_node_struct));
    zz_node_init(p);
    p->gram = flint_malloc(sizeof(fmpz));
    fmpz_init(p->gram);

    /* t contains g(n) and does not contain a zero of the Z function */
    fmpz_set(p->gram, n);
    arf_set(&p->t, arb_midref(t));
    arb_set(&p->v, v);
    p->prec = prec;

    arb_clear(t);
    arb_clear(v);
    acb_clear(z);

    return p;
}

{'fraction_non_alphanumeric': 0.10159557147508955, 'fraction_numerical': 0.009117551286225985, 'mean_word_length': 2.858354271356784, 'pattern_counts': {'":': 0, '<': 18, '<?xml version=': 0, '>': 84, 'https://': 2, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2631035', 'n_tokens_mistral': 4381, 'n_tokens_neox': 4150, 'n_words': 1458}