/*
    Copyright (C) 2010 Juan Arias de Reyna
    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:
 *
 * J. Arias de Reyna, "Programs for Riemann's zeta function", (J. A. J. van
 * 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;
}

/*
 * Count the number of Gram intervals between the Gram point
 * represented by node a and the Gram point represented by node b.
 * Traversing the linked list is not necessary because the Gram indices
 * of nodes a and b can be accessed directly.
 */
static slong
count_gram_intervals(zz_node_srcptr a, zz_node_srcptr b)
{
    slong out = 0;
    if (!a || !b)
    {
        flint_printf("a and b must be non-NULL\n");
        flint_abort();
    }
    if (!zz_node_is_good_gram_node(a) || !zz_node_is_good_gram_node(b))
    {
        flint_printf("both nodes must be good Gram points\n");
        flint_abort();
    }
    else
    {
        fmpz_t m;
        fmpz_init(m);
        fmpz_sub(m, b->gram, a->gram);
        out = fmpz_get_si(m);
        fmpz_clear(m);
    }
    return out;
}

/*
 * Count the observed number of sign changes of Z(t) by traversing
 * a linked list of evaluated points from node a to node b.
 */
static slong
count_sign_changes(zz_node_srcptr a, zz_node_srcptr b)
{
    zz_node_srcptr p, q;
    slong n = 0;
    if (!a || !b)
    {
        flint_printf("a and b must be non-NULL\n");
        flint_abort();
    }
    p = a;
    q = a->next;
    while (p != b)
    {
        if (!q)
        {
            flint_printf("prematurely reached end of list\n");
            flint_abort();
        }
        if (zz_node_sgn(p) != zz_node_sgn(q))
        {
            n++;
        }
        p = q;
        q = q->next;
    }
    return n;
}

/*
 * Modify a linked list that ends with node p,
 * by appending nodes representing Gram points.
 * Continue until a 'good' Gram point is found.
 */
static zz_node_ptr
extend_to_next_good_gram_node(zz_node_t p)
{
    fmpz_t n;
    zz_node_ptr q, r;

    fmpz_init(n);

    if (!zz_node_is_gram_node(p))
    {
        flint_printf("expected to begin at a gram point\n");
        flint_abort();
    }
    if (p->next)
    {
        flint_printf("expected to extend from the end of a list\n");
        flint_abort();
    }
    fmpz_set(n, p->gram);
    q = p;
    while (1)
    {
        fmpz_add_ui(n, n, 1);
        r = create_gram_node(n);
        q->next = r;
        r->prev = q;
        q = r;
        r = NULL;
        if (zz_node_is_good_gram_node(q))
        {
            break;
        }
    }

    fmpz_clear(n);

    return q;
}

/*
 * Modify a linked list that begins with node p,
 * by prepending nodes representing Gram points.
 * Continue until a 'good' Gram point is found.
 */
static zz_node_ptr
extend_to_prev_good_gram_node(zz_node_t p)
{
    fmpz_t n;
    zz_node_ptr q, r;

    fmpz_init(n);

    if (!zz_node_is_gram_node(p))
    {
        flint_printf("expected to begin at a gram point\n");
        flint_abort();
    }
    if (p->prev)
    {
        flint_printf("expected to extend from the start of a list\n");
        flint_abort();
    }
    fmpz_set(n, p->gram);
    q = p;
    while (1)
    {
        fmpz_sub_ui(n, n, 1);
        r = create_gram_node(n);
        q->prev = r;
        r->next = q;
        q = r;
        r = NULL;
        if (zz_node_is_good_gram_node(q))
        {
            break;
        }
    }

    fmpz_clear(n);

    return q;
}

/*
 * TODO: This is probably redundant.
 * Can arb_get_lbound_arf ever give a negative arf given a nonnegative arb?
 * If the answer is no, and it's probably no, then this function
 * can be deleted.
 */
