1 /*
2 Copyright (C) 2010 Juan Arias de Reyna
3 Copyright (C) 2019 D.H.J. Polymath
4
5 This file is part of Arb.
6
7 Arb is free software: you can redistribute it and/or modify it under
8 the terms of the GNU Lesser General Public License (LGPL) as published
9 by the Free Software Foundation; either version 2.1 of the License, or
10 (at your option) any later version. See <http://www.gnu.org/licenses/>.
11 */
12
13 #include "acb_dirichlet.h"
14 #include "arb_calc.h"
15
16 /*
17 * For a detailed explanation of the algorithm implemented in this file, see:
18 *
19 * J. Arias de Reyna, "Programs for Riemann's zeta function", (J. A. J. van
20 * Vonderen, Ed.) Leven met getallen : liber amicorum ter gelegenheid van de
21 * pensionering van Herman te Riele, CWI (2012) 102-112,
22 * https://ir.cwi.nl/pub/19724
23 */
24
25 /*
26 * These constants define the largest n such that for every k <= n
27 * the kth zero of the Hardy Z function is governed by Gram's law
28 * or by Rosser's rule respectively.
29 */
30 static const slong GRAMS_LAW_MAX = 126;
31 static const slong ROSSERS_RULE_MAX = 13999526;
32
33 /*
34 * This structure describes a node of a doubly linked list.
35 * Each node represents a height t at which the Hardy Z-function
36 * Z(t) has been evaluated.
37 *
38 * As it relates to this data structure, the big picture of the algorithm
39 * to isolate the nth zero is roughly:
40 *
41 * 1) Initialize a two-node list consisting of the two points
42 * predicted by Gram's law to surround the nth zero.
43 * 2) Append and prepend as many additional Gram points as are indicated
44 * by the theory behind Turing's method, and occasionally insert points
45 * between existing points for the purpose of certifying blocks of intervals
46 * as 'good'.
47 * 3) Repeatedly insert new points into the list, interleaving existing points,
48 * until the number of sign changes indicated by theory is observed.
49 * This is where the algorithm would enter an infinite loop if it encountered
50 * a counterexample to the Riemann hypothesis.
51 * 4) Traverse the list until we find the pair of adjacent nodes with opposite
52 * signs of Z(t) that corresponds to the nth zero; the t heights of the
53 * two nodes define the endpoints of an isolating interval.
54 * 5) Delete the list.
55 */
56 typedef struct _zz_node_struct
57 {
58 arf_struct t; /* height t where v = Z(t) is evaluated */
59 arb_struct v; /* Z(t) */
60 fmpz *gram; /* Gram point index or NULL if not a Gram point */
61 slong prec; /* precision of evaluation of v */
62 struct _zz_node_struct *prev;
63 struct _zz_node_struct *next;
64 }
65 zz_node_struct;
66
67 typedef zz_node_struct zz_node_t[1];
68 typedef zz_node_struct * zz_node_ptr;
69 typedef const zz_node_struct * zz_node_srcptr;
70
71 static int
zz_node_is_gram_node(const zz_node_t p)72 zz_node_is_gram_node(const zz_node_t p)
73 {
74 return p->gram != NULL;
75 }
76
77 static int
zz_node_sgn(const zz_node_t p)78 zz_node_sgn(const zz_node_t p)
79 {
80 int s = arb_sgn_nonzero(&p->v);
81 if (!s)
82 {
83 flint_printf("unexpectedly imprecise evaluation of Z(t)\n");
84 flint_abort();
85 }
86 return s;
87 }
88
89 /* Good Gram points are Gram points where sgn(Z(g(n)))*(-1)^n > 0. */
90 static int
zz_node_is_good_gram_node(const zz_node_t p)91 zz_node_is_good_gram_node(const zz_node_t p)
92 {
93 if (zz_node_is_gram_node(p))
94 {
95 int s = zz_node_sgn(p);
96 if ((s > 0 && fmpz_is_even(p->gram)) ||
97 (s < 0 && fmpz_is_odd(p->gram)))
98 {
99 return 1;
100 }
101 }
102 return 0;
103 }
104
105 static void
zz_node_init(zz_node_t p)106 zz_node_init(zz_node_t p)
107 {
108 arf_init(&p->t);
109 arb_init(&p->v);
110 arb_indeterminate(&p->v);
111 p->prec = 0;
112 p->gram = NULL;
113 p->prev = NULL;
114 p->next = NULL;
115 }
116
117 static void
zz_node_clear(zz_node_t p)118 zz_node_clear(zz_node_t p)
119 {
120 arf_clear(&p->t);
121 arb_clear(&p->v);
122 if (p->gram)
123 {
124 fmpz_clear(p->gram);
125 flint_free(p->gram);
126 }
127 p->prec = 0;
128 p->gram = NULL;
129 p->prev = NULL;
130 p->next = NULL;
131 }
132
133 /*
134 * The node p represents the evaluation of the Hardy Z-function
135 * at a point t with some precision. This function re-evaluates
136 * the Hardy Z-function at t with greater precision, and it also
137 * guarantees that the precision is high enough to determine the
138 * sign of Z(t).
139 */
140 static int
zz_node_refine(zz_node_t p)141 zz_node_refine(zz_node_t p)
142 {
143 slong default_prec = arf_bits(&p->t) + 8;
144 if (p->prec < default_prec)
145 {
146 p->prec = default_prec;
147 }
148 else
149 {
150 p->prec *= 2;
151 }
152 return _acb_dirichlet_definite_hardy_z(&p->v, &p->t, &p->prec);
153 }
154
155 /*
156 * Create a node representing an evaluation of the Hardy Z-function
157 * at arbitrary point t. Upon creation the sign of Z(t) will
158 * be known with certainty.
159 */
160 static zz_node_ptr
create_non_gram_node(const arf_t t)161 create_non_gram_node(const arf_t t)
162 {
163 zz_node_ptr p = flint_malloc(sizeof(zz_node_struct));
164 zz_node_init(p);
165 arf_set(&p->t, t);
166 zz_node_refine(p);
167 return p;
168 }
169
170 /*
171 * Create a node representing an evaluation of the Hardy Z-function
172 * at the nth Gram point g(n). Upon creation a floating point number t will be
173 * assigned to this node, with the property that there are no zeros of the
174 * Hardy Z-function between t and the actual value of g(n).
175 * The sign of Z(t) will also be known with certainty.
176 */
177 static zz_node_ptr
create_gram_node(const fmpz_t n)178 create_gram_node(const fmpz_t n)
179 {
180 zz_node_ptr p;
181 arb_t t, v;
182 acb_t z;
183 slong prec = fmpz_bits(n) + 8;
184
185 arb_init(t);
186 arb_init(v);
187 acb_init(z);
188
189 while (1)
190 {
191 /* Computing the Gram point to higher precision improves performance
192 significantly. The likely explanation (not verified) is that
193 when evaluating the Z-function at an inexact ball using the
194 Riemann-Siegel formula, error propagation uses a bound for Z'
195 that is far from tight. The extra precision compensates
196 for this lack of tightness. */
197 acb_dirichlet_gram_point(t, n, NULL, NULL, prec + fmpz_bits(n));
198 acb_set_arb(z, t);
199 acb_dirichlet_hardy_z(z, z, NULL, NULL, 1, prec);
200 acb_get_real(v, z);
201 if (!arb_contains_zero(v))
202 {
203 break;
204 }
205 prec *= 2;
206 }
207
208 p = flint_malloc(sizeof(zz_node_struct));
209 zz_node_init(p);
210 p->gram = flint_malloc(sizeof(fmpz));
211 fmpz_init(p->gram);
212
213 /* t contains g(n) and does not contain a zero of the Z function */
214 fmpz_set(p->gram, n);
215 arf_set(&p->t, arb_midref(t));
216 arb_set(&p->v, v);
217 p->prec = prec;
218
219 arb_clear(t);
220 arb_clear(v);
221 acb_clear(z);
222
223 return p;
224 }
225
226 /*
227 * Count the number of Gram intervals between the Gram point
228 * represented by node a and the Gram point represented by node b.
229 * Traversing the linked list is not necessary because the Gram indices
230 * of nodes a and b can be accessed directly.
231 */
232 static slong
count_gram_intervals(zz_node_srcptr a,zz_node_srcptr b)233 count_gram_intervals(zz_node_srcptr a, zz_node_srcptr b)
234 {
235 slong out = 0;
236 if (!a || !b)
237 {
238 flint_printf("a and b must be non-NULL\n");
239 flint_abort();
240 }
241 if (!zz_node_is_good_gram_node(a) || !zz_node_is_good_gram_node(b))
242 {
243 flint_printf("both nodes must be good Gram points\n");
244 flint_abort();
245 }
246 else
247 {
248 fmpz_t m;
249 fmpz_init(m);
250 fmpz_sub(m, b->gram, a->gram);
251 out = fmpz_get_si(m);
252 fmpz_clear(m);
253 }
254 return out;
255 }
256
257 /*
258 * Count the observed number of sign changes of Z(t) by traversing
259 * a linked list of evaluated points from node a to node b.
260 */
261 static slong
count_sign_changes(zz_node_srcptr a,zz_node_srcptr b)262 count_sign_changes(zz_node_srcptr a, zz_node_srcptr b)
263 {
264 zz_node_srcptr p, q;
265 slong n = 0;
266 if (!a || !b)
267 {
268 flint_printf("a and b must be non-NULL\n");
269 flint_abort();
270 }
271 p = a;
272 q = a->next;
273 while (p != b)
274 {
275 if (!q)
276 {
277 flint_printf("prematurely reached end of list\n");
278 flint_abort();
279 }
280 if (zz_node_sgn(p) != zz_node_sgn(q))
281 {
282 n++;
283 }
284 p = q;
285 q = q->next;
286 }
287 return n;
288 }
289
290 /*
291 * Modify a linked list that ends with node p,
292 * by appending nodes representing Gram points.
293 * Continue until a 'good' Gram point is found.
294 */
295 static zz_node_ptr
extend_to_next_good_gram_node(zz_node_t p)296 extend_to_next_good_gram_node(zz_node_t p)
297 {
298 fmpz_t n;
299 zz_node_ptr q, r;
300
301 fmpz_init(n);
302
303 if (!zz_node_is_gram_node(p))
304 {
305 flint_printf("expected to begin at a gram point\n");
306 flint_abort();
307 }
308 if (p->next)
309 {
310 flint_printf("expected to extend from the end of a list\n");
311 flint_abort();
312 }
313 fmpz_set(n, p->gram);
314 q = p;
315 while (1)
316 {
317 fmpz_add_ui(n, n, 1);
318 r = create_gram_node(n);
319 q->next = r;
320 r->prev = q;
321 q = r;
322 r = NULL;
323 if (zz_node_is_good_gram_node(q))
324 {
325 break;
326 }
327 }
328
329 fmpz_clear(n);
330
331 return q;
332 }
333
334 /*
335 * Modify a linked list that begins with node p,
336 * by prepending nodes representing Gram points.
337 * Continue until a 'good' Gram point is found.
338 */
339 static zz_node_ptr
extend_to_prev_good_gram_node(zz_node_t p)340 extend_to_prev_good_gram_node(zz_node_t p)
341 {
342 fmpz_t n;
343 zz_node_ptr q, r;
344
345 fmpz_init(n);
346
347 if (!zz_node_is_gram_node(p))
348 {
349 flint_printf("expected to begin at a gram point\n");
350 flint_abort();
351 }
352 if (p->prev)
353 {
354 flint_printf("expected to extend from the start of a list\n");
355 flint_abort();
356 }
357 fmpz_set(n, p->gram);
358 q = p;
359 while (1)
360 {
361 fmpz_sub_ui(n, n, 1);
362 r = create_gram_node(n);
363 q->prev = r;
364 r->next = q;
365 q = r;
366 r = NULL;
367 if (zz_node_is_good_gram_node(q))
368 {
369 break;
370 }
371 }
372
373 fmpz_clear(n);
374
375 return q;
376 }
377
378 /*
379 * TODO: This is probably redundant.
380 * Can arb_get_lbound_arf ever give a negative arf given a nonnegative arb?
381 * If the answer is no, and it's probably no, then this function
382 * can be deleted.
383 */
384 static void
_arb_get_lbound_arf_nonnegative(arf_t res,const arb_t x,slong prec)385 _arb_get_lbound_arf_nonnegative(arf_t res, const arb_t x, slong prec)
386 {
387 arb_get_lbound_arf(res, x, prec);
388 if (arf_cmp_si(res, 0) < 0)
389 {
390 arf_zero(res);
391 }
392 }
393
394 /*
395 * res = (x1*w1 + x2*w2) / (w1 + w2)
396 * Undefined if weights are not nonnegative.
397 * If w1 and w2 are zero, the resulting interval contains x1 and x2.
398 */
399 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)400 _weighted_arithmetic_mean(arb_t res, const arf_t x1, const arf_t x2,
401 const arb_t w1, const arb_t w2, slong prec)
402 {
403 if (!arb_is_nonnegative(w1) || !arb_is_nonnegative(w2))
404 {
405 arb_indeterminate(res);
406 }
407 else if (arb_is_zero(w1) && arb_is_zero(w2))
408 {
409 arb_set_interval_arf(res, x1, x2, prec);
410 }
411 else if (arb_is_zero(w1))
412 {
413 arb_set_arf(res, x2);
414 }
415 else if (arb_is_zero(w2))
416 {
417 arb_set_arf(res, x1);
418 }
419 else if (arb_is_exact(w1) && arb_is_exact(w2))
420 {
421 arb_t a, b;
422 arb_init(a);
423 arb_init(b);
424 arb_mul_arf(a, w1, x1, prec);
425 arb_addmul_arf(a, w2, x2, prec);
426 arb_add(b, w1, w2, prec);
427 arb_div(res, a, b, prec);
428 arb_clear(a);
429 arb_clear(b);
430 }
431 else
432 {
433 arb_t a, b, r1, r2;
434 arb_init(a);
435 arb_init(b);
436 arb_init(r1);
437 arb_init(r2);
438
439 arb_zero(a);
440 arb_zero(b);
441 _arb_get_lbound_arf_nonnegative(arb_midref(a), w1, prec);
442 arb_get_ubound_arf(arb_midref(b), w2, prec);
443 _weighted_arithmetic_mean(r1, x1, x2, a, b, prec);
444
445 arb_zero(a);
446 arb_zero(b);
447 arb_get_ubound_arf(arb_midref(a), w1, prec);
448 _arb_get_lbound_arf_nonnegative(arb_midref(b), w2, prec);
449 _weighted_arithmetic_mean(r2, x1, x2, a, b, prec);
450
451 arb_union(res, r1, r2, prec);
452
453 arb_clear(a);
454 arb_clear(b);
455 arb_clear(r1);
456 arb_clear(r2);
457 }
458 }
459
460 /*
461 * Split the interval (t1, t2) into the intervals (t1, out) and (out, t2)
462 * in an attempt to increase the number of observed sign changes of Z(t)
463 * between endpoints.
464 * v1 and v2 are the Hardy Z-function values at t1 and t2 respectively.
465 * sign1 and sign2 are the signs of v1 and v2 respectively.
466 */
467 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)468 split_interval(arb_t out,
469 const arf_t t1, const arb_t v1, slong sign1,
470 const arf_t t2, const arb_t v2, slong sign2, slong prec)
471 {
472 if (sign1 == sign2)
473 {
474 /*
475 * out = (sqrt(v2/v1)*t1 + t2) / (sqrt(v2/v1) + 1)
476 * We have f(t1)=v1, f(t2)=v2 where v1 and v2 have the same sign,
477 * and we want to guess t between t1 and t2 so that f(t)
478 * has the opposite sign. Try the vertex of a parabola that would touch
479 * f(t)=0 between t1 and t2 and would pass through (t1,v1) and (t2,v2).
480 */
481 arb_t w1, w2;
482 arb_init(w1);
483 arb_init(w2);
484 arb_abs(w1, v2); /* w1, v2 is deliberate */
485 arb_sqrt(w1, w1, prec);
486 arb_abs(w2, v1); /* w2, v1 is deliberate */
487 arb_sqrt(w2, w2, prec);
488 _weighted_arithmetic_mean(out, t1, t2, w1, w2, prec);
489 arb_clear(w1);
490 arb_clear(w2);
491 }
492 else
493 {
494 /*
495 * out = (t1 + t2) / 2
496 * There is already one sign change in this interval.
497 * To find additional sign changes we would need to evaluate
498 * at least two more points in the interval,
499 * so begin by just splitting the interval in half at the midpoint.
500 */
501 arb_set_arf(out, t1);
502 arb_add_arf(out, out, t2, prec);
503 arb_mul_2exp_si(out, out, -1);
504 }
505 }
506
507 /*
508 * Add a new node between each pair of existing nodes in the linked list
509 * of evaluated values of t, within the sublist demarcated by nodes a and b.
510 */
511 static void
intercalate(zz_node_t a,zz_node_t b)512 intercalate(zz_node_t a, zz_node_t b)
513 {
514 arb_t t;
515 zz_node_ptr q, r, mid_node;
516
517 if (a == NULL || b == NULL)
518 {
519 flint_printf("a and b must be non-NULL\n");
520 flint_abort();
521 }
522 if (!zz_node_is_good_gram_node(a) || !zz_node_is_good_gram_node(b))
523 {
524 flint_printf("a and b must represent good Gram points\n");
525 flint_abort();
526 }
527
528 if (a == b) return;
529
530 arb_init(t);
531
532 q = a;
533 r = a->next;
534 while (q != b)
535 {
536 if (!r)
537 {
538 flint_printf("prematurely reached end of list\n");
539 flint_abort();
540 }
541 while (1)
542 {
543 split_interval(t,
544 &q->t, &q->v, zz_node_sgn(q),
545 &r->t, &r->v, zz_node_sgn(r),
546 FLINT_MIN(q->prec, r->prec));
547 if (!arb_contains_arf(t, &q->t) && !arb_contains_arf(t, &r->t))
548 {
549 break;
550 }
551 zz_node_refine((q->prec < r->prec) ? q : r);
552 }
553 mid_node = create_non_gram_node(arb_midref(t));
554 q->next = mid_node;
555 mid_node->prev = q;
556 mid_node->next = r;
557 r->prev = mid_node;
558 q = r;
559 r = r->next;
560 }
561
562 arb_clear(t);
563 }
564
565 /*
566 * Virtually trim k Gram/Rosser blocks from the head and from the tail
567 * of the sublist (a, b). The resulting sublist is (*A, *B).
568 * No nodes or connections between nodes are modified, added, or deleted.
569 */
570 static void
trim(zz_node_ptr * A,zz_node_ptr * B,zz_node_ptr a,zz_node_ptr b,slong k)571 trim(zz_node_ptr *A, zz_node_ptr *B,
572 zz_node_ptr a, zz_node_ptr b, slong k)
573 {
574 slong n;
575 for (n = 0; n < k; n++)
576 {
577 a = a->next;
578 while (!zz_node_is_good_gram_node(a))
579 {
580 a = a->next;
581 }
582 b = b->prev;
583 while (!zz_node_is_good_gram_node(b))
584 {
585 b = b->prev;
586 }
587 }
588 *A = a;
589 *B = b;
590 }
591
592 /*
593 * Given a linked sublist beginning at U and ending at V defining function
594 * evaluations at points that fully separate zeros of Z(t) in the vicinity
595 * of the nth zero, traverse the list until the nth zero is found.
596 * Continue traversing the list until len consecutive isolating intervals
597 * have been found, or until the end of the sublist is reached.
598 * Return the number of isolated zeros found, starting at the nth zero.
599 */
600 static slong
count_up_separated_zeros(arf_interval_ptr res,zz_node_srcptr U,zz_node_srcptr V,const fmpz_t n,slong len)601 count_up_separated_zeros(arf_interval_ptr res,
602 zz_node_srcptr U, zz_node_srcptr V, const fmpz_t n, slong len)
603 {
604 if (len <= 0)
605 {
606 return 0;
607 }
608 else if (fmpz_sgn(n) < 1)
609 {
610 flint_printf("nonpositive indices of zeros are not supported\n");
611 flint_abort();
612 }
613 else if (U == NULL || V == NULL)
614 {
615 flint_printf("U and V must not be NULL\n");
616 flint_abort();
617 }
618 if (!zz_node_is_good_gram_node(U) || !zz_node_is_good_gram_node(V))
619 {
620 flint_printf("U and V must be good Gram points\n");
621 flint_abort();
622 }
623 else
624 {
625 slong i = 0;
626 zz_node_srcptr p = U;
627 fmpz_t N, k;
628 fmpz_init(N);
629 fmpz_init(k);
630 fmpz_add_ui(N, p->gram, 1);
631 fmpz_set(k, n);
632 while (p != V)
633 {
634 if (!p->next)
635 {
636 flint_printf("prematurely reached end of list\n");
637 flint_abort();
638 }
639 if (zz_node_sgn(p) != zz_node_sgn(p->next))
640 {
641 fmpz_add_ui(N, N, 1);
642 if (fmpz_equal(N, k))
643 {
644 arf_set(&res[i].a, &p->t);
645 arf_set(&res[i].b, &p->next->t);
646 fmpz_add_ui(k, k, 1);
647 i++;
648 if (i == len)
649 break;
650 }
651 }
652 p = p->next;
653 }
654 fmpz_clear(k);
655 return i;
656 }
657 return 0;
658 }
659
660 /*
661 * Find one 'superblock' below n and one 'superblock' above n.
662 * The term 'superblock' in this context means a stretch of
663 * enough consecutive 'good' Rosser/Gram blocks to meet the Turing method bound.
664 * The output *psb is the number of blocks in the superblock.
665 * The output nodes *pu and *pv are the first node of the first superblock
666 * and the last node of the second superblock respectively.
667 */
668 static void
turing_search_near(zz_node_ptr * pu,zz_node_ptr * pv,slong * psb,const fmpz_t n)669 turing_search_near(zz_node_ptr *pu, zz_node_ptr *pv, slong *psb, const fmpz_t n)
670 {
671 zz_node_ptr u, v;
672 slong i;
673 slong zn; /* target number of sign changes */
674 slong sb; /* the number of padding blocks required by Turing's method */
675 slong cgb; /* the number of consecutive good blocks */
676 const slong loopcount = 4;
677 fmpz_t k;
678
679 fmpz_init(k);
680
681 fmpz_sub_ui(k, n, 2);
682 u = create_gram_node(k);
683 fmpz_sub_ui(k, n, 1);
684 v = create_gram_node(k);
685 u->next = v;
686 v->prev = u;
687
688 if (!zz_node_is_good_gram_node(u))
689 u = extend_to_prev_good_gram_node(u);
690 if (!zz_node_is_good_gram_node(v))
691 v = extend_to_next_good_gram_node(v);
692
693 /*
694 * Extend the search to greater heights t.
695 *
696 * Continue adding Gram points until the number of consecutive
697 * 'good' Gram/Rosser blocks reaches the Turing method bound.
698 */
699 sb = 0;
700 cgb = 0;
701 while (1)
702 {
703 zz_node_ptr nv;
704 nv = extend_to_next_good_gram_node(v);
705 zn = count_gram_intervals(v, nv);
706 for (i = 0; i < loopcount && count_sign_changes(v, nv) < zn; i++)
707 {
708 intercalate(v, nv);
709 }
710 if (count_sign_changes(v, nv) >= zn)
711 {
712 cgb++;
713 if (cgb > sb)
714 {
715 sb = cgb;
716 if (acb_dirichlet_turing_method_bound(nv->gram) <= sb)
717 {
718 v = nv;
719 break;
720 }
721 }
722 }
723 else
724 {
725 cgb = 0;
726 }
727 v = nv;
728 }
729
730 /* Extend the search to smaller heights t. */
731 cgb = 0;
732 while (1)
733 {
734 zz_node_ptr pu;
735 pu = extend_to_prev_good_gram_node(u);
736 zn = count_gram_intervals(pu, u);
737 for (i = 0; i < loopcount && count_sign_changes(pu, u) < zn; i++)
738 {
739 intercalate(pu, u);
740 }
741 if (count_sign_changes(pu, u) >= zn)
742 {
743 cgb++;
744 if (cgb == sb)
745 {
746 u = pu;
747 break;
748 }
749 }
750 else
751 {
752 cgb = 0;
753 }
754 u = pu;
755 }
756
757 *pu = u;
758 *pv = v;
759 *psb = sb;
760
761 fmpz_clear(k);
762 }
763
764 /*
765 * Find one 'double superblock' beginning below the point represented
766 * by node u, and find one 'double superblock' ending above the point
767 * represented by node v. The term 'double superblock' in this context
768 * means a stretch of twice as many consecutive 'good' Rosser/Gram blocks
769 * as would meet the Turing method bound.
770 * The output nodes *pu and *pv are the first node of the first double
771 * superblock and the last node of the second double superblock respectively.
772 * The output integer *psb reports one half of the number
773 * of blocks in the double superblock.
774 * The parameter initial_cgb is the number of consecutive good blocks
775 * above and below the points represented by nodes u and v respectively.
776 */
777 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)778 turing_search_far(zz_node_ptr *pu, zz_node_ptr *pv, slong *psb,
779 zz_node_ptr u, zz_node_ptr v, slong initial_cgb)
780 {
781 slong i;
782 slong zn; /* target number of sign changes */
783 slong sb; /* the number of padding blocks required by Turing's method */
784 slong cgb; /* the number of consecutive good blocks */
785 const slong loopcount = 4;
786
787 /*
788 * Extend the search to greater heights t.
789 *
790 * Continue adding Gram points until the number of consecutive
791 * 'good' Gram/Rosser blocks is twice the number required by
792 * the Turing method bound.
793 */
794 sb = 0;
795 cgb = initial_cgb;
796 while (1)
797 {
798 zz_node_ptr nv;
799 nv = extend_to_next_good_gram_node(v);
800 zn = count_gram_intervals(v, nv);
801 for (i = 0; i < loopcount && count_sign_changes(v, nv) < zn; i++)
802 {
803 intercalate(v, nv);
804 }
805 if (count_sign_changes(v, nv) >= zn)
806 {
807 cgb++;
808 if (cgb % 2 == 0 && sb < cgb / 2)
809 {
810 sb = cgb / 2;
811 if (acb_dirichlet_turing_method_bound(nv->gram) <= sb)
812 {
813 v = nv;
814 break;
815 }
816 }
817 }
818 else
819 {
820 cgb = 0;
821 }
822 v = nv;
823 }
824
825 /* Extend the search to smaller heights t. */
826 cgb = initial_cgb;
827 while (1)
828 {
829 zz_node_ptr pu;
830 pu = extend_to_prev_good_gram_node(u);
831 zn = count_gram_intervals(pu, u);
832 for (i = 0; i < loopcount && count_sign_changes(pu, u) < zn; i++)
833 {
834 intercalate(pu, u);
835 }
836 if (count_sign_changes(pu, u) >= zn)
837 {
838 cgb++;
839 if (cgb == sb*2)
840 {
841 u = pu;
842 break;
843 }
844 }
845 else
846 {
847 cgb = 0;
848 }
849 u = pu;
850 }
851
852 *pu = u;
853 *pv = v;
854 *psb = sb;
855 }
856
857 /*
858 * Used for both zero isolation and for N(t).
859 * n is the index of a Hardy Z-function zero of interest.
860 * *pu and *pv are the first and last node of an output list.
861 * *pU and *pV are the first and last nodes of a sublist with
862 * fully separated zeros, within the *pu -- *pv list.
863 */
864 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)865 _separated_turing_list(zz_node_ptr *pU, zz_node_ptr *pV,
866 zz_node_ptr *pu, zz_node_ptr *pv, const fmpz_t n)
867 {
868 zz_node_ptr U, V, u, v;
869 slong i;
870 slong sb_near; /* Turing method bound for near search */
871 slong sb_far; /* Turing method bound for far search */
872 slong zn; /* target number of sign changes */
873 slong variations; /* observed number of sign changes */
874 const slong loopcount = 4;
875 /*
876 * The loopcount controls how hard we try to find zeros in Gram/Rosser
877 * blocks. If the loopcount is too high then when Rosser's rule is violated
878 * we will spend too much time searching for zeros that don't exist.
879 * If the loopcount is too low then we will miss some 'good' blocks,
880 * causing the search to be extended unnecessarily far from the initial
881 * guess before eventually making another pass in which loopcount
882 * is effectively infinite.
883 */
884
885 if (fmpz_cmp_si(n, 2) < 0)
886 {
887 flint_printf("invalid n: "); fmpz_print(n); flint_printf("\n");
888 flint_abort();
889 }
890
891 turing_search_near(&u, &v, &sb_near, n);
892 trim(&U, &V, u, v, sb_near);
893 zn = count_gram_intervals(U, V);
894 for (i = 0; i < loopcount && count_sign_changes(U, V) < zn; i++)
895 {
896 intercalate(U, V);
897 }
898 variations = count_sign_changes(U, V);
899 if (variations > zn)
900 {
901 flint_printf("unexpected number of sign changes\n");
902 flint_abort();
903 }
904 else if (variations < zn)
905 {
906 zz_node_ptr r = U;
907 zz_node_ptr s = V;
908 turing_search_far(&u, &v, &sb_far, u, v, sb_near);
909 trim(&U, &V, u, v, 2*sb_far);
910 zn = count_gram_intervals(U, V);
911 for (i = 0; i < loopcount && count_sign_changes(U, V) < zn; i++)
912 {
913 intercalate(U, r);
914 intercalate(s, V);
915 }
916 variations = count_sign_changes(U, V);
917 if (variations > zn)
918 {
919 flint_printf("unexpected number of sign changes\n");
920 flint_abort();
921 }
922 else if (variations < zn)
923 {
924 trim(&U, &V, u, v, sb_far);
925 zn = count_gram_intervals(U, V);
926 while (count_sign_changes(U, V) < zn)
927 {
928 intercalate(U, V);
929 }
930 if (count_sign_changes(U, V) != zn)
931 {
932 flint_printf("unexpected number of sign changes\n");
933 flint_abort();
934 }
935 }
936 }
937
938 *pu = u;
939 *pv = v;
940 *pU = U;
941 *pV = V;
942 }
943
944 /*
945 * Used for both zero isolation and for N(t).
946 * n is the index of a Hardy Z-function zero of interest.
947 * *pu and *pv are the first and last node of an output list
948 * with fully separated zeros.
949 */
950 static void
_separated_rosser_list(zz_node_ptr * pu,zz_node_ptr * pv,const fmpz_t n)951 _separated_rosser_list(zz_node_ptr *pu, zz_node_ptr *pv, const fmpz_t n)
952 {
953 fmpz_t k;
954 zz_node_ptr u, v;
955
956 if (fmpz_cmp_si(n, 1) < 0 || fmpz_cmp_si(n, ROSSERS_RULE_MAX) > 0)
957 {
958 flint_printf("invalid n: "); fmpz_print(n); flint_printf("\n");
959 flint_abort();
960 }
961
962 fmpz_init(k);
963
964 fmpz_sub_ui(k, n, 2);
965 u = create_gram_node(k);
966 fmpz_sub_ui(k, n, 1);
967 v = create_gram_node(k);
968 u->next = v;
969 v->prev = u;
970
971 if (!zz_node_is_good_gram_node(u))
972 u = extend_to_prev_good_gram_node(u);
973 if (!zz_node_is_good_gram_node(v))
974 v = extend_to_next_good_gram_node(v);
975 while (count_sign_changes(u, v) != count_gram_intervals(u, v))
976 {
977 intercalate(u, v);
978 }
979
980 *pu = u;
981 *pv = v;
982
983 fmpz_clear(k);
984 }
985
986 /*
987 * Used for both zero isolation and for N(t).
988 * n is the index of a Hardy Z-function zero of interest.
989 * *pu and *pv are the first and last node of an output list
990 * with fully separated zeros.
991 */
992 static void
_separated_gram_list(zz_node_ptr * pu,zz_node_ptr * pv,const fmpz_t n)993 _separated_gram_list(zz_node_ptr *pu, zz_node_ptr *pv, const fmpz_t n)
994 {
995 fmpz_t k;
996 zz_node_ptr u, v;
997
998 if (fmpz_cmp_si(n, 1) < 0 || fmpz_cmp_si(n, GRAMS_LAW_MAX) > 0)
999 {
1000 flint_printf("invalid n: "); fmpz_print(n); flint_printf("\n");
1001 flint_abort();
1002 }
1003
1004 fmpz_init(k);
1005
1006 fmpz_sub_ui(k, n, 2);
1007 u = create_gram_node(k);
1008 fmpz_sub_ui(k, n, 1);
1009 v = create_gram_node(k);
1010 u->next = v;
1011 v->prev = u;
1012
1013 *pu = u;
1014 *pv = v;
1015
1016 fmpz_clear(k);
1017 }
1018
1019 /*
1020 * Get a list of points that fully separate zeros of Z.
1021 * Used for both zero isolation and for N(t).
1022 * n is the index of a Hardy Z-function zero of interest.
1023 * *pu and *pv are the first and last node of an output list.
1024 * *pU and *pV are the first and last nodes of a sublist with
1025 * fully separated zeros, within the *pu -- *pv list.
1026 */
1027 static void
_separated_list(zz_node_ptr * pU,zz_node_ptr * pV,zz_node_ptr * pu,zz_node_ptr * pv,const fmpz_t n)1028 _separated_list(zz_node_ptr *pU, zz_node_ptr *pV,
1029 zz_node_ptr *pu, zz_node_ptr *pv, const fmpz_t n)
1030 {
1031 zz_node_ptr U, V, u, v;
1032 if (fmpz_cmp_si(n, GRAMS_LAW_MAX) <= 0)
1033 {
1034 _separated_gram_list(&u, &v, n);
1035 U = u;
1036 V = v;
1037 }
1038 else if (fmpz_cmp_si(n, ROSSERS_RULE_MAX) <= 0)
1039 {
1040 _separated_rosser_list(&u, &v, n);
1041 U = u;
1042 V = v;
1043 }
1044 else
1045 {
1046 _separated_turing_list(&U, &V, &u, &v, n);
1047 }
1048
1049 if (U == NULL || V == NULL)
1050 {
1051 flint_printf("U and V must not be NULL\n");
1052 flint_abort();
1053 }
1054 if (!zz_node_is_good_gram_node(U) || !zz_node_is_good_gram_node(V))
1055 {
1056 flint_printf("U and V must be good Gram points\n");
1057 flint_abort();
1058 }
1059 if (U == V)
1060 {
1061 flint_printf("the list must include at least one interval\n");
1062 flint_abort();
1063 }
1064
1065 *pU = U;
1066 *pV = V;
1067 *pu = u;
1068 *pv = v;
1069 }
1070
1071 /*
1072 * Isolate up to len zeros, starting from the nth zero.
1073 * Return the number of isolated zeros.
1074 */
1075 static slong
_isolate_hardy_z_zeros(arf_interval_ptr res,const fmpz_t n,slong len)1076 _isolate_hardy_z_zeros(arf_interval_ptr res, const fmpz_t n, slong len)
1077 {
1078 zz_node_ptr U, V, u, v;
1079 slong count;
1080
1081 _separated_list(&U, &V, &u, &v, n);
1082 count = count_up_separated_zeros(res, U, V, n, len);
1083
1084 while (u)
1085 {
1086 v = u;
1087 u = u->next;
1088 zz_node_clear(v);
1089 flint_free(v);
1090 }
1091
1092 return count;
1093 }
1094
1095 /* Isolate len zeros, starting from the nth zero. */
1096 void
acb_dirichlet_isolate_hardy_z_zeros(arf_interval_ptr res,const fmpz_t n,slong len)1097 acb_dirichlet_isolate_hardy_z_zeros(arf_interval_ptr res, const fmpz_t n, slong len)
1098 {
1099 if (len <= 0)
1100 {
1101 return;
1102 }
1103 else if (fmpz_sgn(n) < 1)
1104 {
1105 flint_printf("nonpositive indices of zeros are not supported\n");
1106 flint_abort();
1107 }
1108 else
1109 {
1110 slong c = 0;
1111 fmpz_t k;
1112 fmpz_init(k);
1113 while (c < len)
1114 {
1115 fmpz_add_si(k, n, c);
1116 c += _isolate_hardy_z_zeros(res + c, k, len - c);
1117 }
1118 fmpz_clear(k);
1119 }
1120 }
1121
1122 void
acb_dirichlet_isolate_hardy_z_zero(arf_t a,arf_t b,const fmpz_t n)1123 acb_dirichlet_isolate_hardy_z_zero(arf_t a, arf_t b, const fmpz_t n)
1124 {
1125 if (fmpz_sgn(n) < 1)
1126 {
1127 flint_printf("nonpositive indices of zeros are not supported\n");
1128 flint_abort();
1129 }
1130 else
1131 {
1132 arf_interval_t r;
1133 arf_interval_init(r);
1134 _isolate_hardy_z_zeros(r, n, 1);
1135 arf_set(a, &r->a);
1136 arf_set(b, &r->b);
1137 arf_interval_clear(r);
1138 }
1139 }
1140
1141 static void
count_up(arf_t a,arf_t b,zz_node_srcptr U,zz_node_srcptr V,const fmpz_t n)1142 count_up(arf_t a, arf_t b, zz_node_srcptr U, zz_node_srcptr V, const fmpz_t n)
1143 {
1144 arf_interval_t r;
1145 arf_interval_init(r);
1146 count_up_separated_zeros(r, U, V, n, 1);
1147 arf_set(a, &r->a);
1148 arf_set(b, &r->b);
1149 arf_interval_clear(r);
1150 }
1151
1152 void
_acb_dirichlet_isolate_turing_hardy_z_zero(arf_t a,arf_t b,const fmpz_t n)1153 _acb_dirichlet_isolate_turing_hardy_z_zero(arf_t a, arf_t b, const fmpz_t n)
1154 {
1155 zz_node_ptr U, V, u, v;
1156 _separated_turing_list(&U, &V, &u, &v, n);
1157 count_up(a, b, U, V, n);
1158 while (u)
1159 {
1160 v = u;
1161 u = u->next;
1162 zz_node_clear(v);
1163 flint_free(v);
1164 }
1165 }
1166
1167 void
_acb_dirichlet_isolate_rosser_hardy_z_zero(arf_t a,arf_t b,const fmpz_t n)1168 _acb_dirichlet_isolate_rosser_hardy_z_zero(arf_t a, arf_t b, const fmpz_t n)
1169 {
1170 zz_node_ptr u, v;
1171 _separated_rosser_list(&u, &v, n);
1172 count_up(a, b, u, v, n);
1173 while (u)
1174 {
1175 v = u;
1176 u = u->next;
1177 zz_node_clear(v);
1178 flint_free(v);
1179 }
1180 }
1181
1182 void
_acb_dirichlet_isolate_gram_hardy_z_zero(arf_t a,arf_t b,const fmpz_t n)1183 _acb_dirichlet_isolate_gram_hardy_z_zero(arf_t a, arf_t b, const fmpz_t n)
1184 {
1185 zz_node_ptr u, v;
1186 _separated_gram_list(&u, &v, n);
1187 count_up(a, b, u, v, n);
1188 while (u)
1189 {
1190 v = u;
1191 u = u->next;
1192 zz_node_clear(v);
1193 flint_free(v);
1194 }
1195 }
1196
1197 static void
_acb_set_arf(acb_t res,const arf_t t)1198 _acb_set_arf(acb_t res, const arf_t t)
1199 {
1200 acb_zero(res);
1201 arb_set_arf(acb_realref(res), t);
1202 }
1203
1204 /*
1205 * Find the index of the largest Gram point less than t.
1206 * Requires t > 10.
1207 */
1208 static void
gram_index(fmpz_t res,const arf_t t)1209 gram_index(fmpz_t res, const arf_t t)
1210 {
1211 if (arf_cmp_si(t, 10) < 0)
1212 {
1213 flint_printf("gram_index requires t > 10\n");
1214 flint_abort();
1215 }
1216 else
1217 {
1218 acb_t z;
1219 slong prec = arf_abs_bound_lt_2exp_si(t) + 8;
1220 acb_init(z);
1221 while (1)
1222 {
1223 _acb_set_arf(z, t);
1224 acb_dirichlet_hardy_theta(z, z, NULL, NULL, 1, prec);
1225 arb_const_pi(acb_imagref(z), prec);
1226 arb_div(acb_realref(z), acb_realref(z), acb_imagref(z), prec);
1227 arb_floor(acb_realref(z), acb_realref(z), prec);
1228 if (arb_get_unique_fmpz(res, acb_realref(z)))
1229 {
1230 break;
1231 }
1232 prec *= 2;
1233 }
1234 acb_clear(z);
1235 }
1236 }
1237
1238 /*
1239 * Compute nzeros(t) for up to len values of t
1240 * and return the number computed. p must be in increasing order,
1241 * and all entries must be greater than 10.
1242 */
1243 static slong
_exact_zeta_multi_nzeros(fmpz * res,arf_srcptr points,slong len)1244 _exact_zeta_multi_nzeros(fmpz *res, arf_srcptr points, slong len)
1245 {
1246 zz_node_ptr U, V, u, v, p;
1247 arb_t x;
1248 fmpz_t n, N;
1249 slong i;
1250 arf_srcptr t;
1251 int s;
1252 slong prec;
1253
1254 arb_init(x);
1255 fmpz_init(n);
1256 fmpz_init(N);
1257
1258 /*
1259 * The first point is located between two Gram points. Find the unique n
1260 * such that the nth zero is also predicted to be located between these
1261 * Gram points.
1262 */
1263 gram_index(n, points);
1264 fmpz_add_ui(n, n, 2);
1265
1266 /* Get a list of points that fully separate zeros of Z. */
1267 _separated_list(&U, &V, &u, &v, n);
1268
1269 p = U;
1270 fmpz_add_ui(N, U->gram, 1);
1271
1272 /*
1273 * It's possible that one or more points are located between
1274 * the first Gram point and the arf_t representative of that Gram point.
1275 */
1276 for (i = 0, t = points; i < len && arf_cmp(t, &p->t) <= 0; i++, t++)
1277 {
1278 fmpz_set(res + i, N);
1279 }
1280
1281 while (i < len && p != V)
1282 {
1283 if (!p->next)
1284 {
1285 flint_printf("prematurely reached the end of the list\n");
1286 flint_abort();
1287 }
1288 if (zz_node_sgn(p) != zz_node_sgn(p->next))
1289 {
1290 while (i < len && arf_cmp(t, &p->next->t) <= 0)
1291 {
1292 prec = arf_abs_bound_lt_2exp_si(t) + 8;
1293 s = _acb_dirichlet_definite_hardy_z(x, t, &prec);
1294 if (zz_node_sgn(p->next) == s)
1295 {
1296 fmpz_add_ui(res + i, N, 1);
1297 }
1298 else
1299 {
1300 fmpz_set(res + i, N);
1301 }
1302 t++;
1303 i++;
1304 }
1305 fmpz_add_ui(N, N, 1);
1306 }
1307 p = p->next;
1308 }
1309
1310 /*
1311 * It's possible that the first point in the array is located between
1312 * the last Gram point and the arf_t representative of that Gram point.
1313 */
1314 if (i == 0)
1315 {
1316 fmpz_set(res, N);
1317 i++;
1318 }
1319
1320 while (u)
1321 {
1322 v = u;
1323 u = u->next;
1324 zz_node_clear(v);
1325 flint_free(v);
1326 }
1327
1328 arb_clear(x);
1329 fmpz_clear(n);
1330 fmpz_clear(N);
1331
1332 return i;
1333 }
1334
1335 /*
1336 * Compute nzeros for len values of t. The array p must be in increasing order.
1337 */
1338 static void
exact_zeta_multi_nzeros(fmpz * res,arf_srcptr p,slong len)1339 exact_zeta_multi_nzeros(fmpz *res, arf_srcptr p, slong len)
1340 {
1341 slong i, c;
1342 arf_srcptr q;
1343
1344 if (len <= 0)
1345 {
1346 return;
1347 }
1348
1349 for (i = 0, q = p; i < len - 1; i++, q++)
1350 {
1351 if (arf_cmp(q, p) > 0)
1352 {
1353 flint_printf("p must be in increasing order\n");
1354 flint_abort();
1355 }
1356 }
1357
1358 c = 0;
1359 while (c < len)
1360 {
1361 if (arf_cmp_si(p + c, 14) < 0)
1362 {
1363 fmpz_zero(res + c);
1364 c++;
1365 }
1366 else
1367 {
1368 c += _exact_zeta_multi_nzeros(res + c, p + c, len - c);
1369 }
1370 }
1371 }
1372
1373 void
_acb_dirichlet_exact_zeta_nzeros(fmpz_t res,const arf_t t)1374 _acb_dirichlet_exact_zeta_nzeros(fmpz_t res, const arf_t t)
1375 {
1376 exact_zeta_multi_nzeros(res, t, 1);
1377 }
1378
1379 void
acb_dirichlet_hardy_z_zeros(arb_ptr res,const fmpz_t n,slong len,slong prec)1380 acb_dirichlet_hardy_z_zeros(arb_ptr res, const fmpz_t n, slong len, slong prec)
1381 {
1382 if (len <= 0)
1383 {
1384 return;
1385 }
1386 else if (fmpz_sgn(n) < 1)
1387 {
1388 flint_printf("nonpositive indices of zeros are not supported\n");
1389 flint_abort();
1390 }
1391 else
1392 {
1393 slong i;
1394 arf_interval_ptr p = _arf_interval_vec_init(len);
1395 acb_dirichlet_isolate_hardy_z_zeros(p, n, len);
1396 for (i = 0; i < len; i++)
1397 {
1398 _acb_dirichlet_refine_hardy_z_zero(res + i, &p[i].a, &p[i].b, prec);
1399 }
1400 _arf_interval_vec_clear(p, len);
1401 }
1402 }
1403
1404 static void
_arb_set_interval_fmpz(arb_t res,const fmpz_t a,const fmpz_t b)1405 _arb_set_interval_fmpz(arb_t res, const fmpz_t a, const fmpz_t b)
1406 {
1407 fmpz_t n;
1408
1409 fmpz_init(n);
1410
1411 fmpz_add(n, a, b);
1412 arf_set_fmpz(arb_midref(res), n);
1413 fmpz_sub(n, b, a);
1414 mag_set_fmpz(arb_radref(res), n);
1415 arb_mul_2exp_si(res, res, -1);
1416
1417 fmpz_clear(n);
1418 }
1419
1420 void
acb_dirichlet_zeta_nzeros(arb_t res,const arb_t t,slong prec)1421 acb_dirichlet_zeta_nzeros(arb_t res, const arb_t t, slong prec)
1422 {
1423 if (arb_is_exact(t))
1424 {
1425 fmpz_t n;
1426 fmpz_init(n);
1427 _acb_dirichlet_exact_zeta_nzeros(n, arb_midref(t));
1428 arb_set_fmpz(res, n);
1429 fmpz_clear(n);
1430 }
1431 else
1432 {
1433 arf_struct b[2];
1434 fmpz n[2];
1435
1436 arf_init(b);
1437 arf_init(b + 1);
1438 fmpz_init(n);
1439 fmpz_init(n + 1);
1440
1441 arb_get_lbound_arf(b, t, prec);
1442 arb_get_ubound_arf(b + 1, t, prec);
1443 exact_zeta_multi_nzeros(n, b, 2);
1444 _arb_set_interval_fmpz(res, n, n + 1);
1445
1446 arf_clear(b);
1447 arf_clear(b + 1);
1448 fmpz_clear(n);
1449 fmpz_clear(n + 1);
1450 }
1451
1452 arb_set_round(res, res, prec);
1453 }
1454
1455 void
acb_dirichlet_zeta_nzeros_gram(fmpz_t res,const fmpz_t n)1456 acb_dirichlet_zeta_nzeros_gram(fmpz_t res, const fmpz_t n)
1457 {
1458 zz_node_ptr U, V, u, v, p;
1459 fmpz_t k, N;
1460
1461 if (fmpz_cmp_si(n, -1) < 0)
1462 {
1463 flint_printf("n must be >= -1\n");
1464 flint_abort();
1465 }
1466
1467 fmpz_init(k);
1468 fmpz_init(N);
1469
1470 /*
1471 * Get a list of points that fully separate zeros of Z.
1472 * The kth zero is expected to be between the k-2 and the k-1 gram points.
1473 */
1474 fmpz_add_ui(k, n, 2);
1475 _separated_list(&U, &V, &u, &v, k);
1476
1477 p = U;
1478 fmpz_add_ui(N, U->gram, 1);
1479 fmpz_set_si(res, -1);
1480 while (1)
1481 {
1482 if (p == NULL)
1483 {
1484 flint_printf("prematurely reached the end of the list\n");
1485 flint_abort();
1486 }
1487 if (zz_node_is_gram_node(p) && fmpz_equal(n, p->gram))
1488 {
1489 fmpz_set(res, N);
1490 break;
1491 }
1492 if (zz_node_sgn(p) != zz_node_sgn(p->next))
1493 {
1494 fmpz_add_ui(N, N, 1);
1495 }
1496 if (p == V)
1497 {
1498 break;
1499 }
1500 p = p->next;
1501 }
1502
1503 if (fmpz_sgn(res) < 0)
1504 {
1505 flint_printf("failed to find the gram point\n");
1506 flint_abort();
1507 }
1508
1509 while (u)
1510 {
1511 v = u;
1512 u = u->next;
1513 zz_node_clear(v);
1514 flint_free(v);
1515 }
1516
1517 fmpz_clear(k);
1518 fmpz_clear(N);
1519 }
1520