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