static void
_arb_get_lbound_arf_nonnegative(arf_t res, const arb_t x, slong prec)
{
    arb_get_lbound_arf(res, x, prec);
    if (arf_cmp_si(res, 0) < 0)
    {
        arf_zero(res);
    }
}

/*
 * res = (x1*w1 + x2*w2) / (w1 + w2)
 * Undefined if weights are not nonnegative.
 * If w1 and w2 are zero, the resulting interval contains x1 and x2.
 */
static void
_weighted_arithmetic_mean(arb_t res, const arf_t x1, const arf_t x2,
        const arb_t w1, const arb_t w2, slong prec)
{
    if (!arb_is_nonnegative(w1) || !arb_is_nonnegative(w2))
    {
        arb_indeterminate(res);
    }
    else if (arb_is_zero(w1) && arb_is_zero(w2))
    {
        arb_set_interval_arf(res, x1, x2, prec);
    }
    else if (arb_is_zero(w1))
    {
        arb_set_arf(res, x2);
    }
    else if (arb_is_zero(w2))
    {
        arb_set_arf(res, x1);
    }
    else if (arb_is_exact(w1) && arb_is_exact(w2))
    {
        arb_t a, b;
        arb_init(a);
        arb_init(b);
        arb_mul_arf(a, w1, x1, prec);
        arb_addmul_arf(a, w2, x2, prec);
        arb_add(b, w1, w2, prec);
        arb_div(res, a, b, prec);
        arb_clear(a);
        arb_clear(b);
    }
    else
    {
        arb_t a, b, r1, r2;
        arb_init(a);
        arb_init(b);
        arb_init(r1);
        arb_init(r2);

        arb_zero(a);
        arb_zero(b);
        _arb_get_lbound_arf_nonnegative(arb_midref(a), w1, prec);
        arb_get_ubound_arf(arb_midref(b), w2, prec);
        _weighted_arithmetic_mean(r1, x1, x2, a, b, prec);

        arb_zero(a);
        arb_zero(b);
        arb_get_ubound_arf(arb_midref(a), w1, prec);
        _arb_get_lbound_arf_nonnegative(arb_midref(b), w2, prec);
        _weighted_arithmetic_mean(r2, x1, x2, a, b, prec);

        arb_union(res, r1, r2, prec);

        arb_clear(a);
        arb_clear(b);
        arb_clear(r1);
        arb_clear(r2);
    }
}

/*
 * Split the interval (t1, t2) into the intervals (t1, out) and (out, t2)
 * in an attempt to increase the number of observed sign changes of Z(t)
 * between endpoints.
 * v1 and v2 are the Hardy Z-function values at t1 and t2 respectively.
 * sign1 and sign2 are the signs of v1 and v2 respectively.
 */
static void
split_interval(arb_t out,
        const arf_t t1, const arb_t v1, slong sign1,
        const arf_t t2, const arb_t v2, slong sign2, slong prec)
{
    if (sign1 == sign2)
    {
        /*
         * out = (sqrt(v2/v1)*t1 + t2) / (sqrt(v2/v1) + 1)
         * We have f(t1)=v1, f(t2)=v2 where v1 and v2 have the same sign,
         * and we want to guess t between t1 and t2 so that f(t)
         * has the opposite sign. Try the vertex of a parabola that would touch
         * f(t)=0 between t1 and t2 and would pass through (t1,v1) and (t2,v2).
         */
        arb_t w1, w2;
        arb_init(w1);
        arb_init(w2);
        arb_abs(w1, v2); /* w1, v2 is deliberate */
        arb_sqrt(w1, w1, prec);
        arb_abs(w2, v1); /* w2, v1 is deliberate */
        arb_sqrt(w2, w2, prec);
        _weighted_arithmetic_mean(out, t1, t2, w1, w2, prec);
        arb_clear(w1);
        arb_clear(w2);
    }
    else
    {
        /*
         * out = (t1 + t2) / 2
         * There is already one sign change in this interval.
         * To find additional sign changes we would need to evaluate
         * at least two more points in the interval,
         * so begin by just splitting the interval in half at the midpoint.
         */
        arb_set_arf(out, t1);
        arb_add_arf(out, out, t2, prec);
        arb_mul_2exp_si(out, out, -1);
    }
}

