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 2737965 : sieve_init1(long p, ratpoints_sieve_entry *se1, long b1, ratpoints_args *args1)
213 : {
214 2737965 : ratpoints_sieve_entry *se = se1;
215 2737965 : ratpoints_args *args = args1;
216 2737965 : int *isfs = se->is_f_square;
217 2737965 : long b = b1;
218 2737965 : long lmp = BITS_IN_LONG % p;
219 2737965 : long ldp = BITS_IN_LONG / p;
220 2737965 : long p1 = (ldp + 1) * p;
221 2737965 : long diff_shift = p1 & LONG_MASK;
222 2737965 : long diff = BITS_IN_LONG - diff_shift;
223 : ulong help0;
224 : long a;
225 2737965 : long d = se->inverses[b];
226 2737965 : long ab = 0; /* a/b mod p */
227 2737965 : ulong test = 1UL;
228 2737965 : ulong he0 = 0UL;
229 121734264 : for (a = 0; a < p; a++)
230 : {
231 118996299 : if (isfs[ab]) he0 |= test;
232 118996299 : ab += d;
233 118996299 : if (ab >= p) ab -= p;
234 118996299 : test <<= 1;
235 : }
236 2737965 : help0 = he0;
237 : {
238 : ulong help1;
239 : /* repeat bit pattern floor(BITS_IN_LONG/p) times */
240 2737965 : 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 2737965 : help1 = pattern >> lmp;
247 6249223 : for (i = p; i < BITS_IN_LONG; i <<= 1)
248 3511258 : 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 2737965 : ulong *si = (ulong *)args->ba_next;
254 :
255 2737965 : args->ba_next += p + RATPOINTS_CHUNK-1;
256 : /* copy the first chunk into sieve[b][] */
257 2737965 : si[0] = help0;
258 : /* now keep repeating the bit pattern,
259 : rotating it in help0/help1 */
260 118996299 : for (a = 1 ; a < p; a++)
261 : {
262 116258334 : ulong temp = help0 >> diff;
263 116258334 : help0 = help1 | (help0 << diff_shift);
264 116258334 : si[a] = help0;
265 116258334 : help1 = temp;
266 : }
267 314875500 : CODE_INIT_SIEVE_COPY
268 : /* set sieve array */
269 2737965 : se->sieve[b] = (ratpoints_bit_array *)si;
270 2737965 : return (ratpoints_bit_array *)si;
271 : }
272 : }
273 : }
274 :
275 : /* This is for p > BITS_IN_LONG */
276 : static ratpoints_bit_array *
277 9938663 : sieve_init2(long p, ratpoints_sieve_entry *se1, long b1, ratpoints_args *args1)
278 : {
279 9938663 : pari_sp av = avma;
280 9938663 : ratpoints_sieve_entry *se = se1;
281 9938663 : ratpoints_args *args = args1;
282 9938663 : int *isfs = se->is_f_square;
283 9938663 : long b = b1;
284 : /* long ldp = 0; = BITS_IN_LONG / p */
285 : /* long p1 = p; = (ldp + 1) * p; */
286 9938663 : long wp = p >> TWOPOTBITS_IN_LONG;
287 9938663 : long diff_shift = p & LONG_MASK;
288 9938663 : long diff = BITS_IN_LONG - diff_shift;
289 9938663 : ulong *help = (ulong *) new_chunk((p>>TWOPOTBITS_IN_LONG) + 2);
290 :
291 : /* initialize help */
292 : {
293 9938663 : ulong *he = &help[0];
294 9938663 : ulong *he1 = &he[(p>>TWOPOTBITS_IN_LONG) + 2];
295 41996437 : while (he1 != he) { he1--; *he1 = 0UL; }
296 : }
297 9938663 : { ulong work = 0UL;
298 : long a;
299 9938663 : long ab = 0; /* a/b mod p */
300 9938663 : long d = se->inverses[b];
301 9938663 : long n = 0;
302 9938663 : ulong test = 1UL;
303 965864274 : for (a = 0; a < p; )
304 : {
305 955925611 : if (isfs[ab]) work |= test;
306 955925611 : ab += d;
307 955925611 : if (ab >= p) ab -= p;
308 955925611 : test <<= 1;
309 955925611 : a++;
310 955925611 : if ((a & LONG_MASK) == 0)
311 12180448 : { help[n] = work; n++; work = 0UL; test = 1UL; }
312 : }
313 9938663 : 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 9938663 : ulong *si = (ulong *)args->ba_next;
320 : long a1;
321 : long a;
322 :
323 9938663 : args->ba_next += p + RATPOINTS_CHUNK-1;
324 : /* copy the first chunk from help[] into sieve[num][b][] */
325 22119111 : for (a = 0; a < wp; a++) si[a] = help[a];
326 : /* now keep repeating the bit pattern, rotating it in help */
327 953683826 : for (a1 = a ; a < p; a++)
328 : {
329 943745163 : long t = (a1 == wp) ? 0 : a1+1;
330 943745163 : help[a1] |= help[t]<<diff_shift;
331 943745163 : si[a] = help[a1];
332 943745163 : a1 = t;
333 943745163 : help[a1] >>= diff;
334 : }
335 1865556852 : CODE_INIT_SIEVE_COPY
336 : /* set sieve array */
337 9938663 : se->sieve[b] = (ratpoints_bit_array *)si;
338 9938663 : set_avma(av);
339 9938663 : return (ratpoints_bit_array *)si;
340 : }
341 : }
342 :
343 : static GEN
344 12428 : gen_squares(GEN listprime)
345 : {
346 12428 : long nbprime = lg(listprime)-1;
347 12428 : GEN sq = cgetg(nbprime+1, t_VEC);
348 : long n;
349 385268 : for (n = 1; n <= nbprime; n++)
350 : {
351 372840 : ulong i, p = uel(listprime,n);
352 372840 : GEN w = zero_zv(p), work = w+1;
353 372840 : work[0] = 1;
354 : /* record nonzero squares mod p, p odd */
355 10862072 : for (i = 1; i < p; i += 2) work[(i*i) % p] = 1;
356 372840 : gel(sq, n) = w;
357 : }
358 12428 : return sq;
359 : }
360 :
361 : static GEN
362 12428 : gen_offsets(GEN P)
363 : {
364 12428 : long n, l = lg(P);
365 12428 : GEN of = cgetg(l, t_VEC);
366 385268 : for (n = 1; n < l; n++)
367 : {
368 372840 : ulong p = uel(P,n);
369 372840 : uel(of, n) = Fl_inv((2*RBA_LENGTH)%p, p);
370 : }
371 12428 : return of;
372 : }
373 :
374 : static GEN
375 12428 : gen_inverses(GEN P)
376 : {
377 12428 : long n, l = lg(P);
378 12428 : GEN iv = cgetg(l, t_VEC);
379 385268 : for (n = 1; n < l; n++)
380 : {
381 372840 : ulong i, p = uel(P,n);
382 372840 : GEN w = cgetg(p, t_VECSMALL);
383 21351304 : for (i = 1; i < p; i++) uel(w, i) = Fl_inv(i, p);
384 372840 : gel(iv, n) = w;
385 : }
386 12428 : return iv;
387 : }
388 :
389 : static ulong **
390 12428 : gen_sieves0(GEN listprime)
391 : {
392 : long n;
393 12428 : long nbprime = lg(listprime)-1;
394 12428 : ulong ** w = (ulong**) new_chunk(nbprime+1);
395 385268 : for (n = 1; n <= nbprime; n++)
396 : {
397 372840 : ulong a, p = uel(listprime,n);
398 372840 : ulong *si = (ulong *) stack_malloc_align((p+RATPOINTS_CHUNK-1)*RBA_SIZE, RBA_SIZE);
399 21724144 : for (a = 0; a < p; a++) si[a] = ~0UL;
400 22533480 : for (a = 0; a < BITS_IN_LONG; a++)
401 22160640 : uel(si,(p*a)>>TWOPOTBITS_IN_LONG) &= ~(1UL<<((p*a) & LONG_MASK));
402 46524096 : CODE_INIT_SIEVE_COPY
403 372840 : w[n] = si;
404 : }
405 12428 : return w;
406 : }
407 :
408 : static void
409 12428 : gen_sieve(ratpoints_args *args)
410 : {
411 12428 : GEN listprimes = args->listprime;
412 12428 : args->offsets = gen_offsets(listprimes);
413 12428 : args->inverses = gen_inverses(listprimes);
414 12428 : args->sieves0 = gen_sieves0(listprimes);
415 12428 : }
416 :
417 : static GEN
418 12428 : ZX_positive_region(GEN P, long h, long bitprec)
419 : {
420 12428 : long prec = nbits2prec(bitprec);
421 12428 : GEN it = mkvec2(stoi(-h),stoi(h));
422 12428 : GEN R = realroots(P, it, prec);
423 12428 : long nR = lg(R)-1;
424 12428 : long s = signe(ZX_Z_eval(P,gel(it,1)));
425 12428 : long i=1, j;
426 : GEN iv, st, en;
427 12428 : if (s<0 && nR==0) return NULL;
428 11742 : iv = cgetg(((nR+1+(s>=0))>>1)+1, t_VEC);
429 11742 : if (s>=0) st = itor(gel(it,1),prec);
430 5138 : else { st = gel(R,i); i++; }
431 18143 : for (j=1; i<nR; j++)
432 : {
433 6401 : gel(iv, j) = mkvec2(st, gel(R,i));
434 6401 : st = gel(R,i+1);
435 6401 : i+=2;
436 : }
437 11742 : if (i==nR) en = gel(R,i); else en = itor(gel(it,2),prec);
438 11742 : gel(iv,j) = mkvec2(st, en);
439 11742 : return iv;
440 : }
441 :
442 : static long
443 12428 : posint(ratpoints_interval *ivlist, GEN P, long h)
444 : {
445 12428 : GEN R = ZX_positive_region(P, h, 53);
446 12428 : const double eps = 1e-5;
447 : long nR, i;
448 :
449 12428 : if (!R) return 0;
450 11742 : nR = lg(R)-1;
451 11742 : i = 1;
452 29885 : for (i=1; i<=nR; i++)
453 : {
454 18143 : ivlist[i-1].low = rtodbl(gmael(R,i,1))-eps;
455 18143 : ivlist[i-1].up = rtodbl(gmael(R,i,2))+eps;
456 : }
457 11742 : return nR;
458 : }
459 :
460 : static long
461 12428 : ratpoints_compute_sturm(ratpoints_args *args)
462 : {
463 12428 : ratpoints_interval *ivlist = args->domain;
464 12428 : args->num_inter = posint(ivlist, args->cof, (long) ivlist[0].up);
465 12428 : return args->num_inter;
466 : }
467 :
468 : /**************************************************************************
469 : * Try to avoid divisions *
470 : **************************************************************************/
471 : INLINE long
472 888509854 : mod(long a, long b)
473 : {
474 888509854 : long b1 = b << 4; /* b1 = 16*b */
475 :
476 888509854 : if (a < -b1) { a %= b; if (a < 0) { a += b; } return a ; }
477 877811229 : if (a < 0) { a += b1; }
478 355712818 : else { if (a >= b1) { return a % b; } }
479 870405994 : b1 >>= 1; /* b1 = 8*b */
480 870405994 : if (a >= b1) { a -= b1; }
481 870405994 : b1 >>= 1; /* b1 = 4*b */
482 870405994 : if (a >= b1) { a -= b1; }
483 870405994 : b1 >>= 1; /* b1 = 2*b */
484 870405994 : if (a >= b1) { a -= b1; }
485 870405994 : if (a >= b) { a -= b; }
486 870405994 : return a;
487 : }
488 :
489 : static void
490 2326442 : set_bc(long b, ratpoints_args *args)
491 : {
492 2326442 : GEN w0 = gen_1;
493 2326442 : GEN c = args->cof, bc;
494 2326442 : long k, degree = degpol(c);
495 2326442 : bc = cgetg(degree+2, t_POL);
496 11880339 : for (k = degree-1; k >= 0; k--)
497 : {
498 9553897 : w0 = muliu(w0, b);
499 9553897 : gel(bc,k+2) = mulii(gel(c,k+2), w0);
500 : }
501 2326442 : args->bc = bc;
502 2326442 : }
503 :
504 : /**************************************************************************
505 : * Check a `survivor' of the sieve if it really gives a point. *
506 : **************************************************************************/
507 :
508 : static long
509 3282373 : _ratpoints_check_point(long a, long b, ratpoints_args *args, int *quit,
510 : int process(long, long, GEN, void*, int*), void *info)
511 : {
512 3282373 : pari_sp av = avma;
513 3282373 : GEN w0, w2, c = args->cof, bc = args->bc;
514 3282373 : long k, degree = degpol(c);
515 3282373 : 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 3282373 : w2 = gel(c, degree+2);
520 16933659 : for (k = degree-1; k >= 0; k--)
521 : {
522 13651286 : w2 = mulis(w2, a);
523 13651286 : w2 = addii(w2, gel(bc,k+2));
524 : }
525 3282373 : if (odd(degree)) w2 = muliu(w2, b);
526 : /* check if f(x,z) is a square; if so, process the point(s) */
527 3282373 : if (signe(w2) >= 0 && Z_issquareall(w2, &w0))
528 : {
529 44898 : 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 43680 : else (void)process(a, b, w0, info, quit);
535 44898 : 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 44898 : return 1;
546 : }
547 3237475 : set_avma(av);
548 3237475 : return 0;
549 : }
550 :
551 : /**************************************************************************
552 : * The inner loop of the sieving procedure *
553 : **************************************************************************/
554 : static long
555 46550975 : _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 46550975 : long sp1 = args->sp1, sp2 = args->sp2;
561 46550975 : 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 35655390 : ratpoints_bit_array *surv = survivors;
568 :
569 : /* first set the start fields for the first and second phases of sieving */
570 712472724 : for(n = 0; n < sp2; n++)
571 676817334 : 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 243543360 : 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 207887970 : ratpoints_bit_array reg0 = surv[0];
578 207887970 : ratpoints_bit_array reg1 = surv[1];
579 207887970 : ratpoints_bit_array reg2 = surv[2];
580 207887970 : ratpoints_bit_array reg3 = surv[3];
581 207887970 : ratpoints_bit_array reg4 = surv[4];
582 207887970 : ratpoints_bit_array reg5 = surv[5];
583 207887970 : ratpoints_bit_array reg6 = surv[6];
584 207887970 : ratpoints_bit_array reg7 = surv[7];
585 207887970 : ratpoints_bit_array reg8 = surv[8];
586 207887970 : ratpoints_bit_array reg9 = surv[9];
587 207887970 : ratpoints_bit_array reg10 = surv[10];
588 207887970 : ratpoints_bit_array reg11 = surv[11];
589 207887970 : ratpoints_bit_array reg12 = surv[12];
590 207887970 : ratpoints_bit_array reg13 = surv[13];
591 207887970 : ratpoints_bit_array reg14 = surv[14];
592 207887970 : ratpoints_bit_array reg15 = surv[15];
593 :
594 2494655244 : for(n = 0; n < sp1; n++)
595 : { /* retrieve the pointer to the beginning of the relevant bits */
596 2286767274 : ratpoints_bit_array *siv1 = sieves[n].start;
597 2286767274 : reg0 &= *siv1++;
598 2286767274 : reg1 &= *siv1++;
599 2286767274 : reg2 &= *siv1++;
600 2286767274 : reg3 &= *siv1++;
601 2286767274 : reg4 &= *siv1++;
602 2286767274 : reg5 &= *siv1++;
603 2286767274 : reg6 &= *siv1++;
604 2286767274 : reg7 &= *siv1++;
605 2286767274 : reg8 &= *siv1++;
606 2286767274 : reg9 &= *siv1++;
607 2286767274 : reg10 &= *siv1++;
608 2286767274 : reg11 &= *siv1++;
609 2286767274 : reg12 &= *siv1++;
610 2286767274 : reg13 &= *siv1++;
611 2286767274 : reg14 &= *siv1++;
612 2286767274 : reg15 &= *siv1++;
613 :
614 : /* update the pointer for the next round
615 : * (RATPOINTS_CHUNK-1 bit-arrays after sieves[n].end) */
616 3220498164 : while(siv1 >= sieves[n].end) siv1 -= sieves[n].p;
617 2286767274 : sieves[n].start = siv1;
618 : }
619 : /* store the contents of the registers back into memory */
620 207887970 : surv[0] = reg0;
621 207887970 : surv[1] = reg1;
622 207887970 : surv[2] = reg2;
623 207887970 : surv[3] = reg3;
624 207887970 : surv[4] = reg4;
625 207887970 : surv[5] = reg5;
626 207887970 : surv[6] = reg6;
627 207887970 : surv[7] = reg7;
628 207887970 : surv[8] = reg8;
629 207887970 : surv[9] = reg9;
630 207887970 : surv[10] = reg10;
631 207887970 : surv[11] = reg11;
632 207887970 : surv[12] = reg12;
633 207887970 : surv[13] = reg13;
634 207887970 : surv[14] = reg14;
635 207887970 : surv[15] = reg15;
636 : }
637 : #else /* RATPOINTS_CHUNK not between 2 and 16 */
638 10895585 : long range = w_high - w_low;
639 130746954 : for (n = 0; n < sp1; n++)
640 : {
641 119851369 : ratpoints_bit_array *sieve_n = sieves[n].ptr;
642 119851369 : long p = sieves[n].p;
643 119851369 : long r = mod(-w_low-sieves[n].offset, p);
644 119851369 : ratpoints_bit_array *surv = survivors;
645 :
646 119851369 : if (w_high < w_low + r)
647 : { /* if we get here, r > 0, since w_high >= w_low always */
648 6972178 : ratpoints_bit_array *siv1 = &sieve_n[p-r];
649 6972178 : ratpoints_bit_array *siv0 = siv1 + range;
650 :
651 207983792 : while (siv1 != siv0) { *surv &= *siv1++; surv++; }
652 : }
653 : else
654 : {
655 112879191 : ratpoints_bit_array *siv1 = &sieve_n[p-r];
656 112879191 : ratpoints_bit_array *surv_end = &survivors[range - p];
657 : long i;
658 3567574027 : for (i = r; i; i--) { *surv &= *siv1++; surv++; }
659 112879191 : siv1 -= p;
660 572281051 : while (surv <= surv_end)
661 : {
662 15773570368 : for (i = p; i; i--) { *surv &= *siv1++; surv++; }
663 459401860 : siv1 -= p;
664 : }
665 112879191 : surv_end += p;
666 3645060227 : while (surv < surv_end) { *surv &= *siv1++; surv++; }
667 : }
668 : }
669 : /* initialize pointers in sieve for the second phase */
670 97848726 : for(n = sp1; n < sp2; n++)
671 86953141 : 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 46550975 : surv0 = &survivors[0];
676 5418301297 : for (i = w_low; i < w_high; i++, base++)
677 : {
678 5371752079 : ratpoints_bit_array nums = *surv0++;
679 5371752079 : sieve_spec *ssp = &sieves[sp1];
680 : long n;
681 :
682 5870273474 : for (n = sp2-sp1; n && TEST(nums); n--)
683 : {
684 498521395 : ratpoints_bit_array *ptr = (ssp->start) + base;
685 498521395 : long p = ssp->p;
686 1325310826 : while(ptr >= ssp->end) ptr -= p;
687 498521395 : nums &= *ptr;
688 498521395 : ssp++;
689 : }
690 :
691 : /* Check the survivors of the sieve if they really give points */
692 5371752079 : if (TEST(nums))
693 : {
694 : long a0, a, d;
695 10790555 : ulong nums0 = EXT0(nums);
696 : /* a will be the numerator corresponding to the selected bit */
697 10790555 : if (which_bits == num_all)
698 : {
699 6949351 : d = 1; a0 = i * RBA_LENGTH;
700 : }
701 : else
702 : {
703 3841204 : d = 2; a0 = i * 2 * RBA_LENGTH;
704 3841204 : if (which_bits == num_odd) a0++;
705 : }
706 134513795 : for (a = a0; nums0; a += d, nums0 >>= 1)
707 : { /* test one bit */
708 123724301 : if (odd(nums0) && ugcd(labs(a), absb)==1)
709 : {
710 1876711 : if (!args->bc) set_bc(b, args);
711 1876711 : nb += _ratpoints_check_point(a, b, args, quit, process, info);
712 1876711 : if (*quit) return nb;
713 : }
714 : }
715 : #ifdef RBA_USE_VX
716 : {
717 9247380 : ulong k, da = d<<TWOPOTBITS_IN_LONG;
718 18494064 : for (k = 1; k < RBA_PACK; k++)
719 : {
720 9247380 : ulong nums1 = EXT(nums,k);
721 9247380 : a0 += da;
722 112235190 : for (a = a0; nums1; a += d, nums1 >>= 1)
723 : { /* test one bit */
724 102988506 : if (odd(nums1) && ugcd(labs(a), absb)==1)
725 : {
726 1405662 : if (!args->bc) set_bc(b, args);
727 1405662 : nb += _ratpoints_check_point(a, b, args, quit, process, info);
728 1405662 : if (*quit) return nb;
729 : }
730 : }
731 : }
732 : }
733 : #endif
734 : }
735 : }
736 46549218 : 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 12428 : alloc_sieve(long nbprime, long maxp)
750 : {
751 : long i;
752 : ratpoints_sieve_entry * s = (ratpoints_sieve_entry*)
753 12428 : stack_malloc(nbprime*sizeof(ratpoints_sieve_entry));
754 385268 : for (i=0; i<nbprime; i++)
755 372840 : s[i].sieve = (ratpoints_bit_array**) new_chunk(maxp);
756 12428 : return s;
757 : }
758 :
759 : /* NOTE: args->degree must be set */
760 : static void
761 12428 : find_points_init(ratpoints_args *args, long bit_primes)
762 : {
763 12428 : long need = 0, n, nbprime, maxp;
764 12428 : args->listprime = primes_interval_zv(3, 1<<bit_primes);
765 12428 : nbprime = lg(args->listprime)-1;
766 12428 : maxp = args->listprime[nbprime];
767 :
768 : /* allocate space for se_buffer */
769 12428 : args->se_buffer = alloc_sieve(nbprime, maxp);
770 12428 : args->se_next = args->se_buffer;
771 385268 : for (n = 1; n <= nbprime; n++)
772 : {
773 372840 : ulong p = args->listprime[n];
774 372840 : need += p*(p + RATPOINTS_CHUNK-1);
775 : }
776 12428 : args->ba_buffer = (ratpoints_bit_array*)
777 12428 : stack_malloc_align(need*RBA_SIZE,RBA_SIZE);
778 12428 : args->ba_next = args->ba_buffer;
779 :
780 : /* allocate space for int_buffer */
781 12428 : args->int_buffer = (int *) stack_malloc(nbprime*(maxp+1)*sizeof(int));
782 12428 : args->int_next = args->int_buffer;
783 :
784 12428 : args->forb_ba = (forbidden_entry*)
785 12428 : stack_malloc((nbprime + 1)*sizeof(forbidden_entry));
786 12428 : args->forbidden = new_chunk(nbprime + 1);
787 12428 : gen_sieve(args);
788 12428 : 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 51429838 : rpjacobi(long b, GEN lcf)
795 : {
796 : ulong f;
797 51429838 : b >>= vals(b);
798 51429838 : f = umodiu(lcf, b);
799 51429838 : 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 1351 : setup_us1(ratpoints_args *args, GEN w0)
809 : {
810 1351 : GEN F = Z_issmooth_fact(w0, 1000), P, E, S, D;
811 : long i, l;
812 :
813 1351 : if (!F) return;
814 1351 : P = gel(F,1); l = lg(P);
815 1351 : E = gel(F,2);
816 1351 : D = cgetg(1+(1<<(l-1)), t_VECSMALL);
817 : /* factorization is complete, set up array of squarefree divisors */
818 1351 : D[1] = 1;
819 2828 : for (i = 1; i < l; i++)
820 : { /* multiply all divisors known so far by next prime */
821 1477 : long k, n = 1<<(i-1);
822 3080 : for (k=0; k<n; k++) uel(D,1+n+k) = uel(D,1+k) * P[i];
823 : }
824 1351 : S = cgetg(l, t_VECSMALL);
825 : /* set slopes in den_info */
826 2828 : 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 1477 : GEN c = args->cof;
829 1477 : long p = P[i], v = E[i];
830 1477 : long k, n = 1, d = degpol(c);
831 :
832 7042 : for (k = d - 1; k >= 0; k--)
833 : {
834 5565 : long t = 1 + v - Z_lval(gel(c,k+2), p);
835 5565 : long m = CEIL(t, d - k);
836 :
837 5565 : if (m > n) n = m;
838 : }
839 1477 : S[i] = n;
840 : }
841 1351 : args->divisors = D;
842 1351 : args->flags |= RATPOINTS_USE_SQUARES1;
843 1351 : args->den_info = mkvec3(P, E, S);
844 : }
845 :
846 : /************************************************************************
847 : * Consider 2-adic information *
848 : ************************************************************************/
849 :
850 : static bit_selection
851 12428 : get_2adic_info(ratpoints_args *args, ulong *den_bits,
852 : ratpoints_bit_array *num_bits)
853 : {
854 12428 : GEN c = args->cof;
855 12428 : long degree = degpol(c);
856 : int is_f_square16[24];
857 12428 : long *cmp = new_chunk(degree+1);
858 12428 : long npe = 0, npo = 0;
859 : bit_selection result;
860 : long n, a, b;
861 :
862 : /* compute coefficients mod 16 */
863 87175 : for (n = 0; n <= degree; n++) cmp[n] = Mod16(gel(c,n+2));
864 211276 : for (a = 0 ; a < 16; a++)
865 : {
866 198848 : ulong s = cmp[degree];
867 : long n;
868 1195952 : for (n = degree - 1 ; n >= 0 ; n--) s = s*a + cmp[n];
869 198848 : s &= 0xf;
870 198848 : 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 12428 : long np1 = 0, np2 = 0, np3 = 0, np4 = 0;
880 :
881 12428 : if (odd(degree))
882 : {
883 1421 : long a, cf = 4*cmp[degree-1];
884 :
885 1421 : if (degree >= 2) cf += 8*cmp[degree-2];
886 7105 : 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 5684 : long k = 2*a+1;
890 5684 : long s = (2*k*cmp[degree] + cf) & 0xf;
891 5684 : if ((is_f_square16[16+a] = squares16[s])) np1++;
892 : }
893 1421 : if ((is_f_square16[20] = squares16[(4*cmp[degree]) & 0xf])) np2++;
894 1421 : if ((is_f_square16[21] = squares16[(12*cmp[degree]) & 0xf])) np2++;
895 1421 : if ((is_f_square16[22] = squares16[(8*cmp[degree]) & 0xf])) np3++;
896 1421 : is_f_square16[23] = 1; np4++;
897 : }
898 : else
899 : {
900 11007 : long a, cf = (degree >= 2) ? 4*cmp[degree-2] : 0;
901 :
902 11007 : if (degree >= 3) cf += 8*cmp[degree-3];
903 55035 : 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 44028 : long k = 2*a+1;
907 44028 : long s = ((cmp[degree]*k + 2*cmp[degree-1])*k + cf) & 0xf;
908 44028 : if ((is_f_square16[16+a] = squares16[s])) np1++;
909 : }
910 11007 : if ((is_f_square16[20] = squares16[(cmp[degree]+4*cmp[degree-1]) & 0xf]))
911 4715 : np2++;
912 11007 : if ((is_f_square16[21] = squares16[(cmp[degree]+12*cmp[degree-1]) & 0xf]))
913 4687 : np2++;
914 11007 : if ((is_f_square16[22] = squares16[(cmp[degree]+8*cmp[degree-1]) & 0xf]))
915 4561 : np3++;
916 11007 : if ((is_f_square16[23] = squares16[cmp[degree]])) np4++;
917 : }
918 :
919 : /* set den_bits */
920 12428 : { ulong db = 0;
921 : long i;
922 :
923 12428 : if (npe + npo > 0) db |= 0xaaaaUL; /* odd denominators */
924 12428 : if (np1 > 0) db |= 0x4444UL; /* v_2(den) = 1 */
925 12428 : if (np2 > 0) db |= 0x1010UL; /* v_2(den) = 2 */
926 12428 : if (np3 > 0) db |= 0x0100UL; /* v_2(den) = 3 */
927 12428 : if (np4 > 0) db |= 0x0001UL; /* v_2(den) >= 4 */
928 12428 : 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 35012 : for (i = 16; i < BITS_IN_LONG; i <<= 1) db |= db << i;
934 12253 : *den_bits = db;
935 : }
936 12253 : result = (npe == 0) ? (npo == 0) ? num_none : num_odd
937 12253 : : (npo == 0) ? num_even : num_all;
938 : }
939 :
940 : /* set up num_bits[16] */
941 :
942 : /* odd denominators */
943 12253 : switch(result)
944 : {
945 7939 : case num_all:
946 71451 : for (b = 1; b < 16; b += 2)
947 : {
948 63512 : ulong work = 0, bit = 1;
949 63512 : long i, invb = b; /* inverse of b mod 16 */
950 63512 : if (b & 2) invb ^= 8;
951 63512 : if (b & 4) invb ^= 8;
952 1079704 : for (i = 0; i < 16; i++)
953 : {
954 1016192 : if (is_f_square16[(invb*i) & 0xf]) work |= bit;
955 1016192 : bit <<= 1;
956 : }
957 : /* now repeat the 16 bits */
958 181504 : for (i = 16; i < BITS_IN_LONG; i <<= 1) work |= work << i;
959 63512 : num_bits[b] = RBA(work);
960 : }
961 7939 : break;
962 :
963 1842 : case num_odd:
964 16578 : for (b = 1; b < 16; b += 2)
965 : {
966 14736 : ulong work = 0, bit = 1;
967 14736 : long i, invb = b; /* inverse of b mod 16 */
968 14736 : if (b & 2) invb ^= 8;
969 14736 : if (b & 4) invb ^= 8;
970 132624 : for (i = 1; i < 16; i += 2)
971 : {
972 117888 : if (is_f_square16[(invb*i) & 0xf]) work |= bit;
973 117888 : bit <<= 1;
974 : }
975 : /* now repeat the 8 bits */
976 56832 : for (i = 8; i < BITS_IN_LONG; i <<= 1) { work |= work << i; }
977 14736 : num_bits[b] = RBA(work);
978 : }
979 1842 : 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 413 : case num_none:
1000 3717 : for (b = 1; b < 16; b += 2) num_bits[b] = RBA(0UL);
1001 413 : break;
1002 : }
1003 :
1004 : /* v_2(den) = 1 : only odd numerators */
1005 61265 : for (b = 1; b < 8; b += 2)
1006 : {
1007 49012 : ulong work = 0, bit = 1;
1008 : long i;
1009 441108 : for (i = 1; i < 16; i += 2)
1010 : {
1011 392096 : if (is_f_square16[16 + (((b*i)>>1) & 0x3)]) work |= bit;
1012 392096 : bit <<= 1;
1013 : }
1014 : /* now repeat the 8 bits */
1015 189060 : for (i = 8; i < BITS_IN_LONG; i <<= 1) work |= work << i;
1016 49012 : num_bits[2*b] = RBA(work);
1017 : }
1018 :
1019 : /* v_2(den) = 2 : only odd numerators */
1020 36759 : for (b = 1; b < 4; b += 2)
1021 : {
1022 24506 : ulong work = 0, bit = 1;
1023 : long i;
1024 122530 : for (i = 1; i < 8; i += 2)
1025 : {
1026 98024 : if (is_f_square16[20 + (((b*i)>>1) & 0x1)]) work |= bit;
1027 98024 : bit <<= 1;
1028 : }
1029 : /* now repeat the 4 bits */
1030 119036 : for (i = 4; i < BITS_IN_LONG; i <<= 1) work |= work << i;
1031 24506 : num_bits[4*b] = RBA(work);
1032 : }
1033 :
1034 : /* v_2(den) = 3, >= 4 : only odd numerators */
1035 12253 : num_bits[8] = (is_f_square16[22]) ? RBA(~(0UL)) : RBA(0UL);
1036 12253 : num_bits[0] = (is_f_square16[23]) ? RBA(~(0UL)) : RBA(0UL);
1037 :
1038 12253 : 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 1171606 : compare_entries(const void *a, const void *b)
1048 : {
1049 1171606 : double diff = ((entry *)a)->r - ((entry *)b)->r;
1050 1171606 : return (diff > 0) ? 1 : (diff < 0) ? -1 : 0;
1051 : }
1052 :
1053 : /************************************************************************
1054 : * Collect the sieving information *
1055 : ************************************************************************/
1056 :
1057 : static long
1058 12428 : sieving_info(ratpoints_args *args,
1059 : ratpoints_sieve_entry **sieve_list)
1060 : {
1061 12428 : GEN c = args->cof;
1062 12428 : GEN inverses = args->inverses, squares;
1063 12428 : GEN offsets = args->offsets;
1064 12428 : ulong ** sieves0 = args->sieves0;
1065 12428 : long degree = degpol(c);
1066 12428 : long fba = 0, fdc = 0;
1067 12428 : long pn, pnp = 0;
1068 : long n;
1069 12428 : long nbprime = lg(args->listprime)-1;
1070 : long * coeffs_mod_p;
1071 12428 : 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 12428 : forbidden_entry *forb_ba = args->forb_ba;
1076 12428 : long *forbidden = args->forbidden;
1077 12428 : ulong bound = (1UL)<<(BITS_IN_LONG - args->bit_primes);
1078 12428 : pari_sp av = avma;
1079 12428 : squares = gen_squares(args->listprime);
1080 12428 : coeffs_mod_p = (long *) new_chunk(degree+1);
1081 : /* initialize sieve in se_buffer */
1082 381537 : for (pn = 1; pn <= args->num_primes; pn++)
1083 : {
1084 369235 : ulong a, p = args->listprime[pn], np;
1085 : long n;
1086 369235 : int *is_f_square = args->int_next;
1087 :
1088 369235 : args->int_next += p + 1; /* need space for (p+1) int's */
1089 :
1090 : /* compute coefficients mod p */
1091 2586410 : for (n = 0; n <= degree; n++) coeffs_mod_p[n] = umodiu(pel(c,n), p);
1092 :
1093 369235 : np = umael(squares,pn,coeffs_mod_p[0]+1);
1094 369235 : is_f_square[0] = np;
1095 21135487 : for (a = 1 ; a < p; a++)
1096 : {
1097 20766252 : ulong s = coeffs_mod_p[degree];
1098 20766252 : if ((degree+1)*args->bit_primes <= BITS_IN_LONG)
1099 : {
1100 107791784 : for (n = degree - 1 ; n >= 0 ; n--) s = s*a + coeffs_mod_p[n];
1101 : /* here, s < p^(degree+1) <= max. long */
1102 18028248 : s %= p;
1103 : }
1104 : else
1105 : {
1106 16895668 : for (n = degree - 1 ; n >= 0 ; n--)
1107 : {
1108 14157664 : s = s*a + coeffs_mod_p[n];
1109 14157664 : if (s+1 >= bound) s %= p;
1110 : }
1111 2738004 : s %= p;
1112 : }
1113 20766252 : if ((is_f_square[a] = mael(squares,pn,s+1))) np++;
1114 : }
1115 369235 : 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 369235 : if (np == 0 && !is_f_square[p]) return gc_long(av,p);
1119 :
1120 : /* Fill arrays with info for p */
1121 369109 : if (np < p)
1122 : { /* only when there is some information */
1123 : ulong i;
1124 336357 : ratpoints_sieve_entry *se = args->se_next;
1125 862442 : double r = is_f_square[p] ? ((double)(np*(p-1) + p))/((double)(p*p))
1126 336357 : : (double)np/(double)p;
1127 336357 : prec[pnp].r = r;
1128 336357 : args->se_next ++;
1129 336357 : se->p = p;
1130 336357 : se->is_f_square = is_f_square;
1131 336357 : se->inverses = gel(inverses,pn);
1132 336357 : se->offset = offsets[pn];
1133 336357 : se->sieve[0] = (ratpoints_bit_array *)sieves0[pn];
1134 20237033 : for (i = 1; i < p; i++) se->sieve[i] = NULL;
1135 336357 : prec[pnp].ssp = se;
1136 336357 : pnp++;
1137 : }
1138 :
1139 369109 : if ((args->flags & RATPOINTS_CHECK_DENOM)
1140 321859 : && fba + fdc < args->max_forbidden
1141 321859 : && !is_f_square[p])
1142 : { /* record forbidden divisors of the denominator */
1143 148043 : 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 148043 : forb_ba[fba].p = p;
1181 148043 : forb_ba[fba].start = sieves0[pn];
1182 148043 : forb_ba[fba].end = sieves0[pn]+p;
1183 148043 : forb_ba[fba].curr = forb_ba[fba].start;
1184 148043 : fba++;
1185 : }
1186 : }
1187 : } /* end for pn */
1188 :
1189 : /* update sp2 and sp1 if necessary */
1190 12302 : if (args->sp2 > pnp) args->sp2 = pnp;
1191 12302 : if (args->sp1 > args->sp2) args->sp1 = args->sp2;
1192 :
1193 : /* sort the array to get at the best primes */
1194 12302 : qsort(prec, pnp, sizeof(entry), compare_entries);
1195 :
1196 : /* put the sorted entries into sieve_list */
1197 245578 : for (n = 0; n < args->sp2; n++) sieve_list[n] = prec[n].ssp;
1198 :
1199 : /* terminate array of forbidden divisors */
1200 12302 : if (args->flags & RATPOINTS_CHECK_DENOM)
1201 : {
1202 : long n;
1203 :
1204 10727 : for (n = args->num_primes+1;
1205 10727 : 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 10727 : forb_ba[fba].p = 0; /* terminating zero */
1220 10727 : forbidden[fdc] = 0; /* terminating zero */
1221 10727 : args->max_forbidden = fba + fdc; /* note actual number */
1222 : }
1223 :
1224 12302 : if (fba + fdc == 0) args->flags &= ~RATPOINTS_CHECK_DENOM;
1225 12302 : return gc_long(av,0);
1226 : }
1227 :
1228 : /**************************************************************************
1229 : * The sieving procedure itself *
1230 : **************************************************************************/
1231 : static void
1232 31921627 : 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 31921627 : pari_sp av = avma;
1238 31921627 : int do_setup = 1;
1239 31921627 : long k, height = args->height, nb;
1240 31921627 : sieve_spec *ssp = (sieve_spec *) stack_malloc(args->sp2*sizeof(sieve_spec));
1241 :
1242 31921627 : if (odd(b) == 0) which_bits = num_odd; /* even denominator */
1243 :
1244 : /* Note that b is new */
1245 31921627 : args->bc = NULL;
1246 31921627 : nb = 0;
1247 :
1248 75778222 : for (k = 0; k < args->num_inter; k++)
1249 : {
1250 47850070 : ratpoints_interval inter = args->domain[k];
1251 : long low, high;
1252 :
1253 : /* Determine relevant interval [low, high] of numerators. */
1254 47850070 : if (b*inter.low <= -height)
1255 24074646 : low = -height;
1256 : else
1257 : {
1258 23775424 : if (b*inter.low > height) break;
1259 19783706 : low = ceil(b*inter.low);
1260 : }
1261 43858352 : if (b*inter.up >= height)
1262 19682387 : high = height;
1263 : else
1264 : {
1265 24175965 : if (b*inter.up < -height) continue;
1266 20259631 : high = floor(b*inter.up);
1267 : }
1268 :
1269 39942018 : if (do_setup)
1270 : { /* set up the sieve information */
1271 : long n;
1272 :
1273 29574009 : do_setup = 0; /* only do it once for every b */
1274 591086500 : for (n = 0; n < args->sp2; n++)
1275 : {
1276 561512491 : ratpoints_sieve_entry *se = sieve_list[n];
1277 561512491 : long p = se->p;
1278 561512491 : long bp = bp_list[n];
1279 : ratpoints_bit_array *sptr;
1280 :
1281 561512491 : if (which_bits != num_all) /* divide by 2 mod p */
1282 322763170 : bp = odd(bp) ? (bp+p) >> 1 : bp >> 1;
1283 561512491 : sptr = se->sieve[bp];
1284 :
1285 561512491 : ssp[n].p = p;
1286 561512491 : ssp[n].offset = (which_bits == num_odd) ? se->offset : 0;
1287 :
1288 : /* copy if already initialized, else initialize */
1289 574189119 : ssp[n].ptr = sptr ? sptr : (p<BITS_IN_LONG?sieve_init1(p, se, bp, args)
1290 12676628 : :sieve_init2(p, se, bp, args));
1291 561512491 : ssp[n].start = ssp[n].ptr;
1292 561512491 : ssp[n].end = ssp[n].ptr + p;
1293 :
1294 : }
1295 : }
1296 :
1297 39942018 : switch(which_bits)
1298 : {
1299 17083285 : case num_all: break;
1300 0 : case num_none: break;
1301 18217195 : 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 39942018 : high++;
1307 :
1308 39942018 : if (low < high)
1309 : {
1310 39939892 : long w_low, w_high, w_low0, w_high0, range = args->array_size;
1311 :
1312 : /* Now the range of longwords (= bit_arrays) */
1313 39939892 : w_low = low >> RBA_SHIFT;
1314 39939892 : w_high = (high + (long)(RBA_LENGTH-1)) >> RBA_SHIFT;
1315 39939892 : w_low0 = w_low;
1316 39939892 : w_high0 = w_low0 + range;
1317 86489110 : for ( ; w_low0 < w_high; w_low0 = w_high0, w_high0 += range)
1318 : {
1319 : long i;
1320 46550975 : if (w_high0 > w_high) { w_high0 = w_high; range = w_high0 - w_low0; }
1321 : /* initialise the bits */
1322 5205242181 : for (i = range; i; i--) survivors[i-1] = bits16;
1323 : /* boundary words */
1324 46550975 : if (w_low0 == w_low)
1325 39939892 : MASKL(survivors,low - RBA_LENGTH * w_low)
1326 46550975 : if (w_high0 == w_high)
1327 39939748 : MASKU(&survivors[range-1], RBA_LENGTH * w_high - high)
1328 :
1329 : #if (RATPOINTS_CHUNK > 1)
1330 248813166 : while(range%RATPOINTS_CHUNK != 0)
1331 213157776 : { survivors[range] = RBA(0); range++; w_high0++; }
1332 : #endif
1333 46550975 : nb += _ratpoints_sift0(b, w_low0, w_high0, args, which_bits,
1334 : survivors, &ssp[0], quit, process, info);
1335 46550975 : if (*quit) return;
1336 : }
1337 : }
1338 : }
1339 31919870 : 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 12428 : find_points_work(ratpoints_args *args,
1347 : int process(long, long, GEN, void*, int*), void *info)
1348 : {
1349 12428 : int quit = 0;
1350 12428 : GEN c = args->cof;
1351 12428 : long degree = degpol(c);
1352 12428 : long nbprime = lg(args->listprime)-1;
1353 12428 : long height = args->height;
1354 :
1355 12428 : int point_at_infty = 0; /* indicates if there are points at infinity */
1356 12428 : int lcfsq = Z_issquare(pel(c,degree));
1357 :
1358 12428 : forbidden_entry *forb_ba = args->forb_ba;
1359 12428 : 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 12428 : stack_malloc(nbprime*sizeof(ratpoints_sieve_entry*));
1367 12428 : bit_selection which_bits = num_all;
1368 : ulong den_bits;
1369 : ratpoints_bit_array num_bits[16];
1370 :
1371 12428 : args->flags &= RATPOINTS_FLAGS_INPUT_MASK;
1372 12428 : args->flags |= RATPOINTS_CHECK_DENOM;
1373 :
1374 : /* initialize memory management */
1375 12428 : args->se_next = args->se_buffer;
1376 12428 : args->ba_next = args->ba_buffer;
1377 12428 : args->int_next = args->int_buffer;
1378 :
1379 : /* Some sanity checks */
1380 12428 : args->num_inter = 0;
1381 :
1382 12428 : if (args->num_primes > nbprime) args->num_primes = nbprime;
1383 12428 : if (args->sp2 > args->num_primes) args->sp2 = args->num_primes;
1384 12428 : if (args->sp1 > args->sp2) args->sp1 = args->sp2;
1385 :
1386 12428 : if (args->b_low < 1) args->b_low = 1;
1387 12428 : if (args->b_high < 1) args->b_high = height;
1388 12428 : if (args->max_forbidden < 0)
1389 0 : args->max_forbidden = RATPOINTS_DEFAULT_MAX_FORBIDDEN;
1390 12428 : if (args->max_forbidden > nbprime)
1391 0 : args->max_forbidden = nbprime;
1392 12428 : if (args->array_size <= 0) args->array_size = RATPOINTS_ARRAY_SIZE;
1393 : {
1394 12428 : long s = 2*maxss(1,CEIL(height, BITS_IN_LONG));
1395 12428 : if (args->array_size > s) args->array_size = s;
1396 : }
1397 : /* make sure that array size is a multiple of RATPOINTS_CHUNK */
1398 12428 : 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 12428 : 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 12428 : if (!(args->flags & RATPOINTS_NO_REVERSE))
1412 : {
1413 12393 : if (!odd(degree) && degree>0)
1414 : {
1415 11210 : 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 10986 : 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 1183 : 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 12428 : if (args->num_inter == 0)
1444 : { /* default interval (effectively ]-oo,oo[) if none is given */
1445 12428 : args->domain = (ratpoints_interval*) stack_malloc(2*degree*sizeof(ratpoints_interval));
1446 12428 : args->domain[0].low = -height; args->domain[0].up = height;
1447 12428 : args->num_inter = 1;
1448 : }
1449 :
1450 12428 : ratpoints_compute_sturm(args);
1451 :
1452 : /* Point(s) at infinity? */
1453 12428 : if (odd(degree) || lcfsq)
1454 : {
1455 1575 : args->flags &= ~RATPOINTS_CHECK_DENOM;
1456 1575 : point_at_infty = 1;
1457 : }
1458 :
1459 : /* Can use only squares as denoms if degree is odd and poly is +-monic */
1460 12428 : if (odd(degree))
1461 : {
1462 1421 : GEN w1 = pel(c,degree);
1463 1421 : if (is_pm1(w1))
1464 70 : args->flags |= RATPOINTS_USE_SQUARES;
1465 : else /* set up information on divisors of leading coefficient */
1466 1351 : setup_us1(args, absi_shallow(w1));
1467 : }
1468 :
1469 : /* deal with f mod powers of 2 */
1470 12428 : 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 12428 : if (sieving_info(args, sieve_list)) return;
1485 :
1486 : /* deal with point(s) at infinity */
1487 12302 : if (point_at_infty)
1488 : {
1489 1575 : long a = 1, b = 0;
1490 :
1491 1575 : if (args->flags & RATPOINTS_REVERSED) { a = 0; b = 1; }
1492 1575 : if (odd(degree))
1493 1421 : (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 1575 : if (quit) return;
1501 : }
1502 : /* now do the sieving */
1503 : {
1504 : ratpoints_bit_array *survivors = (ratpoints_bit_array *)
1505 12302 : stack_malloc_align((args->array_size)*RBA_SIZE, RBA_SIZE);
1506 12302 : long *bp_list = (long *) new_chunk(args->sp2);
1507 12302 : if (args->flags & (RATPOINTS_USE_SQUARES | RATPOINTS_USE_SQUARES1))
1508 : {
1509 1421 : 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 1351 : GEN den_info = args->den_info;
1541 1351 : GEN divisors = args->divisors;
1542 1351 : long ld = lg(divisors);
1543 : long k;
1544 : long b, bb;
1545 :
1546 4291 : for (k = 1; k < ld; k++)
1547 : {
1548 2947 : long d = divisors[k];
1549 2947 : long last_b = d;
1550 : long n;
1551 2947 : if (d > args->b_high) continue;
1552 58800 : for (n = 0; n < args->sp2; n++)
1553 55860 : bp_list[n] = mod(d, sieve_list[n]->p);
1554 :
1555 279097 : for (b = 1; bb = d*b*b, bb <= args->b_high; b++)
1556 : {
1557 276164 : if (bb >= args->b_low)
1558 : {
1559 276143 : int flag = 1;
1560 276143 : ratpoints_bit_array bits = num_bits[bb & 0xf];
1561 :
1562 276143 : if (EXT0(bits))
1563 : {
1564 227801 : long i, n, l = lg(gel(den_info,1));
1565 227801 : long d = bb - last_b;
1566 :
1567 : /* fill bp_list */
1568 4556020 : for (n = 0; n < args->sp2; n++)
1569 4328219 : bp_list[n] = mod(bp_list[n] + d, sieve_list[n]->p);
1570 227801 : last_b = bb;
1571 :
1572 431893 : for(i = 1; i < l; i++)
1573 : {
1574 256088 : long v = z_lval(bb, mael(den_info,1,i));
1575 256088 : if((v >= mael(den_info,3,i))
1576 124089 : && odd(v + (mael(den_info,2,i)))) { flag = 0; break; }
1577 : }
1578 227801 : if (flag)
1579 : {
1580 175805 : sift(bb, survivors, args, which_bits, bits,
1581 : sieve_list, &bp_list[0], &quit, process, info);
1582 175805 : if (quit) break;
1583 : }
1584 : }
1585 : }
1586 : }
1587 2940 : if (quit) break;
1588 : }
1589 : }
1590 : }
1591 : else
1592 : {
1593 10881 : if (args->flags & RATPOINTS_CHECK_DENOM)
1594 : {
1595 : long *forb;
1596 : long b;
1597 10727 : long last_b = args->b_low;
1598 : ulong b_bits;
1599 : long n;
1600 214211 : for (n = 0; n < args->sp2; n++)
1601 203484 : bp_list[n] = mod(args->b_low, sieve_list[n]->p);
1602 : {
1603 10727 : forbidden_entry *fba = &forb_ba[0];
1604 10727 : long b_low = args->b_low;
1605 10727 : long w_low = (b_low-1) >> TWOPOTBITS_IN_LONG;
1606 :
1607 10727 : b_bits = den_bits;
1608 158756 : while (fba->p)
1609 : {
1610 148029 : fba->curr = fba->start + mod(w_low, fba->p);
1611 148029 : b_bits &= *(fba->curr);
1612 148029 : fba++;
1613 : }
1614 10727 : b_bits >>= (b_low-1) & LONG_MASK;
1615 : }
1616 :
1617 135533359 : for (b = args->b_low; b <= args->b_high; b++)
1618 : {
1619 135524375 : ratpoints_bit_array bits = num_bits[b & 0xf];
1620 :
1621 135524375 : if ((b & LONG_MASK) == 0)
1622 : { /* next b_bits */
1623 2411862 : forbidden_entry *fba = &forb_ba[0];
1624 :
1625 2411862 : b_bits = den_bits;
1626 37831362 : while (fba->p)
1627 : {
1628 35419500 : fba->curr++;
1629 35419500 : if (fba->curr == fba->end) fba->curr = fba->start;
1630 35419500 : b_bits &= *(fba->curr);
1631 35419500 : fba++;
1632 : }
1633 : }
1634 : else
1635 133112513 : b_bits >>= 1;
1636 :
1637 135524375 : if (odd(b_bits) && EXT0(bits))
1638 : { /* check if denominator is excluded */
1639 51429838 : for (forb = &forbidden[0] ; *forb && (b % (*forb)); forb++)
1640 0 : continue;
1641 51429838 : if (*forb == 0 && rpjacobi(b, pel(c,degree)) == 1)
1642 : {
1643 30060768 : long n, d = b - last_b;
1644 :
1645 : /* fill bp_list */
1646 600822079 : for (n = 0; n < args->sp2; n++)
1647 : {
1648 570761311 : long bp = bp_list[n] + d;
1649 570761311 : long p = sieve_list[n]->p;
1650 :
1651 641901231 : while (bp >= p) bp -= p;
1652 570761311 : bp_list[n] = bp;
1653 : }
1654 30060768 : last_b = b;
1655 :
1656 30060768 : sift(b, survivors, args, which_bits, bits,
1657 : sieve_list, &bp_list[0], &quit, process, info);
1658 30060768 : 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 86198 : vec_append_grow(GEN z, long i, GEN x)
1700 : {
1701 86198 : long n = lg(z)-1;
1702 86198 : if (i > n)
1703 : {
1704 1435 : n <<= 1;
1705 1435 : z = vec_lengthen(z,n);
1706 : }
1707 86198 : gel(z,i) = x;
1708 86198 : return z;
1709 : }
1710 :
1711 : struct points
1712 : {
1713 : GEN z;
1714 : long i, f;
1715 : };
1716 :
1717 : static int
1718 89229 : process(long a, long b, GEN y, void *info0, int *quit)
1719 : {
1720 89229 : struct points *p = (struct points *) info0;
1721 89229 : if (b==0 && (p->f&2L)) return 0;
1722 86198 : *quit = (p->f&1);
1723 86198 : p->z = vec_append_grow(p->z, ++p->i, mkvec3(stoi(a),y,stoi(b)));
1724 86198 : return 1;
1725 : }
1726 :
1727 : static int
1728 12435 : args_h(ratpoints_args *args, GEN D)
1729 : {
1730 12435 : long H, h, l = 1;
1731 : GEN L;
1732 12435 : switch(typ(D))
1733 : {
1734 12393 : case t_INT: if (signe(D) <= 0) return 0;
1735 12393 : 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 12428 : args->height = H;
1748 12428 : args->b_low = l;
1749 12428 : args->b_high = h; return 1;
1750 : }
1751 :
1752 : static GEN
1753 12435 : ZX_hyperellratpoints(GEN P, GEN h, long flag)
1754 : {
1755 12435 : pari_sp av = avma;
1756 : ratpoints_args args;
1757 : struct points data;
1758 12435 : ulong flags = 0;
1759 :
1760 12435 : if (!args_h(&args, h))
1761 : {
1762 7 : pari_err_TYPE("hyperellratpoints", h);
1763 : return NULL;/*LCOV_EXCL_LINE*/
1764 : }
1765 12428 : find_points_init(&args, RATPOINTS_DEFAULT_BIT_PRIMES);
1766 :
1767 12428 : args.cof = shallowcopy(P);
1768 12428 : args.num_inter = 0;
1769 12428 : args.sp1 = RATPOINTS_DEFAULT_SP1;
1770 12428 : args.sp2 = RATPOINTS_DEFAULT_SP2;
1771 12428 : args.array_size = RATPOINTS_ARRAY_SIZE;
1772 12428 : args.num_primes = RATPOINTS_DEFAULT_NUM_PRIMES;
1773 12428 : args.bit_primes = RATPOINTS_DEFAULT_BIT_PRIMES;
1774 12428 : args.max_forbidden = RATPOINTS_DEFAULT_MAX_FORBIDDEN;
1775 12428 : args.flags = flags;
1776 12428 : data.z = cgetg(17,t_VEC);
1777 12428 : data.i = 0;
1778 12428 : data.f = flag;
1779 12428 : find_points_work(&args, process, (void *)&data);
1780 :
1781 12428 : setlg(data.z, data.i+1);
1782 12428 : 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 12435 : QX_hyperellratpoints(GEN P, GEN h, long flag, GEN *den)
1789 : {
1790 12435 : GEN Q = Q_remove_denom(P, den);
1791 12435 : if (*den) Q = ZX_Z_mul(Q, *den);
1792 12435 : 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 11539 : hyperell_check(GEN PQ, GEN *P, GEN *Q)
1811 : {
1812 11539 : *P = *Q = NULL;
1813 11539 : if (typ(PQ) == t_POL)
1814 : {
1815 11504 : if (!RgX_is_QX(PQ)) return 0;
1816 11504 : *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 11539 : return 1;
1827 : }
1828 :
1829 : GEN
1830 11539 : hyperellratpoints(GEN PQ, GEN h, long flag)
1831 : {
1832 11539 : pari_sp av = avma;
1833 : GEN P, Q, H, L, den, denQ;
1834 : long i, l, dy, dQ;
1835 :
1836 11539 : if (flag<0 || flag>1) pari_err_FLAG("ellratpoints");
1837 11539 : if (!hyperell_check(PQ, &P, &Q)) pari_err_TYPE("hyperellratpoints",PQ);
1838 11539 : if (!Q)
1839 : {
1840 11518 : L = QX_hyperellratpoints(P, h, flag|2L, &den);
1841 11518 : dy = (degpol(P)+1)>>1;
1842 11518 : l = lg(L);
1843 25434 : 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 11518 : 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 896 : ellratpoints(GEN E, GEN h, long flag)
1871 : {
1872 896 : pari_sp av = avma;
1873 : GEN L, a1, a3, den;
1874 : long i, l;
1875 896 : checkell_Q(E);
1876 896 : if (flag<0 || flag>1) pari_err_FLAG("ellratpoints");
1877 896 : if (!RgV_is_QV(vecslice(E,1,5))) pari_err_TYPE("ellratpoints",E);
1878 896 : a1 = ell_get_a1(E);
1879 896 : a3 = ell_get_a3(E);
1880 896 : L = QX_hyperellratpoints(ec_bmodel(E,0), h, flag|2L, &den);
1881 889 : l = lg(L);
1882 73003 : for (i = 1; i < l; i++)
1883 : {
1884 72114 : GEN P, Li = gel(L,i), x = gel(Li,1), y = gel(Li,2), z = gel(Li,3);
1885 72114 : if (!signe(z))
1886 0 : P = ellinf();
1887 : else
1888 : {
1889 72114 : GEN z2 = sqri(z);
1890 72114 : if (den) y = gdiv(y, den);
1891 72114 : y = gsub(y, gadd(gmul(a1, mulii(x,z)), gmul(a3,z2)));
1892 72114 : P = mkvec2(gdiv(x,z), gdiv(y,shifti(z2,1)));
1893 : }
1894 72114 : gel(L,i) = P;
1895 : }
1896 889 : return gerepilecopy(av, L);
1897 : }
|
