Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - modules - ratpoints.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27783-affec94c65) Lines: 886 915 96.8 %
Date: 2022-07-07 07:34:25 Functions: 33 33 100.0 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2017  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : /* This file is based on ratpoints-2.1.3 by Michael Stoll, see
      16             :  * http://www.mathe2.uni-bayreuth.de/stoll/programs/
      17             :  * Original copyright / license: */
      18             : /***********************************************************************
      19             :  * ratpoints-2.2.1                                                     *
      20             :  * Copyright (C) 2008, 2009, 2022  Michael Stoll                       *
      21             :  *  - A program to find rational points on hyperelliptic curves        *
      22             :  *                                                                     *
      23             :  * This program is free software: you can redistribute it and/or       *
      24             :  * modify it under the terms of the GNU General Public License         *
      25             :  * as published by the Free Software Foundation, either version 2 of   *
      26             :  * the License, or (at your option) any later version.                 *
      27             :  ***********************************************************************/
      28             : 
      29             : #include "paricfg.h"
      30             : #ifdef HAS_AVX
      31             : #include <immintrin.h>
      32             : #elif defined(HAS_SSE2)
      33             : #include <emmintrin.h>
      34             : #endif
      35             : 
      36             : #include "pari.h"
      37             : #include "paripriv.h"
      38             : 
      39             : #define pel(a,b)  gel((a),(b)+2)
      40             : 
      41             : #define RATPOINTS_ARRAY_SIZE 256           /* Array size in longs */
      42             : #define RATPOINTS_DEFAULT_SP1 11           /* Default value for sp1 */
      43             : #define RATPOINTS_DEFAULT_SP2 19           /* Default value for sp2 */
      44             : #define RATPOINTS_DEFAULT_NUM_PRIMES 30    /* Default value for num_primes */
      45             : #define RATPOINTS_DEFAULT_BIT_PRIMES 7     /* Default value for bit_primes */
      46             : #define RATPOINTS_DEFAULT_MAX_FORBIDDEN 30 /* Default value for max_forbidden */
      47             : 
      48             : typedef struct {double low; double up;} ratpoints_interval;
      49             : 
      50             : /* Define the flag bits for the flags component: */
      51             : #define RATPOINTS_NO_REVERSE      0x0004UL
      52             : 
      53             : #define RATPOINTS_FLAGS_INPUT_MASK (RATPOINTS_NO_REVERSE)
      54             : 
      55             : /* Flags bits for internal purposes */
      56             : #define RATPOINTS_REVERSED        0x0100UL
      57             : #define RATPOINTS_CHECK_DENOM     0x0200UL
      58             : #define RATPOINTS_USE_SQUARES     0x0400UL
      59             : #define RATPOINTS_USE_SQUARES1    0x0800UL
      60             : 
      61             : #define LONG_MASK (~(-(1UL<<TWOPOTBITS_IN_LONG)))
      62             : 
      63             : #define CEIL(a,b) (((a) <= 0) ? -(-(a) / (b)) : 1 + ((a)-1) / (b))
      64             : 
      65             : /* define RBA_USE_VX provisionnaly */
      66             : #define RBA_USE_VX
      67             : #ifdef HAS_AVX512
      68             :  /* Use AVX512 512 bit registers for the bit arrays */
      69             : typedef ulong ratpoints_bit_array __attribute__ ((vector_size (64)));
      70             : 
      71             : #define EXT0(a) ((ulong)a[0])
      72             : #define EXT(a,i) ((ulong)a[i])
      73             : #define TEST(a) (  EXT0(a)  || EXT(a,1) || EXT(a,2) ||EXT(a,3)\
      74             :                 || EXT(a,4) || EXT(a,5) || EXT(a,6) ||EXT(a,7) )
      75             : 
      76             : #define RBA(a) ((ratpoints_bit_array) {((ulong) a), ((ulong) a), ((ulong) a), ((ulong) a)\
      77             :                                      , ((ulong) a), ((ulong) a), ((ulong) a), ((ulong) a) })
      78             : #define RBA_SHIFT (9)
      79             : #define MASKL(a,s) { unsigned long *survl = (unsigned long *)(a); long sh = (s); \
      80             :                      long l, qsh = sh>>TWOPOTBITS_IN_LONG, rsh = sh & (BITS_IN_LONG-1); \
      81             :                      for(l = 0; l < qsh; l++) { *survl++ = 0UL; }; *survl &= (~0UL)<<rsh; }
      82             : #define MASKU(a,s) { unsigned long *survl = (unsigned long *)(a); long sh = (s); \
      83             :                      long l, qsh = RBA_PACK-1 - (sh>>TWOPOTBITS_IN_LONG), rsh = sh & (BITS_IN_LONG-1); \
      84             :                      survl += qsh; *survl++ &= (~0UL)>>rsh; \
      85             :                      for(l = qsh+1; l < RBA_PACK; l++) { *survl++ = 0UL; } }
      86             : 
      87             : #elif defined(HAS_AVX)
      88             :  /* Use AVX 256 bit registers for the bit arrays */
      89             : typedef ulong ratpoints_bit_array __attribute__ ((vector_size (32)));
      90             : 
      91             : #define EXT0(a) ((ulong)a[0])
      92             : #define EXT(a,i) ((ulong)a[i])
      93             : 
      94             : #ifdef __AVX2__
      95             : #define TEST(a) ( _mm256_movemask_epi8(_mm256_cmpeq_epi8((__m256i)(a), (__m256i)RBA(0))) != 0xffffffffL )
      96             : #elif defined(__AVX__)
      97             : #define TEST(a) ( !_mm256_testz_si256((__m256i)(a), (__m256i)(a)) )
      98             : #else
      99             : #define TEST(a) (EXT(a,0) || EXT(a,1) || EXT(a,2) || EXT(a,3))
     100             : #endif
     101             : 
     102             : #define RBA(a) ((ratpoints_bit_array){((ulong) a), ((ulong) a), ((ulong) a), ((ulong) a)})
     103             : #define RBA_SHIFT (8)
     104             : #define MASKL(a,s) { unsigned long *survl = (unsigned long *)(a); long sh = (s); \
     105             :                      if(sh >= 2*BITS_IN_LONG) \
     106             :                      { sh -= 2*BITS_IN_LONG; survl[0] = 0UL; survl[1] = 0UL; \
     107             :                        if(sh >= BITS_IN_LONG) \
     108             :                        { survl[2] = 0UL; survl[3] &= (~0UL)<<(sh - BITS_IN_LONG); } \
     109             :                        else { survl[2] &= ~(0UL)<<sh; } } \
     110             :                      else if(sh >= BITS_IN_LONG) { survl[0] = 0UL; survl[1] &= (~0UL)<<(sh - BITS_IN_LONG); } \
     111             :                      else { survl[0] &= ~(0UL)<<sh; } }
     112             : #define MASKU(a,s) { unsigned long *survl = (unsigned long *)(a); long sh = (s); \
     113             :                      if(sh >= 2*BITS_IN_LONG) \
     114             :                      { sh -= 2*BITS_IN_LONG; survl[3] = 0UL; survl[2] = 0UL; \
     115             :                        if(sh >= BITS_IN_LONG) \
     116             :                        { survl[0] &= ~(0UL)>>(sh - BITS_IN_LONG); survl[1] = 0UL; } \
     117             :                        else { survl[1] &= ~(0UL)>>sh; } } \
     118             :                      else if(sh >= BITS_IN_LONG) { survl[2] &= ~(0UL)>>(sh - BITS_IN_LONG); survl[3] = 0UL; } \
     119             :                      else { survl[3] &= ~(0UL)>>sh; } }
     120             : #elif defined(HAS_SSE2) || defined(HAS_NEON)
     121             : 
     122             : #ifdef HAS_SSE2
     123             : /* Use SSE 128 bit registers for the bit arrays */
     124             : typedef __v2di ratpoints_bit_array;
     125             : #define EXT0(a) ((ulong)__builtin_ia32_vec_ext_v2di((__v2di)(a), 0))
     126             : #define EXT(a,i) ((ulong)__builtin_ia32_vec_ext_v2di((__v2di)(a), 1))
     127             : #else
     128             : typedef ulong ratpoints_bit_array __attribute__ ((vector_size (16)));
     129             : #define EXT0(a) ((ulong)a[0])
     130             : #define EXT(a,i) ((ulong)a[i])
     131             : #endif
     132             : 
     133             : #define TEST(a) (EXT0(a) || EXT(a,1))
     134             : #define RBA(a) ((ratpoints_bit_array){((long) a), ((long) a)})
     135             : #define RBA_SHIFT (7)
     136             : #define MASKL(a,s) { unsigned long *survl = (unsigned long *)(a); long sh = (s); \
     137             :                      if(sh >= BITS_IN_LONG) { survl[0] = 0UL; survl[1] &= (~0UL)<<(sh - BITS_IN_LONG); } \
     138             :                      else { survl[0] &= ~(0UL)<<sh; } }
     139             : #define MASKU(a,s) { unsigned long *survl = (unsigned long *)(a); long sh = (s); \
     140             :                      if(sh >= BITS_IN_LONG) { survl[0] &= ~(0UL)>>(sh - BITS_IN_LONG); survl[1] = 0UL; } \
     141             :                      else { survl[1] &= ~(0UL)>>sh; } }
     142             : #else
     143             : 
     144             : /* Use ulong for the bit arrays */
     145             : typedef ulong ratpoints_bit_array;
     146             : 
     147             : #define EXT0(a) (a)
     148             : #define TEST(a) (a)
     149             : #define RBA(a) (a)
     150             : #define RBA_SHIFT TWOPOTBITS_IN_LONG
     151             : #define MASKL(a,s) { *(a) &= ~(0UL)<<(s); }
     152             : #define MASKU(a,s) { *(a) &= ~(0UL)>>(s); }
     153             : #undef RBA_USE_VX
     154             : #endif
     155             : 
     156             : #define RBA_SIZE  (sizeof(ratpoints_bit_array))
     157             : #define RBA_LENGTH  (RBA_SIZE<<3)
     158             : #define RBA_PACK  (RBA_LENGTH>>TWOPOTBITS_IN_LONG)
     159             : 
     160             : #ifdef RBA_USE_VX
     161             : #define RATPOINTS_CHUNK 16
     162             : #define CODE_INIT_SIEVE_COPY \
     163             : { ulong k; \
     164             :       for (a = 0; a < p; a++) \
     165             :         for(k = 1; k < RBA_PACK; k++) \
     166             :           si[a+k*p] = si[a]; \
     167             :       for(a = 0; (ulong)a < (RATPOINTS_CHUNK-1)*RBA_PACK; a++) \
     168             :          si[a+p*RBA_PACK] = si[a];\
     169             : }
     170             : #else
     171             : #define RATPOINTS_CHUNK 1
     172             : #define CODE_INIT_SIEVE_COPY
     173             : #endif
     174             : 
     175             : typedef struct { long p; long offset; ratpoints_bit_array *ptr;
     176             :                  ratpoints_bit_array *start; ratpoints_bit_array *end; } sieve_spec;
     177             : 
     178             : typedef enum { num_all, num_even, num_odd, num_none } bit_selection;
     179             : 
     180             : typedef struct {
     181             :   long p; int *is_f_square;
     182             :   const long *inverses;
     183             :   long offset; ratpoints_bit_array** sieve;
     184             : } ratpoints_sieve_entry;
     185             : 
     186             : typedef struct { long p;
     187             :                  ulong *start;
     188             :                  ulong *end;
     189             :                  ulong *curr; }
     190             :                forbidden_entry;
     191             : 
     192             : typedef struct {
     193             :   GEN cof, listprime;
     194             :   ratpoints_interval *domain;
     195             :   long height, b_low, b_high, sp1, sp2, array_size;
     196             :   long num_inter, num_primes, bit_primes, max_forbidden;
     197             :   ulong flags;
     198             : /* from here: private data */
     199             :   GEN bc;
     200             :   ratpoints_sieve_entry *se_buffer;
     201             :   ratpoints_sieve_entry *se_next;
     202             :   ratpoints_bit_array *ba_buffer;
     203             :   ratpoints_bit_array *ba_next;
     204             :   int *int_buffer, *int_next;
     205             :   forbidden_entry *forb_ba;
     206             :   long *forbidden;
     207             :   GEN inverses, offsets, den_info, divisors;
     208             :   ulong **sieves0;
     209             : } ratpoints_args;
     210             : 
     211             : static ratpoints_bit_array *
     212     2674255 : sieve_init1(long p, ratpoints_sieve_entry *se1, long b1, ratpoints_args *args1)
     213             : {
     214     2674255 :   ratpoints_sieve_entry *se = se1;
     215     2674255 :   ratpoints_args *args = args1;
     216     2674255 :   int *isfs = se->is_f_square;
     217     2674255 :   long b = b1;
     218     2674255 :   long lmp = BITS_IN_LONG % p;
     219     2674255 :   long ldp = BITS_IN_LONG / p;
     220     2674255 :   long p1 = (ldp + 1) * p;
     221     2674255 :   long diff_shift = p1 & LONG_MASK;
     222     2674255 :   long diff = BITS_IN_LONG - diff_shift;
     223             :   ulong help0;
     224             :   long a;
     225     2674255 :   long d = se->inverses[b];
     226     2674255 :   long ab = 0; /* a/b mod p */
     227     2674255 :   ulong test = 1UL;
     228     2674255 :   ulong he0 = 0UL;
     229   118844764 :   for (a = 0; a < p; a++)
     230             :   {
     231   116170509 :     if (isfs[ab]) he0 |= test;
     232   116170509 :     ab += d;
     233   116170509 :     if (ab >= p) ab -= p;
     234   116170509 :     test <<= 1;
     235             :   }
     236     2674255 :   help0 = he0;
     237             :   {
     238             :     ulong help1;
     239             :      /* repeat bit pattern floor(BITS_IN_LONG/p) times */
     240     2674255 :     ulong pattern = help0;
     241             :     long i;
     242             :     /* the p * (floor(BITS_IN_LONG/p) + 1) - BITS_IN_LONG
     243             :             = p - (BITS_IN_LONG mod p)
     244             :        upper bits into help[b][1] :
     245             :        shift away the  BITS_IN_LONG mod p  lower bits */
     246     2674255 :     help1 = pattern >> lmp;
     247     6106045 :     for (i = p; i < BITS_IN_LONG; i <<= 1)
     248     3431790 :       help0 |= help0 << i;
     249             :     { /* fill the bit pattern from help0/help1 into sieve[b][].
     250             :           sieve[b][a0] has the same semantics as help0/help1,
     251             :           but here, a0 runs from 0 to p-1 and all bits are filled. */
     252             :       long a;
     253     2674255 :       ulong *si = (ulong *)args->ba_next;
     254             : 
     255     2674255 :       args->ba_next += p + RATPOINTS_CHUNK-1;
     256             :       /* copy the first chunk into sieve[b][] */
     257     2674255 :       si[0] = help0;
     258             :       /* now keep repeating the bit pattern,
     259             :          rotating it in help0/help1 */
     260   116170509 :       for (a = 1 ; a < p; a++)
     261             :       {
     262   113496254 :         ulong temp = help0 >> diff;
     263   113496254 :         help0 = help1 | (help0 << diff_shift);
     264   113496254 :         si[a] = help0;
     265   113496254 :         help1 = temp;
     266             :       }
     267   307438920 :       CODE_INIT_SIEVE_COPY
     268             :       /* set sieve array */
     269     2674255 :       se->sieve[b] = (ratpoints_bit_array *)si;
     270     2674255 :       return (ratpoints_bit_array *)si;
     271             :     }
     272             :   }
     273             : }
     274             : 
     275             : /* This is for p > BITS_IN_LONG */
     276             : static ratpoints_bit_array *
     277     9743885 : sieve_init2(long p, ratpoints_sieve_entry *se1, long b1, ratpoints_args *args1)
     278     9743885 : {
     279     9743885 :   ratpoints_sieve_entry *se = se1;
     280     9743885 :   ratpoints_args *args = args1;
     281     9743885 :   int *isfs = se->is_f_square;
     282     9743885 :   long b = b1;
     283             :   /* long ldp = 0;  = BITS_IN_LONG / p */
     284             :   /* long p1 = p; = (ldp + 1) * p; */
     285     9743885 :   long wp = p >> TWOPOTBITS_IN_LONG;
     286     9743885 :   long diff_shift = p & LONG_MASK;
     287     9743885 :   long diff = BITS_IN_LONG - diff_shift;
     288     9743885 :   ulong help[(p>>TWOPOTBITS_IN_LONG) + 2];
     289             : 
     290             :   /* initialize help */
     291             :   {
     292     9743885 :     ulong *he = &help[0];
     293     9743885 :     ulong *he1 = &he[(p>>TWOPOTBITS_IN_LONG) + 2];
     294    41179487 :     while (he1 != he) { he1--; *he1 = 0UL; }
     295             :   }
     296     9743885 :   { ulong work = 0UL;
     297             :     long a;
     298     9743885 :     long ab = 0; /* a/b mod p */
     299     9743885 :     long d = se->inverses[b];
     300     9743885 :     long n = 0;
     301     9743885 :     ulong test = 1UL;
     302   946962244 :     for (a = 0; a < p; )
     303             :     {
     304   937218359 :       if (isfs[ab]) work |= test;
     305   937218359 :       ab += d;
     306   937218359 :       if (ab >= p) ab -= p;
     307   937218359 :       test <<= 1;
     308   937218359 :       a++;
     309   937218359 :       if ((a & LONG_MASK) == 0)
     310    11947832 :       { help[n] = work; n++; work = 0UL; test = 1UL; }
     311             :     }
     312     9743885 :     help[n] = work;
     313             :   }
     314             : 
     315             :   { /* fill the bit pattern from help[] into sieve[b][].
     316             :        sieve[b][a0] has the same semantics as help[b][a0],
     317             :        but here, a0 runs from 0 to p-1 and all bits are filled. */
     318     9743885 :     ulong *si = (ulong *)args->ba_next;
     319             :     long a1;
     320             :     long a;
     321             : 
     322     9743885 :     args->ba_next += p + RATPOINTS_CHUNK-1;
     323             :     /* copy the first chunk from help[] into sieve[num][b][] */
     324    21691717 :     for (a = 0; a < wp; a++) si[a] = help[a];
     325             :     /* now keep repeating the bit pattern, rotating it in help */
     326   935014412 :     for (a1 = a ; a < p; a++)
     327             :     {
     328   925270527 :       long t = (a1 == wp) ? 0 : a1+1;
     329   925270527 :       help[a1] |= help[t]<<diff_shift;
     330   925270527 :       si[a] = help[a1];
     331   925270527 :       a1 = t;
     332   925270527 :       help[a1] >>= diff;
     333             :     }
     334  1828338288 :      CODE_INIT_SIEVE_COPY
     335             :     /* set sieve array */
     336     9743885 :     se->sieve[b] = (ratpoints_bit_array *)si;
     337     9743885 :     return (ratpoints_bit_array *)si;
     338             :   }
     339             : }
     340             : 
     341             : static GEN
     342       12194 : gen_squares(GEN listprime)
     343             : {
     344       12194 :   long nbprime = lg(listprime)-1;
     345       12194 :   GEN sq = cgetg(nbprime+1, t_VEC);
     346             :   long n;
     347      378014 :   for (n = 1; n <= nbprime; n++)
     348             :   {
     349      365820 :     ulong i, p = uel(listprime,n);
     350      365820 :     GEN w = zero_zv(p), work = w+1;
     351      365820 :     work[0] = 1;
     352             :     /* record nonzero squares mod p, p odd */
     353    10657556 :     for (i = 1; i < p; i += 2) work[(i*i) % p] = 1;
     354      365820 :     gel(sq, n) = w;
     355             :   }
     356       12194 :   return sq;
     357             : }
     358             : 
     359             : static GEN
     360       12194 : gen_offsets(GEN P)
     361             : {
     362       12194 :   long n, l = lg(P);
     363       12194 :   GEN of = cgetg(l, t_VEC);
     364      378014 :   for (n = 1; n < l; n++)
     365             :   {
     366      365820 :     ulong p = uel(P,n);
     367      365820 :     uel(of, n) = Fl_inv((2*RBA_LENGTH)%p, p);
     368             :   }
     369       12194 :   return of;
     370             : }
     371             : 
     372             : static GEN
     373       12194 : gen_inverses(GEN P)
     374             : {
     375       12194 :   long n, l = lg(P);
     376       12194 :   GEN iv = cgetg(l, t_VEC);
     377      378014 :   for (n = 1; n < l; n++)
     378             :   {
     379      365820 :     ulong i, p = uel(P,n);
     380      365820 :     GEN w = cgetg(p, t_VECSMALL);
     381    20949292 :     for (i = 1; i < p; i++) uel(w, i) = Fl_inv(i, p);
     382      365820 :     gel(iv, n) = w;
     383             :   }
     384       12194 :   return iv;
     385             : }
     386             : 
     387             : static ulong **
     388       12194 : gen_sieves0(GEN listprime)
     389             : {
     390             :   long n;
     391       12194 :   long nbprime = lg(listprime)-1;
     392       12194 :   ulong ** w = (ulong**) new_chunk(nbprime+1);
     393      378014 :   for (n = 1; n <= nbprime; n++)
     394             :   {
     395      365820 :     ulong a, p = uel(listprime,n);
     396      365820 :     ulong *si = (ulong *) stack_malloc_align((p+RATPOINTS_CHUNK-1)*RBA_SIZE, RBA_SIZE);
     397    21315112 :     for (a = 0; a < p; a++) si[a] = ~0UL;
     398    22105980 :     for (a = 0; a < BITS_IN_LONG; a++)
     399    21740160 :       uel(si,(p*a)>>TWOPOTBITS_IN_LONG) &= ~(1UL<<((p*a) & LONG_MASK));
     400    45633432 :     CODE_INIT_SIEVE_COPY
     401      365820 :     w[n] = si;
     402             :   }
     403       12194 :   return w;
     404             : }
     405             : 
     406             : static void
     407       12194 : gen_sieve(ratpoints_args *args)
     408             : {
     409       12194 :   GEN listprimes = args->listprime;
     410       12194 :   args->offsets = gen_offsets(listprimes);
     411       12194 :   args->inverses = gen_inverses(listprimes);
     412       12194 :   args->sieves0 = gen_sieves0(listprimes);
     413       12194 : }
     414             : 
     415             : static GEN
     416       12194 : ZX_positive_region(GEN P, long h, long bitprec)
     417             : {
     418       12194 :   long prec = nbits2prec(bitprec);
     419       12194 :   GEN it = mkvec2(stoi(-h),stoi(h));
     420       12194 :   GEN R = realroots(P, it, prec);
     421       12194 :   long nR = lg(R)-1;
     422       12194 :   long s = signe(ZX_Z_eval(P,gel(it,1)));
     423       12194 :   long i=1, j;
     424             :   GEN iv, st, en;
     425       12194 :   if (s<0 && nR==0) return NULL;
     426       11529 :   iv = cgetg(((nR+1+(s>=0))>>1)+1, t_VEC);
     427       11529 :   if (s>=0) st = itor(gel(it,1),prec);
     428        5033 :   else    { st = gel(R,i); i++; }
     429       17836 :   for (j=1; i<nR; j++)
     430             :   {
     431        6307 :     gel(iv, j) = mkvec2(st, gel(R,i));
     432        6307 :     st = gel(R,i+1);
     433        6307 :     i+=2;
     434             :   }
     435       11529 :   if (i==nR) en = gel(R,i); else en = itor(gel(it,2),prec);
     436       11529 :   gel(iv,j) = mkvec2(st, en);
     437       11529 :   return iv;
     438             : }
     439             : 
     440             : static long
     441       12194 : posint(ratpoints_interval *ivlist, GEN P, long h)
     442             : {
     443       12194 :   GEN R = ZX_positive_region(P, h, 53);
     444       12194 :   const double eps = 1e-5;
     445             :   long nR, i;
     446             : 
     447       12194 :   if (!R) return 0;
     448       11529 :   nR = lg(R)-1;
     449       11529 :   i = 1;
     450       29365 :   for (i=1; i<=nR; i++)
     451             :   {
     452       17836 :     ivlist[i-1].low = rtodbl(gmael(R,i,1))-eps;
     453       17836 :     ivlist[i-1].up  = rtodbl(gmael(R,i,2))+eps;
     454             :   }
     455       11529 :   return nR;
     456             : }
     457             : 
     458             : static long
     459       12194 : ratpoints_compute_sturm(ratpoints_args *args)
     460             : {
     461       12194 :   ratpoints_interval *ivlist = args->domain;
     462       12194 :   args->num_inter = posint(ivlist, args->cof, (long) ivlist[0].up);
     463       12194 :   return args->num_inter;
     464             : }
     465             : 
     466             : /**************************************************************************
     467             :  * Try to avoid divisions                                                 *
     468             :  **************************************************************************/
     469             : INLINE long
     470   811638122 : mod(long a, long b)
     471             : {
     472   811638122 :   long b1 = b << 4; /* b1 = 16*b */
     473             : 
     474   811638122 :   if (a < -b1) { a %= b; if (a < 0) { a += b; } return a ; }
     475   802432872 :   if (a < 0) { a += b1; }
     476   319822442 :   else { if (a >= b1) { return a % b; } }
     477   795740227 :   b1 >>= 1; /* b1 = 8*b */
     478   795740227 :   if (a >= b1) { a -= b1; }
     479   795740227 :   b1 >>= 1; /* b1 = 4*b */
     480   795740227 :   if (a >= b1) { a -= b1; }
     481   795740227 :   b1 >>= 1; /* b1 = 2*b */
     482   795740227 :   if (a >= b1) { a -= b1; }
     483   795740227 :   if (a >= b) { a -= b; }
     484   795740227 :   return a;
     485             : }
     486             : 
     487             : static void
     488     1763251 : set_bc(long b, ratpoints_args *args)
     489             : {
     490     1763251 :   GEN w0 = gen_1;
     491     1763251 :   GEN c = args->cof, bc;
     492     1763251 :   long k, degree = degpol(c);
     493     1763251 :   bc = cgetg(degree+2, t_POL);
     494     9064468 :   for (k = degree-1; k >= 0; k--)
     495             :   {
     496     7301217 :     w0 = muliu(w0, b);
     497     7301217 :     gel(bc,k+2) = mulii(gel(c,k+2), w0);
     498             :   }
     499     1763251 :   args->bc = bc;
     500     1763251 : }
     501             : 
     502             : /**************************************************************************
     503             :  * Check a `survivor' of the sieve if it really gives a point.            *
     504             :  **************************************************************************/
     505             : 
     506             : static long
     507     2553511 : _ratpoints_check_point(long a, long b, ratpoints_args *args, int *quit,
     508             :                  int process(long, long, GEN, void*, int*), void *info)
     509             : {
     510     2553511 :   pari_sp av = avma;
     511     2553511 :   GEN w0, w2, c = args->cof, bc = args->bc;
     512     2553511 :   long k, degree = degpol(c);
     513     2553511 :   int reverse = args->flags & RATPOINTS_REVERSED;
     514             : 
     515             :   /* Compute F(a, b), where F is the homogenized version of f
     516             :      of smallest possible even degree  */
     517     2553511 :   w2 = gel(c, degree+2);
     518    13289720 :   for (k = degree-1; k >= 0; k--)
     519             :   {
     520    10736209 :     w2 = mulis(w2, a);
     521    10736209 :     w2 = addii(w2, gel(bc,k+2));
     522             :   }
     523     2553511 :   if (odd(degree)) w2 = muliu(w2, b);
     524             :   /* check if f(x,z) is a square; if so, process the point(s) */
     525     2553511 :   if (signe(w2) >= 0 && Z_issquareall(w2, &w0))
     526             :   {
     527       44604 :     if (reverse)
     528             :     {
     529        1218 :       if (a >= 0) (void)process(b, a, w0, info, quit);
     530         217 :       else        (void)process(-b, -a, w0, info, quit);
     531             :     }
     532       43386 :     else (void)process(a, b, w0, info, quit);
     533       44604 :     if (!*quit && signe(w0) != 0)
     534             :     {
     535       42357 :       GEN nw0 = negi(w0);
     536       42357 :       if (reverse)
     537             :       {
     538        1155 :         if (a >= 0) (void)process(b, a, nw0, info, quit);
     539         196 :         else        (void)process(-b, -a, nw0, info, quit);
     540             :       }
     541       41202 :       else (void)process(a, b, nw0, info, quit);
     542             :     }
     543       44604 :     return 1;
     544             :   }
     545     2508907 :   set_avma(av);
     546     2508907 :   return 0;
     547             : }
     548             : 
     549             : /**************************************************************************
     550             :  * The inner loop of the sieving procedure                                *
     551             :  **************************************************************************/
     552             : static long
     553    42497850 : _ratpoints_sift0(long b, long w_low, long w_high,
     554             :            ratpoints_args *args, bit_selection which_bits,
     555             :            ratpoints_bit_array *survivors, sieve_spec *sieves, int *quit,
     556             :            int process(long, long, GEN, void*, int*), void *info)
     557             : {
     558    42497850 :   long sp1 = args->sp1, sp2 = args->sp2;
     559    42497850 :   long i, n, nb = 0, absb = labs(b), base = 0;
     560             :   ratpoints_bit_array *surv0;
     561             : 
     562             :   /* now do the sieving (fast!) */
     563             : #if (RATPOINTS_CHUNK == 16)
     564             :   long w_low_new;
     565    32640480 :   ratpoints_bit_array *surv = survivors;
     566             : 
     567             :   /* first set the start fields for the first and second phases of sieving */
     568   652381902 :   for(n = 0; n < sp2; n++)
     569   619741422 :     sieves[n].start = sieves[n].ptr + mod(w_low + sieves[n].offset, sieves[n].p);
     570             :   /* Take RATPOINTS_CHUNK bit-arrays and apply phase 1 to them,
     571             :    * then repeat with the next RATPOINTS_CHUNK bit-arrays. */
     572   216216696 :   for(w_low_new = w_low; w_low_new < w_high; surv += RATPOINTS_CHUNK, w_low_new += RATPOINTS_CHUNK)
     573             :   {
     574             :     /* read data from memory into registers */
     575   183576216 :     ratpoints_bit_array reg0 = surv[0];
     576   183576216 :     ratpoints_bit_array reg1 = surv[1];
     577   183576216 :     ratpoints_bit_array reg2 = surv[2];
     578   183576216 :     ratpoints_bit_array reg3 = surv[3];
     579   183576216 :     ratpoints_bit_array reg4 = surv[4];
     580   183576216 :     ratpoints_bit_array reg5 = surv[5];
     581   183576216 :     ratpoints_bit_array reg6 = surv[6];
     582   183576216 :     ratpoints_bit_array reg7 = surv[7];
     583   183576216 :     ratpoints_bit_array reg8 = surv[8];
     584   183576216 :     ratpoints_bit_array reg9 = surv[9];
     585   183576216 :     ratpoints_bit_array reg10 = surv[10];
     586   183576216 :     ratpoints_bit_array reg11 = surv[11];
     587   183576216 :     ratpoints_bit_array reg12 = surv[12];
     588   183576216 :     ratpoints_bit_array reg13 = surv[13];
     589   183576216 :     ratpoints_bit_array reg14 = surv[14];
     590   183576216 :     ratpoints_bit_array reg15 = surv[15];
     591             : 
     592  2202914592 :     for(n = 0; n < sp1; n++)
     593             :     { /* retrieve the pointer to the beginning of the relevant bits */
     594  2019338376 :       ratpoints_bit_array *siv1 = sieves[n].start;
     595  2019338376 :       reg0 &= *siv1++;
     596  2019338376 :       reg1 &= *siv1++;
     597  2019338376 :       reg2 &= *siv1++;
     598  2019338376 :       reg3 &= *siv1++;
     599  2019338376 :       reg4 &= *siv1++;
     600  2019338376 :       reg5 &= *siv1++;
     601  2019338376 :       reg6 &= *siv1++;
     602  2019338376 :       reg7 &= *siv1++;
     603  2019338376 :       reg8 &= *siv1++;
     604  2019338376 :       reg9 &= *siv1++;
     605  2019338376 :       reg10 &= *siv1++;
     606  2019338376 :       reg11 &= *siv1++;
     607  2019338376 :       reg12 &= *siv1++;
     608  2019338376 :       reg13 &= *siv1++;
     609  2019338376 :       reg14 &= *siv1++;
     610  2019338376 :       reg15 &= *siv1++;
     611             : 
     612             :       /* update the pointer for the next round
     613             :        * (RATPOINTS_CHUNK-1 bit-arrays after sieves[n].end) */
     614  2884644192 :       while(siv1 >= sieves[n].end) siv1 -= sieves[n].p;
     615  2019338376 :       sieves[n].start = siv1;
     616             :     }
     617             :     /* store the contents of the registers back into memory */
     618   183576216 :     surv[0] = reg0;
     619   183576216 :     surv[1] = reg1;
     620   183576216 :     surv[2] = reg2;
     621   183576216 :     surv[3] = reg3;
     622   183576216 :     surv[4] = reg4;
     623   183576216 :     surv[5] = reg5;
     624   183576216 :     surv[6] = reg6;
     625   183576216 :     surv[7] = reg7;
     626   183576216 :     surv[8] = reg8;
     627   183576216 :     surv[9] = reg9;
     628   183576216 :     surv[10] = reg10;
     629   183576216 :     surv[11] = reg11;
     630   183576216 :     surv[12] = reg12;
     631   183576216 :     surv[13] = reg13;
     632   183576216 :     surv[14] = reg14;
     633   183576216 :     surv[15] = reg15;
     634             :   }
     635             : #else /* RATPOINTS_CHUNK not between 2 and 16 */
     636     9857370 :   long range = w_high - w_low;
     637   118288440 :   for (n = 0; n < sp1; n++)
     638             :   {
     639   108431070 :     ratpoints_bit_array *sieve_n = sieves[n].ptr;
     640   108431070 :     long p = sieves[n].p;
     641   108431070 :     long r = mod(-w_low-sieves[n].offset, p);
     642   108431070 :     ratpoints_bit_array *surv = survivors;
     643             : 
     644   108431070 :     if (w_high < w_low + r)
     645             :     { /* if we get here, r > 0, since w_high >= w_low always */
     646     6706840 :       ratpoints_bit_array *siv1 = &sieve_n[p-r];
     647     6706840 :       ratpoints_bit_array *siv0 = siv1 + range;
     648             : 
     649   203595651 :       while (siv1 != siv0) { *surv &= *siv1++; surv++; }
     650             :     }
     651             :     else
     652             :     {
     653   101724230 :       ratpoints_bit_array *siv1 = &sieve_n[p-r];
     654   101724230 :       ratpoints_bit_array *surv_end = &survivors[range - p];
     655             :       long i;
     656  3133365669 :       for (i = r; i; i--) { *surv &= *siv1++; surv++; }
     657   101724230 :       siv1 -= p;
     658   532677976 :       while (surv <= surv_end)
     659             :       {
     660 13998935014 :         for (i = p; i; i--) { *surv &= *siv1++; surv++; }
     661   430953746 :         siv1 -= p;
     662             :       }
     663   101724230 :       surv_end += p;
     664  3216436715 :       while (surv < surv_end) { *surv &= *siv1++; surv++; }
     665             :     }
     666             :   }
     667             :   /* initialize pointers in sieve for the second phase */
     668    88548580 :   for(n = sp1; n < sp2; n++)
     669    78691210 :     sieves[n].start = sieves[n].ptr + mod(w_low + sieves[n].offset, sieves[n].p);
     670             : #endif /* RATPOINTS_CHUNK */
     671             : 
     672             :   /* 2nd phase of the sieve: test each surviving bit array with more primes */
     673    42497850 :   surv0 = &survivors[0];
     674  4789730218 :   for (i = w_low; i < w_high; i++, base++)
     675             :   {
     676  4747234125 :     ratpoints_bit_array nums = *surv0++;
     677  4747234125 :     sieve_spec *ssp = &sieves[sp1];
     678             :     long n;
     679             : 
     680  5137859782 :     for (n = sp2-sp1; n && TEST(nums); n--)
     681             :     {
     682   390625657 :       ratpoints_bit_array *ptr = (ssp->start) + base;
     683   390625657 :       long p = ssp->p;
     684   977695546 :       while(ptr >= ssp->end) ptr -= p;
     685   390625657 :       nums &= *ptr;
     686   390625657 :       ssp++;
     687             :     }
     688             : 
     689             :     /* Check the survivors of the sieve if they really give points */
     690  4747234125 :     if (TEST(nums))
     691             :     {
     692             :       long a0, a, d;
     693     9874619 :       ulong nums0 = EXT0(nums);
     694             :       /* a will be the numerator corresponding to the selected bit */
     695     9874619 :       if (which_bits == num_all)
     696             :       {
     697     6376393 :         d = 1; a0 = i * RBA_LENGTH;
     698             :       }
     699             :       else
     700             :       {
     701     3498226 :         d = 2; a0 = i * 2 * RBA_LENGTH;
     702     3498226 :         if (which_bits == num_odd) a0++;
     703             :       }
     704   118604716 :       for (a = a0; nums0; a += d, nums0 >>= 1)
     705             :       { /* test one bit */
     706   108731098 :         if (odd(nums0) && ugcd(labs(a), absb)==1)
     707             :         {
     708     1460089 :           if (!args->bc) set_bc(b, args);
     709     1460089 :           nb += _ratpoints_check_point(a, b, args, quit, process, info);
     710     1460089 :           if (*quit) return nb;
     711             :         }
     712             :       }
     713             : #ifdef RBA_USE_VX
     714             :       {
     715     8459508 :         ulong k, da = d<<TWOPOTBITS_IN_LONG;
     716    16918260 :         for (k = 1; k < RBA_PACK; k++)
     717             :         {
     718     8459508 :           ulong nums1 = EXT(nums,k);
     719     8459508 :           a0 += da;
     720    98716962 :           for (a = a0; nums1; a += d, nums1 >>= 1)
     721             :           { /* test one bit */
     722    90258210 :             if (odd(nums1) && ugcd(labs(a), absb)==1)
     723             :             {
     724     1093422 :               if (!args->bc) set_bc(b, args);
     725     1093422 :               nb += _ratpoints_check_point(a, b, args, quit, process, info);
     726     1093422 :               if (*quit) return nb;
     727             :             }
     728             :           }
     729             :         }
     730             :       }
     731             : #endif
     732             :     }
     733             :   }
     734    42496093 :   return nb;
     735             : }
     736             : 
     737             : typedef struct { double r; ratpoints_sieve_entry *ssp; } entry;
     738             : 
     739             : static const int squares16[16] = {1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0};
     740             :  /* Says if a is a square mod 16, for a = 0..15 */
     741             : 
     742             : /**************************************************************************
     743             :  * Initialization and cleanup of ratpoints_args structure                 *
     744             :  **************************************************************************/
     745             : 
     746             : static ratpoints_sieve_entry*
     747       12194 : alloc_sieve(long nbprime, long maxprime)
     748             : {
     749             :   long i;
     750             :   ratpoints_sieve_entry * s = (ratpoints_sieve_entry*)
     751       12194 :                         stack_malloc(nbprime*sizeof(ratpoints_sieve_entry));
     752      378014 :   for (i=0; i<nbprime; i++)
     753      365820 :     s[i].sieve = (ratpoints_bit_array**) new_chunk(maxprime);
     754       12194 :   return s;
     755             : }
     756             : 
     757             : /* NOTE: args->degree must be set */
     758             : static void
     759       12194 : find_points_init(ratpoints_args *args, long bit_primes)
     760             : {
     761       12194 :   long need = 0;
     762             :   long n, nbprime,maxprime;
     763       12194 :   args->listprime = primes_interval_zv(3, 1<<bit_primes);
     764       12194 :   nbprime = lg(args->listprime)-1;
     765       12194 :   maxprime = args->listprime[nbprime];
     766             : 
     767             :   /* allocate space for se_buffer */
     768       12194 :   args->se_buffer = alloc_sieve(nbprime, maxprime);
     769       12194 :   args->se_next = args->se_buffer;
     770      378014 :   for (n = 1; n <= nbprime; n++)
     771             :   {
     772      365820 :     ulong p = args->listprime[n];
     773      365820 :     need += p*(p + RATPOINTS_CHUNK-1);
     774             :   }
     775       12194 :   args->ba_buffer = (ratpoints_bit_array*)
     776       12194 :      stack_malloc_align(need*RBA_SIZE,RBA_SIZE);
     777       12194 :   args->ba_next = args->ba_buffer;
     778             : 
     779             :   /* allocate space for int_buffer */
     780       12194 :   args->int_buffer = (int *) stack_malloc(nbprime*(maxprime+1)*sizeof(int));
     781       12194 :   args->int_next = args->int_buffer;
     782             : 
     783       12194 :   args->forb_ba   = (forbidden_entry*)
     784       12194 :     stack_malloc((nbprime + 1)*sizeof(forbidden_entry));
     785       12194 :   args->forbidden = new_chunk(nbprime + 1);
     786       12194 :   gen_sieve(args);
     787       12194 :   return;
     788             : }
     789             : 
     790             : /* f = leading coeff; b = b1*b2, b1 maximal with (b1, 2*f) = 1;
     791             :  * return Jacobi symbol (f, b1) */
     792             : INLINE int
     793    48413876 : rpjacobi(long b, GEN lcf)
     794             : {
     795             :   ulong f;
     796    48413876 :   b >>= vals(b);
     797    48413876 :   f = umodiu(lcf, b);
     798    48413876 :   return krouu(f, u_ppo(b,f));
     799             : }
     800             : 
     801             : /************************************************************************
     802             :  * Set up information on possible denominators                          *
     803             :  * when polynomial is of odd degree with leading coefficient != +-1     *
     804             :  ************************************************************************/
     805             : 
     806             : static void
     807        1316 : setup_us1(ratpoints_args *args, GEN w0)
     808             : {
     809        1316 :   GEN F = Z_issmooth_fact(w0, 1000), P, E, S, D;
     810             :   long i, l;
     811             : 
     812        1316 :   if (!F) return;
     813        1316 :   P = gel(F,1); l = lg(P);
     814        1316 :   E = gel(F,2);
     815        1316 :   D  = cgetg(1+(1<<(l-1)), t_VECSMALL);
     816             :   /* factorization is complete, set up array of squarefree divisors */
     817        1316 :   D[1] = 1;
     818        2758 :   for (i = 1; i < l; i++)
     819             :   { /* multiply all divisors known so far by next prime */
     820        1442 :     long k, n = 1<<(i-1);
     821        3010 :     for (k=0; k<n; k++) uel(D,1+n+k) = uel(D,1+k) * P[i];
     822             :   }
     823        1316 :   S = cgetg(l, t_VECSMALL);
     824             :   /* set slopes in den_info */
     825        2758 :   for (i = 1; i < l; i++)
     826             :   { /* compute min{n : (d-k)*n > v_p(f_d) - v_p(f_k), k = 0,...,d-1} */
     827        1442 :     GEN c = args->cof;
     828        1442 :     long p = P[i], v = E[i];
     829        1442 :     long k, n = 1, d = degpol(c);
     830             : 
     831        6902 :     for (k = d - 1; k >= 0; k--)
     832             :     {
     833        5460 :       long t = 1 + v - Z_lval(gel(c,k+2), p);
     834        5460 :       long m = CEIL(t, d - k);
     835             : 
     836        5460 :       if (m > n) n = m;
     837             :     }
     838        1442 :     S[i] = n;
     839             :   }
     840        1316 :   args->divisors = D;
     841        1316 :   args->flags |= RATPOINTS_USE_SQUARES1;
     842        1316 :   args->den_info = mkvec3(P, E, S);
     843             : }
     844             : 
     845             : /************************************************************************
     846             :  * Consider 2-adic information                                          *
     847             :  ************************************************************************/
     848             : 
     849             : static bit_selection
     850       12194 : get_2adic_info(ratpoints_args *args, ulong *den_bits,
     851             :                ratpoints_bit_array *num_bits)
     852             : {
     853       12194 :   GEN c = args->cof;
     854       12194 :   long degree = degpol(c);
     855             :   int is_f_square16[24];
     856       12194 :   long *cmp = new_chunk(degree+1);
     857       12194 :   long npe = 0, npo = 0;
     858             :   bit_selection result;
     859             :   long n, a, b;
     860             : 
     861             :   /* compute coefficients mod 16 */
     862       85813 :   for (n = 0; n <= degree; n++) cmp[n] = Mod16(gel(c,n+2));
     863      207298 :   for (a = 0 ; a < 16; a++)
     864             :   {
     865      195104 :     ulong s = cmp[degree];
     866             :     long n;
     867     1177904 :     for (n = degree - 1 ; n >= 0 ; n--) s = s*a + cmp[n];
     868      195104 :     s &= 0xf;
     869      195104 :     if ((is_f_square16[a] = squares16[s])) { if (odd(a)) npo++; else npe++; }
     870             :   }
     871             : 
     872             :   /* even denominators:
     873             :      is_f_square16[16+k] says if f((2k+1)/2) is a square, k = 0..3
     874             :      is_f_square16[20+k] says if f((2k+1)/4) is a square, k = 0,1
     875             :      is_f_square16[22]   says if f(odd/8) is a square
     876             :      is_f_square16[23]   says if f(odd/2^n), n >= 4, can be a square */
     877             :   {
     878       12194 :     long np1 = 0, np2 = 0, np3 = 0, np4 = 0;
     879             : 
     880       12194 :     if (odd(degree))
     881             :     {
     882        1379 :       long a, cf = 4*cmp[degree-1];
     883             : 
     884        1379 :       if (degree >= 2) cf += 8*cmp[degree-2];
     885        6895 :       for (a = 0; a < 4; a++)
     886             :       { /* Compute  2 c[d] k^d + 4 c[d-1] k^(d-1) + 8 c[d-2] k^(d-2), k = 2a+1.
     887             :            Note that k^d = k mod 8, k^(d-1) = 1 mod 8. */
     888        5516 :         long k = 2*a+1;
     889        5516 :         long s = (2*k*cmp[degree] + cf) & 0xf;
     890        5516 :         if ((is_f_square16[16+a] = squares16[s])) np1++;
     891             :       }
     892        1379 :       if ((is_f_square16[20] = squares16[(4*cmp[degree])  & 0xf])) np2++;
     893        1379 :       if ((is_f_square16[21] = squares16[(12*cmp[degree]) & 0xf])) np2++;
     894        1379 :       if ((is_f_square16[22] = squares16[(8*cmp[degree])  & 0xf])) np3++;
     895        1379 :       is_f_square16[23] = 1; np4++;
     896             :     }
     897             :     else
     898             :     {
     899       10815 :       long a, cf = (degree >= 2) ? 4*cmp[degree-2] : 0;
     900             : 
     901       10815 :       if (degree >= 3) cf += 8*cmp[degree-3];
     902       54075 :       for (a = 0; a < 4; a++)
     903             :       { /* compute c[d] k^d + 2 c[d-1] k^(d-1) + ... + 8 c[d-3] k^(d-3),
     904             :            k = 2a+1. Note that k^d = k^2 mod 16, k^(d-1) = k mod 8. */
     905       43260 :         long k = 2*a+1;
     906       43260 :         long s = ((cmp[degree]*k + 2*cmp[degree-1])*k + cf) & 0xf;
     907       43260 :         if ((is_f_square16[16+a] = squares16[s])) np1++;
     908             :       }
     909       10815 :       if ((is_f_square16[20] = squares16[(cmp[degree]+4*cmp[degree-1])  & 0xf]))
     910        4585 :         np2++;
     911       10815 :       if ((is_f_square16[21] = squares16[(cmp[degree]+12*cmp[degree-1]) & 0xf]))
     912        4571 :         np2++;
     913       10815 :       if ((is_f_square16[22] = squares16[(cmp[degree]+8*cmp[degree-1])  & 0xf]))
     914        4438 :         np3++;
     915       10815 :       if ((is_f_square16[23] = squares16[cmp[degree]])) np4++;
     916             :     }
     917             : 
     918             :     /* set den_bits */
     919       12194 :     { ulong db = 0;
     920             :       long i;
     921             : 
     922       12194 :       if (npe + npo > 0) db |= 0xaaaaUL; /* odd denominators */
     923       12194 :       if (np1 > 0)       db |= 0x4444UL; /* v_2(den) = 1 */
     924       12194 :       if (np2 > 0)       db |= 0x1010UL; /* v_2(den) = 2 */
     925       12194 :       if (np3 > 0)       db |= 0x0100UL; /* v_2(den) = 3 */
     926       12194 :       if (np4 > 0)       db |= 0x0001UL; /* v_2(den) >= 4 */
     927       12194 :       if (db == 0)
     928             :       {
     929        2975 :         for (i = 0 ; i < 16; i++) num_bits[i] = RBA(0UL);
     930         175 :         *den_bits = 0UL; return num_none;
     931             :       }
     932       34340 :       for (i = 16; i < BITS_IN_LONG; i <<= 1) db |= db << i;
     933       12019 :       *den_bits = db;
     934             :     }
     935       12019 :     result = (npe == 0) ? (npo == 0) ? num_none : num_odd
     936       12019 :                         : (npo == 0) ? num_even : num_all;
     937             :   }
     938             : 
     939             :   /* set up num_bits[16] */
     940             : 
     941             :   /* odd denominators */
     942       12019 :   switch(result)
     943             :   {
     944        7805 :     case num_all:
     945       70245 :       for (b = 1; b < 16; b += 2)
     946             :       {
     947       62440 :         ulong work = 0, bit = 1;
     948       62440 :         long i, invb = b; /* inverse of b mod 16 */
     949       62440 :         if (b & 2) invb ^= 8;
     950       62440 :         if (b & 4) invb ^= 8;
     951     1061480 :         for (i = 0; i < 16; i++)
     952             :         {
     953      999040 :           if (is_f_square16[(invb*i) & 0xf]) work |= bit;
     954      999040 :           bit <<= 1;
     955             :         }
     956             :         /* now repeat the 16 bits */
     957      178400 :         for (i = 16; i < BITS_IN_LONG; i <<= 1) work |= work << i;
     958       62440 :         num_bits[b] = RBA(work);
     959             :       }
     960        7805 :       break;
     961             : 
     962        1820 :     case num_odd:
     963       16380 :       for (b = 1; b < 16; b += 2)
     964             :       {
     965       14560 :         ulong work = 0, bit = 1;
     966       14560 :         long i, invb = b; /* inverse of b mod 16 */
     967       14560 :         if (b & 2) invb ^= 8;
     968       14560 :         if (b & 4) invb ^= 8;
     969      131040 :         for (i = 1; i < 16; i += 2)
     970             :         {
     971      116480 :           if (is_f_square16[(invb*i) & 0xf]) work |= bit;
     972      116480 :           bit <<= 1;
     973             :         }
     974             :         /* now repeat the 8 bits */
     975       56160 :         for (i = 8; i < BITS_IN_LONG; i <<= 1) { work |= work << i; }
     976       14560 :         num_bits[b] = RBA(work);
     977             :       }
     978        1820 :       break;
     979             : 
     980        2002 :     case num_even:
     981       18018 :       for (b = 1; b < 16; b += 2)
     982             :       {
     983       16016 :         ulong work = 0, bit = 1;
     984       16016 :         long i, invb = b; /* inverse of b mod 16 */
     985       16016 :         if (b & 2) invb ^= 8;
     986       16016 :         if (b & 4) invb ^= 8;
     987      144144 :         for (i = 0; i < 16; i += 2)
     988             :         {
     989      128128 :           if (is_f_square16[(invb*i) & 0xf]) work |= bit;
     990      128128 :           bit <<= 1;
     991             :         }
     992             :         /* now repeat the 8 bits */
     993       61776 :         for (i = 8; i < BITS_IN_LONG; i <<= 1) work |= work << i;
     994       16016 :         num_bits[b] = RBA(work);
     995             :       }
     996        2002 :       break;
     997             : 
     998         392 :     case num_none:
     999        3528 :       for (b = 1; b < 16; b += 2) num_bits[b] = RBA(0UL);
    1000         392 :       break;
    1001             :   }
    1002             : 
    1003             :   /* v_2(den) = 1 : only odd numerators */
    1004       60095 :   for (b = 1; b < 8; b += 2)
    1005             :   {
    1006       48076 :     ulong work = 0, bit = 1;
    1007             :     long i;
    1008      432684 :     for (i = 1; i < 16; i += 2)
    1009             :     {
    1010      384608 :       if (is_f_square16[16 + (((b*i)>>1) & 0x3)]) work |= bit;
    1011      384608 :       bit <<= 1;
    1012             :     }
    1013             :     /* now repeat the 8 bits */
    1014      185436 :     for (i = 8; i < BITS_IN_LONG; i <<= 1) work |= work << i;
    1015       48076 :     num_bits[2*b] = RBA(work);
    1016             :   }
    1017             : 
    1018             :   /* v_2(den) = 2 : only odd numerators */
    1019       36057 :   for (b = 1; b < 4; b += 2)
    1020             :   {
    1021       24038 :     ulong work = 0, bit = 1;
    1022             :     long i;
    1023      120190 :     for (i = 1; i < 8; i += 2)
    1024             :     {
    1025       96152 :       if (is_f_square16[20 + (((b*i)>>1) & 0x1)]) work |= bit;
    1026       96152 :       bit <<= 1;
    1027             :     }
    1028             :     /* now repeat the 4 bits */
    1029      116756 :     for (i = 4; i < BITS_IN_LONG; i <<= 1) work |= work << i;
    1030       24038 :     num_bits[4*b] = RBA(work);
    1031             :   }
    1032             : 
    1033             :   /* v_2(den) = 3, >= 4 : only odd numerators */
    1034       12019 :   num_bits[8] = (is_f_square16[22]) ? RBA(~(0UL)) : RBA(0UL);
    1035       12019 :   num_bits[0] = (is_f_square16[23]) ? RBA(~(0UL)) : RBA(0UL);
    1036             : 
    1037       12019 :   return result;
    1038             : }
    1039             : 
    1040             : /**************************************************************************
    1041             :  * This is a comparison function needed for sorting in order to determine *
    1042             :  * the `best' primes for sieving.                                         *
    1043             :  **************************************************************************/
    1044             : 
    1045             : static int
    1046     1151612 : compare_entries(const void *a, const void *b)
    1047             : {
    1048     1151612 :   double diff = ((entry *)a)->r - ((entry *)b)->r;
    1049     1151612 :   return (diff > 0) ? 1 : (diff < 0) ? -1 : 0;
    1050             : }
    1051             : 
    1052             : /************************************************************************
    1053             :  * Collect the sieving information                                      *
    1054             :  ************************************************************************/
    1055             : 
    1056             : static long
    1057       12194 : sieving_info(ratpoints_args *args,
    1058             :              ratpoints_sieve_entry **sieve_list)
    1059             : {
    1060       12194 :   GEN c = args->cof;
    1061       12194 :   GEN inverses = args->inverses, squares;
    1062       12194 :   GEN offsets = args->offsets;
    1063       12194 :   ulong ** sieves0 = args->sieves0;
    1064       12194 :   long degree = degpol(c);
    1065       12194 :   long fba = 0, fdc = 0;
    1066       12194 :   long pn, pnp = 0;
    1067             :   long n;
    1068       12194 :   long nbprime = lg(args->listprime)-1;
    1069       12194 :   entry *prec = (entry*) stack_malloc(nbprime*sizeof(entry));
    1070             :     /* This array is used for sorting in order to
    1071             :        determine the `best' sieving primes. */
    1072             : 
    1073       12194 :   forbidden_entry *forb_ba = args->forb_ba;
    1074       12194 :   long *forbidden = args->forbidden;
    1075       12194 :   ulong bound = (1UL)<<(BITS_IN_LONG - args->bit_primes);
    1076       12194 :   pari_sp av = avma;
    1077       12194 :   squares = gen_squares(args->listprime);
    1078             : 
    1079             :   /* initialize sieve in se_buffer */
    1080      374283 :   for (pn = 1; pn <= args->num_primes; pn++)
    1081      362215 :   {
    1082      362215 :     long coeffs_mod_p[degree+1]; /* The coefficients of f reduced modulo p */
    1083      362215 :     ulong a, p = args->listprime[pn], np;
    1084             :     long n;
    1085      362215 :     int *is_f_square = args->int_next;
    1086             : 
    1087      362215 :     args->int_next += p + 1; /* need space for (p+1) int's */
    1088             : 
    1089             :     /* compute coefficients mod p */
    1090     2545550 :     for (n = 0; n <= degree; n++) coeffs_mod_p[n] = umodiu(pel(c,n), p);
    1091             : 
    1092      362215 :     np = umael(squares,pn,coeffs_mod_p[0]+1);
    1093      362215 :     is_f_square[0] = np;
    1094    20733475 :     for (a = 1 ; a < p; a++)
    1095             :     {
    1096    20371260 :       ulong s = coeffs_mod_p[degree];
    1097    20371260 :       if ((degree+1)*args->bit_primes <= BITS_IN_LONG)
    1098             :       {
    1099   106091968 :         for (n = degree - 1 ; n >= 0 ; n--) s = s*a + coeffs_mod_p[n];
    1100             :         /* here, s < p^(degree+1) <= max. long */
    1101    17673768 :         s %= p;
    1102             :       }
    1103             :       else
    1104             :       {
    1105    16691420 :         for (n = degree - 1 ; n >= 0 ; n--)
    1106             :         {
    1107    13993928 :           s = s*a + coeffs_mod_p[n];
    1108    13993928 :           if (s+1 >= bound) s %= p;
    1109             :         }
    1110     2697492 :         s %= p;
    1111             :       }
    1112    20371260 :       if ((is_f_square[a] = mael(squares,pn,s+1))) np++;
    1113             :     }
    1114      362215 :     is_f_square[p] = odd(degree) || mael(squares,pn,coeffs_mod_p[degree]+1);
    1115             : 
    1116             :     /* check if there are no solutions mod p */
    1117      362215 :     if (np == 0 && !is_f_square[p]) return gc_long(av,p);
    1118             : 
    1119             :     /* Fill arrays with info for p */
    1120      362089 :     if (np < p)
    1121             :     { /* only when there is some information */
    1122             :       ulong i;
    1123      330722 :       ratpoints_sieve_entry *se = args->se_next;
    1124      847476 :       double r = is_f_square[p] ? ((double)(np*(p-1) + p))/((double)(p*p))
    1125      330722 :                                 : (double)np/(double)p;
    1126      330722 :       prec[pnp].r = r;
    1127      330722 :       args->se_next ++;
    1128      330722 :       se->p = p;
    1129      330722 :       se->is_f_square = is_f_square;
    1130      330722 :       se->inverses = gel(inverses,pn);
    1131      330722 :       se->offset = offsets[pn];
    1132      330722 :       se->sieve[0] = (ratpoints_bit_array *)sieves0[pn];
    1133    19885642 :       for (i = 1; i < p; i++) se->sieve[i] = NULL;
    1134      330722 :       prec[pnp].ssp = se;
    1135      330722 :       pnp++;
    1136             :     }
    1137             : 
    1138      362089 :     if ((args->flags & RATPOINTS_CHECK_DENOM)
    1139      316309 :          && fba + fdc < args->max_forbidden
    1140      316309 :          && !is_f_square[p])
    1141             :     { /* record forbidden divisors of the denominator */
    1142      146146 :       if (coeffs_mod_p[degree] == 0)
    1143             :       { /* leading coeff. divisible by p */
    1144             :         GEN r;
    1145           0 :         long v = Z_lvalrem(pel(c,degree), p, &r);
    1146             : 
    1147           0 :         if (odd(v) || !mael(squares,pn, umodiu(r,p)+1))
    1148             :         { /* Can only get something when valuation is odd
    1149             :              or when valuation is even and lcf is not a p-adic square.
    1150             :              Compute smallest n such that if v(den) >= n, the leading
    1151             :              term determines the valuation. Then we must have v(den) < n. */
    1152           0 :           long k, n = 1;
    1153           0 :           for (k = degree-1; k >= 0; k--)
    1154             :           {
    1155           0 :             if (coeffs_mod_p[k] == 0)
    1156             :             {
    1157           0 :               long t = 1 + v - Z_lval(pel(c,k), p);
    1158           0 :               long m = CEIL(t, (degree-k));
    1159           0 :               if (m > n) n = m;
    1160             :             }
    1161             :           }
    1162           0 :           if (n == 1)
    1163             :           {
    1164           0 :             forb_ba[fba].p     = p;
    1165           0 :             forb_ba[fba].start = sieves0[pn];
    1166           0 :             forb_ba[fba].end   = sieves0[pn]+p;
    1167           0 :             forb_ba[fba].curr  = forb_ba[fba].start;
    1168           0 :             fba++;
    1169             :           }
    1170             :           else
    1171             :           {
    1172           0 :             forbidden[fdc] = upowuu(p, n);
    1173           0 :             fdc++;
    1174             :           }
    1175             :         }
    1176             :       }
    1177             :       else /* leading coefficient is a nonsquare mod p */
    1178             :       { /* denominator divisible by p is excluded */
    1179      146146 :         forb_ba[fba].p     = p;
    1180      146146 :         forb_ba[fba].start = sieves0[pn];
    1181      146146 :         forb_ba[fba].end   = sieves0[pn]+p;
    1182      146146 :         forb_ba[fba].curr  = forb_ba[fba].start;
    1183      146146 :         fba++;
    1184             :       }
    1185             :     }
    1186             :   } /* end for pn */
    1187             : 
    1188             :   /* update sp2 and sp1 if necessary */
    1189       12068 :   if (args->sp2 > pnp)       args->sp2 = pnp;
    1190       12068 :   if (args->sp1 > args->sp2) args->sp1 = args->sp2;
    1191             : 
    1192             :   /* sort the array to get at the best primes */
    1193       12068 :   qsort(prec, pnp, sizeof(entry), compare_entries);
    1194             : 
    1195             :   /* put the sorted entries into sieve_list */
    1196      240982 :   for (n = 0; n < args->sp2; n++) sieve_list[n] = prec[n].ssp;
    1197             : 
    1198             :   /* terminate array of forbidden divisors */
    1199       12068 :   if (args->flags & RATPOINTS_CHECK_DENOM)
    1200             :   {
    1201             :     long n;
    1202             : 
    1203       10542 :     for (n = args->num_primes+1;
    1204       10542 :         fba + fdc < args->max_forbidden && n <= nbprime; n++)
    1205             :     {
    1206           0 :       ulong p = args->listprime[n];
    1207             : 
    1208           0 :       if (p*p > (ulong) args->b_high) break;
    1209           0 :       if (kroiu(pel(c,degree), p) == -1)
    1210             :       {
    1211           0 :         forb_ba[fba].p     = p;
    1212           0 :         forb_ba[fba].start = sieves0[n];
    1213           0 :         forb_ba[fba].end   = sieves0[n]+p;
    1214           0 :         forb_ba[fba].curr  = forb_ba[fba].start;
    1215           0 :         fba++;
    1216             :       }
    1217             :     }
    1218       10542 :     forb_ba[fba].p = 0; /* terminating zero */
    1219       10542 :     forbidden[fdc] = 0; /* terminating zero */
    1220       10542 :     args->max_forbidden = fba + fdc; /* note actual number */
    1221             :   }
    1222             : 
    1223       12068 :   if (fba + fdc == 0) args->flags &= ~RATPOINTS_CHECK_DENOM;
    1224       12068 :   return gc_long(av,0);
    1225             : }
    1226             : 
    1227             : /**************************************************************************
    1228             :  * The sieving procedure itself                                           *
    1229             :  **************************************************************************/
    1230             : static void
    1231    29882944 : sift(long b, ratpoints_bit_array *survivors, ratpoints_args *args,
    1232             :      bit_selection which_bits, ratpoints_bit_array bits16,
    1233             :      ratpoints_sieve_entry **sieve_list, long *bp_list, int *quit,
    1234             :      int process(long, long, GEN, void*, int*), void *info)
    1235    29882944 : {
    1236    29882944 :   pari_sp av = avma;
    1237    29882944 :   sieve_spec ssp[args->sp2];
    1238    29882944 :   int do_setup = 1;
    1239    29882944 :   long k, height = args->height, nb;
    1240             : 
    1241    29882944 :   if (odd(b) == 0) which_bits = num_odd; /* even denominator */
    1242             : 
    1243             :   /* Note that b is new */
    1244    29882944 :   args->bc = NULL;
    1245    29882944 :   nb = 0;
    1246             : 
    1247    70818672 :   for (k = 0; k < args->num_inter; k++)
    1248             :   {
    1249    45115498 :     ratpoints_interval inter = args->domain[k];
    1250             :     long low, high;
    1251             : 
    1252             :     /* Determine relevant interval [low, high] of numerators. */
    1253    45115498 :     if (b*inter.low <= -height)
    1254    21912477 :       low = -height;
    1255             :     else
    1256             :     {
    1257    23203021 :       if (b*inter.low > height) break;
    1258    19025008 :       low = ceil(b*inter.low);
    1259             :     }
    1260    40937485 :     if (b*inter.up >= height)
    1261    18063284 :       high = height;
    1262             :     else
    1263             :     {
    1264    22874201 :       if (b*inter.up < -height) continue;
    1265    19299126 :       high = floor(b*inter.up);
    1266             :     }
    1267             : 
    1268    37362410 :     if (do_setup)
    1269             :     { /* set up the sieve information */
    1270             :       long n;
    1271             : 
    1272    27677006 :       do_setup = 0; /* only do it once for every b */
    1273   553150766 :       for (n = 0; n < args->sp2; n++)
    1274             :       {
    1275   525473760 :         ratpoints_sieve_entry *se = sieve_list[n];
    1276   525473760 :         long p = se->p;
    1277   525473760 :         long bp = bp_list[n];
    1278             :         ratpoints_bit_array *sptr;
    1279             : 
    1280   525473760 :         if (which_bits != num_all) /* divide by 2 mod p */
    1281   302566558 :           bp = odd(bp) ? (bp+p) >> 1 : bp >> 1;
    1282   525473760 :         sptr = se->sieve[bp];
    1283             : 
    1284   525473760 :         ssp[n].p = p;
    1285   525473760 :         ssp[n].offset = (which_bits == num_odd) ? se->offset : 0;
    1286             : 
    1287             :         /* copy if already initialized, else initialize */
    1288   537891900 :         ssp[n].ptr = sptr ? sptr : (p<BITS_IN_LONG?sieve_init1(p, se, bp, args)
    1289    12418140 :                                                   :sieve_init2(p, se, bp, args));
    1290   525473760 :         ssp[n].start = ssp[n].ptr;
    1291   525473760 :         ssp[n].end = ssp[n].ptr + p;
    1292             : 
    1293             :       }
    1294             :     }
    1295             : 
    1296    37362410 :     switch(which_bits)
    1297             :     {
    1298    15961156 :       case num_all: break;
    1299           0 :       case num_none: break;
    1300    16775066 :       case num_odd: low >>= 1; high--; high >>= 1; break;
    1301     4626188 :       case num_even: low++; low >>= 1; high >>= 1; break;
    1302             :     }
    1303             : 
    1304             :     /* now turn the bit interval into [low, high[ */
    1305    37362410 :     high++;
    1306             : 
    1307    37362410 :     if (low < high)
    1308             :     {
    1309    37361143 :       long w_low, w_high, w_low0, w_high0, range = args->array_size;
    1310             : 
    1311             :       /* Now the range of longwords (= bit_arrays) */
    1312    37361143 :       w_low = low >> RBA_SHIFT;
    1313    37361143 :       w_high = (high + (long)(RBA_LENGTH-1)) >> RBA_SHIFT;
    1314    37361143 :       w_low0 = w_low;
    1315    37361143 :       w_high0 = w_low0 + range;
    1316    79857236 :       for ( ; w_low0 < w_high; w_low0 = w_high0, w_high0 += range)
    1317             :       {
    1318             :         long i;
    1319    42497850 :         if (w_high0 > w_high) { w_high0 = w_high; range = w_high0 - w_low0; }
    1320             :         /* initialise the bits */
    1321  4592675083 :         for (i = range; i; i--) survivors[i-1] = bits16;
    1322             :         /* boundary words */
    1323    42497850 :         if (w_low0 == w_low)
    1324    37361143 :           MASKL(survivors,low - RBA_LENGTH * w_low)
    1325    42497850 :         if (w_high0 == w_high)
    1326    37360996 :           MASKU(&survivors[range-1], RBA_LENGTH * w_high - high)
    1327             : 
    1328             : #if (RATPOINTS_CHUNK > 1)
    1329   229793976 :         while(range%RATPOINTS_CHUNK != 0)
    1330   197153496 :           { survivors[range] = RBA(0); range++; w_high0++; }
    1331             : #endif
    1332    42497850 :         nb += _ratpoints_sift0(b, w_low0, w_high0, args, which_bits,
    1333             :                          survivors, &ssp[0], quit, process, info);
    1334    42497850 :         if (*quit) return;
    1335             :       }
    1336             :     }
    1337             :   }
    1338    29881187 :   if (nb==0) set_avma(av);
    1339             : }
    1340             : 
    1341             : /**************************************************************************
    1342             :  * Find points by looping over the denominators and sieving numerators    *
    1343             :  **************************************************************************/
    1344             : static void
    1345       12194 : find_points_work(ratpoints_args *args,
    1346             :                  int process(long, long, GEN, void*, int*), void *info)
    1347             : {
    1348       12194 :   int quit = 0;
    1349       12194 :   GEN c = args->cof;
    1350       12194 :   long degree = degpol(c);
    1351       12194 :   long nbprime = lg(args->listprime)-1;
    1352       12194 :   long height = args->height;
    1353             : 
    1354       12194 :   int point_at_infty = 0; /* indicates if there are points at infinity */
    1355       12194 :   int lcfsq = Z_issquare(pel(c,degree));
    1356             : 
    1357       12194 :   forbidden_entry *forb_ba = args->forb_ba;
    1358       12194 :   long *forbidden = args->forbidden;
    1359             :     /* The forbidden divisors, a zero-terminated array.
    1360             :        Used when degree is even and leading coefficient is not a square */
    1361             : 
    1362             :     /* These are used when degree is odd and leading coeff. is not +-1 */
    1363             : 
    1364             :   ratpoints_sieve_entry **sieve_list = (ratpoints_sieve_entry**)
    1365       12194 :      stack_malloc(nbprime*sizeof(ratpoints_sieve_entry*));
    1366       12194 :   bit_selection which_bits = num_all;
    1367             :   ulong den_bits;
    1368             :   ratpoints_bit_array num_bits[16];
    1369             : 
    1370       12194 :   args->flags &= RATPOINTS_FLAGS_INPUT_MASK;
    1371       12194 :   args->flags |= RATPOINTS_CHECK_DENOM;
    1372             : 
    1373             :   /* initialize memory management */
    1374       12194 :   args->se_next = args->se_buffer;
    1375       12194 :   args->ba_next = args->ba_buffer;
    1376       12194 :   args->int_next = args->int_buffer;
    1377             : 
    1378             :   /* Some sanity checks */
    1379       12194 :   args->num_inter = 0;
    1380             : 
    1381       12194 :   if (args->num_primes > nbprime) args->num_primes = nbprime;
    1382       12194 :   if (args->sp2 > args->num_primes) args->sp2 = args->num_primes;
    1383       12194 :   if (args->sp1 > args->sp2)        args->sp1 = args->sp2;
    1384             : 
    1385       12194 :   if (args->b_low < 1)  args->b_low = 1;
    1386       12194 :   if (args->b_high < 1) args->b_high = height;
    1387       12194 :   if (args->max_forbidden < 0)
    1388           0 :     args->max_forbidden = RATPOINTS_DEFAULT_MAX_FORBIDDEN;
    1389       12194 :   if (args->max_forbidden > nbprime)
    1390           0 :     args->max_forbidden = nbprime;
    1391       12194 :   if (args->array_size <= 0) args->array_size = RATPOINTS_ARRAY_SIZE;
    1392             :   {
    1393       12194 :     long s = 2*maxss(1,CEIL(height, BITS_IN_LONG));
    1394       12194 :     if (args->array_size > s) args->array_size = s;
    1395             :   }
    1396             :   /* make sure that array size is a multiple of RATPOINTS_CHUNK */
    1397       12194 :   args->array_size = CEIL(args->array_size, RATPOINTS_CHUNK)*RATPOINTS_CHUNK;
    1398             : 
    1399             :   /* Don't reverse if intervals are specified or limits for the denominator
    1400             :      are given */
    1401       12194 :   if (args->num_inter > 0 || args->b_low > 1 || args->b_high != height)
    1402          35 :     args->flags |= RATPOINTS_NO_REVERSE;
    1403             : 
    1404             :   /* Check if reversal of polynomial might be better:
    1405             :    * case 1: degree is even, but trailing coefficient is zero
    1406             :    * case 2: degree is even, leading coefficient is a square, but
    1407             :    *         trailing coefficient is not
    1408             :    * case 3: degree is odd, |leading coefficient| > 1,
    1409             :    *         trailing coefficient is zero, |coeff. of x| = 1 */
    1410       12194 :   if (!(args->flags & RATPOINTS_NO_REVERSE))
    1411             :   {
    1412       12159 :     if (!odd(degree))
    1413             :     {
    1414       11018 :       if (signe(pel(c,0)) == 0)
    1415             :       { /* case 1 */
    1416             :         long n;
    1417         224 :         args->flags |= RATPOINTS_REVERSED;
    1418         896 :         for (n = 0; n < degree>>1; n++) swap(pel(c,n), pel(c,degree-n));
    1419         224 :         degree--;
    1420         224 :         setlg(c,degree+3);
    1421             :       }
    1422       10794 :       else if (lcfsq && !Z_issquare(pel(c,0)))
    1423             :       { /* case 2 */
    1424             :         long n;
    1425         735 :         args->flags |= RATPOINTS_REVERSED;
    1426        2940 :         for (n = 0; n < degree>>1; n++) swap(pel(c,n), pel(c,degree-n));
    1427         735 :         lcfsq = 0;
    1428             :       }
    1429             :     }
    1430             :     else
    1431             :     { /* odd degree, case 3*/
    1432        1141 :       if (!is_pm1(pel(c,degree)) && !signe(pel(c,0)) && is_pm1(pel(c,1)))
    1433             :       {
    1434             :         long n;
    1435           7 :         args->flags |= RATPOINTS_REVERSED;
    1436          14 :         for (n = 1; n < degree>>1; n++) swap(pel(c,n),pel(c,degree+1-n));
    1437             :       }
    1438             :     }
    1439             :   }
    1440             : 
    1441             :   /* Deal with the intervals */
    1442       12194 :   if (args->num_inter == 0)
    1443             :   { /* default interval (effectively ]-oo,oo[) if none is given */
    1444       12194 :     args->domain = (ratpoints_interval*) stack_malloc(2*degree*sizeof(ratpoints_interval));
    1445       12194 :     args->domain[0].low = -height; args->domain[0].up = height;
    1446       12194 :     args->num_inter = 1;
    1447             :   }
    1448             : 
    1449       12194 :   ratpoints_compute_sturm(args);
    1450             : 
    1451             :   /* Point(s) at infinity? */
    1452       12194 :   if (odd(degree) || lcfsq)
    1453             :   {
    1454        1526 :     args->flags &= ~RATPOINTS_CHECK_DENOM;
    1455        1526 :     point_at_infty = 1;
    1456             :   }
    1457             : 
    1458             :   /* Can use only squares as denoms if degree is odd and poly is +-monic */
    1459       12194 :   if (odd(degree))
    1460             :   {
    1461        1379 :     GEN w1 = pel(c,degree);
    1462        1379 :     if (is_pm1(w1))
    1463          63 :       args->flags |= RATPOINTS_USE_SQUARES;
    1464             :     else /* set up information on divisors of leading coefficient */
    1465        1316 :       setup_us1(args, absi_shallow(w1));
    1466             :   }
    1467             : 
    1468             :   /* deal with f mod powers of 2 */
    1469       12194 :   which_bits = get_2adic_info(args, &den_bits, &num_bits[0]);
    1470             :   /* which_bits says whether to consider even and/or odd numerators
    1471             :    * when the denominator is odd.
    1472             :    *
    1473             :    * Bit k in den_bits is 0 if b congruent to k mod BITS_IN_LONG need
    1474             :    * not be considered as a denominator.
    1475             :    *
    1476             :    * Bit k in num_bits[b] is 0 is numerators congruent to
    1477             :    *  k (which_bits = den_all)
    1478             :    *  2k (which_bits = den_even)
    1479             :    *  2k+1 (which_bits = den_odd)
    1480             :    * need not be considered for denominators congruent to b mod 16. */
    1481             : 
    1482             :   /* set up the sieve data structure */
    1483       12194 :   if (sieving_info(args, sieve_list)) return;
    1484             : 
    1485             :   /* deal with point(s) at infinity */
    1486       12068 :   if (point_at_infty)
    1487             :   {
    1488        1526 :     long a = 1, b = 0;
    1489             : 
    1490        1526 :     if (args->flags & RATPOINTS_REVERSED) { a = 0; b = 1; }
    1491        1526 :     if (odd(degree))
    1492        1379 :       (void)process(a, b, gen_0, info, &quit);
    1493             :     else
    1494             :     {
    1495         147 :       GEN w0 = sqrti(pel(c,degree));
    1496         147 :       (void)process(a, b, w0, info, &quit);
    1497         147 :       (void)process(a, b, negi(w0), info, &quit);
    1498             :     }
    1499        1526 :     if (quit) return;
    1500             :   }
    1501             :   /* now do the sieving */
    1502             :   {
    1503             :     ratpoints_bit_array *survivors = (ratpoints_bit_array *)
    1504       12068 :       stack_malloc_align((args->array_size)*RBA_SIZE, RBA_SIZE);
    1505       12068 :     if (args->flags & (RATPOINTS_USE_SQUARES | RATPOINTS_USE_SQUARES1))
    1506             :     {
    1507        1379 :       if (args->flags & RATPOINTS_USE_SQUARES)
    1508             :       /* need only take squares as denoms */
    1509          63 :       {
    1510             :         long b, bb;
    1511          63 :         long bp_list[args->sp2];
    1512          63 :         long last_b = args->b_low;
    1513             :         long n;
    1514        1260 :         for (n = 0; n < args->sp2; n++)
    1515        1197 :           bp_list[n] = mod(args->b_low, sieve_list[n]->p);
    1516             : 
    1517        7875 :         for (b = 1; bb = b*b, bb <= args->b_high; b++)
    1518        7812 :           if (bb >= args->b_low)
    1519             :           {
    1520        7812 :             ratpoints_bit_array bits = num_bits[bb & 0xf];
    1521        7812 :             if (TEST(bits))
    1522             :             {
    1523             :               long n;
    1524        6916 :               long d = bb - last_b;
    1525             : 
    1526             :               /* fill bp_list */
    1527      138320 :               for (n = 0; n < args->sp2; n++)
    1528      131404 :                 bp_list[n] = mod(bp_list[n] + d, sieve_list[n]->p);
    1529        6916 :               last_b = bb;
    1530             : 
    1531        6916 :               sift(bb, survivors, args, which_bits, bits,
    1532             :                    sieve_list, &bp_list[0], &quit, process, info);
    1533        6916 :               if (quit) break;
    1534             :             }
    1535             :           }
    1536             :       }
    1537             :       else /* args->flags & RATPOINTS_USE_SQUARES1 */
    1538        1316 :       {
    1539        1316 :         GEN den_info = args->den_info;
    1540        1316 :         GEN divisors = args->divisors;
    1541        1316 :         long ld = lg(divisors);
    1542             :         long k;
    1543             :         long b, bb;
    1544        1316 :         long bp_list[args->sp2];
    1545             : 
    1546        4186 :         for (k = 1; k < ld; k++)
    1547             :         {
    1548        2877 :           long d = divisors[k];
    1549        2877 :           long last_b = d;
    1550             :           long n;
    1551        2877 :           if (d > args->b_high) continue;
    1552       57400 :           for (n = 0; n < args->sp2; n++)
    1553       54530 :             bp_list[n] = mod(d, sieve_list[n]->p);
    1554             : 
    1555      273077 :           for (b = 1; bb = d*b*b, bb <= args->b_high; b++)
    1556             :           {
    1557      270214 :             if (bb >= args->b_low)
    1558             :             {
    1559      270193 :               int flag = 1;
    1560      270193 :               ratpoints_bit_array bits = num_bits[bb & 0xf];
    1561             : 
    1562      270193 :               if (EXT0(bits))
    1563             :               {
    1564      223076 :                 long i, n, l = lg(gel(den_info,1));
    1565      223076 :                 long d = bb - last_b;
    1566             : 
    1567             :                 /* fill bp_list */
    1568     4461520 :                 for (n = 0; n < args->sp2; n++)
    1569     4238444 :                   bp_list[n] = mod(bp_list[n] + d, sieve_list[n]->p);
    1570      223076 :                 last_b = bb;
    1571             : 
    1572      423668 :                 for(i = 1; i < l; i++)
    1573             :                 {
    1574      251363 :                   long v = z_lval(bb, mael(den_info,1,i));
    1575      251363 :                   if((v >= mael(den_info,3,i))
    1576      121114 :                       && odd(v + (mael(den_info,2,i)))) { flag = 0; break; }
    1577             :                 }
    1578      223076 :                 if (flag)
    1579             :                 {
    1580      172305 :                   sift(bb, survivors, args, which_bits, bits,
    1581             :                        sieve_list, &bp_list[0], &quit, process, info);
    1582      172305 :                   if (quit) break;
    1583             :                 }
    1584             :               }
    1585             :             }
    1586             :           }
    1587        2870 :           if (quit) break;
    1588             :         }
    1589             :       }
    1590             :     }
    1591             :     else
    1592             :     {
    1593       10689 :       if (args->flags & RATPOINTS_CHECK_DENOM)
    1594       10542 :       {
    1595             :         long *forb;
    1596             :         long b;
    1597       10542 :         long bp_list[args->sp2];
    1598       10542 :         long last_b = args->b_low;
    1599             :         ulong b_bits;
    1600             :         long n;
    1601      210462 :         for (n = 0; n < args->sp2; n++)
    1602      199920 :           bp_list[n] = mod(args->b_low, sieve_list[n]->p);
    1603             :         {
    1604       10542 :           forbidden_entry *fba = &forb_ba[0];
    1605       10542 :           long b_low = args->b_low;
    1606       10542 :           long w_low = (b_low-1) >> TWOPOTBITS_IN_LONG;
    1607             : 
    1608       10542 :           b_bits = den_bits;
    1609      156674 :           while (fba->p)
    1610             :           {
    1611      146132 :             fba->curr = fba->start + mod(w_low, fba->p);
    1612      146132 :             b_bits &= *(fba->curr);
    1613      146132 :             fba++;
    1614             :           }
    1615       10542 :           b_bits >>= (b_low-1) & LONG_MASK;
    1616             :         }
    1617             : 
    1618   131962383 :         for (b = args->b_low; b <= args->b_high; b++)
    1619             :         {
    1620   131953584 :           ratpoints_bit_array bits = num_bits[b & 0xf];
    1621             : 
    1622   131953584 :           if ((b & LONG_MASK) == 0)
    1623             :           { /* next b_bits */
    1624     2348804 :             forbidden_entry *fba = &forb_ba[0];
    1625             : 
    1626     2348804 :             b_bits = den_bits;
    1627    37228876 :             while (fba->p)
    1628             :             {
    1629    34880072 :               fba->curr++;
    1630    34880072 :               if (fba->curr == fba->end) fba->curr = fba->start;
    1631    34880072 :               b_bits &= *(fba->curr);
    1632    34880072 :               fba++;
    1633             :             }
    1634             :           }
    1635             :           else
    1636   129604780 :             b_bits >>= 1;
    1637             : 
    1638   131953584 :           if (odd(b_bits) && EXT0(bits))
    1639             :           { /* check if denominator is excluded */
    1640    48413876 :             for (forb = &forbidden[0] ; *forb && (b % (*forb)); forb++)
    1641           0 :               continue;
    1642    48413876 :             if (*forb == 0 && rpjacobi(b, pel(c,degree)) == 1)
    1643             :             {
    1644    28026495 :               long n, d = b - last_b;
    1645             : 
    1646             :               /* fill bp_list */
    1647   560140546 :               for (n = 0; n < args->sp2; n++)
    1648             :               {
    1649   532114051 :                 long bp = bp_list[n] + d;
    1650   532114051 :                 long p = sieve_list[n]->p;
    1651             : 
    1652   601685595 :                 while (bp >= p) bp -= p;
    1653   532114051 :                 bp_list[n] = bp;
    1654             :               }
    1655    28026495 :               last_b = b;
    1656             : 
    1657    28026495 :               sift(b, survivors, args, which_bits, bits,
    1658             :                    sieve_list, &bp_list[0], &quit, process, info);
    1659    28026495 :               if (quit) break;
    1660             :             }
    1661             :           }
    1662             :         }
    1663             :       } /* if (args->flags & RATPOINTS_CHECK_DENOM) */
    1664             :       else
    1665         147 :       {
    1666             :         long b, n;
    1667         147 :         long bp_list[args->sp2];
    1668         147 :         long last_b = args->b_low;
    1669        2940 :         for (n = 0; n < args->sp2; n++)
    1670        2793 :           bp_list[n] = mod(args->b_low, sieve_list[n]->p);
    1671     2179156 :         for (b = args->b_low; b <= args->b_high; b++)
    1672             :         {
    1673     2179016 :           ratpoints_bit_array bits = num_bits[b & 0xf];
    1674     2179016 :           if (EXT0(bits))
    1675             :           {
    1676             :             long n;
    1677     1677228 :             long d = b - last_b;
    1678             : 
    1679             :             /* fill bp_list */
    1680    33544560 :             for (n = 0; n < args->sp2; n++)
    1681             :             {
    1682    31867332 :               long bp = bp_list[n] + d;
    1683    31867332 :               long p = sieve_list[n]->p;
    1684             : 
    1685    32980773 :               while (bp >= p) bp -= p;
    1686    31867332 :               bp_list[n] = bp;
    1687             :             }
    1688     1677228 :             last_b = b;
    1689             : 
    1690     1677228 :             sift(b, survivors, args, which_bits, bits,
    1691             :                  sieve_list, &bp_list[0], &quit, process, info);
    1692     1677228 :             if (quit) break;
    1693             :           }
    1694             :         }
    1695             :       }
    1696             :     }
    1697             :   }
    1698             : }
    1699             : 
    1700             : static GEN
    1701       85659 : vec_append_grow(GEN z, long i, GEN x)
    1702             : {
    1703       85659 :   long n = lg(z)-1;
    1704       85659 :   if (i > n)
    1705             :   {
    1706        1414 :     n <<= 1;
    1707        1414 :     z = vec_lengthen(z,n);
    1708             :   }
    1709       85659 :   gel(z,i) = x;
    1710       85659 :   return z;
    1711             : }
    1712             : 
    1713             : struct points
    1714             : {
    1715             :   GEN z;
    1716             :   long i, f;
    1717             : };
    1718             : 
    1719             : static int
    1720       88634 : process(long a, long b, GEN y, void *info0, int *quit)
    1721             : {
    1722       88634 :   struct points *p = (struct points *) info0;
    1723       88634 :   if (b==0 && (p->f&2L)) return 0;
    1724       85659 :   *quit = (p->f&1);
    1725       85659 :   p->z = vec_append_grow(p->z, ++p->i, mkvec3(stoi(a),y,stoi(b)));
    1726       85659 :   return 1;
    1727             : }
    1728             : 
    1729             : static int
    1730       12201 : args_h(ratpoints_args *args, GEN D)
    1731             : {
    1732       12201 :   long H, h, l = 1;
    1733             :   GEN L;
    1734       12201 :   switch(typ(D))
    1735             :   {
    1736       12159 :     case t_INT: if (signe(D) <= 0) return 0;
    1737       12159 :       H = h = itos(D); break;
    1738          42 :     case t_VEC: if (lg(D) != 3) return 0;
    1739          42 :       H = gtos(gel(D,1));
    1740          42 :       L = gel(D,2);
    1741          42 :       if (typ(L)==t_INT) { h = itos(L); break; }
    1742          14 :       if (typ(L)==t_VEC && lg(L)==3)
    1743             :       {
    1744           7 :         l = gtos(gel(L,1));
    1745           7 :         h = gtos(gel(L,2)); break;
    1746             :       }
    1747           7 :     default: return 0;
    1748             :   }
    1749       12194 :   args->height = H;
    1750       12194 :   args->b_low  = l;
    1751       12194 :   args->b_high = h; return 1;
    1752             : }
    1753             : 
    1754             : static GEN
    1755       12201 : ZX_hyperellratpoints(GEN P, GEN h, long flag)
    1756             : {
    1757       12201 :   pari_sp av = avma;
    1758             :   ratpoints_args args;
    1759             :   struct points data;
    1760       12201 :   ulong flags = 0;
    1761             : 
    1762       12201 :   if (!args_h(&args, h))
    1763             :   {
    1764           7 :     pari_err_TYPE("hyperellratpoints", h);
    1765             :     return NULL;/*LCOV_EXCL_LINE*/
    1766             :   }
    1767       12194 :   find_points_init(&args, RATPOINTS_DEFAULT_BIT_PRIMES);
    1768             : 
    1769       12194 :   args.cof           = shallowcopy(P);
    1770       12194 :   args.num_inter     = 0;
    1771       12194 :   args.sp1           = RATPOINTS_DEFAULT_SP1;
    1772       12194 :   args.sp2           = RATPOINTS_DEFAULT_SP2;
    1773       12194 :   args.array_size    = RATPOINTS_ARRAY_SIZE;
    1774       12194 :   args.num_primes    = RATPOINTS_DEFAULT_NUM_PRIMES;
    1775       12194 :   args.bit_primes    = RATPOINTS_DEFAULT_BIT_PRIMES;
    1776       12194 :   args.max_forbidden = RATPOINTS_DEFAULT_MAX_FORBIDDEN;
    1777       12194 :   args.flags         = flags;
    1778       12194 :   data.z = cgetg(17,t_VEC);
    1779       12194 :   data.i = 0;
    1780       12194 :   data.f = flag;
    1781       12194 :   find_points_work(&args, process, (void *)&data);
    1782             : 
    1783       12194 :   setlg(data.z, data.i+1);
    1784       12194 :   return gerepilecopy(av, data.z);
    1785             : }
    1786             : 
    1787             : /* The ordinates of the points returned need to be divided by den
    1788             :  * by the caller to get the actual solutions */
    1789             : static GEN
    1790       12201 : QX_hyperellratpoints(GEN P, GEN h, long flag, GEN *den)
    1791             : {
    1792       12201 :   GEN Q = Q_remove_denom(P, den);
    1793       12201 :   if (*den) Q = ZX_Z_mul(Q, *den);
    1794       12201 :   return ZX_hyperellratpoints(Q, h, flag);
    1795             : }
    1796             : 
    1797             : static GEN
    1798         168 : QX_homogenous_evalpow(GEN P, GEN A, GEN B)
    1799             : {
    1800         168 :   pari_sp av = avma;
    1801         168 :   GEN den, Q = Q_remove_denom(P, &den);
    1802         168 :   long i, d = degpol(Q);
    1803         168 :   GEN s = gel(Q, d+2);
    1804         672 :   for (i = d-1; i >= 0; i--)
    1805         504 :     s = addii(mulii(s, A), mulii(gel(B,d+1-i), gel(Q,i+2)));
    1806         168 :   return gerepileupto(av, den ? gdiv(s,den): s);
    1807             : }
    1808             : 
    1809             : static GEN
    1810          70 : to_RgX(GEN a, long v) { return typ(a)==t_POL? a: scalarpol(a,v); }
    1811             : 
    1812             : static int
    1813       11340 : hyperell_check(GEN PQ, GEN *P, GEN *Q)
    1814             : {
    1815       11340 :   *P = *Q = NULL;
    1816       11340 :   if (typ(PQ) == t_POL)
    1817             :   {
    1818       11305 :     if (!RgX_is_QX(PQ)) return 0;
    1819       11305 :     *P = PQ;
    1820             :   }
    1821             :   else
    1822             :   {
    1823          35 :     long v = gvar(PQ);
    1824          35 :     if (v == NO_VARIABLE || typ(PQ) != t_VEC || lg(PQ) != 3) return 0;
    1825          35 :     *P = to_RgX(gel(PQ,1), v); if (!RgX_is_QX(*P)) return 0;
    1826          35 :     *Q = to_RgX(gel(PQ,2), v); if (!RgX_is_QX(*Q)) return 0;
    1827          35 :     if (!signe(*Q)) *Q = NULL;
    1828             :   }
    1829       11340 :   return 1;
    1830             : }
    1831             : 
    1832             : GEN
    1833       11340 : hyperellratpoints(GEN PQ, GEN h, long flag)
    1834             : {
    1835       11340 :   pari_sp av = avma;
    1836             :   GEN P, Q, H, L, den;
    1837             :   long i, l, dy, dQ;
    1838             : 
    1839       11340 :   if (flag<0 || flag>1) pari_err_FLAG("ellratpoints");
    1840       11340 :   if (!hyperell_check(PQ, &P, &Q)) pari_err_TYPE("hyperellratpoints",PQ);
    1841       11340 :   if (!Q)
    1842             :   {
    1843       11319 :     L = QX_hyperellratpoints(P, h, flag|2L, &den);
    1844       11319 :     dy = (degpol(P)+1)>>1;
    1845       11319 :     l = lg(L);
    1846       25004 :     for (i = 1; i < l; i++)
    1847             :     {
    1848       13685 :       GEN Li = gel(L,i), x = gel(Li,1), y = gel(Li,2), z = gel(Li,3);
    1849       13685 :       GEN zdy = powiu(z, dy);
    1850       13685 :       if (den) zdy = mulii(zdy, den);
    1851       13685 :       gel(L,i) = mkvec2(gdiv(x,z), gdiv(y, zdy));
    1852             :     }
    1853       11319 :     return gerepilecopy(av, L);
    1854             :   }
    1855          21 :   H = RgX_add(RgX_muls(P,4), RgX_sqr(Q));
    1856          21 :   dy = (degpol(H)+1)>>1; dQ = degpol(Q);
    1857          21 :   L = QX_hyperellratpoints(H, h, flag|2L, &den);
    1858          21 :   l = lg(L);
    1859         189 :   for (i = 1; i < l; i++)
    1860             :   {
    1861         168 :     GEN Li = gel(L,i), x = gel(Li,1), y = gel(Li,2), z = gel(Li,3);
    1862         168 :     GEN B = gpowers(z, dQ);
    1863         168 :     GEN Qx = gdiv(QX_homogenous_evalpow(Q, x, B), gel(B, dQ+1));
    1864         168 :     GEN zdy = powiu(z, dy);
    1865         168 :     if (den) zdy = mulii(zdy, den);
    1866         168 :     gel(L,i) = mkvec2(gdiv(x,z), gmul2n(gsub(gdiv(y,zdy),Qx),-1));
    1867             :   }
    1868          21 :   return gerepilecopy(av, L);
    1869             : }
    1870             : 
    1871             : GEN
    1872         861 : ellratpoints(GEN E, GEN h, long flag)
    1873             : {
    1874         861 :   pari_sp av = avma;
    1875             :   GEN L, a1, a3, den;
    1876             :   long i, l;
    1877         861 :   checkell_Q(E);
    1878         861 :   if (flag<0 || flag>1) pari_err_FLAG("ellratpoints");
    1879         861 :   if (!RgV_is_QV(vecslice(E,1,5))) pari_err_TYPE("ellratpoints",E);
    1880         861 :   a1 = ell_get_a1(E);
    1881         861 :   a3 = ell_get_a3(E);
    1882         861 :   L = QX_hyperellratpoints(ec_bmodel(E), h, flag|2L, &den);
    1883         854 :   l = lg(L);
    1884       72660 :   for (i = 1; i < l; i++)
    1885             :   {
    1886       71806 :     GEN P, Li = gel(L,i), x = gel(Li,1), y = gel(Li,2), z = gel(Li,3);
    1887       71806 :     if (!signe(z))
    1888           0 :       P = ellinf();
    1889             :     else
    1890             :     {
    1891       71806 :       GEN z2 = sqri(z);
    1892       71806 :       if (den) y = gdiv(y, den);
    1893       71806 :       y = gsub(y, gadd(gmul(a1, mulii(x,z)), gmul(a3,z2)));
    1894       71806 :       P = mkvec2(gdiv(x,z), gdiv(y,shifti(z2,1)));
    1895             :     }
    1896       71806 :     gel(L,i) = P;
    1897             :   }
    1898         854 :   return gerepilecopy(av, L);
    1899             : }

Generated by: LCOV version 1.13