Code coverage tests

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

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

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

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

Generated by: LCOV version 1.13