/* Copyright (C) 2010 Juan Arias de Reyna Copyright (C) 2019 D.H.J. Polymath Copyright (C) 2019 Fredrik Johansson 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 . */ #include "acb_dirichlet.h" #include "arb_calc.h" static const slong LOOPCOUNT = 4; /* * For a detailed explanation of the variant of Turing's method * 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 */ /* * This structure describes the local context in which Platt's scaled Lambda * function f may be estimated by Gaussian-windowed Whittaker-Shannon * interpolation. */ typedef struct { /* grid location and shape */ fmpz T; /* midpoint */ slong A; /* resolution */ slong B; /* width */ /* interpolation tuning parameters */ slong Ns_max; /* max number of support points per side */ arb_struct H; /* standard deviation of Gaussian window */ slong sigma; arb_ptr p; /* f evaluated at N = A*B points on the grid */ acb_dirichlet_platt_ws_precomp_struct pre; /* precomp interpolation stuff */ } platt_ctx_struct; typedef platt_ctx_struct platt_ctx_t[1]; typedef platt_ctx_struct * platt_ctx_ptr; typedef const platt_ctx_struct * platt_ctx_srcptr; static void platt_ctx_init(platt_ctx_t ctx, 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) { fmpz_init(&ctx->T); arb_init(&ctx->H); ctx->p = _arb_vec_init(A*B); ctx->A = A; ctx->B = B; ctx->Ns_max = Ns_max; ctx->sigma = sigma_interp; fmpz_set(&ctx->T, T); arb_set(&ctx->H, H); acb_dirichlet_platt_ws_precomp_init(&ctx->pre, A, H, sigma_interp, prec); acb_dirichlet_platt_multieval(ctx->p, T, A, B, h, J, K, sigma_grid, prec); } static void platt_ctx_clear(platt_ctx_t ctx) { slong N = ctx->A * ctx->B; fmpz_clear(&ctx->T); arb_clear(&ctx->H); _arb_vec_clear(ctx->p, N); acb_dirichlet_platt_ws_precomp_clear(&ctx->pre); } static void platt_ctx_interpolate(arb_t res, arf_t deriv, const platt_ctx_t ctx, const arb_t t0, slong prec) { acb_dirichlet_platt_ws_interpolation_precomp(res, deriv, &ctx->pre, t0, ctx->p, &ctx->T, ctx->A, ctx->B, ctx->Ns_max, &ctx->H, ctx->sigma, prec); } static void platt_ctx_interpolate_arf(arb_t res, arf_t deriv, const platt_ctx_t ctx, const arf_t t0, slong prec) { arb_t t; arb_init(t); arb_set_arf(t, t0); platt_ctx_interpolate(res, deriv, ctx, t, prec); arb_clear(t); } /* * This structure describes a node of a doubly linked list. * Each node represents a height t at which Platt's scaled Lambda function * f(t) has been evaluated by interpolation on a grid. */ typedef struct _zz_node_struct { arf_struct t; /* height t where v = f(t) is evaluated */ arb_struct v; /* f(t) */ fmpz *gram; /* Gram point index or NULL if not a Gram point */ 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 f(t)\n"); flint_abort(); } return s; } /* Good Gram points are Gram points where sgn(f(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->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->gram = NULL; p->prev = NULL; p->next = NULL; } static void delete_list_to(zz_node_ptr head, zz_node_srcptr target) { zz_node_ptr u, v; if (head) { if (head->prev) { flint_printf("expected the first node in the list\n"); flint_abort(); } } u = head; while (u != target) { if (u == NULL) { flint_printf("failed to find target within list\n"); flint_abort(); } v = u; u = u->next; zz_node_clear(v); flint_free(v); } if (u != NULL) { u->prev = NULL; } } static void delete_list(zz_node_ptr head) { delete_list_to(head, NULL); } /* * Create a node representing an evaluation of the scaled Lambda function * at arbitrary point t. Upon creation the sign of f(t) will * be known with certainty. If this is not possible then NULL is returned. */ static zz_node_ptr create_non_gram_node(const arf_t t, const platt_ctx_t ctx, slong prec) { zz_node_ptr p = flint_malloc(sizeof(zz_node_struct)); zz_node_init(p); arf_set(&p->t, t); platt_ctx_interpolate_arf(&p->v, NULL, ctx, t, prec); if (arb_contains_zero(&p->v)) { zz_node_clear(p); p = NULL; } return p; } /* * Create a node representing an evaluation of the scaled Lambda 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 * scaled Lambda function between t and the actual value of g(n). * The sign of f(t) will also be known with certainty, otherwise NULL * is returned. */ static zz_node_ptr create_gram_node(const fmpz_t n, const platt_ctx_t ctx, slong prec) { zz_node_ptr p = NULL; arb_t t, v; acb_t z; arb_init(t); arb_init(v); acb_init(z); acb_dirichlet_gram_point(t, n, NULL, NULL, prec + fmpz_sizeinbase(n, 2)); acb_set_arb(z, t); platt_ctx_interpolate(v, NULL, ctx, t, prec); if (!arb_contains_zero(v)) { /* t contains g(n) and does not contain a zero of the f function */ p = flint_malloc(sizeof(zz_node_struct)); zz_node_init(p); p->gram = flint_malloc(sizeof(fmpz)); fmpz_init(p->gram); fmpz_set(p->gram, n); arf_set(&p->t, arb_midref(t)); arb_set(&p->v, v); } 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 f(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. * Returns nonzero on success. * Sets *out to the last point in the list. */ static int extend_to_next_good_gram_node( zz_node_ptr *out, zz_node_t p, const platt_ctx_t ctx, slong prec) { fmpz_t n; zz_node_ptr q, r; int result = 1; 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, ctx, prec); if (r) { q->next = r; r->prev = q; q = r; r = NULL; if (zz_node_is_good_gram_node(q)) { break; } } else { result = 0; break; } } fmpz_clear(n); *out = q; return result; } /* * Modify a linked list that begins with node p, * by prepending nodes representing Gram points. * Continue until a 'good' Gram point is found. * Returns nonzero on success. * Sets *out to the first point in the list. */ static int extend_to_prev_good_gram_node(zz_node_ptr *out, zz_node_t p, const platt_ctx_t ctx, slong prec) { fmpz_t n; zz_node_ptr q, r; int result = 1; 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, ctx, prec); if (r) { q->prev = r; r->next = q; q = r; r = NULL; if (zz_node_is_good_gram_node(q)) { break; } } else { result = 0; break; } } fmpz_clear(n); *out = q; return result; } static zz_node_ptr _scan_to_prev_good_gram_node(zz_node_ptr p) { zz_node_ptr u = p->prev; while (u) { if (zz_node_is_good_gram_node(u)) { return u; } u = u->prev; } return NULL; } static zz_node_ptr scan_to_prev_good_gram_node(zz_node_ptr p, slong count) { slong i; zz_node_ptr u = p; for (i = 0; i < count; i++) { if ((u = _scan_to_prev_good_gram_node(u)) == NULL) { return NULL; } } return u; } static zz_node_ptr _scan_to_next_good_gram_node(zz_node_ptr p) { zz_node_ptr u = p->next; while (u) { if (zz_node_is_good_gram_node(u)) { return u; } u = u->next; } return NULL; } static zz_node_ptr scan_to_next_good_gram_node(zz_node_ptr p, slong count) { slong i; zz_node_ptr u = p; for (i = 0; i < count; i++) { if ((u = _scan_to_next_good_gram_node(u)) == NULL) { return NULL; } } return u; } /* * 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(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(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 f(t) * between endpoints. * v1 and v2 are the scaled Lambda 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. * Returns nonzero on success. */ static int intercalate(const platt_ctx_t ctx, zz_node_t a, zz_node_t b, slong prec) { arb_t t; zz_node_ptr q, r, mid_node; int result = 1; 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 result; arb_init(t); q = a; r = a->next; while (q != b) { if (!r) { flint_printf("prematurely reached end of list\n"); flint_abort(); } split_interval(t, &q->t, &q->v, zz_node_sgn(q), &r->t, &r->v, zz_node_sgn(r), prec); if (arb_contains_arf(t, &q->t) || arb_contains_arf(t, &r->t)) { result = 0; break; } mid_node = create_non_gram_node(arb_midref(t), ctx, prec); if (!mid_node) { result = 0; break; } q->next = mid_node; mid_node->prev = q; mid_node->next = r; r->prev = mid_node; q = r; r = r->next; } arb_clear(t); return result; } /* * Given a linked sublist beginning at U and ending at V defining function * evaluations at points that fully separate zeros of f(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; } /* * Create a small linked list defining the Gram block that is expected * to contain the nth zero according to the Gram heuristic. * Returns 0 if unable to create the Gram block. * The output node *p is the first node in the Gram block on success * or NULL on failure. * The output node *q is the last node in the Gram block on success * or NULL on failure. * Success does not necessarily mean that the Gram block contains the nth zero. */ static int create_initial_gram_block(zz_node_ptr *p, zz_node_ptr *q, const platt_ctx_t ctx, const fmpz_t n, slong prec) { zz_node_ptr u, v; fmpz_t k; slong result = 1; fmpz_init(k); *p = NULL; *q = NULL; fmpz_sub_ui(k, n, 2); u = create_gram_node(k, ctx, prec); if (!u) { result = 0; goto finish; } fmpz_sub_ui(k, n, 1); v = create_gram_node(k, ctx, prec); if (!v) { result = 0; goto finish; } u->next = v; v->prev = u; if (!zz_node_is_good_gram_node(u)) { if (!extend_to_prev_good_gram_node(&u, u, ctx, prec)) { result = 0; goto finish; } } if (!zz_node_is_good_gram_node(v)) { if (!extend_to_next_good_gram_node(&v, v, ctx, prec)) { result = 0; goto finish; } } finish: if (result) { *p = u; *q = v; } else { delete_list(u); } return result; } /* * On failure returns 0 and output variables are NULL or zero. * The output variable *pu is the head of the output list. * The output variable *pv is the tail of the output list. * The first *pbound Gram blocks in the output list are certified as 'good', * and the list probably contains a few more trailing Gram blocks (not * necessarily certified as 'good'). */ static int create_initial_double_superblock(zz_node_ptr *pu, zz_node_ptr *pv, slong *pbound, const platt_ctx_t ctx, const fmpz_t n, slong prec) { zz_node_ptr p, q, u, v; slong i, k, bound, zn; slong good_block_count; slong result = 1; *pu = NULL; *pv = NULL; *pbound = 0; if (!create_initial_gram_block(&p, &q, ctx, n, prec)) { result = 0; goto finish; } /* * Add blocks in the forward direction until we have at least k Gram blocks, * where k is the turing method bound. * Note that these blocks are not necessarily good Gram blocks * and will not necessarily belong to the initial double superblock. * The bound may increase as the list is extended. */ for (k = 1; k < acb_dirichlet_turing_method_bound(q->gram); k++) { if (!extend_to_next_good_gram_node(&q, q, ctx, prec)) { result = 0; goto finish; } } bound = k; /* * Scan the list backwards, attempting to certify blocks as 'good' * and tracking the current number of consecutive good Gram blocks. */ good_block_count = 0; v = q; while ((u = _scan_to_prev_good_gram_node(v)) != NULL) { zn = count_gram_intervals(u, v); for (i = 0; i < LOOPCOUNT && count_sign_changes(u, v) < zn; i++) { if (!intercalate(ctx, u, v, prec)) { result = 0; goto finish; } } if (count_sign_changes(u, v) >= zn) { good_block_count++; } else { good_block_count = 0; } v = u; } if (v != p) { flint_printf("unexpected endpoint of backwards scan\n"); flint_abort(); } /* * Add blocks in the backwards direction until the number * of consecutive good Gram blocks is twice the computed bound. */ while (good_block_count < 2*bound) { if (!extend_to_prev_good_gram_node(&p, v, ctx, prec)) { result = 0; goto finish; } zn = count_gram_intervals(p, v); for (i = 0; i < LOOPCOUNT && count_sign_changes(p, v) < zn; i++) { if (!intercalate(ctx, p, v, prec)) { result = 0; goto finish; } } if (count_sign_changes(p, v) >= zn) { good_block_count++; } else { good_block_count = 0; } v = p; } finish: if (result) { *pu = p; *pv = q; *pbound = bound; } else { delete_list(p); } return result; } static slong _isolate_zeros(arf_interval_ptr res, const platt_ctx_t ctx, const fmpz_t n, slong len, slong prec) { zz_node_ptr x, y; /* Anchor nodes where N(t) is known */ zz_node_ptr u, v; /* For certifying Gram nodes as 'good' */ zz_node_ptr p, q; fmpz_t nnext; slong i, k, bound, zn, zc, zeros_count; fmpz_init(nnext); fmpz_set(nnext, n); p = NULL; zeros_count = 0; if (!create_initial_double_superblock(&p, &q, &bound, ctx, n, prec)) { goto finish; } /* * Set the anchor to the central good Gram node of the * initial double superblock. * Delete the nodes in the list before the anchor. */ x = scan_to_next_good_gram_node(p, bound); if (x == NULL) { flint_printf("missing or incomplete initial block\n"); flint_abort(); } delete_list_to(p, x); p = x; /* * Set v to the forward-most node in the list, * and track the number of consecutive good Gram blocks at that point. */ v = scan_to_next_good_gram_node(p, bound); if (v == NULL) { flint_printf("missing or incomplete initial block\n"); flint_abort(); } k = 2*bound; u = v; while ((v = _scan_to_next_good_gram_node(v)) != NULL) { zn = count_gram_intervals(u, v); if (count_sign_changes(u, v) >= zn) { k++; } else { k = 0; } u = v; } if (u != q) { flint_printf("failed to scan the initial list\n"); flint_abort(); } v = u; /* * Iterate through Gram blocks. The central good Gram point in each * stretch of 2*bound consecutive good Gram blocks is an 'anchor' point * where the number of zeros less than that point is known. * Therefore the number of zeros between each pair of anchor points * is known. As anchor points are certified, isolate the zeros falling * between each pair. */ while (1) { u = v; if (!extend_to_next_good_gram_node(&v, v, ctx, prec)) { goto finish; } zn = count_gram_intervals(u, v); for (i = 0; i < LOOPCOUNT && count_sign_changes(u, v) < zn; i++) { if (!intercalate(ctx, u, v, prec)) { goto finish; } } if (count_sign_changes(u, v) >= zn) { k++; } else { k = 0; } bound = acb_dirichlet_turing_method_bound(v->gram); if (k >= 2*bound && fmpz_cmp(x->gram, v->gram) < 0) { /* There are exactly zn zeros between the anchor points x and y. */ y = scan_to_prev_good_gram_node(v, bound); if (!y) { flint_printf("failed to scan backwards to anchor point\n"); flint_abort(); } zn = count_gram_intervals(x, y); while (count_sign_changes(x, y) < zn) { if (!intercalate(ctx, x, y, prec)) { goto finish; } } zc = count_up_separated_zeros(res + zeros_count, x, y, nnext, len - zeros_count); if (zc < 0 || zc > len - zeros_count) { flint_printf("unexpected number of isolated zeros\n"); flint_abort(); } zeros_count += zc; if (zeros_count == len) { goto finish; } fmpz_add_ui(nnext, nnext, zc); x = y; delete_list_to(p, x); p = x; } } finish: fmpz_clear(nnext); delete_list(p); return zeros_count; } slong _acb_dirichlet_platt_isolate_local_hardy_z_zeros( arf_interval_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) { slong zeros_count; platt_ctx_t ctx; platt_ctx_init(ctx, T, A, B, h, J, K, sigma_grid, Ns_max, H, sigma_interp, prec); zeros_count = _isolate_zeros(res, ctx, n, len, prec); platt_ctx_clear(ctx); return zeros_count; } static void _refine_local_hardy_z_zero_illinois(arb_t res, const platt_ctx_t ctx, const arf_t ra, const arf_t rb, slong prec) { arf_t a, b, fa, fb, c, fc, t; arb_t z; slong k, nmag, abs_tol, wp; int asign, bsign, csign; arf_init(a); arf_init(b); arf_init(c); arf_init(fa); arf_init(fb); arf_init(fc); arf_init(t); arb_init(z); arf_set(a, ra); arf_set(b, rb); nmag = arf_abs_bound_lt_2exp_si(b); abs_tol = nmag - prec - 4; wp = prec + nmag + 8; platt_ctx_interpolate_arf(z, NULL, ctx, a, wp); asign = arb_sgn_nonzero(z); arf_set(fa, arb_midref(z)); platt_ctx_interpolate_arf(z, NULL, ctx, b, wp); bsign = arb_sgn_nonzero(z); arf_set(fb, arb_midref(z)); if (!asign || !bsign) { flint_printf("the function evaluations at the endpoints of the initial " "interval must not contain zero\n"); flint_abort(); } if (asign == bsign) { flint_printf("isolate a zero before bisecting the interval\n"); flint_abort(); } for (k = 0; k < 40; k++) { /* c = a - fa * (b - a) / (fb - fa) */ arf_sub(c, b, a, wp, ARF_RND_NEAR); arf_sub(t, fb, fa, wp, ARF_RND_NEAR); arf_div(c, c, t, wp, ARF_RND_NEAR); arf_mul(c, c, fa, wp, ARF_RND_NEAR); arf_sub(c, a, c, wp, ARF_RND_NEAR); /* if c is not sandwiched between a and b, fall back to one bisection step */ if (!arf_is_finite(c) || !((arf_cmp(a, c) < 0 && arf_cmp(c, b) < 0) || (arf_cmp(b, c) < 0 && arf_cmp(c, a) < 0))) { /* flint_printf("no sandwich (k = %wd)\n", k); */ arf_add(c, a, b, ARF_PREC_EXACT, ARF_RND_DOWN); arf_mul_2exp_si(c, c, -1); } platt_ctx_interpolate_arf(z, NULL, ctx, c, wp); csign = arb_sgn_nonzero(z); /* If the guess is close enough to a zero that the sign * cannot be determined, then use the derivative to * make an appropriately small interval around the guess. */ if (!csign) { arf_t deriv, aprime, bprime, faprime, fbprime, err, delta; slong i, aprimesign, bprimesign; arf_init(deriv); arf_init(aprime); arf_init(bprime); arf_init(faprime); arf_init(fbprime); arf_init(err); arf_init(delta); arf_set_mag(err, arb_radref(z)); platt_ctx_interpolate_arf(NULL, deriv, ctx, c, wp); arf_div(delta, err, deriv, wp, ARF_RND_NEAR); arf_mul_si(delta, delta, 3, wp, ARF_RND_NEAR); arf_mul_2exp_si(delta, delta, -1); arf_set(aprime, c); arf_set(bprime, c); /* When the context allows the interval endpoints to * be evaluated to relatively high precision, * this should not require more than one or two iterations. */ for (i = 0; i < 5; i++) { arf_sub(aprime, aprime, delta, wp, ARF_RND_DOWN); arf_add(bprime, bprime, delta, wp, ARF_RND_UP); if (arf_cmp(a, b) < 0) { if (arf_cmp(aprime, a) < 0) arf_set(aprime, a); if (arf_cmp(b, bprime) < 0) arf_set(bprime, b); } else { if (arf_cmp(aprime, b) < 0) arf_set(aprime, b); if (arf_cmp(a, bprime) < 0) arf_set(bprime, a); } platt_ctx_interpolate_arf(z, NULL, ctx, aprime, wp); arf_set(faprime, arb_midref(z)); aprimesign = arb_sgn_nonzero(z); platt_ctx_interpolate_arf(z, NULL, ctx, bprime, wp); arf_set(fbprime, arb_midref(z)); bprimesign = arb_sgn_nonzero(z); if (aprimesign && bprimesign && aprimesign != bprimesign) { arf_set(a, aprime); arf_set(b, bprime); arf_set(fa, faprime); arf_set(fb, fbprime); break; } } arf_clear(deriv); arf_clear(aprime); arf_clear(bprime); arf_clear(faprime); arf_clear(fbprime); arf_clear(err); arf_clear(delta); break; } arf_set(fc, arb_midref(z)); if (csign != bsign) { arf_set(a, b); arf_set(fa, fb); asign = bsign; arf_set(b, c); arf_set(fb, fc); bsign = csign; } else { arf_set(b, c); arf_set(fb, fc); bsign = csign; arf_mul_2exp_si(fa, fa, -1); } arf_sub(t, a, b, wp, ARF_RND_DOWN); arf_abs(t, t); if (arf_cmpabs_2exp_si(t, abs_tol) < 0) break; } /* a and b may have changed places */ if (arf_cmp(a, b) > 0) arf_swap(a, b); arb_set_interval_arf(res, a, b, prec); arf_clear(a); arf_clear(b); arf_clear(c); arf_clear(fa); arf_clear(fb); arf_clear(fc); arf_clear(t); arb_clear(z); } 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) { slong zeros_count, i; arf_interval_ptr p; platt_ctx_t ctx; platt_ctx_init( ctx, T, A, B, h, J, K, sigma_grid, Ns_max, H, sigma_interp, prec); p = _arf_interval_vec_init(len); zeros_count = _isolate_zeros(p, ctx, n, len, prec); for (i = 0; i < zeros_count; i++) { _refine_local_hardy_z_zero_illinois(res+i, ctx, &p[i].a, &p[i].b, prec); } platt_ctx_clear(ctx); _arf_interval_vec_clear(p, len); return zeros_count; } static void _arb_get_lbound_fmpz(fmpz_t z, const arb_t x, slong prec) { arf_t u; arf_init(u); arb_get_lbound_arf(u, x, prec); arf_get_fmpz(z, u, ARF_RND_DOWN); arf_clear(u); } /* Compares f to g=a*10^b. * Returns a negative value if f < g, positive value if g < f, otherwise 0. */ static int _fmpz_cmp_a_10exp_b(const fmpz_t f, slong a, slong b) { int result; fmpz_t g; fmpz_init(g); fmpz_set_ui(g, 10); fmpz_pow_ui(g, g, b); fmpz_mul_si(g, g, a); result = fmpz_cmp(f, g); fmpz_clear(g); return result; } static platt_ctx_ptr _create_heuristic_context(const fmpz_t n, slong prec) { platt_ctx_ptr p = NULL; slong J, K, A, B, Ns_max, sigma_grid, sigma_interp; slong kbits; fmpz_t T, k; arb_t g, h, H, logT; double dlogJ, dK, dgrid, dh, dH, dinterp; double x, x2, x3, x4; fmpz_init(T); fmpz_init(k); arb_init(g); arb_init(h); arb_init(H); arb_init(logT); /* Estimate the height of the nth zero using gram points -- * it's predicted to fall between g(n-2) and g(n-1). */ fmpz_sub_ui(k, n, 2); kbits = fmpz_sizeinbase(k, 2); acb_dirichlet_gram_point(g, k, NULL, NULL, prec + kbits); /* Let T be the integer at the center of the evaluation grid. */ _arb_get_lbound_fmpz(T, g, prec + kbits); arb_log_fmpz(logT, T, prec); x = arf_get_d(arb_midref(logT), ARF_RND_NEAR); x2 = x*x; x3 = x2*x; x4 = x2*x2; if (_fmpz_cmp_a_10exp_b(n, 1, 4) < 0) { goto finish; } else if (_fmpz_cmp_a_10exp_b(n, 1, 5) < 0) { /* interpolated for n in [1e4, 1e5] */ A = 4; B = 64; Ns_max = 100; dinterp = 25; dK = 28; dgrid = 31; dlogJ = 8.4398 + -0.40306*x + 0.029866*x2 + -2.2858e-05*x3; dh = 1.0844 + 0.25524*x + -0.0046997*x2 + -6.3447e-05*x3; dH = -11.882 + 3.9521*x + -0.38654*x2 + 0.012728*x3; } else if (_fmpz_cmp_a_10exp_b(n, 1, 7) < 0) { /* interpolated for n in [1e4, 1e7] */ A = 8; B = 4096; Ns_max = 200; dinterp = 25; dlogJ = 0.88323 + 0.21392*x + 0.020846*x2 + -0.00053151*x3; dK = 137.27 + -15.609*x + 1.0778*x2 + -0.025927*x3; dgrid = -1711.1 + 701.03*x + -48.424*x2 + 1.2075*x3; dh = 448.2 + -84.087*x + 6.2089*x2 + -0.14565*x3; dH = 0.94123 + 0.021136*x + -0.00093042*x2 + 3.1007e-05*x3; } else if (_fmpz_cmp_a_10exp_b(n, 2, 17) < 0) { /* interpolated for n in [1e7, 5e22] */ A = 8; B = 4096; Ns_max = 200; dlogJ = -0.4035 + 0.49086*x + 0.00016299*x2 + -3.6139e-06*x3 + 2.9323e-08*x4; dK = 79.032 + -1.781*x + 0.039243*x2 + -0.00094859*x3 + 7.3149e-06*x4; dgrid = 1186.9 + 130.17*x + -7.4059*x2 + 0.17895*x3 + -0.001602*x4; dinterp = -24.252 + 7.3231*x + -0.38971*x2 + 0.0088745*x3 + -7.4331e-05*x4; dh = 178.66 + -15.127*x + 0.93132*x2 + -0.02311*x3 + 0.00022146*x4; dH = 2.5499 + -0.24402*x + 0.014953*x2 + -0.00037347*x3 + 3.5596e-06*x4; } else if (_fmpz_cmp_a_10exp_b(n, 1, 37) < 0) { /* interpolated for n in [1e7, 1e37] */ A = 16; B = 8192; Ns_max = 300; dlogJ = -0.50566 + 0.49723*x + 1.7964e-05*x2 + -2.3664e-07*x3 + 1.1234e-09*x4; dK = 100.97 + -0.709*x + -0.0020664*x2 + 3.1633e-05*x3 + -2.2912e-07*x4; dgrid = 3998.1 + 6.68*x + -0.3202*x2 + 0.0051782*x3 + -3.3829e-05*x4; dinterp = 21.203 + -0.2797*x + 0.01191*x2 + -0.00019769*x3 + 1.0395e-06*x4; dh = 137.6 + -0.16471*x + 0.039086*x2 + -0.00063299*x3 + 4.9674e-06*x4; dH = 0.64172 + -0.0017413*x + 0.0002195*x2 + -3.5247e-06*x3 + 2.6633e-08*x4; } else { goto finish; } arb_set_d(h, dh); arb_set_d(H, dH); J = (slong) exp(dlogJ); K = (slong) dK; sigma_grid = ((slong) (dgrid/2))*2 + 1; sigma_interp = ((slong) (dinterp/2))*2 + 1; p = malloc(sizeof(platt_ctx_struct)); platt_ctx_init(p, T, A, B, h, J, K, sigma_grid, Ns_max, H, sigma_interp, prec); finish: fmpz_clear(T); fmpz_clear(k); arb_clear(g); arb_clear(h); arb_clear(H); arb_clear(logT); return p; } /* Returns the number of zeros found. */ slong acb_dirichlet_platt_isolate_local_hardy_z_zeros( arf_interval_ptr res, const fmpz_t n, slong len, slong prec) { if (len <= 0 || fmpz_sizeinbase(n, 10) < 5) { return 0; } else if (fmpz_sgn(n) < 1) { flint_printf("Nonpositive indices of Hardy Z zeros are not supported.\n"); flint_abort(); } else { slong zeros_count = 0; platt_ctx_ptr ctx = _create_heuristic_context(n, prec); if (ctx) { zeros_count = _isolate_zeros(res, ctx, n, len, prec); platt_ctx_clear(ctx); free(ctx); } return zeros_count; } return 0; } /* Returns the number of zeros found. */ slong acb_dirichlet_platt_local_hardy_z_zeros( arb_ptr res, const fmpz_t n, slong len, slong prec) { if (len <= 0 || fmpz_sizeinbase(n, 10) < 5) { return 0; } else if (fmpz_sgn(n) < 1) { flint_printf("Nonpositive indices of Hardy Z zeros are not supported.\n"); flint_abort(); } else { slong zeros_count = 0; platt_ctx_ptr ctx; ctx = _create_heuristic_context(n, prec); if (ctx) { slong i; arf_interval_ptr p = _arf_interval_vec_init(len); zeros_count = _isolate_zeros(p, ctx, n, len, prec); for (i = 0; i < zeros_count; i++) { _refine_local_hardy_z_zero_illinois( res+i, ctx, &p[i].a, &p[i].b, prec); } _arf_interval_vec_clear(p, len); platt_ctx_clear(ctx); free(ctx); } return zeros_count; } return 0; }