/*
 * Add a new node between each pair of existing nodes in the linked list
 * of evaluated values of t, within the sublist demarcated by nodes a and b.
 */
static void
intercalate(zz_node_t a, zz_node_t b)
{
    arb_t t;
    zz_node_ptr q, r, mid_node;

    if (a == NULL || b == NULL)
    {
        flint_printf("a and b must be non-NULL\n");
        flint_abort();
    }
    if (!zz_node_is_good_gram_node(a) || !zz_node_is_good_gram_node(b))
    {
        flint_printf("a and b must represent good Gram points\n");
        flint_abort();
    }

    if (a == b) return;

    arb_init(t);

    q = a;
    r = a->next;
    while (q != b)
    {
        if (!r)
        {
            flint_printf("prematurely reached end of list\n");
            flint_abort();
        }
        while (1)
        {
            split_interval(t,
                    &q->t, &q->v, zz_node_sgn(q),
                    &r->t, &r->v, zz_node_sgn(r),
                    FLINT_MIN(q->prec, r->prec));
            if (!arb_contains_arf(t, &q->t) && !arb_contains_arf(t, &r->t))
            {
                break;
            }
            zz_node_refine((q->prec < r->prec) ? q : r);
        }
        mid_node = create_non_gram_node(arb_midref(t));
        q->next = mid_node;
        mid_node->prev = q;
        mid_node->next = r;
        r->prev = mid_node;
        q = r;
        r = r->next;
    }

    arb_clear(t);
}

/*
 * Virtually trim k Gram/Rosser blocks from the head and from the tail
 * of the sublist (a, b). The resulting sublist is (*A, *B).
 * No nodes or connections between nodes are modified, added, or deleted.
 */
static void
trim(zz_node_ptr *A, zz_node_ptr *B,
        zz_node_ptr a, zz_node_ptr b, slong k)
{
    slong n;
    for (n = 0; n < k; n++)
    {
        a = a->next;
        while (!zz_node_is_good_gram_node(a))
        {
            a = a->next;
        }
        b = b->prev;
        while (!zz_node_is_good_gram_node(b))
        {
            b = b->prev;
        }
    }
    *A = a;
    *B = b;
}

/*
 * Given a linked sublist beginning at U and ending at V defining function
 * evaluations at points that fully separate zeros of Z(t) in the vicinity
 * of the nth zero, traverse the list until the nth zero is found.
 * Continue traversing the list until len consecutive isolating intervals
 * have been found, or until the end of the sublist is reached.
 * Return the number of isolated zeros found, starting at the nth zero.
 */
static slong
count_up_separated_zeros(arf_interval_ptr res,
        zz_node_srcptr U, zz_node_srcptr V, const fmpz_t n, slong len)
{
    if (len <= 0)
    {
        return 0;
    }
    else if (fmpz_sgn(n) < 1)
    {
        flint_printf("nonpositive indices of zeros are not supported\n");
        flint_abort();
    }
    else 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();
    }
    else
    {
        slong i = 0;
        zz_node_srcptr p = U;
        fmpz_t N, k;
        fmpz_init(N);
        fmpz_init(k);
        fmpz_add_ui(N, p->gram, 1);
        fmpz_set(k, n);
        while (p != V)
        {
            if (!p->next)
            {
                flint_printf("prematurely reached end of list\n");
                flint_abort();
            }
            if (zz_node_sgn(p) != zz_node_sgn(p->next))
            {
                fmpz_add_ui(N, N, 1);
                if (fmpz_equal(N, k))
                {
                    arf_set(&res[i].a, &p->t);
                    arf_set(&res[i].b, &p->next->t);
                    fmpz_add_ui(k, k, 1);
                    i++;
                    if (i == len)
                        break;
                }
            }
            p = p->next;
        }
        fmpz_clear(k);
        return i;
    }
    return 0;
}

