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