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.16.2 lcov report (development 29401-b333f5c7d5) Lines: 879 908 96.8 %
Date: 2024-06-19 09:03:24 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.2.1 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     2725899 : sieve_init1(long p, ratpoints_sieve_entry *se1, long b1, ratpoints_args *args1)
     213             : {
     214     2725899 :   ratpoints_sieve_entry *se = se1;
     215     2725899 :   ratpoints_args *args = args1;
     216     2725899 :   int *isfs = se->is_f_square;
     217     2725899 :   long b = b1;
     218     2725899 :   long lmp = BITS_IN_LONG % p;
     219     2725899 :   long ldp = BITS_IN_LONG / p;
     220     2725899 :   long p1 = (ldp + 1) * p;
     221     2725899 :   long diff_shift = p1 & LONG_MASK;
     222     2725899 :   long diff = BITS_IN_LONG - diff_shift;
     223             :   ulong help0;
     224             :   long a;
     225     2725899 :   long d = se->inverses[b];
     226     2725899 :   long ab = 0; /* a/b mod p */
     227     2725899 :   ulong test = 1UL;
     228     2725899 :   ulong he0 = 0UL;
     229   121243572 :   for (a = 0; a < p; a++)
     230             :   {
     231   118517673 :     if (isfs[ab]) he0 |= test;
     232   118517673 :     ab += d;
     233   118517673 :     if (ab >= p) ab -= p;
     234   118517673 :     test <<= 1;
     235             :   }
     236     2725899 :   help0 = he0;
     237             :   {
     238             :     ulong help1;
     239             :      /* repeat bit pattern floor(BITS_IN_LONG/p) times */
     240     2725899 :     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     2725899 :     help1 = pattern >> lmp;
     247     6219973 :     for (i = p; i < BITS_IN_LONG; i <<= 1)
     248     3494074 :       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     2725899 :       ulong *si = (ulong *)args->ba_next;
     254             : 
     255     2725899 :       args->ba_next += p + RATPOINTS_CHUNK-1;
     256             :       /* copy the first chunk into sieve[b][] */
     257     2725899 :       si[0] = help0;
     258             :       /* now keep repeating the bit pattern,
     259             :          rotating it in help0/help1 */
     260   118517673 :       for (a = 1 ; a < p; a++)
     261             :       {
     262   115791774 :         ulong temp = help0 >> diff;
     263   115791774 :         help0 = help1 | (help0 << diff_shift);
     264   115791774 :         si[a] = help0;
     265   115791774 :         help1 = temp;
     266             :       }
     267   313601562 :       CODE_INIT_SIEVE_COPY
     268             :       /* set sieve array */
     269     2725899 :       se->sieve[b] = (ratpoints_bit_array *)si;
     270     2725899 :       return (ratpoints_bit_array *)si;
     271             :     }
     272             :   }
     273             : }
     274             : 
     275             : /* This is for p > BITS_IN_LONG */
     276             : static ratpoints_bit_array *
     277     9886497 : sieve_init2(long p, ratpoints_sieve_entry *se1, long b1, ratpoints_args *args1)
     278             : {
     279     9886497 :   pari_sp av = avma;
     280     9886497 :   ratpoints_sieve_entry *se = se1;
     281     9886497 :   ratpoints_args *args = args1;
     282     9886497 :   int *isfs = se->is_f_square;
     283     9886497 :   long b = b1;
     284             :   /* long ldp = 0;  = BITS_IN_LONG / p */
     285             :   /* long p1 = p; = (ldp + 1) * p; */
     286     9886497 :   long wp = p >> TWOPOTBITS_IN_LONG;
     287     9886497 :   long diff_shift = p & LONG_MASK;
     288     9886497 :   long diff = BITS_IN_LONG - diff_shift;
     289     9886497 :   ulong *help = (ulong *) new_chunk((p>>TWOPOTBITS_IN_LONG) + 2);
     290             : 
     291             :   /* initialize help */
     292             :   {
     293     9886497 :     ulong *he = &help[0];
     294     9886497 :     ulong *he1 = &he[(p>>TWOPOTBITS_IN_LONG) + 2];
     295    41776494 :     while (he1 != he) { he1--; *he1 = 0UL; }
     296             :   }
     297     9886497 :   { ulong work = 0UL;
     298             :     long a;
     299     9886497 :     long ab = 0; /* a/b mod p */
     300     9886497 :     long d = se->inverses[b];
     301     9886497 :     long n = 0;
     302     9886497 :     ulong test = 1UL;
     303   960877354 :     for (a = 0; a < p; )
     304             :     {
     305   950990857 :       if (isfs[ab]) work |= test;
     306   950990857 :       ab += d;
     307   950990857 :       if (ab >= p) ab -= p;
     308   950990857 :       test <<= 1;
     309   950990857 :       a++;
     310   950990857 :       if ((a & LONG_MASK) == 0)
     311    12117003 :       { help[n] = work; n++; work = 0UL; test = 1UL; }
     312             :     }
     313     9886497 :     help[n] = work;
     314             :   }
     315             : 
     316             :   { /* fill the bit pattern from help[] into sieve[b][].
     317             :        sieve[b][a0] has the same semantics as help[b][a0],
     318             :        but here, a0 runs from 0 to p-1 and all bits are filled. */
     319     9886497 :     ulong *si = (ulong *)args->ba_next;
     320             :     long a1;
     321             :     long a;
     322             : 
     323     9886497 :     args->ba_next += p + RATPOINTS_CHUNK-1;
     324             :     /* copy the first chunk from help[] into sieve[num][b][] */
     325    22003500 :     for (a = 0; a < wp; a++) si[a] = help[a];
     326             :     /* now keep repeating the bit pattern, rotating it in help */
     327   948760351 :     for (a1 = a ; a < p; a++)
     328             :     {
     329   938873854 :       long t = (a1 == wp) ? 0 : a1+1;
     330   938873854 :       help[a1] |= help[t]<<diff_shift;
     331   938873854 :       si[a] = help[a1];
     332   938873854 :       a1 = t;
     333   938873854 :       help[a1] >>= diff;
     334             :     }
     335  1855843974 :      CODE_INIT_SIEVE_COPY
     336             :     /* set sieve array */
     337     9886497 :     se->sieve[b] = (ratpoints_bit_array *)si;
     338     9886497 :     set_avma(av);
     339     9886497 :     return (ratpoints_bit_array *)si;
     340             :   }
     341             : }
     342             : 
     343             : static GEN
     344       12358 : gen_squares(GEN listprime)
     345             : {
     346       12358 :   long nbprime = lg(listprime)-1;
     347       12358 :   GEN sq = cgetg(nbprime+1, t_VEC);
     348             :   long n;
     349      383098 :   for (n = 1; n <= nbprime; n++)
     350             :   {
     351      370740 :     ulong i, p = uel(listprime,n);
     352      370740 :     GEN w = zero_zv(p), work = w+1;
     353      370740 :     work[0] = 1;
     354             :     /* record nonzero squares mod p, p odd */
     355    10800892 :     for (i = 1; i < p; i += 2) work[(i*i) % p] = 1;
     356      370740 :     gel(sq, n) = w;
     357             :   }
     358       12358 :   return sq;
     359             : }
     360             : 
     361             : static GEN
     362       12358 : gen_offsets(GEN P)
     363             : {
     364       12358 :   long n, l = lg(P);
     365       12358 :   GEN of = cgetg(l, t_VEC);
     366      383098 :   for (n = 1; n < l; n++)
     367             :   {
     368      370740 :     ulong p = uel(P,n);
     369      370740 :     uel(of, n) = Fl_inv((2*RBA_LENGTH)%p, p);
     370             :   }
     371       12358 :   return of;
     372             : }
     373             : 
     374             : static GEN
     375       12358 : gen_inverses(GEN P)
     376             : {
     377       12358 :   long n, l = lg(P);
     378       12358 :   GEN iv = cgetg(l, t_VEC);
     379      383098 :   for (n = 1; n < l; n++)
     380             :   {
     381      370740 :     ulong i, p = uel(P,n);
     382      370740 :     GEN w = cgetg(p, t_VECSMALL);
     383    21231044 :     for (i = 1; i < p; i++) uel(w, i) = Fl_inv(i, p);
     384      370740 :     gel(iv, n) = w;
     385             :   }
     386       12358 :   return iv;
     387             : }
     388             : 
     389             : static ulong **
     390       12358 : gen_sieves0(GEN listprime)
     391             : {
     392             :   long n;
     393       12358 :   long nbprime = lg(listprime)-1;
     394       12358 :   ulong ** w = (ulong**) new_chunk(nbprime+1);
     395      383098 :   for (n = 1; n <= nbprime; n++)
     396             :   {
     397      370740 :     ulong a, p = uel(listprime,n);
     398      370740 :     ulong *si = (ulong *) stack_malloc_align((p+RATPOINTS_CHUNK-1)*RBA_SIZE, RBA_SIZE);
     399    21601784 :     for (a = 0; a < p; a++) si[a] = ~0UL;
     400    22406580 :     for (a = 0; a < BITS_IN_LONG; a++)
     401    22035840 :       uel(si,(p*a)>>TWOPOTBITS_IN_LONG) &= ~(1UL<<((p*a) & LONG_MASK));
     402    46262136 :     CODE_INIT_SIEVE_COPY
     403      370740 :     w[n] = si;
     404             :   }
     405       12358 :   return w;
     406             : }
     407             : 
     408             : static void
     409       12358 : gen_sieve(ratpoints_args *args)
     410             : {
     411       12358 :   GEN listprimes = args->listprime;
     412       12358 :   args->offsets = gen_offsets(listprimes);
     413       12358 :   args->inverses = gen_inverses(listprimes);
     414       12358 :   args->sieves0 = gen_sieves0(listprimes);
     415       12358 : }
     416             : 
     417             : static GEN
     418       12358 : ZX_positive_region(GEN P, long h, long bitprec)
     419             : {
     420       12358 :   long prec = nbits2prec(bitprec);
     421       12358 :   GEN it = mkvec2(stoi(-h),stoi(h));
     422       12358 :   GEN R = realroots(P, it, prec);
     423       12358 :   long nR = lg(R)-1;
     424       12358 :   long s = signe(ZX_Z_eval(P,gel(it,1)));
     425       12358 :   long i=1, j;
     426             :   GEN iv, st, en;
     427       12358 :   if (s<0 && nR==0) return NULL;
     428       11686 :   iv = cgetg(((nR+1+(s>=0))>>1)+1, t_VEC);
     429       11686 :   if (s>=0) st = itor(gel(it,1),prec);
     430        5089 :   else    { st = gel(R,i); i++; }
     431       18080 :   for (j=1; i<nR; j++)
     432             :   {
     433        6394 :     gel(iv, j) = mkvec2(st, gel(R,i));
     434        6394 :     st = gel(R,i+1);
     435        6394 :     i+=2;
     436             :   }
     437       11686 :   if (i==nR) en = gel(R,i); else en = itor(gel(it,2),prec);
     438       11686 :   gel(iv,j) = mkvec2(st, en);
     439       11686 :   return iv;
     440             : }
     441             : 
     442             : static long
     443       12358 : posint(ratpoints_interval *ivlist, GEN P, long h)
     444             : {
     445       12358 :   GEN R = ZX_positive_region(P, h, 53);
     446       12358 :   const double eps = 1e-5;
     447             :   long nR, i;
     448             : 
     449       12358 :   if (!R) return 0;
     450       11686 :   nR = lg(R)-1;
     451       11686 :   i = 1;
     452       29766 :   for (i=1; i<=nR; i++)
     453             :   {
     454       18080 :     ivlist[i-1].low = rtodbl(gmael(R,i,1))-eps;
     455       18080 :     ivlist[i-1].up  = rtodbl(gmael(R,i,2))+eps;
     456             :   }
     457       11686 :   return nR;
     458             : }
     459             : 
     460             : static long
     461       12358 : ratpoints_compute_sturm(ratpoints_args *args)
     462             : {
     463       12358 :   ratpoints_interval *ivlist = args->domain;
     464       12358 :   args->num_inter = posint(ivlist, args->cof, (long) ivlist[0].up);
     465       12358 :   return args->num_inter;
     466             : }
     467             : 
     468             : /**************************************************************************
     469             :  * Try to avoid divisions                                                 *
     470             :  **************************************************************************/
     471             : INLINE long
     472   886184974 : mod(long a, long b)
     473             : {
     474   886184974 :   long b1 = b << 4; /* b1 = 16*b */
     475             : 
     476   886184974 :   if (a < -b1) { a %= b; if (a < 0) { a += b; } return a ; }
     477   875491073 :   if (a < 0) { a += b1; }
     478   354346373 :   else { if (a >= b1) { return a % b; } }
     479   868087807 :   b1 >>= 1; /* b1 = 8*b */
     480   868087807 :   if (a >= b1) { a -= b1; }
     481   868087807 :   b1 >>= 1; /* b1 = 4*b */
     482   868087807 :   if (a >= b1) { a -= b1; }
     483   868087807 :   b1 >>= 1; /* b1 = 2*b */
     484   868087807 :   if (a >= b1) { a -= b1; }
     485   868087807 :   if (a >= b) { a -= b; }
     486   868087807 :   return a;
     487             : }
     488             : 
     489             : static void
     490     2325903 : set_bc(long b, ratpoints_args *args)
     491             : {
     492     2325903 :   GEN w0 = gen_1;
     493     2325903 :   GEN c = args->cof, bc;
     494     2325903 :   long k, degree = degpol(c);
     495     2325903 :   bc = cgetg(degree+2, t_POL);
     496    11877658 :   for (k = degree-1; k >= 0; k--)
     497             :   {
     498     9551755 :     w0 = muliu(w0, b);
     499     9551755 :     gel(bc,k+2) = mulii(gel(c,k+2), w0);
     500             :   }
     501     2325903 :   args->bc = bc;
     502     2325903 : }
     503             : 
     504             : /**************************************************************************
     505             :  * Check a `survivor' of the sieve if it really gives a point.            *
     506             :  **************************************************************************/
     507             : 
     508             : static long
     509     3281834 : _ratpoints_check_point(long a, long b, ratpoints_args *args, int *quit,
     510             :                  int process(long, long, GEN, void*, int*), void *info)
     511             : {
     512     3281834 :   pari_sp av = avma;
     513     3281834 :   GEN w0, w2, c = args->cof, bc = args->bc;
     514     3281834 :   long k, degree = degpol(c);
     515     3281834 :   int reverse = args->flags & RATPOINTS_REVERSED;
     516             : 
     517             :   /* Compute F(a, b), where F is the homogenized version of f
     518             :      of smallest possible even degree  */
     519     3281834 :   w2 = gel(c, degree+2);
     520    16930978 :   for (k = degree-1; k >= 0; k--)
     521             :   {
     522    13649144 :     w2 = mulis(w2, a);
     523    13649144 :     w2 = addii(w2, gel(bc,k+2));
     524             :   }
     525     3281834 :   if (odd(degree)) w2 = muliu(w2, b);
     526             :   /* check if f(x,z) is a square; if so, process the point(s) */
     527     3281834 :   if (signe(w2) >= 0 && Z_issquareall(w2, &w0))
     528             :   {
     529       44891 :     if (reverse)
     530             :     {
     531        1218 :       if (a >= 0) (void)process(b, a, w0, info, quit);
     532         217 :       else        (void)process(-b, -a, w0, info, quit);
     533             :     }
     534       43673 :     else (void)process(a, b, w0, info, quit);
     535       44891 :     if (!*quit && signe(w0) != 0)
     536             :     {
     537       42602 :       GEN nw0 = negi(w0);
     538       42602 :       if (reverse)
     539             :       {
     540        1155 :         if (a >= 0) (void)process(b, a, nw0, info, quit);
     541         196 :         else        (void)process(-b, -a, nw0, info, quit);
     542             :       }
     543       41447 :       else (void)process(a, b, nw0, info, quit);
     544             :     }
     545       44891 :     return 1;
     546             :   }
     547     3236943 :   set_avma(av);
     548     3236943 :   return 0;
     549             : }
     550             : 
     551             : /**************************************************************************
     552             :  * The inner loop of the sieving procedure                                *
     553             :  **************************************************************************/
     554             : static long
     555    46430623 : _ratpoints_sift0(long b, long w_low, long w_high,
     556             :            ratpoints_args *args, bit_selection which_bits,
     557             :            ratpoints_bit_array *survivors, sieve_spec *sieves, int *quit,
     558             :            int process(long, long, GEN, void*, int*), void *info)
     559             : {
     560    46430623 :   long sp1 = args->sp1, sp2 = args->sp2;
     561    46430623 :   long i, n, nb = 0, absb = labs(b), base = 0;
     562             :   ratpoints_bit_array *surv0;
     563             : 
     564             :   /* now do the sieving (fast!) */
     565             : #if (RATPOINTS_CHUNK == 16)
     566             :   long w_low_new;
     567    35552238 :   ratpoints_bit_array *surv = survivors;
     568             : 
     569             :   /* first set the start fields for the first and second phases of sieving */
     570   710409684 :   for(n = 0; n < sp2; n++)
     571   674857446 :     sieves[n].start = sieves[n].ptr + mod(w_low + sieves[n].offset, sieves[n].p);
     572             :   /* Take RATPOINTS_CHUNK bit-arrays and apply phase 1 to them,
     573             :    * then repeat with the next RATPOINTS_CHUNK bit-arrays. */
     574   243270858 :   for(w_low_new = w_low; w_low_new < w_high; surv += RATPOINTS_CHUNK, w_low_new += RATPOINTS_CHUNK)
     575             :   {
     576             :     /* read data from memory into registers */
     577   207718620 :     ratpoints_bit_array reg0 = surv[0];
     578   207718620 :     ratpoints_bit_array reg1 = surv[1];
     579   207718620 :     ratpoints_bit_array reg2 = surv[2];
     580   207718620 :     ratpoints_bit_array reg3 = surv[3];
     581   207718620 :     ratpoints_bit_array reg4 = surv[4];
     582   207718620 :     ratpoints_bit_array reg5 = surv[5];
     583   207718620 :     ratpoints_bit_array reg6 = surv[6];
     584   207718620 :     ratpoints_bit_array reg7 = surv[7];
     585   207718620 :     ratpoints_bit_array reg8 = surv[8];
     586   207718620 :     ratpoints_bit_array reg9 = surv[9];
     587   207718620 :     ratpoints_bit_array reg10 = surv[10];
     588   207718620 :     ratpoints_bit_array reg11 = surv[11];
     589   207718620 :     ratpoints_bit_array reg12 = surv[12];
     590   207718620 :     ratpoints_bit_array reg13 = surv[13];
     591   207718620 :     ratpoints_bit_array reg14 = surv[14];
     592   207718620 :     ratpoints_bit_array reg15 = surv[15];
     593             : 
     594  2492623044 :     for(n = 0; n < sp1; n++)
     595             :     { /* retrieve the pointer to the beginning of the relevant bits */
     596  2284904424 :       ratpoints_bit_array *siv1 = sieves[n].start;
     597  2284904424 :       reg0 &= *siv1++;
     598  2284904424 :       reg1 &= *siv1++;
     599  2284904424 :       reg2 &= *siv1++;
     600  2284904424 :       reg3 &= *siv1++;
     601  2284904424 :       reg4 &= *siv1++;
     602  2284904424 :       reg5 &= *siv1++;
     603  2284904424 :       reg6 &= *siv1++;
     604  2284904424 :       reg7 &= *siv1++;
     605  2284904424 :       reg8 &= *siv1++;
     606  2284904424 :       reg9 &= *siv1++;
     607  2284904424 :       reg10 &= *siv1++;
     608  2284904424 :       reg11 &= *siv1++;
     609  2284904424 :       reg12 &= *siv1++;
     610  2284904424 :       reg13 &= *siv1++;
     611  2284904424 :       reg14 &= *siv1++;
     612  2284904424 :       reg15 &= *siv1++;
     613             : 
     614             :       /* update the pointer for the next round
     615             :        * (RATPOINTS_CHUNK-1 bit-arrays after sieves[n].end) */
     616  3218102898 :       while(siv1 >= sieves[n].end) siv1 -= sieves[n].p;
     617  2284904424 :       sieves[n].start = siv1;
     618             :     }
     619             :     /* store the contents of the registers back into memory */
     620   207718620 :     surv[0] = reg0;
     621   207718620 :     surv[1] = reg1;
     622   207718620 :     surv[2] = reg2;
     623   207718620 :     surv[3] = reg3;
     624   207718620 :     surv[4] = reg4;
     625   207718620 :     surv[5] = reg5;
     626   207718620 :     surv[6] = reg6;
     627   207718620 :     surv[7] = reg7;
     628   207718620 :     surv[8] = reg8;
     629   207718620 :     surv[9] = reg9;
     630   207718620 :     surv[10] = reg10;
     631   207718620 :     surv[11] = reg11;
     632   207718620 :     surv[12] = reg12;
     633   207718620 :     surv[13] = reg13;
     634   207718620 :     surv[14] = reg14;
     635   207718620 :     surv[15] = reg15;
     636             :   }
     637             : #else /* RATPOINTS_CHUNK not between 2 and 16 */
     638    10878385 :   long range = w_high - w_low;
     639   130540554 :   for (n = 0; n < sp1; n++)
     640             :   {
     641   119662169 :     ratpoints_bit_array *sieve_n = sieves[n].ptr;
     642   119662169 :     long p = sieves[n].p;
     643   119662169 :     long r = mod(-w_low-sieves[n].offset, p);
     644   119662169 :     ratpoints_bit_array *surv = survivors;
     645             : 
     646   119662169 :     if (w_high < w_low + r)
     647             :     { /* if we get here, r > 0, since w_high >= w_low always */
     648     6919017 :       ratpoints_bit_array *siv1 = &sieve_n[p-r];
     649     6919017 :       ratpoints_bit_array *siv0 = siv1 + range;
     650             : 
     651   206350409 :       while (siv1 != siv0) { *surv &= *siv1++; surv++; }
     652             :     }
     653             :     else
     654             :     {
     655   112743152 :       ratpoints_bit_array *siv1 = &sieve_n[p-r];
     656   112743152 :       ratpoints_bit_array *surv_end = &survivors[range - p];
     657             :       long i;
     658  3564181740 :       for (i = r; i; i--) { *surv &= *siv1++; surv++; }
     659   112743152 :       siv1 -= p;
     660   572026248 :       while (surv <= surv_end)
     661             :       {
     662 15769155204 :         for (i = p; i; i--) { *surv &= *siv1++; surv++; }
     663   459283096 :         siv1 -= p;
     664             :       }
     665   112743152 :       surv_end += p;
     666  3641042738 :       while (surv < surv_end) { *surv &= *siv1++; surv++; }
     667             :     }
     668             :   }
     669             :   /* initialize pointers in sieve for the second phase */
     670    97693926 :   for(n = sp1; n < sp2; n++)
     671    86815541 :     sieves[n].start = sieves[n].ptr + mod(w_low + sieves[n].offset, sieves[n].p);
     672             : #endif /* RATPOINTS_CHUNK */
     673             : 
     674             :   /* 2nd phase of the sieve: test each surviving bit array with more primes */
     675    46430623 :   surv0 = &survivors[0];
     676  5414288225 :   for (i = w_low; i < w_high; i++, base++)
     677             :   {
     678  5367859359 :     ratpoints_bit_array nums = *surv0++;
     679  5367859359 :     sieve_spec *ssp = &sieves[sp1];
     680             :     long n;
     681             : 
     682  5866235905 :     for (n = sp2-sp1; n && TEST(nums); n--)
     683             :     {
     684   498376546 :       ratpoints_bit_array *ptr = (ssp->start) + base;
     685   498376546 :       long p = ssp->p;
     686  1324996669 :       while(ptr >= ssp->end) ptr -= p;
     687   498376546 :       nums &= *ptr;
     688   498376546 :       ssp++;
     689             :     }
     690             : 
     691             :     /* Check the survivors of the sieve if they really give points */
     692  5367859359 :     if (TEST(nums))
     693             :     {
     694             :       long a0, a, d;
     695    10789757 :       ulong nums0 = EXT0(nums);
     696             :       /* a will be the numerator corresponding to the selected bit */
     697    10789757 :       if (which_bits == num_all)
     698             :       {
     699     6949092 :         d = 1; a0 = i * RBA_LENGTH;
     700             :       }
     701             :       else
     702             :       {
     703     3840665 :         d = 2; a0 = i * 2 * RBA_LENGTH;
     704     3840665 :         if (which_bits == num_odd) a0++;
     705             :       }
     706   134499944 :       for (a = a0; nums0; a += d, nums0 >>= 1)
     707             :       { /* test one bit */
     708   123711248 :         if (odd(nums0) && ugcd(labs(a), absb)==1)
     709             :         {
     710     1876370 :           if (!args->bc) set_bc(b, args);
     711     1876370 :           nb += _ratpoints_check_point(a, b, args, quit, process, info);
     712     1876370 :           if (*quit) return nb;
     713             :         }
     714             :       }
     715             : #ifdef RBA_USE_VX
     716             :       {
     717     9246696 :         ulong k, da = d<<TWOPOTBITS_IN_LONG;
     718    18492696 :         for (k = 1; k < RBA_PACK; k++)
     719             :         {
     720     9246696 :           ulong nums1 = EXT(nums,k);
     721     9246696 :           a0 += da;
     722   112226682 :           for (a = a0; nums1; a += d, nums1 >>= 1)
     723             :           { /* test one bit */
     724   102980682 :             if (odd(nums1) && ugcd(labs(a), absb)==1)
     725             :             {
     726     1405464 :               if (!args->bc) set_bc(b, args);
     727     1405464 :               nb += _ratpoints_check_point(a, b, args, quit, process, info);
     728     1405464 :               if (*quit) return nb;
     729             :             }
     730             :           }
     731             :         }
     732             :       }
     733             : #endif
     734             :     }
     735             :   }
     736    46428866 :   return nb;
     737             : }
     738             : 
     739             : typedef struct { double r; ratpoints_sieve_entry *ssp; } entry;
     740             : 
     741             : static const int squares16[16] = {1,1,0,0,1,0,0,0,0,1,0,0,0,0,0,0};
     742             :  /* Says if a is a square mod 16, for a = 0..15 */
     743             : 
     744             : /**************************************************************************
     745             :  * Initialization and cleanup of ratpoints_args structure                 *
     746             :  **************************************************************************/
     747             : 
     748             : static ratpoints_sieve_entry*
     749       12358 : alloc_sieve(long nbprime, long maxp)
     750             : {
     751             :   long i;
     752             :   ratpoints_sieve_entry * s = (ratpoints_sieve_entry*)
     753       12358 :                         stack_malloc(nbprime*sizeof(ratpoints_sieve_entry));
     754      383098 :   for (i=0; i<nbprime; i++)
     755      370740 :     s[i].sieve = (ratpoints_bit_array**) new_chunk(maxp);
     756       12358 :   return s;
     757             : }
     758             : 
     759             : /* NOTE: args->degree must be set */
     760             : static void
     761       12358 : find_points_init(ratpoints_args *args, long bit_primes)
     762             : {
     763       12358 :   long need = 0, n, nbprime, maxp;
     764       12358 :   args->listprime = primes_interval_zv(3, 1<<bit_primes);
     765       12358 :   nbprime = lg(args->listprime)-1;
     766       12358 :   maxp = args->listprime[nbprime];
     767             : 
     768             :   /* allocate space for se_buffer */
     769       12358 :   args->se_buffer = alloc_sieve(nbprime, maxp);
     770       12358 :   args->se_next = args->se_buffer;
     771      383098 :   for (n = 1; n <= nbprime; n++)
     772             :   {
     773      370740 :     ulong p = args->listprime[n];
     774      370740 :     need += p*(p + RATPOINTS_CHUNK-1);
     775             :   }
     776       12358 :   args->ba_buffer = (ratpoints_bit_array*)
     777       12358 :      stack_malloc_align(need*RBA_SIZE,RBA_SIZE);
     778       12358 :   args->ba_next = args->ba_buffer;
     779             : 
     780             :   /* allocate space for int_buffer */
     781       12358 :   args->int_buffer = (int *) stack_malloc(nbprime*(maxp+1)*sizeof(int));
     782       12358 :   args->int_next = args->int_buffer;
     783             : 
     784       12358 :   args->forb_ba   = (forbidden_entry*)
     785       12358 :     stack_malloc((nbprime + 1)*sizeof(forbidden_entry));
     786       12358 :   args->forbidden = new_chunk(nbprime + 1);
     787       12358 :   gen_sieve(args);
     788       12358 :   return;
     789             : }
     790             : 
     791             : /* f = leading coeff; b = b1*b2, b1 maximal with (b1, 2*f) = 1;
     792             :  * return Jacobi symbol (f, b1) */
     793             : INLINE int
     794    51190767 : rpjacobi(long b, GEN lcf)
     795             : {
     796             :   ulong f;
     797    51190767 :   b >>= vals(b);
     798    51190767 :   f = umodiu(lcf, b);
     799    51190767 :   return krouu(f, u_ppo(b,f));
     800             : }
     801             : 
     802             : /************************************************************************
     803             :  * Set up information on possible denominators                          *
     804             :  * when polynomial is of odd degree with leading coefficient != +-1     *
     805             :  ************************************************************************/
     806             : 
     807             : static void
     808        1337 : setup_us1(ratpoints_args *args, GEN w0)
     809             : {
     810        1337 :   GEN F = Z_issmooth_fact(w0, 1000), P, E, S, D;
     811             :   long i, l;
     812             : 
     813        1337 :   if (!F) return;
     814        1337 :   P = gel(F,1); l = lg(P);
     815        1337 :   E = gel(F,2);
     816        1337 :   D  = cgetg(1+(1<<(l-1)), t_VECSMALL);
     817             :   /* factorization is complete, set up array of squarefree divisors */
     818        1337 :   D[1] = 1;
     819        2800 :   for (i = 1; i < l; i++)
     820             :   { /* multiply all divisors known so far by next prime */
     821        1463 :     long k, n = 1<<(i-1);
     822        3052 :     for (k=0; k<n; k++) uel(D,1+n+k) = uel(D,1+k) * P[i];
     823             :   }
     824        1337 :   S = cgetg(l, t_VECSMALL);
     825             :   /* set slopes in den_info */
     826        2800 :   for (i = 1; i < l; i++)
     827             :   { /* compute min{n : (d-k)*n > v_p(f_d) - v_p(f_k), k = 0,...,d-1} */
     828        1463 :     GEN c = args->cof;
     829        1463 :     long p = P[i], v = E[i];
     830        1463 :     long k, n = 1, d = degpol(c);
     831             : 
     832        6986 :     for (k = d - 1; k >= 0; k--)
     833             :     {
     834        5523 :       long t = 1 + v - Z_lval(gel(c,k+2), p);
     835        5523 :       long m = CEIL(t, d - k);
     836             : 
     837        5523 :       if (m > n) n = m;
     838             :     }
     839        1463 :     S[i] = n;
     840             :   }
     841        1337 :   args->divisors = D;
     842        1337 :   args->flags |= RATPOINTS_USE_SQUARES1;
     843        1337 :   args->den_info = mkvec3(P, E, S);
     844             : }
     845             : 
     846             : /************************************************************************
     847             :  * Consider 2-adic information                                          *
     848             :  ************************************************************************/
     849             : 
     850             : static bit_selection
     851       12358 : get_2adic_info(ratpoints_args *args, ulong *den_bits,
     852             :                ratpoints_bit_array *num_bits)
     853             : {
     854       12358 :   GEN c = args->cof;
     855       12358 :   long degree = degpol(c);
     856             :   int is_f_square16[24];
     857       12358 :   long *cmp = new_chunk(degree+1);
     858       12358 :   long npe = 0, npo = 0;
     859             :   bit_selection result;
     860             :   long n, a, b;
     861             : 
     862             :   /* compute coefficients mod 16 */
     863       86769 :   for (n = 0; n <= degree; n++) cmp[n] = Mod16(gel(c,n+2));
     864      210086 :   for (a = 0 ; a < 16; a++)
     865             :   {
     866      197728 :     ulong s = cmp[degree];
     867             :     long n;
     868     1190576 :     for (n = degree - 1 ; n >= 0 ; n--) s = s*a + cmp[n];
     869      197728 :     s &= 0xf;
     870      197728 :     if ((is_f_square16[a] = squares16[s])) { if (odd(a)) npo++; else npe++; }
     871             :   }
     872             : 
     873             :   /* even denominators:
     874             :      is_f_square16[16+k] says if f((2k+1)/2) is a square, k = 0..3
     875             :      is_f_square16[20+k] says if f((2k+1)/4) is a square, k = 0,1
     876             :      is_f_square16[22]   says if f(odd/8) is a square
     877             :      is_f_square16[23]   says if f(odd/2^n), n >= 4, can be a square */
     878             :   {
     879       12358 :     long np1 = 0, np2 = 0, np3 = 0, np4 = 0;
     880             : 
     881       12358 :     if (odd(degree))
     882             :     {
     883        1407 :       long a, cf = 4*cmp[degree-1];
     884             : 
     885        1407 :       if (degree >= 2) cf += 8*cmp[degree-2];
     886        7035 :       for (a = 0; a < 4; a++)
     887             :       { /* Compute  2 c[d] k^d + 4 c[d-1] k^(d-1) + 8 c[d-2] k^(d-2), k = 2a+1.
     888             :            Note that k^d = k mod 8, k^(d-1) = 1 mod 8. */
     889        5628 :         long k = 2*a+1;
     890        5628 :         long s = (2*k*cmp[degree] + cf) & 0xf;
     891        5628 :         if ((is_f_square16[16+a] = squares16[s])) np1++;
     892             :       }
     893        1407 :       if ((is_f_square16[20] = squares16[(4*cmp[degree])  & 0xf])) np2++;
     894        1407 :       if ((is_f_square16[21] = squares16[(12*cmp[degree]) & 0xf])) np2++;
     895        1407 :       if ((is_f_square16[22] = squares16[(8*cmp[degree])  & 0xf])) np3++;
     896        1407 :       is_f_square16[23] = 1; np4++;
     897             :     }
     898             :     else
     899             :     {
     900       10951 :       long a, cf = (degree >= 2) ? 4*cmp[degree-2] : 0;
     901             : 
     902       10951 :       if (degree >= 3) cf += 8*cmp[degree-3];
     903       54755 :       for (a = 0; a < 4; a++)
     904             :       { /* compute c[d] k^d + 2 c[d-1] k^(d-1) + ... + 8 c[d-3] k^(d-3),
     905             :            k = 2a+1. Note that k^d = k^2 mod 16, k^(d-1) = k mod 8. */
     906       43804 :         long k = 2*a+1;
     907       43804 :         long s = ((cmp[degree]*k + 2*cmp[degree-1])*k + cf) & 0xf;
     908       43804 :         if ((is_f_square16[16+a] = squares16[s])) np1++;
     909             :       }
     910       10951 :       if ((is_f_square16[20] = squares16[(cmp[degree]+4*cmp[degree-1])  & 0xf]))
     911        4694 :         np2++;
     912       10951 :       if ((is_f_square16[21] = squares16[(cmp[degree]+12*cmp[degree-1]) & 0xf]))
     913        4666 :         np2++;
     914       10951 :       if ((is_f_square16[22] = squares16[(cmp[degree]+8*cmp[degree-1])  & 0xf]))
     915        4540 :         np3++;
     916       10951 :       if ((is_f_square16[23] = squares16[cmp[degree]])) np4++;
     917             :     }
     918             : 
     919             :     /* set den_bits */
     920       12358 :     { ulong db = 0;
     921             :       long i;
     922             : 
     923       12358 :       if (npe + npo > 0) db |= 0xaaaaUL; /* odd denominators */
     924       12358 :       if (np1 > 0)       db |= 0x4444UL; /* v_2(den) = 1 */
     925       12358 :       if (np2 > 0)       db |= 0x1010UL; /* v_2(den) = 2 */
     926       12358 :       if (np3 > 0)       db |= 0x0100UL; /* v_2(den) = 3 */
     927       12358 :       if (np4 > 0)       db |= 0x0001UL; /* v_2(den) >= 4 */
     928       12358 :       if (db == 0)
     929             :       {
     930        2975 :         for (i = 0 ; i < 16; i++) num_bits[i] = RBA(0UL);
     931         175 :         *den_bits = 0UL; return num_none;
     932             :       }
     933       34812 :       for (i = 16; i < BITS_IN_LONG; i <<= 1) db |= db << i;
     934       12183 :       *den_bits = db;
     935             :     }
     936       12183 :     result = (npe == 0) ? (npo == 0) ? num_none : num_odd
     937       12183 :                         : (npo == 0) ? num_even : num_all;
     938             :   }
     939             : 
     940             :   /* set up num_bits[16] */
     941             : 
     942             :   /* odd denominators */
     943       12183 :   switch(result)
     944             :   {
     945        7918 :     case num_all:
     946       71262 :       for (b = 1; b < 16; b += 2)
     947             :       {
     948       63344 :         ulong work = 0, bit = 1;
     949       63344 :         long i, invb = b; /* inverse of b mod 16 */
     950       63344 :         if (b & 2) invb ^= 8;
     951       63344 :         if (b & 4) invb ^= 8;
     952     1076848 :         for (i = 0; i < 16; i++)
     953             :         {
     954     1013504 :           if (is_f_square16[(invb*i) & 0xf]) work |= bit;
     955     1013504 :           bit <<= 1;
     956             :         }
     957             :         /* now repeat the 16 bits */
     958      181024 :         for (i = 16; i < BITS_IN_LONG; i <<= 1) work |= work << i;
     959       63344 :         num_bits[b] = RBA(work);
     960             :       }
     961        7918 :       break;
     962             : 
     963        1814 :     case num_odd:
     964       16326 :       for (b = 1; b < 16; b += 2)
     965             :       {
     966       14512 :         ulong work = 0, bit = 1;
     967       14512 :         long i, invb = b; /* inverse of b mod 16 */
     968       14512 :         if (b & 2) invb ^= 8;
     969       14512 :         if (b & 4) invb ^= 8;
     970      130608 :         for (i = 1; i < 16; i += 2)
     971             :         {
     972      116096 :           if (is_f_square16[(invb*i) & 0xf]) work |= bit;
     973      116096 :           bit <<= 1;
     974             :         }
     975             :         /* now repeat the 8 bits */
     976       55968 :         for (i = 8; i < BITS_IN_LONG; i <<= 1) { work |= work << i; }
     977       14512 :         num_bits[b] = RBA(work);
     978             :       }
     979        1814 :       break;
     980             : 
     981        2059 :     case num_even:
     982       18531 :       for (b = 1; b < 16; b += 2)
     983             :       {
     984       16472 :         ulong work = 0, bit = 1;
     985       16472 :         long i, invb = b; /* inverse of b mod 16 */
     986       16472 :         if (b & 2) invb ^= 8;
     987       16472 :         if (b & 4) invb ^= 8;
     988      148248 :         for (i = 0; i < 16; i += 2)
     989             :         {
     990      131776 :           if (is_f_square16[(invb*i) & 0xf]) work |= bit;
     991      131776 :           bit <<= 1;
     992             :         }
     993             :         /* now repeat the 8 bits */
     994       63528 :         for (i = 8; i < BITS_IN_LONG; i <<= 1) work |= work << i;
     995       16472 :         num_bits[b] = RBA(work);
     996             :       }
     997        2059 :       break;
     998             : 
     999         392 :     case num_none:
    1000        3528 :       for (b = 1; b < 16; b += 2) num_bits[b] = RBA(0UL);
    1001         392 :       break;
    1002             :   }
    1003             : 
    1004             :   /* v_2(den) = 1 : only odd numerators */
    1005       60915 :   for (b = 1; b < 8; b += 2)
    1006             :   {
    1007       48732 :     ulong work = 0, bit = 1;
    1008             :     long i;
    1009      438588 :     for (i = 1; i < 16; i += 2)
    1010             :     {
    1011      389856 :       if (is_f_square16[16 + (((b*i)>>1) & 0x3)]) work |= bit;
    1012      389856 :       bit <<= 1;
    1013             :     }
    1014             :     /* now repeat the 8 bits */
    1015      187980 :     for (i = 8; i < BITS_IN_LONG; i <<= 1) work |= work << i;
    1016       48732 :     num_bits[2*b] = RBA(work);
    1017             :   }
    1018             : 
    1019             :   /* v_2(den) = 2 : only odd numerators */
    1020       36549 :   for (b = 1; b < 4; b += 2)
    1021             :   {
    1022       24366 :     ulong work = 0, bit = 1;
    1023             :     long i;
    1024      121830 :     for (i = 1; i < 8; i += 2)
    1025             :     {
    1026       97464 :       if (is_f_square16[20 + (((b*i)>>1) & 0x1)]) work |= bit;
    1027       97464 :       bit <<= 1;
    1028             :     }
    1029             :     /* now repeat the 4 bits */
    1030      118356 :     for (i = 4; i < BITS_IN_LONG; i <<= 1) work |= work << i;
    1031       24366 :     num_bits[4*b] = RBA(work);
    1032             :   }
    1033             : 
    1034             :   /* v_2(den) = 3, >= 4 : only odd numerators */
    1035       12183 :   num_bits[8] = (is_f_square16[22]) ? RBA(~(0UL)) : RBA(0UL);
    1036       12183 :   num_bits[0] = (is_f_square16[23]) ? RBA(~(0UL)) : RBA(0UL);
    1037             : 
    1038       12183 :   return result;
    1039             : }
    1040             : 
    1041             : /**************************************************************************
    1042             :  * This is a comparison function needed for sorting in order to determine *
    1043             :  * the `best' primes for sieving.                                         *
    1044             :  **************************************************************************/
    1045             : 
    1046             : static int
    1047     1165117 : compare_entries(const void *a, const void *b)
    1048             : {
    1049     1165117 :   double diff = ((entry *)a)->r - ((entry *)b)->r;
    1050     1165117 :   return (diff > 0) ? 1 : (diff < 0) ? -1 : 0;
    1051             : }
    1052             : 
    1053             : /************************************************************************
    1054             :  * Collect the sieving information                                      *
    1055             :  ************************************************************************/
    1056             : 
    1057             : static long
    1058       12358 : sieving_info(ratpoints_args *args,
    1059             :              ratpoints_sieve_entry **sieve_list)
    1060             : {
    1061       12358 :   GEN c = args->cof;
    1062       12358 :   GEN inverses = args->inverses, squares;
    1063       12358 :   GEN offsets = args->offsets;
    1064       12358 :   ulong ** sieves0 = args->sieves0;
    1065       12358 :   long degree = degpol(c);
    1066       12358 :   long fba = 0, fdc = 0;
    1067       12358 :   long pn, pnp = 0;
    1068             :   long n;
    1069       12358 :   long nbprime = lg(args->listprime)-1;
    1070             :   long * coeffs_mod_p;
    1071       12358 :   entry *prec = (entry*) stack_malloc(nbprime*sizeof(entry));
    1072             :     /* This array is used for sorting in order to
    1073             :        determine the `best' sieving primes. */
    1074             : 
    1075       12358 :   forbidden_entry *forb_ba = args->forb_ba;
    1076       12358 :   long *forbidden = args->forbidden;
    1077       12358 :   ulong bound = (1UL)<<(BITS_IN_LONG - args->bit_primes);
    1078       12358 :   pari_sp av = avma;
    1079       12358 :   squares = gen_squares(args->listprime);
    1080       12358 :   coeffs_mod_p = (long *) new_chunk(degree+1);
    1081             :   /* initialize sieve in se_buffer */
    1082      379367 :   for (pn = 1; pn <= args->num_primes; pn++)
    1083             :   {
    1084      367135 :     ulong a, p = args->listprime[pn], np;
    1085             :     long n;
    1086      367135 :     int *is_f_square = args->int_next;
    1087             : 
    1088      367135 :     args->int_next += p + 1; /* need space for (p+1) int's */
    1089             : 
    1090             :     /* compute coefficients mod p */
    1091     2574230 :     for (n = 0; n <= degree; n++) coeffs_mod_p[n] = umodiu(pel(c,n), p);
    1092             : 
    1093      367135 :     np = umael(squares,pn,coeffs_mod_p[0]+1);
    1094      367135 :     is_f_square[0] = np;
    1095    21015227 :     for (a = 1 ; a < p; a++)
    1096             :     {
    1097    20648092 :       ulong s = coeffs_mod_p[degree];
    1098    20648092 :       if ((degree+1)*args->bit_primes <= BITS_IN_LONG)
    1099             :       {
    1100   107292136 :         for (n = degree - 1 ; n >= 0 ; n--) s = s*a + coeffs_mod_p[n];
    1101             :         /* here, s < p^(degree+1) <= max. long */
    1102    17923592 :         s %= p;
    1103             :       }
    1104             :       else
    1105             :       {
    1106    16828148 :         for (n = degree - 1 ; n >= 0 ; n--)
    1107             :         {
    1108    14103648 :           s = s*a + coeffs_mod_p[n];
    1109    14103648 :           if (s+1 >= bound) s %= p;
    1110             :         }
    1111     2724500 :         s %= p;
    1112             :       }
    1113    20648092 :       if ((is_f_square[a] = mael(squares,pn,s+1))) np++;
    1114             :     }
    1115      367135 :     is_f_square[p] = odd(degree) || mael(squares,pn,coeffs_mod_p[degree]+1);
    1116             : 
    1117             :     /* check if there are no solutions mod p */
    1118      367135 :     if (np == 0 && !is_f_square[p]) return gc_long(av,p);
    1119             : 
    1120             :     /* Fill arrays with info for p */
    1121      367009 :     if (np < p)
    1122             :     { /* only when there is some information */
    1123             :       ulong i;
    1124      334509 :       ratpoints_sieve_entry *se = args->se_next;
    1125      857577 :       double r = is_f_square[p] ? ((double)(np*(p-1) + p))/((double)(p*p))
    1126      334509 :                                 : (double)np/(double)p;
    1127      334509 :       prec[pnp].r = r;
    1128      334509 :       args->se_next ++;
    1129      334509 :       se->p = p;
    1130      334509 :       se->is_f_square = is_f_square;
    1131      334509 :       se->inverses = gel(inverses,pn);
    1132      334509 :       se->offset = offsets[pn];
    1133      334509 :       se->sieve[0] = (ratpoints_bit_array *)sieves0[pn];
    1134    20126685 :       for (i = 1; i < p; i++) se->sieve[i] = NULL;
    1135      334509 :       prec[pnp].ssp = se;
    1136      334509 :       pnp++;
    1137             :     }
    1138             : 
    1139      367009 :     if ((args->flags & RATPOINTS_CHECK_DENOM)
    1140      320179 :          && fba + fdc < args->max_forbidden
    1141      320179 :          && !is_f_square[p])
    1142             :     { /* record forbidden divisors of the denominator */
    1143      147357 :       if (coeffs_mod_p[degree] == 0)
    1144             :       { /* leading coeff. divisible by p */
    1145             :         GEN r;
    1146           0 :         long v = Z_lvalrem(pel(c,degree), p, &r);
    1147             : 
    1148           0 :         if (odd(v) || !mael(squares,pn, umodiu(r,p)+1))
    1149             :         { /* Can only get something when valuation is odd
    1150             :              or when valuation is even and lcf is not a p-adic square.
    1151             :              Compute smallest n such that if v(den) >= n, the leading
    1152             :              term determines the valuation. Then we must have v(den) < n. */
    1153           0 :           long k, n = 1;
    1154           0 :           for (k = degree-1; k >= 0; k--)
    1155             :           {
    1156           0 :             if (coeffs_mod_p[k] == 0)
    1157             :             {
    1158           0 :               long t = 1 + v - Z_lval(pel(c,k), p);
    1159           0 :               long m = CEIL(t, (degree-k));
    1160           0 :               if (m > n) n = m;
    1161             :             }
    1162             :           }
    1163           0 :           if (n == 1)
    1164             :           {
    1165           0 :             forb_ba[fba].p     = p;
    1166           0 :             forb_ba[fba].start = sieves0[pn];
    1167           0 :             forb_ba[fba].end   = sieves0[pn]+p;
    1168           0 :             forb_ba[fba].curr  = forb_ba[fba].start;
    1169           0 :             fba++;
    1170             :           }
    1171             :           else
    1172             :           {
    1173           0 :             forbidden[fdc] = upowuu(p, n);
    1174           0 :             fdc++;
    1175             :           }
    1176             :         }
    1177             :       }
    1178             :       else /* leading coefficient is a nonsquare mod p */
    1179             :       { /* denominator divisible by p is excluded */
    1180      147357 :         forb_ba[fba].p     = p;
    1181      147357 :         forb_ba[fba].start = sieves0[pn];
    1182      147357 :         forb_ba[fba].end   = sieves0[pn]+p;
    1183      147357 :         forb_ba[fba].curr  = forb_ba[fba].start;
    1184      147357 :         fba++;
    1185             :       }
    1186             :     }
    1187             :   } /* end for pn */
    1188             : 
    1189             :   /* update sp2 and sp1 if necessary */
    1190       12232 :   if (args->sp2 > pnp)       args->sp2 = pnp;
    1191       12232 :   if (args->sp1 > args->sp2) args->sp1 = args->sp2;
    1192             : 
    1193             :   /* sort the array to get at the best primes */
    1194       12232 :   qsort(prec, pnp, sizeof(entry), compare_entries);
    1195             : 
    1196             :   /* put the sorted entries into sieve_list */
    1197      244178 :   for (n = 0; n < args->sp2; n++) sieve_list[n] = prec[n].ssp;
    1198             : 
    1199             :   /* terminate array of forbidden divisors */
    1200       12232 :   if (args->flags & RATPOINTS_CHECK_DENOM)
    1201             :   {
    1202             :     long n;
    1203             : 
    1204       10671 :     for (n = args->num_primes+1;
    1205       10671 :         fba + fdc < args->max_forbidden && n <= nbprime; n++)
    1206             :     {
    1207           0 :       ulong p = args->listprime[n];
    1208             : 
    1209           0 :       if (p*p > (ulong) args->b_high) break;
    1210           0 :       if (kroiu(pel(c,degree), p) == -1)
    1211             :       {
    1212           0 :         forb_ba[fba].p     = p;
    1213           0 :         forb_ba[fba].start = sieves0[n];
    1214           0 :         forb_ba[fba].end   = sieves0[n]+p;
    1215           0 :         forb_ba[fba].curr  = forb_ba[fba].start;
    1216           0 :         fba++;
    1217             :       }
    1218             :     }
    1219       10671 :     forb_ba[fba].p = 0; /* terminating zero */
    1220       10671 :     forbidden[fdc] = 0; /* terminating zero */
    1221       10671 :     args->max_forbidden = fba + fdc; /* note actual number */
    1222             :   }
    1223             : 
    1224       12232 :   if (fba + fdc == 0) args->flags &= ~RATPOINTS_CHECK_DENOM;
    1225       12232 :   return gc_long(av,0);
    1226             : }
    1227             : 
    1228             : /**************************************************************************
    1229             :  * The sieving procedure itself                                           *
    1230             :  **************************************************************************/
    1231             : static void
    1232    31760858 : sift(long b, ratpoints_bit_array *survivors, ratpoints_args *args,
    1233             :      bit_selection which_bits, ratpoints_bit_array bits16,
    1234             :      ratpoints_sieve_entry **sieve_list, long *bp_list, int *quit,
    1235             :      int process(long, long, GEN, void*, int*), void *info)
    1236             : {
    1237    31760858 :   pari_sp av = avma;
    1238    31760858 :   int do_setup = 1;
    1239    31760858 :   long k, height = args->height, nb;
    1240    31760858 :   sieve_spec *ssp = (sieve_spec *) stack_malloc(args->sp2*sizeof(sieve_spec));
    1241             : 
    1242    31760858 :   if (odd(b) == 0) which_bits = num_odd; /* even denominator */
    1243             : 
    1244             :   /* Note that b is new */
    1245    31760858 :   args->bc = NULL;
    1246    31760858 :   nb = 0;
    1247             : 
    1248    75483508 :   for (k = 0; k < args->num_inter; k++)
    1249             :   {
    1250    47692437 :     ratpoints_interval inter = args->domain[k];
    1251             :     long low, high;
    1252             : 
    1253             :     /* Determine relevant interval [low, high] of numerators. */
    1254    47692437 :     if (b*inter.low <= -height)
    1255    24024561 :       low = -height;
    1256             :     else
    1257             :     {
    1258    23667876 :       if (b*inter.low > height) break;
    1259    19699846 :       low = ceil(b*inter.low);
    1260             :     }
    1261    43724407 :     if (b*inter.up >= height)
    1262    19667260 :       high = height;
    1263             :     else
    1264             :     {
    1265    24057147 :       if (b*inter.up < -height) continue;
    1266    20154393 :       high = floor(b*inter.up);
    1267             :     }
    1268             : 
    1269    39821653 :     if (do_setup)
    1270             :     { /* set up the sieve information */
    1271             :       long n;
    1272             : 
    1273    29453868 :       do_setup = 0; /* only do it once for every b */
    1274   588683680 :       for (n = 0; n < args->sp2; n++)
    1275             :       {
    1276   559229812 :         ratpoints_sieve_entry *se = sieve_list[n];
    1277   559229812 :         long p = se->p;
    1278   559229812 :         long bp = bp_list[n];
    1279             :         ratpoints_bit_array *sptr;
    1280             : 
    1281   559229812 :         if (which_bits != num_all) /* divide by 2 mod p */
    1282   320662036 :           bp = odd(bp) ? (bp+p) >> 1 : bp >> 1;
    1283   559229812 :         sptr = se->sieve[bp];
    1284             : 
    1285   559229812 :         ssp[n].p = p;
    1286   559229812 :         ssp[n].offset = (which_bits == num_odd) ? se->offset : 0;
    1287             : 
    1288             :         /* copy if already initialized, else initialize */
    1289   571842208 :         ssp[n].ptr = sptr ? sptr : (p<BITS_IN_LONG?sieve_init1(p, se, bp, args)
    1290    12612396 :                                                   :sieve_init2(p, se, bp, args));
    1291   559229812 :         ssp[n].start = ssp[n].ptr;
    1292   559229812 :         ssp[n].end = ssp[n].ptr + p;
    1293             : 
    1294             :       }
    1295             :     }
    1296             : 
    1297    39821653 :     switch(which_bits)
    1298             :     {
    1299    17073730 :       case num_all: break;
    1300           0 :       case num_none: break;
    1301    18106385 :       case num_odd: low >>= 1; high--; high >>= 1; break;
    1302     4641538 :       case num_even: low++; low >>= 1; high >>= 1; break;
    1303             :     }
    1304             : 
    1305             :     /* now turn the bit interval into [low, high[ */
    1306    39821653 :     high++;
    1307             : 
    1308    39821653 :     if (low < high)
    1309             :     {
    1310    39819548 :       long w_low, w_high, w_low0, w_high0, range = args->array_size;
    1311             : 
    1312             :       /* Now the range of longwords (= bit_arrays) */
    1313    39819548 :       w_low = low >> RBA_SHIFT;
    1314    39819548 :       w_high = (high + (long)(RBA_LENGTH-1)) >> RBA_SHIFT;
    1315    39819548 :       w_low0 = w_low;
    1316    39819548 :       w_high0 = w_low0 + range;
    1317    86248414 :       for ( ; w_low0 < w_high; w_low0 = w_high0, w_high0 += range)
    1318             :       {
    1319             :         long i;
    1320    46430623 :         if (w_high0 > w_high) { w_high0 = w_high; range = w_high0 - w_low0; }
    1321             :         /* initialise the bits */
    1322  5202069841 :         for (i = range; i; i--) survivors[i-1] = bits16;
    1323             :         /* boundary words */
    1324    46430623 :         if (w_low0 == w_low)
    1325    39819548 :           MASKL(survivors,low - RBA_LENGTH * w_low)
    1326    46430623 :         if (w_high0 == w_high)
    1327    39819404 :           MASKU(&survivors[range-1], RBA_LENGTH * w_high - high)
    1328             : 
    1329             : #if (RATPOINTS_CHUNK > 1)
    1330   247869282 :         while(range%RATPOINTS_CHUNK != 0)
    1331   212317044 :           { survivors[range] = RBA(0); range++; w_high0++; }
    1332             : #endif
    1333    46430623 :         nb += _ratpoints_sift0(b, w_low0, w_high0, args, which_bits,
    1334             :                          survivors, &ssp[0], quit, process, info);
    1335    46430623 :         if (*quit) return;
    1336             :       }
    1337             :     }
    1338             :   }
    1339    31759101 :   if (nb==0) set_avma(av);
    1340             : }
    1341             : 
    1342             : /**************************************************************************
    1343             :  * Find points by looping over the denominators and sieving numerators    *
    1344             :  **************************************************************************/
    1345             : static void
    1346       12358 : find_points_work(ratpoints_args *args,
    1347             :                  int process(long, long, GEN, void*, int*), void *info)
    1348             : {
    1349       12358 :   int quit = 0;
    1350       12358 :   GEN c = args->cof;
    1351       12358 :   long degree = degpol(c);
    1352       12358 :   long nbprime = lg(args->listprime)-1;
    1353       12358 :   long height = args->height;
    1354             : 
    1355       12358 :   int point_at_infty = 0; /* indicates if there are points at infinity */
    1356       12358 :   int lcfsq = Z_issquare(pel(c,degree));
    1357             : 
    1358       12358 :   forbidden_entry *forb_ba = args->forb_ba;
    1359       12358 :   long *forbidden = args->forbidden;
    1360             :     /* The forbidden divisors, a zero-terminated array.
    1361             :        Used when degree is even and leading coefficient is not a square */
    1362             : 
    1363             :     /* These are used when degree is odd and leading coeff. is not +-1 */
    1364             : 
    1365             :   ratpoints_sieve_entry **sieve_list = (ratpoints_sieve_entry**)
    1366       12358 :      stack_malloc(nbprime*sizeof(ratpoints_sieve_entry*));
    1367       12358 :   bit_selection which_bits = num_all;
    1368             :   ulong den_bits;
    1369             :   ratpoints_bit_array num_bits[16];
    1370             : 
    1371       12358 :   args->flags &= RATPOINTS_FLAGS_INPUT_MASK;
    1372       12358 :   args->flags |= RATPOINTS_CHECK_DENOM;
    1373             : 
    1374             :   /* initialize memory management */
    1375       12358 :   args->se_next = args->se_buffer;
    1376       12358 :   args->ba_next = args->ba_buffer;
    1377       12358 :   args->int_next = args->int_buffer;
    1378             : 
    1379             :   /* Some sanity checks */
    1380       12358 :   args->num_inter = 0;
    1381             : 
    1382       12358 :   if (args->num_primes > nbprime) args->num_primes = nbprime;
    1383       12358 :   if (args->sp2 > args->num_primes) args->sp2 = args->num_primes;
    1384       12358 :   if (args->sp1 > args->sp2)        args->sp1 = args->sp2;
    1385             : 
    1386       12358 :   if (args->b_low < 1)  args->b_low = 1;
    1387       12358 :   if (args->b_high < 1) args->b_high = height;
    1388       12358 :   if (args->max_forbidden < 0)
    1389           0 :     args->max_forbidden = RATPOINTS_DEFAULT_MAX_FORBIDDEN;
    1390       12358 :   if (args->max_forbidden > nbprime)
    1391           0 :     args->max_forbidden = nbprime;
    1392       12358 :   if (args->array_size <= 0) args->array_size = RATPOINTS_ARRAY_SIZE;
    1393             :   {
    1394       12358 :     long s = 2*maxss(1,CEIL(height, BITS_IN_LONG));
    1395       12358 :     if (args->array_size > s) args->array_size = s;
    1396             :   }
    1397             :   /* make sure that array size is a multiple of RATPOINTS_CHUNK */
    1398       12358 :   args->array_size = CEIL(args->array_size, RATPOINTS_CHUNK)*RATPOINTS_CHUNK;
    1399             : 
    1400             :   /* Don't reverse if intervals are specified or limits for the denominator
    1401             :      are given */
    1402       12358 :   if (args->num_inter > 0 || args->b_low > 1 || args->b_high != height)
    1403          35 :     args->flags |= RATPOINTS_NO_REVERSE;
    1404             : 
    1405             :   /* Check if reversal of polynomial might be better:
    1406             :    * case 1: degree is even, but trailing coefficient is zero
    1407             :    * case 2: degree is even, leading coefficient is a square, but
    1408             :    *         trailing coefficient is not
    1409             :    * case 3: degree is odd, |leading coefficient| > 1,
    1410             :    *         trailing coefficient is zero, |coeff. of x| = 1 */
    1411       12358 :   if (!(args->flags & RATPOINTS_NO_REVERSE))
    1412             :   {
    1413       12323 :     if (!odd(degree) && degree>0)
    1414             :     {
    1415       11154 :       if (signe(pel(c,0)) == 0 && signe(pel(c,1))!=0)
    1416         224 :       { /* case 1 */
    1417             :         long n;
    1418         224 :         args->flags |= RATPOINTS_REVERSED;
    1419         896 :         for (n = 0; n < degree>>1; n++) swap(pel(c,n), pel(c,degree-n));
    1420         224 :         degree--;
    1421         224 :         setlg(c,degree+3);
    1422             :       }
    1423       10930 :       else if (lcfsq && !Z_issquare(pel(c,0)))
    1424             :       { /* case 2 */
    1425             :         long n;
    1426         735 :         args->flags |= RATPOINTS_REVERSED;
    1427        2940 :         for (n = 0; n < degree>>1; n++) swap(pel(c,n), pel(c,degree-n));
    1428         735 :         lcfsq = 0;
    1429             :       }
    1430             :     }
    1431             :     else
    1432             :     { /* odd degree, case 3*/
    1433        1169 :       if (!is_pm1(pel(c,degree)) && !signe(pel(c,0)) && is_pm1(pel(c,1)))
    1434             :       {
    1435             :         long n;
    1436           7 :         args->flags |= RATPOINTS_REVERSED;
    1437          21 :         for (n = 1; n <= degree>>1; n++) swap(pel(c,n),pel(c,degree+1-n));
    1438             :       }
    1439             :     }
    1440             :   }
    1441             : 
    1442             :   /* Deal with the intervals */
    1443       12358 :   if (args->num_inter == 0)
    1444             :   { /* default interval (effectively ]-oo,oo[) if none is given */
    1445       12358 :     args->domain = (ratpoints_interval*) stack_malloc(2*degree*sizeof(ratpoints_interval));
    1446       12358 :     args->domain[0].low = -height; args->domain[0].up = height;
    1447       12358 :     args->num_inter = 1;
    1448             :   }
    1449             : 
    1450       12358 :   ratpoints_compute_sturm(args);
    1451             : 
    1452             :   /* Point(s) at infinity? */
    1453       12358 :   if (odd(degree) || lcfsq)
    1454             :   {
    1455        1561 :     args->flags &= ~RATPOINTS_CHECK_DENOM;
    1456        1561 :     point_at_infty = 1;
    1457             :   }
    1458             : 
    1459             :   /* Can use only squares as denoms if degree is odd and poly is +-monic */
    1460       12358 :   if (odd(degree))
    1461             :   {
    1462        1407 :     GEN w1 = pel(c,degree);
    1463        1407 :     if (is_pm1(w1))
    1464          70 :       args->flags |= RATPOINTS_USE_SQUARES;
    1465             :     else /* set up information on divisors of leading coefficient */
    1466        1337 :       setup_us1(args, absi_shallow(w1));
    1467             :   }
    1468             : 
    1469             :   /* deal with f mod powers of 2 */
    1470       12358 :   which_bits = get_2adic_info(args, &den_bits, &num_bits[0]);
    1471             :   /* which_bits says whether to consider even and/or odd numerators
    1472             :    * when the denominator is odd.
    1473             :    *
    1474             :    * Bit k in den_bits is 0 if b congruent to k mod BITS_IN_LONG need
    1475             :    * not be considered as a denominator.
    1476             :    *
    1477             :    * Bit k in num_bits[b] is 0 is numerators congruent to
    1478             :    *  k (which_bits = den_all)
    1479             :    *  2k (which_bits = den_even)
    1480             :    *  2k+1 (which_bits = den_odd)
    1481             :    * need not be considered for denominators congruent to b mod 16. */
    1482             : 
    1483             :   /* set up the sieve data structure */
    1484       12358 :   if (sieving_info(args, sieve_list)) return;
    1485             : 
    1486             :   /* deal with point(s) at infinity */
    1487       12232 :   if (point_at_infty)
    1488             :   {
    1489        1561 :     long a = 1, b = 0;
    1490             : 
    1491        1561 :     if (args->flags & RATPOINTS_REVERSED) { a = 0; b = 1; }
    1492        1561 :     if (odd(degree))
    1493        1407 :       (void)process(a, b, gen_0, info, &quit);
    1494             :     else
    1495             :     {
    1496         154 :       GEN w0 = sqrti(pel(c,degree));
    1497         154 :       (void)process(a, b, w0, info, &quit);
    1498         154 :       (void)process(a, b, negi(w0), info, &quit);
    1499             :     }
    1500        1561 :     if (quit) return;
    1501             :   }
    1502             :   /* now do the sieving */
    1503             :   {
    1504             :     ratpoints_bit_array *survivors = (ratpoints_bit_array *)
    1505       12232 :       stack_malloc_align((args->array_size)*RBA_SIZE, RBA_SIZE);
    1506       12232 :     long *bp_list = (long *) new_chunk(args->sp2);
    1507       12232 :     if (args->flags & (RATPOINTS_USE_SQUARES | RATPOINTS_USE_SQUARES1))
    1508             :     {
    1509        1407 :       if (args->flags & RATPOINTS_USE_SQUARES)
    1510             :       /* need only take squares as denoms */
    1511             :       {
    1512             :         long b, bb;
    1513          70 :         long last_b = args->b_low;
    1514             :         long n;
    1515        1400 :         for (n = 0; n < args->sp2; n++)
    1516        1330 :           bp_list[n] = mod(args->b_low, sieve_list[n]->p);
    1517             : 
    1518        8771 :         for (b = 1; bb = b*b, bb <= args->b_high; b++)
    1519        8701 :           if (bb >= args->b_low)
    1520             :           {
    1521        8701 :             ratpoints_bit_array bits = num_bits[bb & 0xf];
    1522        8701 :             if (TEST(bits))
    1523             :             {
    1524             :               long n;
    1525        7805 :               long d = bb - last_b;
    1526             : 
    1527             :               /* fill bp_list */
    1528      156100 :               for (n = 0; n < args->sp2; n++)
    1529      148295 :                 bp_list[n] = mod(bp_list[n] + d, sieve_list[n]->p);
    1530        7805 :               last_b = bb;
    1531             : 
    1532        7805 :               sift(bb, survivors, args, which_bits, bits,
    1533             :                    sieve_list, &bp_list[0], &quit, process, info);
    1534        7805 :               if (quit) break;
    1535             :             }
    1536             :           }
    1537             :       }
    1538             :       else /* args->flags & RATPOINTS_USE_SQUARES1 */
    1539             :       {
    1540        1337 :         GEN den_info = args->den_info;
    1541        1337 :         GEN divisors = args->divisors;
    1542        1337 :         long ld = lg(divisors);
    1543             :         long k;
    1544             :         long b, bb;
    1545             : 
    1546        4249 :         for (k = 1; k < ld; k++)
    1547             :         {
    1548        2919 :           long d = divisors[k];
    1549        2919 :           long last_b = d;
    1550             :           long n;
    1551        2919 :           if (d > args->b_high) continue;
    1552       58240 :           for (n = 0; n < args->sp2; n++)
    1553       55328 :             bp_list[n] = mod(d, sieve_list[n]->p);
    1554             : 
    1555      276689 :           for (b = 1; bb = d*b*b, bb <= args->b_high; b++)
    1556             :           {
    1557      273784 :             if (bb >= args->b_low)
    1558             :             {
    1559      273763 :               int flag = 1;
    1560      273763 :               ratpoints_bit_array bits = num_bits[bb & 0xf];
    1561             : 
    1562      273763 :               if (EXT0(bits))
    1563             :               {
    1564      225911 :                 long i, n, l = lg(gel(den_info,1));
    1565      225911 :                 long d = bb - last_b;
    1566             : 
    1567             :                 /* fill bp_list */
    1568     4518220 :                 for (n = 0; n < args->sp2; n++)
    1569     4292309 :                   bp_list[n] = mod(bp_list[n] + d, sieve_list[n]->p);
    1570      225911 :                 last_b = bb;
    1571             : 
    1572      428603 :                 for(i = 1; i < l; i++)
    1573             :                 {
    1574      254198 :                   long v = z_lval(bb, mael(den_info,1,i));
    1575      254198 :                   if((v >= mael(den_info,3,i))
    1576      122899 :                       && odd(v + (mael(den_info,2,i)))) { flag = 0; break; }
    1577             :                 }
    1578      225911 :                 if (flag)
    1579             :                 {
    1580      174405 :                   sift(bb, survivors, args, which_bits, bits,
    1581             :                        sieve_list, &bp_list[0], &quit, process, info);
    1582      174405 :                   if (quit) break;
    1583             :                 }
    1584             :               }
    1585             :             }
    1586             :           }
    1587        2912 :           if (quit) break;
    1588             :         }
    1589             :       }
    1590             :     }
    1591             :     else
    1592             :     {
    1593       10825 :       if (args->flags & RATPOINTS_CHECK_DENOM)
    1594             :       {
    1595             :         long *forb;
    1596             :         long b;
    1597       10671 :         long last_b = args->b_low;
    1598             :         ulong b_bits;
    1599             :         long n;
    1600      213091 :         for (n = 0; n < args->sp2; n++)
    1601      202420 :           bp_list[n] = mod(args->b_low, sieve_list[n]->p);
    1602             :         {
    1603       10671 :           forbidden_entry *fba = &forb_ba[0];
    1604       10671 :           long b_low = args->b_low;
    1605       10671 :           long w_low = (b_low-1) >> TWOPOTBITS_IN_LONG;
    1606             : 
    1607       10671 :           b_bits = den_bits;
    1608      158014 :           while (fba->p)
    1609             :           {
    1610      147343 :             fba->curr = fba->start + mod(w_low, fba->p);
    1611      147343 :             b_bits &= *(fba->curr);
    1612      147343 :             fba++;
    1613             :           }
    1614       10671 :           b_bits >>= (b_low-1) & LONG_MASK;
    1615             :         }
    1616             : 
    1617   134973303 :         for (b = args->b_low; b <= args->b_high; b++)
    1618             :         {
    1619   134964375 :           ratpoints_bit_array bits = num_bits[b & 0xf];
    1620             : 
    1621   134964375 :           if ((b & LONG_MASK) == 0)
    1622             :           { /* next b_bits */
    1623     2401878 :             forbidden_entry *fba = &forb_ba[0];
    1624             : 
    1625     2401878 :             b_bits = den_bits;
    1626    37699074 :             while (fba->p)
    1627             :             {
    1628    35297196 :               fba->curr++;
    1629    35297196 :               if (fba->curr == fba->end) fba->curr = fba->start;
    1630    35297196 :               b_bits &= *(fba->curr);
    1631    35297196 :               fba++;
    1632             :             }
    1633             :           }
    1634             :           else
    1635   132562497 :             b_bits >>= 1;
    1636             : 
    1637   134964375 :           if (odd(b_bits) && EXT0(bits))
    1638             :           { /* check if denominator is excluded */
    1639    51190767 :             for (forb = &forbidden[0] ; *forb && (b % (*forb)); forb++)
    1640           0 :               continue;
    1641    51190767 :             if (*forb == 0 && rpjacobi(b, pel(c,degree)) == 1)
    1642             :             {
    1643    29901399 :               long n, d = b - last_b;
    1644             : 
    1645             :               /* fill bp_list */
    1646   597634699 :               for (n = 0; n < args->sp2; n++)
    1647             :               {
    1648   567733300 :                 long bp = bp_list[n] + d;
    1649   567733300 :                 long p = sieve_list[n]->p;
    1650             : 
    1651   638576595 :                 while (bp >= p) bp -= p;
    1652   567733300 :                 bp_list[n] = bp;
    1653             :               }
    1654    29901399 :               last_b = b;
    1655             : 
    1656    29901399 :               sift(b, survivors, args, which_bits, bits,
    1657             :                    sieve_list, &bp_list[0], &quit, process, info);
    1658    29901399 :               if (quit) break;
    1659             :             }
    1660             :           }
    1661             :         }
    1662             :       } /* if (args->flags & RATPOINTS_CHECK_DENOM) */
    1663             :       else
    1664             :       {
    1665             :         long b, n;
    1666         154 :         long last_b = args->b_low;
    1667        2947 :         for (n = 0; n < args->sp2; n++)
    1668        2793 :           bp_list[n] = mod(args->b_low, sieve_list[n]->p);
    1669     2179184 :         for (b = args->b_low; b <= args->b_high; b++)
    1670             :         {
    1671     2179037 :           ratpoints_bit_array bits = num_bits[b & 0xf];
    1672     2179037 :           if (EXT0(bits))
    1673             :           {
    1674             :             long n;
    1675     1677249 :             long d = b - last_b;
    1676             : 
    1677             :             /* fill bp_list */
    1678    33544581 :             for (n = 0; n < args->sp2; n++)
    1679             :             {
    1680    31867332 :               long bp = bp_list[n] + d;
    1681    31867332 :               long p = sieve_list[n]->p;
    1682             : 
    1683    32980773 :               while (bp >= p) bp -= p;
    1684    31867332 :               bp_list[n] = bp;
    1685             :             }
    1686     1677249 :             last_b = b;
    1687             : 
    1688     1677249 :             sift(b, survivors, args, which_bits, bits,
    1689             :                  sieve_list, &bp_list[0], &quit, process, info);
    1690     1677249 :             if (quit) break;
    1691             :           }
    1692             :         }
    1693             :       }
    1694             :     }
    1695             :   }
    1696             : }
    1697             : 
    1698             : static GEN
    1699       86191 : vec_append_grow(GEN z, long i, GEN x)
    1700             : {
    1701       86191 :   long n = lg(z)-1;
    1702       86191 :   if (i > n)
    1703             :   {
    1704        1435 :     n <<= 1;
    1705        1435 :     z = vec_lengthen(z,n);
    1706             :   }
    1707       86191 :   gel(z,i) = x;
    1708       86191 :   return z;
    1709             : }
    1710             : 
    1711             : struct points
    1712             : {
    1713             :   GEN z;
    1714             :   long i, f;
    1715             : };
    1716             : 
    1717             : static int
    1718       89208 : process(long a, long b, GEN y, void *info0, int *quit)
    1719             : {
    1720       89208 :   struct points *p = (struct points *) info0;
    1721       89208 :   if (b==0 && (p->f&2L)) return 0;
    1722       86191 :   *quit = (p->f&1);
    1723       86191 :   p->z = vec_append_grow(p->z, ++p->i, mkvec3(stoi(a),y,stoi(b)));
    1724       86191 :   return 1;
    1725             : }
    1726             : 
    1727             : static int
    1728       12365 : args_h(ratpoints_args *args, GEN D)
    1729             : {
    1730       12365 :   long H, h, l = 1;
    1731             :   GEN L;
    1732       12365 :   switch(typ(D))
    1733             :   {
    1734       12323 :     case t_INT: if (signe(D) <= 0) return 0;
    1735       12323 :       H = h = itos(D); break;
    1736          42 :     case t_VEC: if (lg(D) != 3) return 0;
    1737          42 :       H = gtos(gel(D,1));
    1738          42 :       L = gel(D,2);
    1739          42 :       if (typ(L)==t_INT) { h = itos(L); break; }
    1740          14 :       if (typ(L)==t_VEC && lg(L)==3)
    1741             :       {
    1742           7 :         l = gtos(gel(L,1));
    1743           7 :         h = gtos(gel(L,2)); break;
    1744             :       }
    1745           7 :     default: return 0;
    1746             :   }
    1747       12358 :   args->height = H;
    1748       12358 :   args->b_low  = l;
    1749       12358 :   args->b_high = h; return 1;
    1750             : }
    1751             : 
    1752             : static GEN
    1753       12365 : ZX_hyperellratpoints(GEN P, GEN h, long flag)
    1754             : {
    1755       12365 :   pari_sp av = avma;
    1756             :   ratpoints_args args;
    1757             :   struct points data;
    1758       12365 :   ulong flags = 0;
    1759             : 
    1760       12365 :   if (!args_h(&args, h))
    1761             :   {
    1762           7 :     pari_err_TYPE("hyperellratpoints", h);
    1763             :     return NULL;/*LCOV_EXCL_LINE*/
    1764             :   }
    1765       12358 :   find_points_init(&args, RATPOINTS_DEFAULT_BIT_PRIMES);
    1766             : 
    1767       12358 :   args.cof           = shallowcopy(P);
    1768       12358 :   args.num_inter     = 0;
    1769       12358 :   args.sp1           = RATPOINTS_DEFAULT_SP1;
    1770       12358 :   args.sp2           = RATPOINTS_DEFAULT_SP2;
    1771       12358 :   args.array_size    = RATPOINTS_ARRAY_SIZE;
    1772       12358 :   args.num_primes    = RATPOINTS_DEFAULT_NUM_PRIMES;
    1773       12358 :   args.bit_primes    = RATPOINTS_DEFAULT_BIT_PRIMES;
    1774       12358 :   args.max_forbidden = RATPOINTS_DEFAULT_MAX_FORBIDDEN;
    1775       12358 :   args.flags         = flags;
    1776       12358 :   data.z = cgetg(17,t_VEC);
    1777       12358 :   data.i = 0;
    1778       12358 :   data.f = flag;
    1779       12358 :   find_points_work(&args, process, (void *)&data);
    1780             : 
    1781       12358 :   setlg(data.z, data.i+1);
    1782       12358 :   return gerepilecopy(av, data.z);
    1783             : }
    1784             : 
    1785             : /* The ordinates of the points returned need to be divided by den
    1786             :  * by the caller to get the actual solutions */
    1787             : static GEN
    1788       12365 : QX_hyperellratpoints(GEN P, GEN h, long flag, GEN *den)
    1789             : {
    1790       12365 :   GEN Q = Q_remove_denom(P, den);
    1791       12365 :   if (*den) Q = ZX_Z_mul(Q, *den);
    1792       12365 :   return ZX_hyperellratpoints(Q, h, flag);
    1793             : }
    1794             : 
    1795             : static GEN
    1796         168 : ZX_homogenous_evalpow(GEN Q, GEN A, GEN B)
    1797             : {
    1798         168 :   pari_sp av = avma;
    1799         168 :   long i, d = degpol(Q);
    1800         168 :   GEN s = gel(Q, d+2);
    1801         672 :   for (i = d-1; i >= 0; i--)
    1802         504 :     s = addii(mulii(s, A), mulii(gel(B,d+1-i), gel(Q,i+2)));
    1803         168 :   return d? gerepileupto(av, s): s;
    1804             : }
    1805             : 
    1806             : static GEN
    1807          70 : to_RgX(GEN a, long v) { return typ(a)==t_POL? a: scalarpol(a,v); }
    1808             : 
    1809             : static int
    1810       11483 : hyperell_check(GEN PQ, GEN *P, GEN *Q)
    1811             : {
    1812       11483 :   *P = *Q = NULL;
    1813       11483 :   if (typ(PQ) == t_POL)
    1814             :   {
    1815       11448 :     if (!RgX_is_QX(PQ)) return 0;
    1816       11448 :     *P = PQ;
    1817             :   }
    1818             :   else
    1819             :   {
    1820          35 :     long v = gvar(PQ);
    1821          35 :     if (v == NO_VARIABLE || typ(PQ) != t_VEC || lg(PQ) != 3) return 0;
    1822          35 :     *P = to_RgX(gel(PQ,1), v); if (!RgX_is_QX(*P)) return 0;
    1823          35 :     *Q = to_RgX(gel(PQ,2), v); if (!RgX_is_QX(*Q)) return 0;
    1824          35 :     if (!signe(*Q)) *Q = NULL;
    1825             :   }
    1826       11483 :   return 1;
    1827             : }
    1828             : 
    1829             : GEN
    1830       11483 : hyperellratpoints(GEN PQ, GEN h, long flag)
    1831             : {
    1832       11483 :   pari_sp av = avma;
    1833             :   GEN P, Q, H, L, den, denQ;
    1834             :   long i, l, dy, dQ;
    1835             : 
    1836       11483 :   if (flag<0 || flag>1) pari_err_FLAG("ellratpoints");
    1837       11483 :   if (!hyperell_check(PQ, &P, &Q)) pari_err_TYPE("hyperellratpoints",PQ);
    1838       11483 :   if (!Q)
    1839             :   {
    1840       11462 :     L = QX_hyperellratpoints(P, h, flag|2L, &den);
    1841       11462 :     dy = (degpol(P)+1)>>1;
    1842       11462 :     l = lg(L);
    1843       25378 :     for (i = 1; i < l; i++)
    1844             :     {
    1845       13916 :       GEN Li = gel(L,i), x = gel(Li,1), y = gel(Li,2), z = gel(Li,3);
    1846       13916 :       GEN zdy = powiu(z, dy);
    1847       13916 :       if (den) zdy = mulii(zdy, den);
    1848       13916 :       gel(L,i) = mkvec2(gdiv(x,z), gdiv(y, zdy));
    1849             :     }
    1850       11462 :     return gerepilecopy(av, L);
    1851             :   }
    1852          21 :   H = RgX_add(RgX_muls(P,4), RgX_sqr(Q));
    1853          21 :   dy = (degpol(H)+1)>>1; dQ = degpol(Q);
    1854          21 :   L = QX_hyperellratpoints(H, h, flag|2L, &den);
    1855          21 :   Q = Q_remove_denom(Q, &denQ);
    1856          21 :   l = lg(L);
    1857         189 :   for (i = 1; i < l; i++)
    1858             :   {
    1859         168 :     GEN Li = gel(L,i), x = gel(Li,1), y = gel(Li,2), z = gel(Li,3);
    1860         168 :     GEN Qx, B = gpowers(z, dQ), zdy = powiu(z, dy), dQx = gel(B, dQ+1);
    1861         168 :     if (denQ) dQx = mulii(dQx, denQ);
    1862         168 :     Qx = gdiv(ZX_homogenous_evalpow(Q, x, B), dQx);
    1863         168 :     if (den) zdy = mulii(zdy, den);
    1864         168 :     gel(L,i) = mkvec2(gdiv(x,z), gmul2n(gsub(gdiv(y,zdy),Qx),-1));
    1865             :   }
    1866          21 :   return gerepilecopy(av, L);
    1867             : }
    1868             : 
    1869             : GEN
    1870         882 : ellratpoints(GEN E, GEN h, long flag)
    1871             : {
    1872         882 :   pari_sp av = avma;
    1873             :   GEN L, a1, a3, den;
    1874             :   long i, l;
    1875         882 :   checkell_Q(E);
    1876         882 :   if (flag<0 || flag>1) pari_err_FLAG("ellratpoints");
    1877         882 :   if (!RgV_is_QV(vecslice(E,1,5))) pari_err_TYPE("ellratpoints",E);
    1878         882 :   a1 = ell_get_a1(E);
    1879         882 :   a3 = ell_get_a3(E);
    1880         882 :   L = QX_hyperellratpoints(ec_bmodel(E,0), h, flag|2L, &den);
    1881         875 :   l = lg(L);
    1882       72982 :   for (i = 1; i < l; i++)
    1883             :   {
    1884       72107 :     GEN P, Li = gel(L,i), x = gel(Li,1), y = gel(Li,2), z = gel(Li,3);
    1885       72107 :     if (!signe(z))
    1886           0 :       P = ellinf();
    1887             :     else
    1888             :     {
    1889       72107 :       GEN z2 = sqri(z);
    1890       72107 :       if (den) y = gdiv(y, den);
    1891       72107 :       y = gsub(y, gadd(gmul(a1, mulii(x,z)), gmul(a3,z2)));
    1892       72107 :       P = mkvec2(gdiv(x,z), gdiv(y,shifti(z2,1)));
    1893             :     }
    1894       72107 :     gel(L,i) = P;
    1895             :   }
    1896         875 :   return gerepilecopy(av, L);
    1897             : }

Generated by: LCOV version 1.16