/*
 * Find one 'superblock' below n and one 'superblock' above n.
 * The term 'superblock' in this context means a stretch of
 * enough consecutive 'good' Rosser/Gram blocks to meet the Turing method bound.
 * The output *psb is the number of blocks in the superblock.
 * The output nodes *pu and *pv are the first node of the first superblock
 * and the last node of the second superblock respectively.
 */
static void
turing_search_near(zz_node_ptr *pu, zz_node_ptr *pv, slong *psb, const fmpz_t n)
{
    zz_node_ptr u, v;
    slong i;
    slong zn; /* target number of sign changes */
    slong sb; /* the number of padding blocks required by Turing's method */
    slong cgb; /* the number of consecutive good blocks */
    const slong loopcount = 4;
    fmpz_t k;

    fmpz_init(k);

    fmpz_sub_ui(k, n, 2);
    u = create_gram_node(k);
    fmpz_sub_ui(k, n, 1);
    v = create_gram_node(k);
    u->next = v;
    v->prev = u;

    if (!zz_node_is_good_gram_node(u))
        u = extend_to_prev_good_gram_node(u);
    if (!zz_node_is_good_gram_node(v))
        v = extend_to_next_good_gram_node(v);

    /*
     * Extend the search to greater heights t.
     * 
     * Continue adding Gram points until the number of consecutive
     * 'good' Gram/Rosser blocks reaches the Turing method bound.
     */
    sb = 0;
    cgb = 0;
    while (1)
    {
        zz_node_ptr nv;
        nv = extend_to_next_good_gram_node(v);
        zn = count_gram_intervals(v, nv);
        for (i = 0; i < loopcount && count_sign_changes(v, nv) < zn; i++)
        {
            intercalate(v, nv);
        }
        if (count_sign_changes(v, nv) >= zn)
        {
            cgb++;
            if (cgb > sb)
            {
                sb = cgb;
                if (acb_dirichlet_turing_method_bound(nv->gram) <= sb)
                {
                    v = nv;
                    break;
                }
            }
        }
        else
        {
            cgb = 0;
        }
        v = nv;
    }

    /* Extend the search to smaller heights t. */
    cgb = 0;
    while (1)
    {
        zz_node_ptr pu;
        pu = extend_to_prev_good_gram_node(u);
        zn = count_gram_intervals(pu, u);
        for (i = 0; i < loopcount && count_sign_changes(pu, u) < zn; i++)
        {
            intercalate(pu, u);
        }
        if (count_sign_changes(pu, u) >= zn)
        {
            cgb++;
            if (cgb == sb)
            {
                u = pu;
                break;
            }
        }
        else
        {
            cgb = 0;
        }
        u = pu;
    }

    *pu = u;
    *pv = v;
    *psb = sb;

    fmpz_clear(k);
}

/*
 * Find one 'double superblock' beginning below the point represented
 * by node u, and find one 'double superblock' ending above the point
 * represented by node v. The term 'double superblock' in this context
 * means a stretch of twice as many consecutive 'good' Rosser/Gram blocks
 * as would meet the Turing method bound.
 * The output nodes *pu and *pv are the first node of the first double
 * superblock and the last node of the second double superblock respectively.
 * The output integer *psb reports one half of the number
 * of blocks in the double superblock.
 * The parameter initial_cgb is the number of consecutive good blocks
 * above and below the points represented by nodes u and v respectively.
 */
