Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - modules - ratpoints.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27102-e8476a37b5) Lines: 839 861 97.4 %
Date: 2021-12-02 07:05:14 Functions: 33 33 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2017  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /* This file is based on ratpoints-2.1.3 by Michael Stoll, see
      16             :  * http://www.mathe2.uni-bayreuth.de/stoll/programs/
      17             :  * Original copyright / license: */
      18             : /***********************************************************************
      19             :  * ratpoints-2.1.3                                                     *
      20             :  * Copyright (C) 2008, 2009  Michael Stoll                             *
      21             :  *  - A program to find rational points on hyperelliptic curves        *
      22             :  *                                                                     *
      23             :  * This program is free software: you can redistribute it and/or       *
      24             :  * modify it under the terms of the GNU General Public License         *
      25             :  * as published by the Free Software Foundation, either version 2 of   *
      26             :  * the License, or (at your option) any later version.                 *
      27             :  ***********************************************************************/
      28             : 
      29             : #include "pari.h"
      30             : #include "paripriv.h"
      31             : 
      32             : #define pel(a,b)  gel((a),(b)+2)
      33             : 
      34             : #define RATPOINTS_ARRAY_SIZE 256           /* Array size in longs */
      35             : #define RATPOINTS_DEFAULT_SP1 9            /* Default value for sp1 */
      36             : #define RATPOINTS_DEFAULT_SP2 16           /* Default value for sp2 */
      37             : #define RATPOINTS_DEFAULT_NUM_PRIMES 28    /* Default value for num_primes */
      38             : #define RATPOINTS_DEFAULT_BIT_PRIMES 7     /* Default value for bit_primes */
      39             : #define RATPOINTS_DEFAULT_MAX_FORBIDDEN 30 /* Default value for max_forbidden */
      40             : 
      41             : typedef struct {double low; double up;} ratpoints_interval;
      42             : 
      43             : /* Define the flag bits for the flags component: */
      44             : #define RATPOINTS_NO_REVERSE      0x0004UL
      45             : 
      46             : #define RATPOINTS_FLAGS_INPUT_MASK (RATPOINTS_NO_REVERSE)
      47             : 
      48             : /* Flags bits for internal purposes */
      49             : #define RATPOINTS_REVERSED        0x0100UL
      50             : #define RATPOINTS_CHECK_DENOM     0x0200UL
      51             : #define RATPOINTS_USE_SQUARES     0x0400UL
      52             : #define RATPOINTS_USE_SQUARES1    0x0800UL
      53             : 
      54             : #define LONG_MASK (~(-(1UL<<TWOPOTBITS_IN_LONG)))
      55             : 
      56             : #define CEIL(a,b) (((a) <= 0) ? -(-(a) / (b)) : 1 + ((a)-1) / (b))
      57             : 
      58             : /* Some Macros for working with SSE registers */
      59             : #ifdef HAS_SSE2
      60             : #include <emmintrin.h>
      61             : 
      62             : #define AND(a,b) ((ratpoints_bit_array)__builtin_ia32_andps((__v4sf)(a), (__v4sf)(b)))
      63             : #define EXT0(a) ((ulong)__builtin_ia32_vec_ext_v2di((__v2di)(a), 0))
      64             : #define EXT1(a) ((ulong)__builtin_ia32_vec_ext_v2di((__v2di)(a), 1))
      65             : #define TEST(a) (EXT0(a) || EXT1(a))
      66             : #define RBA(a) ((__v2di){((long) a), ((long) a)})
      67             : 
      68             : /* Use SSE 128 bit registers for the bit arrays */
      69             : typedef __v2di ratpoints_bit_array;
      70             : 
      71             : #define RBA_LENGTH (128)
      72             : #define RBA_SHIFT (7)
      73             : #define RBA_ALIGN  (sizeof(ratpoints_bit_array))
      74             : #define RBA_MASK (~(-(1UL<<RBA_SHIFT)))
      75             : 
      76             : #else
      77             : /* Use ulong for the bit arrays */
      78             : typedef ulong ratpoints_bit_array;
      79             : #define AND(a,b) ((a)&(b))
      80             : #define EXT0(a) (a)
      81             : #define TEST(a) (a)
      82             : #define RBA(a) (a)
      83             : 
      84             : #define RBA_LENGTH BITS_IN_LONG
      85             : #define RBA_SHIFT TWOPOTBITS_IN_LONG
      86             : #define RBA_ALIGN  (sizeof(long))
      87             : #define RBA_MASK LONG_MASK
      88             : 
      89             : #endif
      90             : 
      91             : typedef struct { long p; long offset; ratpoints_bit_array *ptr; } sieve_spec;
      92             : 
      93             : typedef enum { num_all, num_even, num_odd, num_none } bit_selection;
      94             : 
      95             : typedef struct {
      96             :   long p; int *is_f_square;
      97             :   const long *inverses;
      98             :   long offset; ratpoints_bit_array** sieve;
      99             : } ratpoints_sieve_entry;
     100             : 
     101             : typedef struct { long p;
     102             :                  ulong *start;
     103             :                  ulong *end;
     104             :                  ulong *curr; }
     105             :                forbidden_entry;
     106             : 
     107             : typedef struct {
     108             :   GEN cof, listprime;
     109             :   ratpoints_interval *domain;
     110             :   long height, b_low, b_high, sp1, sp2, array_size;
     111             :   long num_inter, num_primes, bit_primes, max_forbidden;
     112             :   ulong flags;
     113             : /* from here: private data */
     114             :   GEN bc;
     115             :   ratpoints_sieve_entry *se_buffer;
     116             :   ratpoints_sieve_entry *se_next;
     117             :   ratpoints_bit_array *ba_buffer;
     118             :   ratpoints_bit_array *ba_next;
     119             :   int *int_buffer, *int_next;
     120             :   forbidden_entry *forb_ba;
     121             :   long *forbidden;
     122             :   GEN inverses, offsets, den_info, divisors;
     123             :   ulong **sieves0;
     124             : } ratpoints_args;
     125             : 
     126             : #ifdef HAS_SSE2
     127             : #define CODE_INIT_SIEVE_COPY for (a = 0; a < p; a++) si[a+p] = si[a];
     128             : #else
     129             : #define CODE_INIT_SIEVE_COPY
     130             : #endif
     131             : 
     132             : static ratpoints_bit_array *
     133     2470231 : sieve_init1(long p, ratpoints_sieve_entry *se1, long b1, ratpoints_args *args1)
     134             : {
     135     2470231 :   ratpoints_sieve_entry *se = se1;
     136     2470231 :   ratpoints_args *args = args1;
     137     2470231 :   int *isfs = se->is_f_square;
     138     2470231 :   long b = b1;
     139     2470231 :   long lmp = BITS_IN_LONG % p;
     140     2470231 :   long ldp = BITS_IN_LONG / p;
     141     2470231 :   long p1 = (ldp + 1) * p;
     142     2470231 :   long diff_shift = p1 & LONG_MASK;
     143     2470231 :   long diff = BITS_IN_LONG - diff_shift;
     144             :   ulong help0;
     145             :   long a;
     146     2470231 :   long d = se->inverses[b];
     147     2470231 :   long ab = 0; /* a/b mod p */
     148     2470231 :   ulong test = 1UL;
     149     2470231 :   ulong he0 = 0UL;
     150   110298124 :   for (a = 0; a < p; a++)
     151             :   {
     152   107827893 :     if (isfs[ab]) he0 |= test;
     153   107827893 :     ab += d;
     154   107827893 :     if (ab >= p) ab -= p;
     155   107827893 :     test <<= 1;
     156             :   }
     157     2470231 :   help0 = he0;
     158             :   {
     159             :     ulong help1;
     160             :      /* repeat bit pattern floor(BITS_IN_LONG/p) times */
     161     2470231 :     ulong pattern = help0;
     162             :     long i;
     163             :     /* the p * (floor(BITS_IN_LONG/p) + 1) - BITS_IN_LONG
     164             :             = p - (BITS_IN_LONG mod p)
     165             :        upper bits into help[b][1] :
     166             :        shift away the  BITS_IN_LONG mod p  lower bits */
     167     2470231 :     help1 = pattern >> lmp;
     168     5630383 :     for (i = p; i < BITS_IN_LONG; i <<= 1)
     169     3160152 :       help0 |= help0 << i;
     170             :     { /* fill the bit pattern from help0/help1 into sieve[b][].
     171             :           sieve[b][a0] has the same semantics as help0/help1,
     172             :           but here, a0 runs from 0 to p-1 and all bits are filled. */
     173             :       long a;
     174     2470231 :       ulong *si = (ulong *)args->ba_next;
     175             : 
     176     2470231 :       args->ba_next += p;
     177             :       /* copy the first chunk into sieve[b][] */
     178     2470231 :       si[0] = help0;
     179             :       /* now keep repeating the bit pattern,
     180             :          rotating it in help0/help1 */
     181   107827893 :       for (a = 1 ; a < p; a++)
     182             :       {
     183   105357662 :         ulong temp = help0 >> diff;
     184   105357662 :         help0 = help1 | (help0 << diff_shift);
     185   105357662 :         si[a] = help0;
     186   105357662 :         help1 = temp;
     187             :       }
     188   108088488 :       CODE_INIT_SIEVE_COPY
     189             :       /* set sieve array */
     190     2470231 :       se->sieve[b] = (ratpoints_bit_array *)si;
     191     2470231 :       return (ratpoints_bit_array *)si;
     192             :     }
     193             :   }
     194             : }
     195             : 
     196             : /* This is for p > BITS_IN_LONG */
     197             : static ratpoints_bit_array *
     198     7347583 : sieve_init2(long p, ratpoints_sieve_entry *se1, long b1, ratpoints_args *args1)
     199     7347583 : {
     200     7347583 :   ratpoints_sieve_entry *se = se1;
     201     7347583 :   ratpoints_args *args = args1;
     202     7347583 :   int *isfs = se->is_f_square;
     203     7347583 :   long b = b1;
     204             :   /* long ldp = 0;  = BITS_IN_LONG / p */
     205             :   /* long p1 = p; = (ldp + 1) * p; */
     206     7347583 :   long wp = p >> TWOPOTBITS_IN_LONG;
     207     7347583 :   long diff_shift = p & LONG_MASK;
     208     7347583 :   long diff = BITS_IN_LONG - diff_shift;
     209     7347583 :   ulong help[(p>>TWOPOTBITS_IN_LONG) + 2];
     210             : 
     211             :   /* initialize help */
     212             :   {
     213     7347583 :     ulong *he = &help[0];
     214     7347583 :     ulong *he1 = &he[(p>>TWOPOTBITS_IN_LONG) + 2];
     215    30948872 :     while (he1 != he) { he1--; *he1 = 0UL; }
     216             :   }
     217     7347583 :   { ulong work = 0UL;
     218             :     long a;
     219     7347583 :     long ab = 0; /* a/b mod p */
     220     7347583 :     long d = se->inverses[b];
     221     7347583 :     long n = 0;
     222     7347583 :     ulong test = 1UL;
     223   670719292 :     for (a = 0; a < p; )
     224             :     {
     225   663371709 :       if (isfs[ab]) work |= test;
     226   663371709 :       ab += d;
     227   663371709 :       if (ab >= p) ab -= p;
     228   663371709 :       test <<= 1;
     229   663371709 :       a++;
     230   663371709 :       if ((a & LONG_MASK) == 0)
     231     8906123 :       { help[n] = work; n++; work = 0UL; test = 1UL; }
     232             :     }
     233     7347583 :     help[n] = work;
     234             :   }
     235             : 
     236             :   { /* fill the bit pattern from help[] into sieve[b][].
     237             :        sieve[b][a0] has the same semantics as help[b][a0],
     238             :        but here, a0 runs from 0 to p-1 and all bits are filled. */
     239     7347583 :     ulong *si = (ulong *)args->ba_next;
     240             :     long a1;
     241             :     long a;
     242             : 
     243     7347583 :     args->ba_next += p;
     244             :     /* copy the first chunk from help[] into sieve[num][b][] */
     245    16253706 :     for (a = 0; a < wp; a++) si[a] = help[a];
     246             :     /* now keep repeating the bit pattern, rotating it in help */
     247   661813169 :     for (a1 = a ; a < p; a++)
     248             :     {
     249   654465586 :       long t = (a1 == wp) ? 0 : a1+1;
     250   654465586 :       help[a1] |= help[t]<<diff_shift;
     251   654465586 :       si[a] = help[a1];
     252   654465586 :       a1 = t;
     253   654465586 :       help[a1] >>= diff;
     254             :     }
     255   561329448 :      CODE_INIT_SIEVE_COPY
     256             :     /* set sieve array */
     257     7347583 :     se->sieve[b] = (ratpoints_bit_array *)si;
     258     7347583 :     return (ratpoints_bit_array *)si;
     259             :   }
     260             : }
     261             : 
     262             : static GEN
     263       12243 : gen_squares(GEN listprime)
     264             : {
     265       12243 :   long nbprime = lg(listprime)-1;
     266       12243 :   GEN sq = cgetg(nbprime+1, t_VEC);
     267             :   long n;
     268      379533 :   for (n = 1; n <= nbprime; n++)
     269             :   {
     270      367290 :     ulong i, p = uel(listprime,n);
     271      367290 :     GEN w = zero_zv(p), work = w+1;
     272      367290 :     work[0] = 1;
     273             :     /* record nonzero squares mod p, p odd */
     274    10700382 :     for (i = 1; i < p; i += 2) work[(i*i) % p] = 1;
     275      367290 :     gel(sq, n) = w;
     276             :   }
     277       12243 :   return sq;
     278             : }
     279             : 
     280             : static GEN
     281       12243 : gen_offsets(GEN P)
     282             : {
     283       12243 :   long n, l = lg(P);
     284       12243 :   GEN of = cgetg(l, t_VEC);
     285      379533 :   for (n = 1; n < l; n++)
     286             :   {
     287      367290 :     ulong p = uel(P,n);
     288      367290 :     uel(of, n) = Fl_inv((2*RBA_LENGTH)%p, p);
     289             :   }
     290       12243 :   return of;
     291             : }
     292             : 
     293             : static GEN
     294       12243 : gen_inverses(GEN P)
     295             : {
     296       12243 :   long n, l = lg(P);
     297       12243 :   GEN iv = cgetg(l, t_VEC);
     298      379533 :   for (n = 1; n < l; n++)
     299             :   {
     300      367290 :     ulong i, p = uel(P,n);
     301      367290 :     GEN w = cgetg(p, t_VECSMALL);
     302    21033474 :     for (i = 1; i < p; i++) uel(w, i) = Fl_inv(i, p);
     303      367290 :     gel(iv, n) = w;
     304             :   }
     305       12243 :   return iv;
     306             : }
     307             : 
     308             : static ulong **
     309       12243 : gen_sieves0(GEN listprime)
     310             : {
     311             :   long n;
     312       12243 :   long nbprime = lg(listprime)-1;
     313       12243 :   ulong ** si = (ulong**) new_chunk(nbprime+1);
     314      379533 :   for (n = 1; n <= nbprime; n++)
     315             :   {
     316      367290 :     ulong i, p = uel(listprime,n);
     317      367290 :     ulong *w = (ulong *) stack_malloc_align(2*p*sizeof(ulong), RBA_ALIGN);
     318    21400764 :     for (i = 0; i < p; i++) uel(w,i) = ~0UL;
     319    22194810 :     for (i = 0; i < BITS_IN_LONG; i++)
     320    21827520 :       uel(w,(p*i)>>TWOPOTBITS_IN_LONG) &= ~(1UL<<((p*i) & LONG_MASK));
     321    21400764 :     for (i = 0; i < p; i++) uel(w,i+p) = uel(w,i);
     322      367290 :     si[n] = w;
     323             :   }
     324       12243 :   return si;
     325             : }
     326             : 
     327             : static void
     328       12243 : gen_sieve(ratpoints_args *args)
     329             : {
     330       12243 :   GEN listprimes = args->listprime;
     331       12243 :   args->offsets = gen_offsets(listprimes);
     332       12243 :   args->inverses = gen_inverses(listprimes);
     333       12243 :   args->sieves0 = gen_sieves0(listprimes);
     334       12243 : }
     335             : 
     336             : static GEN
     337       12243 : ZX_positive_region(GEN P, long h, long bitprec)
     338             : {
     339       12243 :   long prec = nbits2prec(bitprec);
     340       12243 :   GEN it = mkvec2(stoi(-h),stoi(h));
     341       12243 :   GEN R = ZX_Uspensky(P, it, 1, bitprec);
     342       12243 :   long nR = lg(R)-1;
     343       12243 :   long s = signe(ZX_Z_eval(P,gel(it,1)));
     344       12243 :   long i=1, j;
     345             :   GEN iv, st, en;
     346       12243 :   if (s<0 && nR==0) return NULL;
     347       11578 :   iv = cgetg(((nR+1+(s>=0))>>1)+1, t_VEC);
     348       11578 :   if (s>=0) st = itor(gel(it,1),prec);
     349        5025 :   else    { st = gel(R,i); i++; }
     350       18019 :   for (j=1; i<nR; j++)
     351             :   {
     352        6441 :     gel(iv, j) = mkvec2(st, gel(R,i));
     353        6441 :     st = gel(R,i+1);
     354        6441 :     i+=2;
     355             :   }
     356       11578 :   if (i==nR) en = gel(R,i); else en = itor(gel(it,2),prec);
     357       11578 :   gel(iv,j) = mkvec2(st, en);
     358       11578 :   return iv;
     359             : }
     360             : 
     361             : static long
     362       12243 : posint(ratpoints_interval *ivlist, GEN P, long h)
     363             : {
     364       12243 :   GEN R = ZX_positive_region(P, h, 53);
     365       12243 :   const double eps = 1e-5;
     366             :   long nR, i;
     367             : 
     368       12243 :   if (!R) return 0;
     369       11578 :   nR = lg(R)-1;
     370       11578 :   i = 1;
     371       29597 :   for (i=1; i<=nR; i++)
     372             :   {
     373       18019 :     ivlist[i-1].low = rtodbl(gmael(R,i,1))-eps;
     374       18019 :     ivlist[i-1].up  = rtodbl(gmael(R,i,2))+eps;
     375             :   }
     376       11578 :   return nR;
     377             : }
     378             : 
     379             : static long
     380       12243 : ratpoints_compute_sturm(ratpoints_args *args)
     381             : {
     382       12243 :   ratpoints_interval *ivlist = args->domain;
     383       12243 :   args->num_inter = posint(ivlist, args->cof, (long) ivlist[0].up);
     384       12243 :   return args->num_inter;
     385             : }
     386             : 
     387             : /**************************************************************************
     388             :  * Try to avoid divisions                                                 *
     389             :  **************************************************************************/
     390             : INLINE long
     391  1636302303 : mod(long a, long b)
     392             : {
     393  1636302303 :   long b1 = b << 4; /* b1 = 16*b */
     394             : 
     395  1636302303 :   if (a < -b1) { a %= b; if (a < 0) { a += b; } return a ; }
     396  1630797127 :   if (a < 0) { a += b1; }
     397   991029655 :   else { if (a >= b1) { return a % b; } }
     398  1620420501 :   b1 >>= 1; /* b1 = 8*b */
     399  1620420501 :   if (a >= b1) { a -= b1; }
     400  1620420501 :   b1 >>= 1; /* b1 = 4*b */
     401  1620420501 :   if (a >= b1) { a -= b1; }
     402  1620420501 :   b1 >>= 1; /* b1 = 2*b */
     403  1620420501 :   if (a >= b1) { a -= b1; }
     404  1620420501 :   if (a >= b) { a -= b; }
     405  1620420501 :   return a;
     406             : }
     407             : 
     408             : static void
     409     6082133 : set_bc(long b, ratpoints_args *args)
     410             : {
     411     6082133 :   GEN w0 = gen_1;
     412     6082133 :   GEN c = args->cof, bc;
     413     6082133 :   long k, degree = degpol(c);
     414     6082133 :   bc = cgetg(degree+2, t_POL);
     415    31716627 :   for (k = degree-1; k >= 0; k--)
     416             :   {
     417    25634494 :     w0 = muliu(w0, b);
     418    25634494 :     gel(bc,k+2) = mulii(gel(c,k+2), w0);
     419             :   }
     420     6082133 :   args->bc = bc;
     421     6082133 : }
     422             : 
     423             : /**************************************************************************
     424             :  * Check a `survivor' of the sieve if it really gives a point.            *
     425             :  **************************************************************************/
     426             : 
     427             : static long
     428    12915871 : _ratpoints_check_point(long a, long b, ratpoints_args *args, int *quit,
     429             :                  int process(long, long, GEN, void*, int*), void *info)
     430             : {
     431    12915871 :   pari_sp av = avma;
     432    12915871 :   GEN w0, w2, c = args->cof, bc = args->bc;
     433    12915871 :   long k, degree = degpol(c);
     434    12915871 :   int reverse = args->flags & RATPOINTS_REVERSED;
     435             : 
     436             :   /* Compute F(a, b), where F is the homogenized version of f
     437             :      of smallest possible even degree  */
     438    12915871 :   w2 = gel(c, degree+2);
     439    67053050 :   for (k = degree-1; k >= 0; k--)
     440             :   {
     441    54137179 :     w2 = mulis(w2, a);
     442    54137179 :     w2 = addii(w2, gel(bc,k+2));
     443             :   }
     444    12915871 :   if (odd(degree)) w2 = muliu(w2, b);
     445             :   /* check if f(x,z) is a square; if so, process the point(s) */
     446    12915871 :   if (signe(w2) >= 0 && Z_issquareall(w2, &w0))
     447             :   {
     448       44576 :     if (reverse)
     449             :     {
     450        1218 :       if (a >= 0) (void)process(b, a, w0, info, quit);
     451         217 :       else        (void)process(-b, -a, w0, info, quit);
     452             :     }
     453       43358 :     else (void)process(a, b, w0, info, quit);
     454       44576 :     if (!*quit && signe(w0) != 0)
     455             :     {
     456       42329 :       GEN nw0 = negi(w0);
     457       42329 :       if (reverse)
     458             :       {
     459        1155 :         if (a >= 0) (void)process(b, a, nw0, info, quit);
     460         196 :         else        (void)process(-b, -a, nw0, info, quit);
     461             :       }
     462       41174 :       else (void)process(a, b, nw0, info, quit);
     463             :     }
     464       44576 :     return 1;
     465             :   }
     466    12871295 :   set_avma(av);
     467    12871295 :   return 0;
     468             : }
     469             : 
     470             : /**************************************************************************
     471             :  * The inner loop of the sieving procedure                                *
     472             :  **************************************************************************/
     473             : static long
     474    44636132 : _ratpoints_sift0(long b, long w_low, long w_high,
     475             :            ratpoints_args *args, bit_selection which_bits,
     476             :            ratpoints_bit_array *survivors, sieve_spec *sieves, int *quit,
     477             :            int process(long, long, GEN, void*, int*), void *info)
     478             : {
     479    44636132 :   long range = w_high - w_low;
     480    44636132 :   long sp1 = args->sp1;
     481    44636132 :   long sp2 = args->sp2;
     482    44636132 :   long i, n, nb = 0, absb = labs(b);
     483             :   ratpoints_bit_array *surv0;
     484             : 
     485             :   /* now do the sieving (fast!) */
     486   446361320 :   for (n = 0; n < sp1; n++)
     487             :   {
     488   401725188 :     ratpoints_bit_array *sieve_n = sieves[n].ptr;
     489   401725188 :     long p = sieves[n].p;
     490   401725188 :     long r = mod(-w_low-sieves[n].offset, p);
     491   401725188 :     ratpoints_bit_array *surv = survivors;
     492             : 
     493   401725188 :     if (w_high < w_low + r)
     494             :     { /* if we get here, r > 0, since w_high >= w_low always */
     495    64055035 :       ratpoints_bit_array *siv1 = &sieve_n[p-r];
     496    64055035 :       ratpoints_bit_array *siv0 = siv1 + range;
     497             : 
     498  1943216567 :       while (siv1 != siv0) { *surv = AND(*surv, *siv1++); surv++; }
     499             :     }
     500             :     else
     501             :     {
     502   337670153 :       ratpoints_bit_array *siv1 = &sieve_n[p-r];
     503   337670153 :       ratpoints_bit_array *surv_end = &survivors[range - p];
     504             :       long i;
     505  8803867120 :       for (i = r; i; i--) { *surv = AND(*surv, *siv1++); surv++; }
     506   337670153 :       siv1 -= p;
     507  1265024160 :       while (surv <= surv_end)
     508             :       {
     509 25035530538 :         for (i = p; i; i--) { *surv = AND(*surv, *siv1++); surv++; }
     510   927354007 :         siv1 -= p;
     511             :       }
     512   337670153 :       surv_end += p;
     513  8894518722 :       while (surv < surv_end) { *surv = AND(*surv, *siv1++); surv++; }
     514             :     }
     515             :   } /* for n */
     516             :   /* 2nd phase of the sieve: test each surviving bit array with more primes */
     517    44636132 :   surv0 = &survivors[0];
     518  4823487783 :   for (i = w_low; i < w_high; i++)
     519             :   {
     520  4778853408 :     ratpoints_bit_array nums = *surv0++;
     521  4778853408 :     sieve_spec *ssp = &sieves[sp1];
     522             :     long n;
     523             : 
     524  6009405170 :     for (n = sp2-sp1; n && TEST(nums); n--)
     525             :     {
     526  1230551762 :       long p = ssp->p;
     527  1230551762 :       nums = AND(nums, ssp->ptr[mod(i + ssp->offset, p)]);
     528  1230551762 :       ssp++;
     529             :     }
     530             : 
     531             :     /* Check the survivors of the sieve if they really give points */
     532  4778853408 :     if (TEST(nums))
     533             :     {
     534             :       long a0, a, d;
     535             : #ifdef HAS_SSE2
     536    19255248 :       long da = (which_bits == num_all) ? RBA_LENGTH/2 : RBA_LENGTH;
     537             : #endif
     538             :            /* a will be the numerator corresponding to the selected bit */
     539    22489112 :       if (which_bits == num_all)
     540             :       {
     541    13605455 :         d = 1; a0 = i * RBA_LENGTH;
     542             :       }
     543             :       else
     544             :       {
     545     8883657 :         d = 2; a0 = i * 2 * RBA_LENGTH;
     546     8883657 :         if (which_bits == num_odd) a0++;
     547             :       }
     548             :       {
     549             : #ifdef HAS_SSE2
     550    19255248 :         ulong nums0 = EXT0(nums);
     551    19255248 :         ulong nums1 = EXT1(nums);
     552             : #else
     553     3233864 :         ulong nums0 = nums;
     554             : #endif
     555             : 
     556   336384062 :         for (a = a0; nums0; a += d, nums0 >>= 1)
     557             :         { /* test one bit */
     558   313895963 :           if (odd(nums0) && ugcd(labs(a), absb)==1)
     559             :           {
     560     7379899 :             if (!args->bc) set_bc(b, args);
     561     7379899 :             nb += _ratpoints_check_point(a, b, args, quit, process, info);
     562     7379899 :             if (*quit) return nb;
     563             :           }
     564             :         }
     565             : #ifdef HAS_SSE2
     566   285338316 :         for (a = a0 + da; nums1; a += d, nums1 >>= 1)
     567             :         { /* test one bit */
     568   266084574 :           if (odd(nums1) && ugcd(labs(a), absb)==1)
     569             :           {
     570     5535972 :             if (!args->bc) set_bc(b, args);
     571     5535972 :             nb += _ratpoints_check_point(a, b, args, quit, process, info);
     572     5535972 :             if (*quit) return nb;
     573             :           }
     574             :         }
     575             : #endif
     576             :       }
     577             :     }
     578             :   }
     579    44634375 :   return nb;
     580             : }
     581             : 
     582             : typedef struct { double r; ratpoints_sieve_entry *ssp; } entry;
     583             : 
     584             : static const int squares16[16] = {1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0};
     585             :  /* Says if a is a square mod 16, for a = 0..15 */
     586             : 
     587             : /**************************************************************************
     588             :  * Initialization and cleanup of ratpoints_args structure                 *
     589             :  **************************************************************************/
     590             : 
     591             : static ratpoints_sieve_entry*
     592       12243 : alloc_sieve(long nbprime, long maxprime)
     593             : {
     594             :   long i;
     595             :   ratpoints_sieve_entry * s = (ratpoints_sieve_entry*)
     596       12243 :                         stack_malloc(nbprime*sizeof(ratpoints_sieve_entry));
     597      379533 :   for (i=0; i<nbprime; i++)
     598      367290 :     s[i].sieve = (ratpoints_bit_array**) new_chunk(maxprime);
     599       12243 :   return s;
     600             : }
     601             : 
     602             : /* NOTE: args->degree must be set */
     603             : static void
     604       12243 : find_points_init(ratpoints_args *args, long bit_primes)
     605             : {
     606       12243 :   long need = 0;
     607             :   long n, nbprime,maxprime;
     608       12243 :   args->listprime = primes_interval_zv(3, 1<<bit_primes);
     609       12243 :   nbprime = lg(args->listprime)-1;
     610       12243 :   maxprime = args->listprime[nbprime];
     611             : 
     612             :   /* allocate space for se_buffer */
     613       12243 :   args->se_buffer = alloc_sieve(nbprime, maxprime);
     614       12243 :   args->se_next = args->se_buffer;
     615      379533 :   for (n = 1; n <= nbprime; n++)
     616             :   {
     617      367290 :     ulong p = args->listprime[n];
     618      367290 :     need += p*p;
     619             :   }
     620       12243 :   args->ba_buffer = (ratpoints_bit_array*)
     621       12243 :      stack_malloc_align(need*sizeof(ratpoints_bit_array),RBA_ALIGN);
     622       12243 :   args->ba_next = args->ba_buffer;
     623             : 
     624             :   /* allocate space for int_buffer */
     625       12243 :   args->int_buffer = (int *) stack_malloc(nbprime*(maxprime+1)*sizeof(int));
     626       12243 :   args->int_next = args->int_buffer;
     627             : 
     628       12243 :   args->forb_ba   = (forbidden_entry*)
     629       12243 :     stack_malloc((nbprime + 1)*sizeof(forbidden_entry));
     630       12243 :   args->forbidden = new_chunk(nbprime + 1);
     631       12243 :   gen_sieve(args);
     632       12243 :   return;
     633             : }
     634             : 
     635             : /* f = leading coeff; b = b1*b2, b1 maximal with (b1, 2*f) = 1;
     636             :  * return Jacobi symbol (f, b1) */
     637             : INLINE int
     638    49874164 : rpjacobi(long b, GEN lcf)
     639             : {
     640             :   ulong f;
     641    49874164 :   b >>= vals(b);
     642    49874164 :   f = umodiu(lcf, b);
     643    49874164 :   return krouu(f, u_ppo(b,f));
     644             : }
     645             : 
     646             : /************************************************************************
     647             :  * Set up information on possible denominators                          *
     648             :  * when polynomial is of odd degree with leading coefficient != +-1     *
     649             :  ************************************************************************/
     650             : 
     651             : static void
     652        1309 : setup_us1(ratpoints_args *args, GEN w0)
     653             : {
     654        1309 :   GEN F = Z_issmooth_fact(w0, 1000), P, E, S, D;
     655             :   long i, l;
     656             : 
     657        1309 :   if (!F) return;
     658        1309 :   P = gel(F,1); l = lg(P);
     659        1309 :   E = gel(F,2);
     660        1309 :   D  = cgetg(1+(1<<(l-1)), t_VECSMALL);
     661             :   /* factorization is complete, set up array of squarefree divisors */
     662        1309 :   D[1] = 1;
     663        2744 :   for (i = 1; i < l; i++)
     664             :   { /* multiply all divisors known so far by next prime */
     665        1435 :     long k, n = 1<<(i-1);
     666        2996 :     for (k=0; k<n; k++) uel(D,1+n+k) = uel(D,1+k) * P[i];
     667             :   }
     668        1309 :   S = cgetg(l, t_VECSMALL);
     669             :   /* set slopes in den_info */
     670        2744 :   for (i = 1; i < l; i++)
     671             :   { /* compute min{n : (d-k)*n > v_p(f_d) - v_p(f_k), k = 0,...,d-1} */
     672        1435 :     GEN c = args->cof;
     673        1435 :     long p = P[i], v = E[i];
     674        1435 :     long k, n = 1, d = degpol(c);
     675             : 
     676        6874 :     for (k = d - 1; k >= 0; k--)
     677             :     {
     678        5439 :       long t = 1 + v - Z_lval(gel(c,k+2), p);
     679        5439 :       long m = CEIL(t, d - k);
     680             : 
     681        5439 :       if (m > n) n = m;
     682             :     }
     683        1435 :     S[i] = n;
     684             :   }
     685        1309 :   args->divisors = D;
     686        1309 :   args->flags |= RATPOINTS_USE_SQUARES1;
     687        1309 :   args->den_info = mkvec3(P, E, S);
     688             : }
     689             : 
     690             : /************************************************************************
     691             :  * Consider 2-adic information                                          *
     692             :  ************************************************************************/
     693             : 
     694             : static bit_selection
     695       12243 : get_2adic_info(ratpoints_args *args, ulong *den_bits,
     696             :                ratpoints_bit_array *num_bits)
     697             : {
     698       12243 :   GEN c = args->cof;
     699       12243 :   long degree = degpol(c);
     700             :   int is_f_square16[24];
     701       12243 :   long *cmp = new_chunk(degree+1);
     702       12243 :   long npe = 0, npo = 0;
     703             :   bit_selection result;
     704             :   long n, a, b;
     705             : 
     706             :   /* compute coefficients mod 16 */
     707       86114 :   for (n = 0; n <= degree; n++) cmp[n] = Mod16(gel(c,n+2));
     708      208131 :   for (a = 0 ; a < 16; a++)
     709             :   {
     710      195888 :     ulong s = cmp[degree];
     711             :     long n;
     712     1181936 :     for (n = degree - 1 ; n >= 0 ; n--) s = s*a + cmp[n];
     713      195888 :     s &= 0xf;
     714      195888 :     if ((is_f_square16[a] = squares16[s])) { if (odd(a)) npo++; else npe++; }
     715             :   }
     716             : 
     717             :   /* even denominators:
     718             :      is_f_square16[16+k] says if f((2k+1)/2) is a square, k = 0..3
     719             :      is_f_square16[20+k] says if f((2k+1)/4) is a square, k = 0,1
     720             :      is_f_square16[22]   says if f(odd/8) is a square
     721             :      is_f_square16[23]   says if f(odd/2^n), n >= 4, can be a square */
     722             :   {
     723       12243 :     long np1 = 0, np2 = 0, np3 = 0, np4 = 0;
     724             : 
     725       12243 :     if (odd(degree))
     726             :     {
     727        1372 :       long a, cf = 4*cmp[degree-1];
     728             : 
     729        1372 :       if (degree >= 2) cf += 8*cmp[degree-2];
     730        6860 :       for (a = 0; a < 4; a++)
     731             :       { /* Compute  2 c[d] k^d + 4 c[d-1] k^(d-1) + 8 c[d-2] k^(d-2), k = 2a+1.
     732             :            Note that k^d = k mod 8, k^(d-1) = 1 mod 8. */
     733        5488 :         long k = 2*a+1;
     734        5488 :         long s = (2*k*cmp[degree] + cf) & 0xf;
     735        5488 :         if ((is_f_square16[16+a] = squares16[s])) np1++;
     736             :       }
     737        1372 :       if ((is_f_square16[20] = squares16[(4*cmp[degree])  & 0xf])) np2++;
     738        1372 :       if ((is_f_square16[21] = squares16[(12*cmp[degree]) & 0xf])) np2++;
     739        1372 :       if ((is_f_square16[22] = squares16[(8*cmp[degree])  & 0xf])) np3++;
     740        1372 :       is_f_square16[23] = 1; np4++;
     741             :     }
     742             :     else
     743             :     {
     744       10871 :       long a, cf = (degree >= 2) ? 4*cmp[degree-2] : 0;
     745             : 
     746       10871 :       if (degree >= 3) cf += 8*cmp[degree-3];
     747       54355 :       for (a = 0; a < 4; a++)
     748             :       { /* compute c[d] k^d + 2 c[d-1] k^(d-1) + ... + 8 c[d-3] k^(d-3),
     749             :            k = 2a+1. Note that k^d = k^2 mod 16, k^(d-1) = k mod 8. */
     750       43484 :         long k = 2*a+1;
     751       43484 :         long s = ((cmp[degree]*k + 2*cmp[degree-1])*k + cf) & 0xf;
     752       43484 :         if ((is_f_square16[16+a] = squares16[s])) np1++;
     753             :       }
     754       10871 :       if ((is_f_square16[20] = squares16[(cmp[degree]+4*cmp[degree-1])  & 0xf]))
     755        4473 :         np2++;
     756       10871 :       if ((is_f_square16[21] = squares16[(cmp[degree]+12*cmp[degree-1]) & 0xf]))
     757        4466 :         np2++;
     758       10871 :       if ((is_f_square16[22] = squares16[(cmp[degree]+8*cmp[degree-1])  & 0xf]))
     759        4270 :         np3++;
     760       10871 :       if ((is_f_square16[23] = squares16[cmp[degree]])) np4++;
     761             :     }
     762             : 
     763             :     /* set den_bits */
     764       12243 :     { ulong db = 0;
     765             :       long i;
     766             : 
     767       12243 :       if (npe + npo > 0) db |= 0xaaaaUL; /* odd denominators */
     768       12243 :       if (np1 > 0)       db |= 0x4444UL; /* v_2(den) = 1 */
     769       12243 :       if (np2 > 0)       db |= 0x1010UL; /* v_2(den) = 2 */
     770       12243 :       if (np3 > 0)       db |= 0x0100UL; /* v_2(den) = 3 */
     771       12243 :       if (np4 > 0)       db |= 0x0001UL; /* v_2(den) >= 4 */
     772       12243 :       if (db == 0) { *den_bits = 0UL; return num_none; }
     773             : 
     774       34480 :       for (i = 16; i < BITS_IN_LONG; i <<= 1) db |= db << i;
     775       12068 :       *den_bits = db;
     776             :     }
     777       12068 :     result = (npe == 0) ? (npo == 0) ? num_none : num_odd
     778       12068 :                         : (npo == 0) ? num_even : num_all;
     779             :   }
     780             : 
     781             :   /* set up num_bits[16] */
     782             : 
     783             :   /* odd denominators */
     784       12068 :   switch(result)
     785             :   {
     786        7609 :     case num_all:
     787       68481 :       for (b = 1; b < 16; b += 2)
     788             :       {
     789       60872 :         ulong work = 0, bit = 1;
     790       60872 :         long i, invb = b; /* inverse of b mod 16 */
     791       60872 :         if (b & 2) invb ^= 8;
     792       60872 :         if (b & 4) invb ^= 8;
     793     1034824 :         for (i = 0; i < 16; i++)
     794             :         {
     795      973952 :           if (is_f_square16[(invb*i) & 0xf]) work |= bit;
     796      973952 :           bit <<= 1;
     797             :         }
     798             :         /* now repeat the 16 bits */
     799      173920 :         for (i = 16; i < BITS_IN_LONG; i <<= 1) work |= work << i;
     800       60872 :         num_bits[b] = RBA(work);
     801             :       }
     802        7609 :       break;
     803             : 
     804        1960 :     case num_odd:
     805       17640 :       for (b = 1; b < 16; b += 2)
     806             :       {
     807       15680 :         ulong work = 0, bit = 1;
     808       15680 :         long i, invb = b; /* inverse of b mod 16 */
     809       15680 :         if (b & 2) invb ^= 8;
     810       15680 :         if (b & 4) invb ^= 8;
     811      141120 :         for (i = 1; i < 16; i += 2)
     812             :         {
     813      125440 :           if (is_f_square16[(invb*i) & 0xf]) work |= bit;
     814      125440 :           bit <<= 1;
     815             :         }
     816             :         /* now repeat the 8 bits */
     817       60480 :         for (i = 8; i < BITS_IN_LONG; i <<= 1) { work |= work << i; }
     818       15680 :         num_bits[b] = RBA(work);
     819             :       }
     820        1960 :       break;
     821             : 
     822        2114 :     case num_even:
     823       19026 :       for (b = 1; b < 16; b += 2)
     824             :       {
     825       16912 :         ulong work = 0, bit = 1;
     826       16912 :         long i, invb = b; /* inverse of b mod 16 */
     827       16912 :         if (b & 2) invb ^= 8;
     828       16912 :         if (b & 4) invb ^= 8;
     829      152208 :         for (i = 0; i < 16; i += 2)
     830             :         {
     831      135296 :           if (is_f_square16[(invb*i) & 0xf]) work |= bit;
     832      135296 :           bit <<= 1;
     833             :         }
     834             :         /* now repeat the 8 bits */
     835       65232 :         for (i = 8; i < BITS_IN_LONG; i <<= 1) work |= work << i;
     836       16912 :         num_bits[b] = RBA(work);
     837             :       }
     838        2114 :       break;
     839             : 
     840         385 :     case num_none:
     841        3465 :       for (b = 1; b < 16; b += 2) num_bits[b] = RBA(0UL);
     842         385 :       break;
     843             :   }
     844             : 
     845             :   /* v_2(den) = 1 : only odd numerators */
     846       60340 :   for (b = 1; b < 8; b += 2)
     847             :   {
     848       48272 :     ulong work = 0, bit = 1;
     849             :     long i;
     850      434448 :     for (i = 1; i < 16; i += 2)
     851             :     {
     852      386176 :       if (is_f_square16[16 + (((b*i)>>1) & 0x3)]) work |= bit;
     853      386176 :       bit <<= 1;
     854             :     }
     855             :     /* now repeat the 8 bits */
     856      186192 :     for (i = 8; i < BITS_IN_LONG; i <<= 1) work |= work << i;
     857       48272 :     num_bits[2*b] = RBA(work);
     858             :   }
     859             : 
     860             :   /* v_2(den) = 2 : only odd numerators */
     861       36204 :   for (b = 1; b < 4; b += 2)
     862             :   {
     863       24136 :     ulong work = 0, bit = 1;
     864             :     long i;
     865      120680 :     for (i = 1; i < 8; i += 2)
     866             :     {
     867       96544 :       if (is_f_square16[20 + (((b*i)>>1) & 0x1)]) work |= bit;
     868       96544 :       bit <<= 1;
     869             :     }
     870             :     /* now repeat the 4 bits */
     871      117232 :     for (i = 4; i < BITS_IN_LONG; i <<= 1) work |= work << i;
     872       24136 :     num_bits[4*b] = RBA(work);
     873             :   }
     874             : 
     875             :   /* v_2(den) = 3, >= 4 : only odd numerators */
     876       12068 :   num_bits[8] = (is_f_square16[22]) ? RBA(~(0UL)) : RBA(0UL);
     877       12068 :   num_bits[0] = (is_f_square16[23]) ? RBA(~(0UL)) : RBA(0UL);
     878             : 
     879       12068 :   return result;
     880             : }
     881             : 
     882             : /**************************************************************************
     883             :  * This is a comparison function needed for sorting in order to determine *
     884             :  * the `best' primes for sieving.                                         *
     885             :  **************************************************************************/
     886             : 
     887             : static int
     888     1037321 : compare_entries(const void *a, const void *b)
     889             : {
     890     1037321 :   double diff = ((entry *)a)->r - ((entry *)b)->r;
     891     1037321 :   return (diff > 0) ? 1 : (diff < 0) ? -1 : 0;
     892             : }
     893             : 
     894             : /************************************************************************
     895             :  * Collect the sieving information                                      *
     896             :  ************************************************************************/
     897             : 
     898             : static long
     899       12243 : sieving_info(ratpoints_args *args,
     900             :              ratpoints_sieve_entry **sieve_list)
     901             : {
     902       12243 :   GEN c = args->cof;
     903       12243 :   GEN inverses = args->inverses, squares;
     904       12243 :   GEN offsets = args->offsets;
     905       12243 :   ulong ** sieves0 = args->sieves0;
     906       12243 :   long degree = degpol(c);
     907       12243 :   long fba = 0, fdc = 0;
     908       12243 :   long pn, pnp = 0;
     909             :   long n;
     910       12243 :   long nbprime = lg(args->listprime)-1;
     911       12243 :   entry *prec = (entry*) stack_malloc(nbprime*sizeof(entry));
     912             :     /* This array is used for sorting in order to
     913             :        determine the `best' sieving primes. */
     914             : 
     915       12243 :   forbidden_entry *forb_ba = args->forb_ba;
     916       12243 :   long *forbidden = args->forbidden;
     917       12243 :   ulong bound = (1UL)<<(BITS_IN_LONG - args->bit_primes);
     918       12243 :   pari_sp av = avma;
     919       12243 :   squares = gen_squares(args->listprime);
     920             : 
     921             :   /* initialize sieve in se_buffer */
     922      351568 :   for (pn = 1; pn <= args->num_primes; pn++)
     923      339451 :   {
     924      339451 :     long coeffs_mod_p[degree+1]; /* The coefficients of f reduced modulo p */
     925      339451 :     ulong a, p = args->listprime[pn], np;
     926             :     long n;
     927      339451 :     int *is_f_square = args->int_next;
     928             : 
     929      339451 :     args->int_next += p + 1; /* need space for (p+1) int's */
     930             : 
     931             :     /* compute coefficients mod p */
     932     2384368 :     for (n = 0; n <= degree; n++) coeffs_mod_p[n] = umodiu(pel(c,n), p);
     933             : 
     934      339451 :     np = umael(squares,pn,coeffs_mod_p[0]+1);
     935      339451 :     is_f_square[0] = np;
     936    17909577 :     for (a = 1 ; a < p; a++)
     937             :     {
     938    17570126 :       ulong s = coeffs_mod_p[degree];
     939    17570126 :       if ((degree+1)*args->bit_primes <= BITS_IN_LONG)
     940             :       {
     941    91441306 :         for (n = degree - 1 ; n >= 0 ; n--) s = s*a + coeffs_mod_p[n];
     942             :         /* here, s < p^(degree+1) <= max. long */
     943    15241358 :         s %= p;
     944             :       }
     945             :       else
     946             :       {
     947    14396076 :         for (n = degree - 1 ; n >= 0 ; n--)
     948             :         {
     949    12067308 :           s = s*a + coeffs_mod_p[n];
     950    12067308 :           if (s+1 >= bound) s %= p;
     951             :         }
     952     2328768 :         s %= p;
     953             :       }
     954    17570126 :       if ((is_f_square[a] = mael(squares,pn,s+1))) np++;
     955             :     }
     956      339451 :     is_f_square[p] = odd(degree) || mael(squares,pn,coeffs_mod_p[degree]+1);
     957             : 
     958             :     /* check if there are no solutions mod p */
     959      339451 :     if (np == 0 && !is_f_square[p]) return gc_long(av,p);
     960             : 
     961             :     /* Fill arrays with info for p */
     962      339325 :     if (np < p)
     963             :     { /* only when there is some information */
     964             :       ulong i;
     965      307412 :       ratpoints_sieve_entry *se = args->se_next;
     966      788421 :       double r = is_f_square[p] ? ((double)(np*(p-1) + p))/((double)(p*p))
     967      307412 :                                 : (double)np/(double)p;
     968      307412 :       prec[pnp].r = r;
     969      307412 :       args->se_next ++;
     970      307412 :       se->p = p;
     971      307412 :       se->is_f_square = is_f_square;
     972      307412 :       se->inverses = gel(inverses,pn);
     973      307412 :       se->offset = offsets[pn];
     974      307412 :       se->sieve[0] = (ratpoints_bit_array *)sieves0[pn];
     975    17056018 :       for (i = 1; i < p; i++) se->sieve[i] = NULL;
     976      307412 :       prec[pnp].ssp = se;
     977      307412 :       pnp++;
     978             :     }
     979             : 
     980      339325 :     if ((args->flags & RATPOINTS_CHECK_DENOM)
     981      296793 :          && fba + fdc < args->max_forbidden
     982      296793 :          && !is_f_square[p])
     983             :     { /* record forbidden divisors of the denominator */
     984      135208 :       if (coeffs_mod_p[degree] == 0)
     985             :       { /* leading coeff. divisible by p */
     986             :         GEN r;
     987           0 :         long v = Z_lvalrem(pel(c,degree), p, &r);
     988             : 
     989           0 :         if (odd(v) || !mael(squares,pn, umodiu(r,p)+1))
     990             :         { /* Can only get something when valuation is odd
     991             :              or when valuation is even and lcf is not a p-adic square.
     992             :              Compute smallest n such that if v(den) >= n, the leading
     993             :              term determines the valuation. Then we must have v(den) < n. */
     994           0 :           long k, n = 1;
     995           0 :           for (k = degree-1; k >= 0; k--)
     996             :           {
     997           0 :             if (coeffs_mod_p[k] == 0)
     998             :             {
     999           0 :               long t = 1 + v - Z_lval(pel(c,k), p);
    1000           0 :               long m = CEIL(t, (degree-k));
    1001           0 :               if (m > n) n = m;
    1002             :             }
    1003             :           }
    1004           0 :           if (n == 1)
    1005             :           {
    1006           0 :             forb_ba[fba].p     = p;
    1007           0 :             forb_ba[fba].start = sieves0[pn];
    1008           0 :             forb_ba[fba].end   = sieves0[pn]+p;
    1009           0 :             forb_ba[fba].curr  = forb_ba[fba].start;
    1010           0 :             fba++;
    1011             :           }
    1012             :           else
    1013             :           {
    1014           0 :             forbidden[fdc] = upowuu(p, n);
    1015           0 :             fdc++;
    1016             :           }
    1017             :         }
    1018             :       }
    1019             :       else /* leading coefficient is a nonsquare mod p */
    1020             :       { /* denominator divisible by p is excluded */
    1021      135208 :         forb_ba[fba].p     = p;
    1022      135208 :         forb_ba[fba].start = sieves0[pn];
    1023      135208 :         forb_ba[fba].end   = sieves0[pn]+p;
    1024      135208 :         forb_ba[fba].curr  = forb_ba[fba].start;
    1025      135208 :         fba++;
    1026             :       }
    1027             :     }
    1028             :   } /* end for pn */
    1029             : 
    1030             :   /* update sp2 and sp1 if necessary */
    1031       12117 :   if (args->sp2 > pnp)       args->sp2 = pnp;
    1032       12117 :   if (args->sp1 > args->sp2) args->sp1 = args->sp2;
    1033             : 
    1034             :   /* sort the array to get at the best primes */
    1035       12117 :   qsort(prec, pnp, sizeof(entry), compare_entries);
    1036             : 
    1037             :   /* put the sorted entries into sieve_list */
    1038      205835 :   for (n = 0; n < args->sp2; n++) sieve_list[n] = prec[n].ssp;
    1039             : 
    1040             :   /* terminate array of forbidden divisors */
    1041       12117 :   if (args->flags & RATPOINTS_CHECK_DENOM)
    1042             :   {
    1043             :     long n;
    1044             : 
    1045       10598 :     for (n = args->num_primes+1;
    1046       24052 :         fba + fdc < args->max_forbidden && n <= nbprime; n++)
    1047             :     {
    1048       17325 :       ulong p = args->listprime[n];
    1049             : 
    1050       17325 :       if (p*p > (ulong) args->b_high) break;
    1051       13454 :       if (kroiu(pel(c,degree), p) == -1)
    1052             :       {
    1053        6937 :         forb_ba[fba].p     = p;
    1054        6937 :         forb_ba[fba].start = sieves0[n];
    1055        6937 :         forb_ba[fba].end   = sieves0[n]+p;
    1056        6937 :         forb_ba[fba].curr  = forb_ba[fba].start;
    1057        6937 :         fba++;
    1058             :       }
    1059             :     }
    1060       10598 :     forb_ba[fba].p = 0; /* terminating zero */
    1061       10598 :     forbidden[fdc] = 0; /* terminating zero */
    1062       10598 :     args->max_forbidden = fba + fdc; /* note actual number */
    1063             :   }
    1064             : 
    1065       12117 :   if (fba + fdc == 0) args->flags &= ~RATPOINTS_CHECK_DENOM;
    1066       12117 :   return gc_long(av,0);
    1067             : }
    1068             : 
    1069             : /**************************************************************************
    1070             :  * The sieving procedure itself                                           *
    1071             :  **************************************************************************/
    1072             : static void
    1073    30737590 : sift(long b, ratpoints_bit_array *survivors, ratpoints_args *args,
    1074             :      bit_selection which_bits, ratpoints_bit_array bits16,
    1075             :      ratpoints_sieve_entry **sieve_list, long *bp_list, int *quit,
    1076             :      int process(long, long, GEN, void*, int*), void *info)
    1077    30737590 : {
    1078    30737590 :   pari_sp av = avma;
    1079    30737590 :   sieve_spec ssp[args->sp2];
    1080    30737590 :   int do_setup = 1;
    1081    30737590 :   long k, height = args->height, nb;
    1082             : 
    1083    30737590 :   if (odd(b) == 0) which_bits = num_odd; /* even denominator */
    1084             : 
    1085             :   /* Note that b is new */
    1086    30737590 :   args->bc = NULL;
    1087    30737590 :   nb = 0;
    1088             : 
    1089    74116849 :   for (k = 0; k < args->num_inter; k++)
    1090             :   {
    1091    47567804 :     ratpoints_interval inter = args->domain[k];
    1092             :     long low, high;
    1093             : 
    1094             :     /* Determine relevant interval [low, high] of numerators. */
    1095    47567804 :     if (b*inter.low <= -height)
    1096    23259181 :       low = -height;
    1097             :     else
    1098             :     {
    1099    24308623 :       if (b*inter.low > height) break;
    1100    20121835 :       low = ceil(b*inter.low);
    1101             :     }
    1102    43381016 :     if (b*inter.up >= height)
    1103    18844176 :       high = height;
    1104             :     else
    1105             :     {
    1106    24536840 :       if (b*inter.up < -height) continue;
    1107    20401813 :       high = floor(b*inter.up);
    1108             :     }
    1109             : 
    1110    39245989 :     if (do_setup)
    1111             :     { /* set up the sieve information */
    1112             :       long n;
    1113             : 
    1114    28550314 :       do_setup = 0; /* only do it once for every b */
    1115   485244584 :       for (n = 0; n < args->sp2; n++)
    1116             :       {
    1117   456694270 :         ratpoints_sieve_entry *se = sieve_list[n];
    1118   456694270 :         long p = se->p;
    1119   456694270 :         long bp = bp_list[n];
    1120             :         ratpoints_bit_array *sptr;
    1121             : 
    1122   456694270 :         if (which_bits != num_all) /* divide by 2 mod p */
    1123   268468115 :           bp = odd(bp) ? (bp+p) >> 1 : bp >> 1;
    1124   456694270 :         sptr = se->sieve[bp];
    1125             : 
    1126   456694270 :         ssp[n].p = p;
    1127   456694270 :         ssp[n].offset = (which_bits == num_odd) ? se->offset : 0;
    1128             : 
    1129             :         /* copy if already initialized, else initialize */
    1130   466512084 :         ssp[n].ptr = sptr ? sptr : (p<BITS_IN_LONG?sieve_init1(p, se, bp, args)
    1131     9817814 :                                                   :sieve_init2(p, se, bp, args));
    1132             :       }
    1133             :     }
    1134             : 
    1135    39245989 :     switch(which_bits)
    1136             :     {
    1137    16327768 :       case num_all: break;
    1138           0 :       case num_none: break;
    1139    17907968 :       case num_odd: low >>= 1; high--; high >>= 1; break;
    1140     5010253 :       case num_even: low++; low >>= 1; high >>= 1; break;
    1141             :     }
    1142             : 
    1143             :     /* now turn the bit interval into [low, high[ */
    1144    39245989 :     high++;
    1145             : 
    1146    39245989 :     if (low < high)
    1147             :     {
    1148    39244715 :       long w_low, w_high, w_low0, w_high0, range = args->array_size;
    1149             : 
    1150             :       /* Now the range of longwords (= bit_arrays) */
    1151    39244715 :       w_low = low >> RBA_SHIFT;
    1152    39244715 :       w_high = (high + (long)(RBA_LENGTH-1)) >> RBA_SHIFT;
    1153    39244715 :       w_low0 = w_low;
    1154    39244715 :       w_high0 = w_low0 + range;
    1155    83879090 :       for ( ; w_low0 < w_high; w_low0 = w_high0, w_high0 += range)
    1156             :       {
    1157             :         long i;
    1158    44636132 :         if (w_high0 > w_high) { w_high0 = w_high; range = w_high0 - w_low0; }
    1159             :         /* initialise the bits */
    1160  4823567643 :         for (i = range; i; i--) survivors[i-1] = bits16;
    1161             :         /* boundary words */
    1162    44636132 :         if (w_low0 == w_low)
    1163             :         {
    1164    39244715 :           long sh = low - RBA_LENGTH * w_low;
    1165    39244715 :           ulong *survl = (ulong *)survivors;
    1166             : 
    1167             : #ifdef HAS_SSE2
    1168    33637350 :           if (sh >= BITS_IN_LONG)
    1169             :           {
    1170    17825478 :             survl[0] = 0UL;
    1171    17825478 :             survl[1] &= (~0UL)<<(sh - BITS_IN_LONG);
    1172             :           }
    1173             :           else
    1174    15811872 :             survl[0] &= ~(0UL)<<sh;
    1175             : #else
    1176     5607365 :           survl[0] &= ~(0UL)<<sh;
    1177             : #endif
    1178             :         }
    1179    44636132 :         if (w_high0 == w_high)
    1180             :         {
    1181    39244591 :           long sh = RBA_LENGTH * w_high - high;
    1182    39244591 :           ulong *survl = (ulong *)&survivors[range-1];
    1183             : 
    1184             : #ifdef HAS_SSE2
    1185    33637344 :           if (sh >= BITS_IN_LONG)
    1186             :           {
    1187    17959932 :             survl[0] &= ~(0UL)>>(sh - BITS_IN_LONG);
    1188    17959932 :             survl[1] = 0UL;
    1189             :           }
    1190             :           else
    1191    15677412 :             survl[1] &= ~(0UL)>>sh;
    1192             : #else
    1193     5607247 :           survl[0] &= ~(0UL)>>sh;
    1194             : #endif
    1195             :         }
    1196    44636132 :         nb += _ratpoints_sift0(b, w_low0, w_high0, args, which_bits,
    1197             :                          survivors, &ssp[0], quit, process, info);
    1198    44636132 :         if (*quit) return;
    1199             :       }
    1200             :     }
    1201             :   }
    1202    30735833 :   if (nb==0) set_avma(av);
    1203             : }
    1204             : 
    1205             : /**************************************************************************
    1206             :  * Find points by looping over the denominators and sieving numerators    *
    1207             :  **************************************************************************/
    1208             : static void
    1209       12243 : find_points_work(ratpoints_args *args,
    1210             :                  int process(long, long, GEN, void*, int*), void *info)
    1211             : {
    1212       12243 :   int quit = 0;
    1213       12243 :   GEN c = args->cof;
    1214       12243 :   long degree = degpol(c);
    1215       12243 :   long nbprime = lg(args->listprime)-1;
    1216       12243 :   long height = args->height;
    1217             : 
    1218       12243 :   int point_at_infty = 0; /* indicates if there are points at infinity */
    1219       12243 :   int lcfsq = Z_issquare(pel(c,degree));
    1220             : 
    1221       12243 :   forbidden_entry *forb_ba = args->forb_ba;
    1222       12243 :   long *forbidden = args->forbidden;
    1223             :     /* The forbidden divisors, a zero-terminated array.
    1224             :        Used when degree is even and leading coefficient is not a square */
    1225             : 
    1226             :     /* These are used when degree is odd and leading coeff. is not +-1 */
    1227             : 
    1228             :   ratpoints_sieve_entry **sieve_list = (ratpoints_sieve_entry**)
    1229       12243 :      stack_malloc(nbprime*sizeof(ratpoints_sieve_entry*));
    1230       12243 :   bit_selection which_bits = num_all;
    1231             :   ulong den_bits;
    1232             :   ratpoints_bit_array num_bits[16];
    1233             : 
    1234       12243 :   args->flags &= RATPOINTS_FLAGS_INPUT_MASK;
    1235       12243 :   args->flags |= RATPOINTS_CHECK_DENOM;
    1236             : 
    1237             :   /* initialize memory management */
    1238       12243 :   args->se_next = args->se_buffer;
    1239       12243 :   args->ba_next = args->ba_buffer;
    1240       12243 :   args->int_next = args->int_buffer;
    1241             : 
    1242             :   /* Some sanity checks */
    1243       12243 :   args->num_inter = 0;
    1244             : 
    1245       12243 :   if (args->num_primes > nbprime) args->num_primes = nbprime;
    1246       12243 :   if (args->sp2 > args->num_primes) args->sp2 = args->num_primes;
    1247       12243 :   if (args->sp1 > args->sp2)        args->sp1 = args->sp2;
    1248             : 
    1249       12243 :   if (args->b_low < 1)  args->b_low = 1;
    1250       12243 :   if (args->b_high < 1) args->b_high = height;
    1251       12243 :   if (args->max_forbidden < 0)
    1252           0 :     args->max_forbidden = RATPOINTS_DEFAULT_MAX_FORBIDDEN;
    1253       12243 :   if (args->max_forbidden > nbprime)
    1254           0 :     args->max_forbidden = nbprime;
    1255       12243 :   if (args->array_size <= 0) args->array_size = RATPOINTS_ARRAY_SIZE;
    1256             :   {
    1257       12243 :     long s = 2*maxss(1,CEIL(height, BITS_IN_LONG));
    1258       12243 :     if (args->array_size > s) args->array_size = s;
    1259             :   }
    1260             :   /* Don't reverse if intervals are specified or limits for the denominator
    1261             :      are given */
    1262       12243 :   if (args->num_inter > 0 || args->b_low > 1 || args->b_high != height)
    1263          35 :     args->flags |= RATPOINTS_NO_REVERSE;
    1264             : 
    1265             :   /* Check if reversal of polynomial might be better:
    1266             :    * case 1: degree is even, but trailing coefficient is zero
    1267             :    * case 2: degree is even, leading coefficient is a square, but
    1268             :    *         trailing coefficient is not
    1269             :    * case 3: degree is odd, |leading coefficient| > 1,
    1270             :    *         trailing coefficient is zero, |coeff. of x| = 1 */
    1271       12243 :   if (!(args->flags & RATPOINTS_NO_REVERSE))
    1272             :   {
    1273       12208 :     if (!odd(degree))
    1274             :     {
    1275       11074 :       if (signe(pel(c,0)) == 0)
    1276             :       { /* case 1 */
    1277             :         long n;
    1278         224 :         args->flags |= RATPOINTS_REVERSED;
    1279         896 :         for (n = 0; n < degree>>1; n++) swap(pel(c,n), pel(c,degree-n));
    1280         224 :         degree--;
    1281         224 :         setlg(c,degree+3);
    1282             :       }
    1283       10850 :       else if (lcfsq && !Z_issquare(pel(c,0)))
    1284             :       { /* case 2 */
    1285             :         long n;
    1286         735 :         args->flags |= RATPOINTS_REVERSED;
    1287        2940 :         for (n = 0; n < degree>>1; n++) swap(pel(c,n), pel(c,degree-n));
    1288         735 :         lcfsq = 0;
    1289             :       }
    1290             :     }
    1291             :     else
    1292             :     { /* odd degree, case 3*/
    1293        1134 :       if (!is_pm1(pel(c,degree)) && !signe(pel(c,0)) && is_pm1(pel(c,1)))
    1294             :       {
    1295             :         long n;
    1296           7 :         args->flags |= RATPOINTS_REVERSED;
    1297          14 :         for (n = 1; n < degree>>1; n++) swap(pel(c,n),pel(c,degree+1-n));
    1298             :       }
    1299             :     }
    1300             :   }
    1301             : 
    1302             :   /* Deal with the intervals */
    1303       12243 :   if (args->num_inter == 0)
    1304             :   { /* default interval (effectively ]-oo,oo[) if none is given */
    1305       12243 :     args->domain = (ratpoints_interval*) stack_malloc(2*degree*sizeof(ratpoints_interval));
    1306       12243 :     args->domain[0].low = -height; args->domain[0].up = height;
    1307       12243 :     args->num_inter = 1;
    1308             :   }
    1309             : 
    1310       12243 :   ratpoints_compute_sturm(args);
    1311             : 
    1312             :   /* Point(s) at infinity? */
    1313       12243 :   if (odd(degree) || lcfsq)
    1314             :   {
    1315        1519 :     args->flags &= ~RATPOINTS_CHECK_DENOM;
    1316        1519 :     point_at_infty = 1;
    1317             :   }
    1318             : 
    1319             :   /* Can use only squares as denoms if degree is odd and poly is +-monic */
    1320       12243 :   if (odd(degree))
    1321             :   {
    1322        1372 :     GEN w1 = pel(c,degree);
    1323        1372 :     if (is_pm1(w1))
    1324          63 :       args->flags |= RATPOINTS_USE_SQUARES;
    1325             :     else /* set up information on divisors of leading coefficient */
    1326        1309 :       setup_us1(args, absi_shallow(w1));
    1327             :   }
    1328             : 
    1329             :   /* deal with f mod powers of 2 */
    1330       12243 :   which_bits = get_2adic_info(args, &den_bits, &num_bits[0]);
    1331             :   /* which_bits says whether to consider even and/or odd numerators
    1332             :    * when the denominator is odd.
    1333             :    *
    1334             :    * Bit k in den_bits is 0 if b congruent to k mod BITS_IN_LONG need
    1335             :    * not be considered as a denominator.
    1336             :    *
    1337             :    * Bit k in num_bits[b] is 0 is numerators congruent to
    1338             :    *  k (which_bits = den_all)
    1339             :    *  2k (which_bits = den_even)
    1340             :    *  2k+1 (which_bits = den_odd)
    1341             :    * need not be considered for denominators congruent to b mod 16. */
    1342             : 
    1343             :   /* set up the sieve data structure */
    1344       12243 :   if (sieving_info(args, sieve_list)) return;
    1345             : 
    1346             :   /* deal with point(s) at infinity */
    1347       12117 :   if (point_at_infty)
    1348             :   {
    1349        1519 :     long a = 1, b = 0;
    1350             : 
    1351        1519 :     if (args->flags & RATPOINTS_REVERSED) { a = 0; b = 1; }
    1352        1519 :     if (odd(degree))
    1353        1372 :       (void)process(a, b, gen_0, info, &quit);
    1354             :     else
    1355             :     {
    1356         147 :       GEN w0 = sqrti(pel(c,degree));
    1357         147 :       (void)process(a, b, w0, info, &quit);
    1358         147 :       (void)process(a, b, negi(w0), info, &quit);
    1359             :     }
    1360        1519 :     if (quit) return;
    1361             :   }
    1362             :   /* now do the sieving */
    1363             :   {
    1364             :     ratpoints_bit_array *survivors = (ratpoints_bit_array *)
    1365       12117 :       stack_malloc_align((args->array_size)*sizeof(ratpoints_bit_array), RBA_ALIGN);
    1366       12117 :     if (args->flags & (RATPOINTS_USE_SQUARES | RATPOINTS_USE_SQUARES1))
    1367             :     {
    1368        1372 :       if (args->flags & RATPOINTS_USE_SQUARES)
    1369             :       /* need only take squares as denoms */
    1370          63 :       {
    1371             :         long b, bb;
    1372          63 :         long bp_list[args->sp2];
    1373          63 :         long last_b = args->b_low;
    1374             :         long n;
    1375        1071 :         for (n = 0; n < args->sp2; n++)
    1376        1008 :           bp_list[n] = mod(args->b_low, sieve_list[n]->p);
    1377             : 
    1378        7875 :         for (b = 1; bb = b*b, bb <= args->b_high; b++)
    1379        7812 :           if (bb >= args->b_low)
    1380             :           {
    1381        7812 :             ratpoints_bit_array bits = num_bits[bb & 0xf];
    1382        7812 :             if (TEST(bits))
    1383             :             {
    1384             :               long n;
    1385        6916 :               long d = bb - last_b;
    1386             : 
    1387             :               /* fill bp_list */
    1388      117572 :               for (n = 0; n < args->sp2; n++)
    1389      110656 :                 bp_list[n] = mod(bp_list[n] + d, sieve_list[n]->p);
    1390        6916 :               last_b = bb;
    1391             : 
    1392        6916 :               sift(bb, survivors, args, which_bits, bits,
    1393             :                    sieve_list, &bp_list[0], &quit, process, info);
    1394        6916 :               if (quit) break;
    1395             :             }
    1396             :           }
    1397             :       }
    1398             :       else /* args->flags & RATPOINTS_USE_SQUARES1 */
    1399        1309 :       {
    1400        1309 :         GEN den_info = args->den_info;
    1401        1309 :         GEN divisors = args->divisors;
    1402        1309 :         long ld = lg(divisors);
    1403             :         long k;
    1404             :         long b, bb;
    1405        1309 :         long bp_list[args->sp2];
    1406             : 
    1407        4165 :         for (k = 1; k < ld; k++)
    1408             :         {
    1409        2863 :           long d = divisors[k];
    1410        2863 :           long last_b = d;
    1411             :           long n;
    1412        2863 :           if (d > args->b_high) continue;
    1413       48552 :           for (n = 0; n < args->sp2; n++)
    1414       45696 :             bp_list[n] = mod(d, sieve_list[n]->p);
    1415             : 
    1416      271873 :           for (b = 1; bb = d*b*b, bb <= args->b_high; b++)
    1417             :           {
    1418      269024 :             if (bb >= args->b_low)
    1419             :             {
    1420      269003 :               int flag = 1;
    1421      269003 :               ratpoints_bit_array bits = num_bits[bb & 0xf];
    1422             : 
    1423      269003 :               if (EXT0(bits))
    1424             :               {
    1425      222131 :                 long i, n, l = lg(gel(den_info,1));
    1426      222131 :                 long d = bb - last_b;
    1427             : 
    1428             :                 /* fill bp_list */
    1429     3776227 :                 for (n = 0; n < args->sp2; n++)
    1430     3554096 :                   bp_list[n] = mod(bp_list[n] + d, sieve_list[n]->p);
    1431      222131 :                 last_b = bb;
    1432             : 
    1433      421897 :                 for(i = 1; i < l; i++)
    1434             :                 {
    1435      250418 :                   long v = z_lval(bb, mael(den_info,1,i));
    1436      250418 :                   if((v >= mael(den_info,3,i))
    1437      120820 :                       && odd(v + (mael(den_info,2,i)))) { flag = 0; break; }
    1438             :                 }
    1439      222131 :                 if (flag)
    1440             :                 {
    1441      171479 :                   sift(bb, survivors, args, which_bits, bits,
    1442             :                        sieve_list, &bp_list[0], &quit, process, info);
    1443      171479 :                   if (quit) break;
    1444             :                 }
    1445             :               }
    1446             :             }
    1447             :           }
    1448        2856 :           if (quit) break;
    1449             :         }
    1450             :       }
    1451             :     }
    1452             :     else
    1453             :     {
    1454       10745 :       if (args->flags & RATPOINTS_CHECK_DENOM)
    1455       10598 :       {
    1456             :         long *forb;
    1457             :         long b;
    1458       10598 :         long bp_list[args->sp2];
    1459       10598 :         long last_b = args->b_low;
    1460             :         ulong b_bits;
    1461             :         long n;
    1462      180012 :         for (n = 0; n < args->sp2; n++)
    1463      169414 :           bp_list[n] = mod(args->b_low, sieve_list[n]->p);
    1464             :         {
    1465       10598 :           forbidden_entry *fba = &forb_ba[0];
    1466       10598 :           long b_low = args->b_low;
    1467       10598 :           long w_low = (b_low-1) >> TWOPOTBITS_IN_LONG;
    1468             : 
    1469       10598 :           b_bits = den_bits;
    1470      152729 :           while (fba->p)
    1471             :           {
    1472      142131 :             fba->curr = fba->start + mod(w_low, fba->p);
    1473      142131 :             b_bits &= *(fba->curr);
    1474      142131 :             fba++;
    1475             :           }
    1476       10598 :           b_bits >>= (b_low-1) & LONG_MASK;
    1477             :         }
    1478             : 
    1479   133780346 :         for (b = args->b_low; b <= args->b_high; b++)
    1480             :         {
    1481   133771491 :           ratpoints_bit_array bits = num_bits[b & 0xf];
    1482             : 
    1483   133771491 :           if ((b & LONG_MASK) == 0)
    1484             :           { /* next b_bits */
    1485     2381246 :             forbidden_entry *fba = &forb_ba[0];
    1486             : 
    1487     2381246 :             b_bits = den_bits;
    1488    36969088 :             while (fba->p)
    1489             :             {
    1490    34587842 :               fba->curr++;
    1491    34587842 :               if (fba->curr == fba->end) fba->curr = fba->start;
    1492    34587842 :               b_bits &= *(fba->curr);
    1493    34587842 :               fba++;
    1494             :             }
    1495             :           }
    1496             :           else
    1497   131390245 :             b_bits >>= 1;
    1498             : 
    1499   133771491 :           if (odd(b_bits) && EXT0(bits))
    1500             :           { /* check if denominator is excluded */
    1501    49874164 :             for (forb = &forbidden[0] ; *forb && (b % (*forb)); forb++)
    1502           0 :               continue;
    1503    49874164 :             if (*forb == 0 && rpjacobi(b, pel(c,degree)) == 1)
    1504             :             {
    1505    28881967 :               long n, d = b - last_b;
    1506             : 
    1507             :               /* fill bp_list */
    1508   490882685 :               for (n = 0; n < args->sp2; n++)
    1509             :               {
    1510   462000718 :                 long bp = bp_list[n] + d;
    1511   462000718 :                 long p = sieve_list[n]->p;
    1512             : 
    1513   524751148 :                 while (bp >= p) bp -= p;
    1514   462000718 :                 bp_list[n] = bp;
    1515             :               }
    1516    28881967 :               last_b = b;
    1517             : 
    1518    28881967 :               sift(b, survivors, args, which_bits, bits,
    1519             :                    sieve_list, &bp_list[0], &quit, process, info);
    1520    28881967 :               if (quit) break;
    1521             :             }
    1522             :           }
    1523             :         }
    1524             :       } /* if (args->flags & RATPOINTS_CHECK_DENOM) */
    1525             :       else
    1526         147 :       {
    1527             :         long b, n;
    1528         147 :         long bp_list[args->sp2];
    1529         147 :         long last_b = args->b_low;
    1530        2499 :         for (n = 0; n < args->sp2; n++)
    1531        2352 :           bp_list[n] = mod(args->b_low, sieve_list[n]->p);
    1532     2179156 :         for (b = args->b_low; b <= args->b_high; b++)
    1533             :         {
    1534     2179016 :           ratpoints_bit_array bits = num_bits[b & 0xf];
    1535     2179016 :           if (EXT0(bits))
    1536             :           {
    1537             :             long n;
    1538     1677228 :             long d = b - last_b;
    1539             : 
    1540             :             /* fill bp_list */
    1541    28512876 :             for (n = 0; n < args->sp2; n++)
    1542             :             {
    1543    26835648 :               long bp = bp_list[n] + d;
    1544    26835648 :               long p = sieve_list[n]->p;
    1545             : 
    1546    27855807 :               while (bp >= p) bp -= p;
    1547    26835648 :               bp_list[n] = bp;
    1548             :             }
    1549     1677228 :             last_b = b;
    1550             : 
    1551     1677228 :             sift(b, survivors, args, which_bits, bits,
    1552             :                  sieve_list, &bp_list[0], &quit, process, info);
    1553     1677228 :             if (quit) break;
    1554             :           }
    1555             :         }
    1556             :       }
    1557             :     }
    1558             :   }
    1559             : }
    1560             : 
    1561             : static GEN
    1562       85603 : vec_append_grow(GEN z, long i, GEN x)
    1563             : {
    1564       85603 :   long n = lg(z)-1;
    1565       85603 :   if (i > n)
    1566             :   {
    1567        1414 :     n <<= 1;
    1568        1414 :     z = vec_lengthen(z,n);
    1569             :   }
    1570       85603 :   gel(z,i) = x;
    1571       85603 :   return z;
    1572             : }
    1573             : 
    1574             : struct points
    1575             : {
    1576             :   GEN z;
    1577             :   long i, f;
    1578             : };
    1579             : 
    1580             : static int
    1581       88571 : process(long a, long b, GEN y, void *info0, int *quit)
    1582             : {
    1583       88571 :   struct points *p = (struct points *) info0;
    1584       88571 :   if (b==0 && (p->f&2L)) return 0;
    1585       85603 :   *quit = (p->f&1);
    1586       85603 :   p->z = vec_append_grow(p->z, ++p->i, mkvec3(stoi(a),y,stoi(b)));
    1587       85603 :   return 1;
    1588             : }
    1589             : 
    1590             : static int
    1591       12250 : args_h(ratpoints_args *args, GEN D)
    1592             : {
    1593       12250 :   long H, h, l = 1;
    1594             :   GEN L;
    1595       12250 :   switch(typ(D))
    1596             :   {
    1597       12208 :     case t_INT: if (signe(D) <= 0) return 0;
    1598       12208 :       H = h = itos(D); break;
    1599          42 :     case t_VEC: if (lg(D) != 3) return 0;
    1600          42 :       H = gtos(gel(D,1));
    1601          42 :       L = gel(D,2);
    1602          42 :       if (typ(L)==t_INT) { h = itos(L); break; }
    1603          14 :       if (typ(L)==t_VEC && lg(L)==3)
    1604             :       {
    1605           7 :         l = gtos(gel(L,1));
    1606           7 :         h = gtos(gel(L,2)); break;
    1607             :       }
    1608           7 :     default: return 0;
    1609             :   }
    1610       12243 :   args->height = H;
    1611       12243 :   args->b_low  = l;
    1612       12243 :   args->b_high = h; return 1;
    1613             : }
    1614             : 
    1615             : static GEN
    1616       12250 : ZX_hyperellratpoints(GEN P, GEN h, long flag)
    1617             : {
    1618       12250 :   pari_sp av = avma;
    1619             :   ratpoints_args args;
    1620             :   struct points data;
    1621       12250 :   ulong flags = 0;
    1622             : 
    1623       12250 :   if (!ZX_is_squarefree(P))
    1624           0 :     pari_err_DOMAIN("hyperellratpoints","issquarefree(pol)","=",gen_0, P);
    1625       12250 :   if (!args_h(&args, h))
    1626             :   {
    1627           7 :     pari_err_TYPE("hyperellratpoints", h);
    1628             :     return NULL;/*LCOV_EXCL_LINE*/
    1629             :   }
    1630       12243 :   find_points_init(&args, RATPOINTS_DEFAULT_BIT_PRIMES);
    1631             : 
    1632       12243 :   args.cof           = shallowcopy(P);
    1633       12243 :   args.num_inter     = 0;
    1634       12243 :   args.sp1           = RATPOINTS_DEFAULT_SP1;
    1635       12243 :   args.sp2           = RATPOINTS_DEFAULT_SP2;
    1636       12243 :   args.array_size    = RATPOINTS_ARRAY_SIZE;
    1637       12243 :   args.num_primes    = RATPOINTS_DEFAULT_NUM_PRIMES;
    1638       12243 :   args.bit_primes    = RATPOINTS_DEFAULT_BIT_PRIMES;
    1639       12243 :   args.max_forbidden = RATPOINTS_DEFAULT_MAX_FORBIDDEN;
    1640       12243 :   args.flags         = flags;
    1641       12243 :   data.z = cgetg(17,t_VEC);
    1642       12243 :   data.i = 0;
    1643       12243 :   data.f = flag;
    1644       12243 :   find_points_work(&args, process, (void *)&data);
    1645             : 
    1646       12243 :   setlg(data.z, data.i+1);
    1647       12243 :   return gerepilecopy(av, data.z);
    1648             : }
    1649             : 
    1650             : /* The ordinates of the points returned need to be divided by den
    1651             :  * by the caller to get the actual solutions */
    1652             : static GEN
    1653       12250 : QX_hyperellratpoints(GEN P, GEN h, long flag, GEN *den)
    1654             : {
    1655       12250 :   GEN Q = Q_remove_denom(P, den);
    1656       12250 :   if (*den) Q = ZX_Z_mul(Q, *den);
    1657       12250 :   return ZX_hyperellratpoints(Q, h, flag);
    1658             : }
    1659             : 
    1660             : static GEN
    1661         112 : ZX_homogenous_evalpow(GEN Q, GEN A, GEN B)
    1662             : {
    1663         112 :   pari_sp av = avma;
    1664         112 :   long i, d = degpol(Q);
    1665         112 :   GEN s = gel(Q, d+2);
    1666         448 :   for (i = d-1; i >= 0; i--)
    1667         336 :     s = addii(mulii(s, A), mulii(gel(B,d+1-i), gel(Q,i+2)));
    1668         112 :   return gerepileuptoint(av, s);
    1669             : }
    1670             : 
    1671             : static GEN
    1672          56 : to_RgX(GEN a, long v) { return typ(a)==t_POL? a: scalarpol(a,v); }
    1673             : 
    1674             : static int
    1675       11389 : hyperell_check(GEN PQ, GEN *P, GEN *Q)
    1676             : {
    1677       11389 :   *P = *Q = NULL;
    1678       11389 :   if (typ(PQ) == t_POL)
    1679             :   {
    1680       11361 :     if (!RgX_is_QX(PQ)) return 0;
    1681       11361 :     *P = PQ;
    1682             :   }
    1683             :   else
    1684             :   {
    1685          28 :     long v = gvar(PQ);
    1686          28 :     if (v == NO_VARIABLE || typ(PQ) != t_VEC || lg(PQ) != 3) return 0;
    1687          28 :     *P = to_RgX(gel(PQ,1), v); if (!RgX_is_QX(*P)) return 0;
    1688          28 :     *Q = to_RgX(gel(PQ,2), v); if (!RgX_is_QX(*Q)) return 0;
    1689          28 :     if (!signe(*Q)) *Q = NULL;
    1690             :   }
    1691       11389 :   return 1;
    1692             : }
    1693             : 
    1694             : GEN
    1695       11389 : hyperellratpoints(GEN PQ, GEN h, long flag)
    1696             : {
    1697       11389 :   pari_sp av = avma;
    1698             :   GEN P, Q, H, L, den;
    1699             :   long i, l, dy, dQ;
    1700             : 
    1701       11389 :   if (flag<0 || flag>1) pari_err_FLAG("ellratpoints");
    1702       11389 :   if (!hyperell_check(PQ, &P, &Q)) pari_err_TYPE("hyperellratpoints",PQ);
    1703       11389 :   if (!Q)
    1704             :   {
    1705       11375 :     L = QX_hyperellratpoints(P, h, flag|2L, &den);
    1706       11375 :     dy = (degpol(P)+1)>>1;
    1707       11375 :     l = lg(L);
    1708       25060 :     for (i = 1; i < l; i++)
    1709             :     {
    1710       13685 :       GEN Li = gel(L,i), x = gel(Li,1), y = gel(Li,2), z = gel(Li,3);
    1711       13685 :       GEN zdy = powiu(z, dy);
    1712       13685 :       if (den) zdy = mulii(zdy, den);
    1713       13685 :       gel(L,i) = mkvec2(gdiv(x,z), gdiv(y, zdy));
    1714             :     }
    1715       11375 :     return gerepilecopy(av, L);
    1716             :   }
    1717          14 :   H = RgX_add(RgX_muls(P,4), RgX_sqr(Q));
    1718          14 :   dy = (degpol(H)+1)>>1; dQ = degpol(Q);
    1719          14 :   L = QX_hyperellratpoints(H, h, flag|2L, &den);
    1720          14 :   l = lg(L);
    1721         126 :   for (i = 1; i < l; i++)
    1722             :   {
    1723         112 :     GEN Li = gel(L,i), x = gel(Li,1), y = gel(Li,2), z = gel(Li,3);
    1724         112 :     GEN B = gpowers(z, dQ);
    1725         112 :     GEN Qx = gdiv(ZX_homogenous_evalpow(Q, x, B), gel(B, dQ+1));
    1726         112 :     GEN zdy = powiu(z, dy);
    1727         112 :     if (den) zdy = mulii(zdy, den);
    1728         112 :     gel(L,i) = mkvec2(gdiv(x,z), gmul2n(gsub(gdiv(y,zdy),Qx),-1));
    1729             :   }
    1730          14 :   return gerepilecopy(av, L);
    1731             : }
    1732             : 
    1733             : GEN
    1734         861 : ellratpoints(GEN E, GEN h, long flag)
    1735             : {
    1736         861 :   pari_sp av = avma;
    1737             :   GEN L, a1, a3, den;
    1738             :   long i, l;
    1739         861 :   checkell_Q(E);
    1740         861 :   if (flag<0 || flag>1) pari_err_FLAG("ellratpoints");
    1741         861 :   if (!RgV_is_QV(vecslice(E,1,5))) pari_err_TYPE("ellratpoints",E);
    1742         861 :   a1 = ell_get_a1(E);
    1743         861 :   a3 = ell_get_a3(E);
    1744         861 :   L = QX_hyperellratpoints(ec_bmodel(E), h, flag|2L, &den);
    1745         854 :   l = lg(L);
    1746       72660 :   for (i = 1; i < l; i++)
    1747             :   {
    1748       71806 :     GEN P, Li = gel(L,i), x = gel(Li,1), y = gel(Li,2), z = gel(Li,3);
    1749       71806 :     if (!signe(z))
    1750           0 :       P = ellinf();
    1751             :     else
    1752             :     {
    1753       71806 :       GEN z2 = sqri(z);
    1754       71806 :       if (den) y = gdiv(y, den);
    1755       71806 :       y = gsub(y, gadd(gmul(a1, mulii(x,z)), gmul(a3,z2)));
    1756       71806 :       P = mkvec2(gdiv(x,z), gdiv(y,shifti(z2,1)));
    1757             :     }
    1758       71806 :     gel(L,i) = P;
    1759             :   }
    1760         854 :   return gerepilecopy(av, L);
    1761             : }

Generated by: LCOV version 1.13