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

Generated by: LCOV version 1.11