static void
turing_search_far(zz_node_ptr *pu, zz_node_ptr *pv, slong *psb,
        zz_node_ptr u, zz_node_ptr v, slong initial_cgb)
{
    slong i;
    slong zn; /* target number of sign changes */
    slong sb; /* the number of padding blocks required by Turing's method */
    slong cgb; /* the number of consecutive good blocks */
    const slong loopcount = 4;

    /*
     * Extend the search to greater heights t.
     * 
     * Continue adding Gram points until the number of consecutive
     * 'good' Gram/Rosser blocks is twice the number required by
     * the Turing method bound.
     */
    sb = 0;
    cgb = initial_cgb;
    while (1)
    {
        zz_node_ptr nv;
        nv = extend_to_next_good_gram_node(v);
        zn = count_gram_intervals(v, nv);
        for (i = 0; i < loopcount && count_sign_changes(v, nv) < zn; i++)
        {
            intercalate(v, nv);
        }
        if (count_sign_changes(v, nv) >= zn)
        {
            cgb++;
            if (cgb % 2 == 0 && sb < cgb / 2)
            {
                sb = cgb / 2;
                if (acb_dirichlet_turing_method_bound(nv->gram) <= sb)
                {
                    v = nv;
                    break;
                }
            }
        }
        else
        {
            cgb = 0;
        }
        v = nv;
    }

    /* Extend the search to smaller heights t. */
    cgb = initial_cgb;
    while (1)
    {
        zz_node_ptr pu;
        pu = extend_to_prev_good_gram_node(u);
        zn = count_gram_intervals(pu, u);
        for (i = 0; i < loopcount && count_sign_changes(pu, u) < zn; i++)
        {
            intercalate(pu, u);
        }
        if (count_sign_changes(pu, u) >= zn)
        {
            cgb++;
            if (cgb == sb*2)
            {
                u = pu;
                break;
            }
        }
        else
        {
            cgb = 0;
        }
        u = pu;
    }

    *pu = u;
    *pv = v;
    *psb = sb;
}

/*
 * 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 void
_separated_turing_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 i;
    slong sb_near; /* Turing method bound for near search */
    slong sb_far; /* Turing method bound for far search */
    slong zn; /* target number of sign changes */
    slong variations; /* observed number of sign changes */
    const slong loopcount = 4;
    /*
     * The loopcount controls how hard we try to find zeros in Gram/Rosser
     * blocks. If the loopcount is too high then when Rosser's rule is violated
     * we will spend too much time searching for zeros that don't exist.
     * If the loopcount is too low then we will miss some 'good' blocks,
     * causing the search to be extended unnecessarily far from the initial
     * guess before eventually making another pass in which loopcount
     * is effectively infinite.
     */

    if (fmpz_cmp_si(n, 2) < 0)
    {
        flint_printf("invalid n: "); fmpz_print(n); flint_printf("\n");
        flint_abort();
    }

    turing_search_near(&u, &v, &sb_near, n);
    trim(&U, &V, u, v, sb_near);
    zn = count_gram_intervals(U, V);
    for (i = 0; i < loopcount && count_sign_changes(U, V) < zn; i++)
    {
        intercalate(U, V);
    }
    variations = count_sign_changes(U, V);
    if (variations > zn)
    {
        flint_printf("unexpected number of sign changes\n");
        flint_abort();
    }
    else if (variations < zn)
    {
        zz_node_ptr r = U;
        zz_node_ptr s = V;
        turing_search_far(&u, &v, &sb_far, u, v, sb_near);
        trim(&U, &V, u, v, 2*sb_far);
        zn = count_gram_intervals(U, V);
        for (i = 0; i < loopcount && count_sign_changes(U, V) < zn; i++)
        {
            intercalate(U, r);
            intercalate(s, V);
        }
        variations = count_sign_changes(U, V);
        if (variations > zn)
        {
            flint_printf("unexpected number of sign changes\n");
            flint_abort();
        }
        else if (variations < zn)
        {
            trim(&U, &V, u, v, sb_far);
            zn = count_gram_intervals(U, V);
            while (count_sign_changes(U, V) < zn)
            {
                intercalate(U, V);
            }
            if (count_sign_changes(U, V) != zn)
            {
                flint_printf("unexpected number of sign changes\n");
                flint_abort();
            }
        }
    }

    *pu = u;
    *pv = v;
    *pU = U;
    *pV = V;
}

/*
 * 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
 * with fully separated zeros.
 */
static void
_separated_rosser_list(zz_node_ptr *pu, zz_node_ptr *pv, const fmpz_t n)
{
    fmpz_t k;
    zz_node_ptr u, v;

    if (fmpz_cmp_si(n, 1) < 0 || fmpz_cmp_si(n, ROSSERS_RULE_MAX) > 0)
    {
        flint_printf("invalid n: "); fmpz_print(n); flint_printf("\n");
        flint_abort();
    }

    fmpz_init(k);

    fmpz_sub_ui(k, n, 2);
    u = create_gram_node(k);
    fmpz_sub_ui(k, n, 1);
    v = create_gram_node(k);
    u->next = v;
    v->prev = u;

    if (!zz_node_is_good_gram_node(u))
        u = extend_to_prev_good_gram_node(u);
    if (!zz_node_is_good_gram_node(v))
        v = extend_to_next_good_gram_node(v);
    while (count_sign_changes(u, v) != count_gram_intervals(u, v))
    {
        intercalate(u, v);
    }

    *pu = u;
    *pv = v;

    fmpz_clear(k);
}

/*
 * 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
 * with fully separated zeros.
 */
static void
_separated_gram_list(zz_node_ptr *pu, zz_node_ptr *pv, const fmpz_t n)
{
    fmpz_t k;
    zz_node_ptr u, v;

    if (fmpz_cmp_si(n, 1) < 0 || fmpz_cmp_si(n, GRAMS_LAW_MAX) > 0)
    {
        flint_printf("invalid n: "); fmpz_print(n); flint_printf("\n");
        flint_abort();
    }

    fmpz_init(k);

    fmpz_sub_ui(k, n, 2);
    u = create_gram_node(k);
    fmpz_sub_ui(k, n, 1);
    v = create_gram_node(k);
    u->next = v;
    v->prev = u;

    *pu = u;
    *pv = v;

    fmpz_clear(k);
}

/*
 * Get a list of points that fully separate zeros of Z.
 * 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 void
_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;
    if (fmpz_cmp_si(n, GRAMS_LAW_MAX) <= 0)
    {
        _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();
    }

    *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)
    {
        v = u;
        u = u->next;
        zz_node_clear(v);
        flint_free(v);
    }

    return count;
}

/* Isolate len zeros, starting from the nth zero. */
void
acb_dirichlet_isolate_hardy_z_zeros(arf_interval_ptr res, const fmpz_t n, slong len)
{
    if (len <= 0)
    {
        return;
    }
    else if (fmpz_sgn(n) < 1)
    {
        flint_printf("nonpositive indices of zeros are not supported\n");
        flint_abort();
    }
    else
    {
        slong c = 0;
        fmpz_t k;
        fmpz_init(k);
        while (c < len)
        {
            fmpz_add_si(k, n, c);
            c += _isolate_hardy_z_zeros(res + c, k, len - c);
        }
        fmpz_clear(k);
    }
}

void
acb_dirichlet_isolate_hardy_z_zero(arf_t a, arf_t b, const fmpz_t n)
{
    if (fmpz_sgn(n) < 1)
    {
        flint_printf("nonpositive indices of zeros are not supported\n");
        flint_abort();
    }
    else
    {
        arf_interval_t r;
        arf_interval_init(r);
        _isolate_hardy_z_zeros(r, n, 1);
        arf_set(a, &r->a);
        arf_set(b, &r->b);
        arf_interval_clear(r);
    }
}

static void
count_up(arf_t a, arf_t b, zz_node_srcptr U, zz_node_srcptr V, const fmpz_t n)
{
    arf_interval_t r;
    arf_interval_init(r);
    count_up_separated_zeros(r, U, V, n, 1);
    arf_set(a, &r->a);
    arf_set(b, &r->b);
    arf_interval_clear(r);
}

void
_acb_dirichlet_isolate_turing_hardy_z_zero(arf_t a, arf_t b, const fmpz_t n)
{
    zz_node_ptr U, V, u, v;
    _separated_turing_list(&U, &V, &u, &v, n);
    count_up(a, b, U, V, n);
    while (u)
    {
        v = u;
        u = u->next;
        zz_node_clear(v);
        flint_free(v);
    }
}

void
_acb_dirichlet_isolate_rosser_hardy_z_zero(arf_t a, arf_t b, const fmpz_t n)
{
    zz_node_ptr u, v;
    _separated_rosser_list(&u, &v, n);
    count_up(a, b, u, v, n);
    while (u)
    {
        v = u;
        u = u->next;
        zz_node_clear(v);
        flint_free(v);
    }
}

void
_acb_dirichlet_isolate_gram_hardy_z_zero(arf_t a, arf_t b, const fmpz_t n)
{
    zz_node_ptr u, v;
    _separated_gram_list(&u, &v, n);
    count_up(a, b, u, v, n);
    while (u)
    {
        v = u;
        u = u->next;
        zz_node_clear(v);
        flint_free(v);
    }
}

static void
_acb_set_arf(acb_t res, const arf_t t)
{
    acb_zero(res);
    arb_set_arf(acb_realref(res), t);
}

/*
 * Find the index of the largest Gram point less than t.
 * Requires t > 10.
 */
static void
gram_index(fmpz_t res, const arf_t t)
{
    if (arf_cmp_si(t, 10) < 0)
    {
        flint_printf("gram_index requires t > 10\n");
        flint_abort();
    }
    else
    {
        acb_t z;
        slong prec = arf_abs_bound_lt_2exp_si(t) + 8;
        acb_init(z);
        while (1)
        {
            _acb_set_arf(z, t);
            acb_dirichlet_hardy_theta(z, z, NULL, NULL, 1, prec);
            arb_const_pi(acb_imagref(z), prec);
            arb_div(acb_realref(z), acb_realref(z), acb_imagref(z), prec);
            arb_floor(acb_realref(z), acb_realref(z), prec);
            if (arb_get_unique_fmpz(res, acb_realref(z)))
            {
                break;
            }
            prec *= 2;
        }
        acb_clear(z);
    }
}

/*
 * Compute nzeros(t) for up to len values of t
 * and return the number computed. p must be in increasing order,
 * and all entries must be greater than 10.
 */
static slong
_exact_zeta_multi_nzeros(fmpz *res, arf_srcptr points, slong len)
{
    zz_node_ptr U, V, u, v, p;
    arb_t x;
    fmpz_t n, N;
    slong i;
    arf_srcptr t;
    int s;
    slong prec;

    arb_init(x);
    fmpz_init(n);
    fmpz_init(N);

    /*
     * The first point is located between two Gram points. Find the unique n
     * such that the nth zero is also predicted to be located between these
     * Gram points.
     */
    gram_index(n, points);
    fmpz_add_ui(n, n, 2);

    /* Get a list of points that fully separate zeros of Z. */
    _separated_list(&U, &V, &u, &v, n);

    p = U;
    fmpz_add_ui(N, U->gram, 1);

    /*
     * It's possible that one or more points are located between
     * the first Gram point and the arf_t representative of that Gram point.
     */
    for (i = 0, t = points; i < len && arf_cmp(t, &p->t) <= 0; i++, t++)
    {
        fmpz_set(res + i, N);
    }

    while (i < len && p != V)
    {
        if (!p->next)
        {
            flint_printf("prematurely reached the end of the list\n");
            flint_abort();
        }
        if (zz_node_sgn(p) != zz_node_sgn(p->next))
        {
            while (i < len && arf_cmp(t, &p->next->t) <= 0)
            {
                prec = arf_abs_bound_lt_2exp_si(t) + 8;
                s = _acb_dirichlet_definite_hardy_z(x, t, &prec);
                if (zz_node_sgn(p->next) == s)
                {
                    fmpz_add_ui(res + i, N, 1);
                }
                else
                {
                    fmpz_set(res + i, N);
                }
                t++;
                i++;
            }
            fmpz_add_ui(N, N, 1);
        }
        p = p->next;
    }

    /*
     * It's possible that the first point in the array is located between
     * the last Gram point and the arf_t representative of that Gram point.
     */
    if (i == 0)
    {
        fmpz_set(res, N);
        i++;
    }

    while (u)
    {
        v = u;
        u = u->next;
        zz_node_clear(v);
        flint_free(v);
    }

    arb_clear(x);
    fmpz_clear(n);
    fmpz_clear(N);

    return i;
}

/*
 * Compute nzeros for len values of t. The array p must be in increasing order.
 */
static void
exact_zeta_multi_nzeros(fmpz *res, arf_srcptr p, slong len)
{
    slong i, c;
    arf_srcptr q;

    if (len <= 0)
    {
        return;
    }

    for (i = 0, q = p; i < len - 1; i++, q++)
    {
        if (arf_cmp(q, p) > 0)
        {
            flint_printf("p must be in increasing order\n");
            flint_abort();
        }
    }

    c = 0;
    while (c < len)
    {
        if (arf_cmp_si(p + c, 14) < 0)
        {
            fmpz_zero(res + c);
            c++;
        }
        else
        {
            c += _exact_zeta_multi_nzeros(res + c, p + c, len - c);
        }
    }
}

void
_acb_dirichlet_exact_zeta_nzeros(fmpz_t res, const arf_t t)
{
    exact_zeta_multi_nzeros(res, t, 1);
}

void
acb_dirichlet_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, slong prec)
{
    if (len <= 0)
    {
        return;
    }
    else if (fmpz_sgn(n) < 1)
    {
        flint_printf("nonpositive indices of zeros are not supported\n");
        flint_abort();
    }
    else
    {
        slong i;
        arf_interval_ptr p = _arf_interval_vec_init(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);
        }
        _arf_interval_vec_clear(p, len);
    }
}

static void
_arb_set_interval_fmpz(arb_t res, const fmpz_t a, const fmpz_t b)
{
    fmpz_t n;

    fmpz_init(n);

    fmpz_add(n, a, b);
    arf_set_fmpz(arb_midref(res), n);
    fmpz_sub(n, b, a);
    mag_set_fmpz(arb_radref(res), n);
    arb_mul_2exp_si(res, res, -1);

    fmpz_clear(n);
}

void
acb_dirichlet_zeta_nzeros(arb_t res, const arb_t t, slong prec)
{
    if (arb_is_exact(t))
    {
        fmpz_t n;
        fmpz_init(n);
        _acb_dirichlet_exact_zeta_nzeros(n, arb_midref(t));
        arb_set_fmpz(res, n);
        fmpz_clear(n);
    }
    else
    {
        arf_struct b[2];
        fmpz n[2];
        
        arf_init(b);
        arf_init(b + 1);
        fmpz_init(n);
        fmpz_init(n + 1);

        arb_get_lbound_arf(b, t, prec);
        arb_get_ubound_arf(b + 1, t, prec);
        exact_zeta_multi_nzeros(n, b, 2);
        _arb_set_interval_fmpz(res, n, n + 1);

        arf_clear(b);
        arf_clear(b + 1);
        fmpz_clear(n);
        fmpz_clear(n + 1);
    }

    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);
}