Line data Source code
1 : /* Copyright (C) 2014 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. It is distributed in the hope that it will be useful, but WITHOUT
8 : ANY WARRANTY WHATSOEVER.
9 :
10 : Check the License for details. You should have received a copy of it, along
11 : with the package; see the file 'COPYING'. If not, write to the Free Software
12 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
13 :
14 : #include "pari.h"
15 : #include "paripriv.h"
16 :
17 : #define dbg_printf(lvl) if (DEBUGLEVEL >= (lvl) + 3) err_printf
18 :
19 : /**
20 : * START Code from AVSs "class_inv.h"
21 : */
22 :
23 : /* actually just returns the square-free part of the level, which is
24 : * all we care about */
25 : long
26 30814 : modinv_level(long inv)
27 : {
28 30814 : switch (inv) {
29 : case INV_J:
30 22943 : return 1;
31 : case INV_G2:
32 : case INV_W3W3E2:
33 833 : return 3;
34 : case INV_F:
35 : case INV_F2:
36 : case INV_F4:
37 : case INV_F8:
38 1224 : return 6;
39 : case INV_F3:
40 70 : return 2;
41 : case INV_W3W3:
42 343 : return 6;
43 : case INV_W2W7E2:
44 : case INV_W2W7:
45 1421 : return 14;
46 : case INV_W3W5:
47 269 : return 15;
48 : case INV_W2W3E2:
49 : case INV_W2W3:
50 378 : return 6;
51 : case INV_W2W5E2:
52 : case INV_W2W5:
53 630 : return 30;
54 : case INV_W2W13:
55 336 : return 26;
56 : case INV_W3W7:
57 1809 : return 42;
58 : case INV_W5W7:
59 502 : return 35;
60 : case INV_W3W13:
61 56 : return 39;
62 : }
63 0 : pari_err_BUG("modinv_level");
64 0 : return 0;
65 : }
66 :
67 : /* Where applicable, returns N=p1*p2 (possibly p2=1) s.t. two j's
68 : * related to the same f are N-isogenous, and 0 otherwise. This is
69 : * often (but not necessarily) equal to the level. */
70 : long
71 8846144 : modinv_degree(long *P1, long *P2, long inv)
72 : {
73 : long *p1, *p2, ignored;
74 :
75 8846144 : p1 = P1 ? P1 : &ignored;
76 8846144 : p2 = P2 ? P2 : &ignored;
77 8846144 : *p1 = 1;
78 8846144 : *p2 = 1;
79 :
80 8846144 : switch (inv) {
81 : case INV_W3W5:
82 361311 : return (*p1 = 3) * (*p2 = 5);
83 : case INV_W2W3E2:
84 : case INV_W2W3:
85 533134 : return (*p1 = 2) * (*p2 = 3);
86 : case INV_W2W5E2:
87 : case INV_W2W5:
88 1903832 : return (*p1 = 2) * (*p2 = 5);
89 : case INV_W2W7E2:
90 : case INV_W2W7:
91 1125099 : return (*p1 = 2) * (*p2 = 7);
92 : case INV_W2W13:
93 1801107 : return (*p1 = 2) * (*p2 = 13);
94 : case INV_W3W7:
95 628896 : return (*p1 = 3) * (*p2 = 7);
96 : case INV_W3W3E2:
97 : case INV_W3W3:
98 942071 : return (*p1 = 3) * (*p2 = 3);
99 : case INV_W5W7:
100 415703 : return (*p1 = 5) * (*p2 = 7);
101 : case INV_W3W13:
102 242067 : return (*p1 = 3) * (*p2 = 13);
103 : }
104 :
105 892924 : return 0;
106 : }
107 :
108 : /* Certain invariants require that D not have 2 in it's conductor, but
109 : * this doesn't apply to every invariant with even level so we handle
110 : * it separately */
111 : INLINE int
112 430367 : modinv_odd_conductor(long inv)
113 : {
114 430367 : switch (inv) {
115 : case INV_F:
116 : case INV_W3W3:
117 : case INV_W3W7:
118 64586 : return 1;
119 : }
120 365781 : return 0;
121 : }
122 :
123 : long
124 73169192 : modinv_height_factor(long inv)
125 : {
126 73169192 : switch (inv) {
127 : case INV_J:
128 6207910 : return 1;
129 : case INV_G2:
130 5285217 : return 3;
131 : case INV_F:
132 13805171 : return 72;
133 : case INV_F2:
134 28 : return 36;
135 : case INV_F3:
136 1491651 : return 24;
137 : case INV_F4:
138 49 : return 18;
139 : case INV_F8:
140 49 : return 9;
141 : case INV_W2W3:
142 70 : return 72;
143 : case INV_W3W3:
144 7181909 : return 36;
145 : case INV_W2W5:
146 11538289 : return 54;
147 : case INV_W2W7:
148 4080147 : return 48;
149 : case INV_W3W5:
150 1512 : return 36;
151 : case INV_W2W13:
152 11215582 : return 42;
153 : case INV_W3W7:
154 2905819 : return 32;
155 : case INV_W2W3E2:
156 2940112 : return 36;
157 : case INV_W2W5E2:
158 571368 : return 27;
159 : case INV_W2W7E2:
160 3255616 : return 24;
161 : case INV_W3W3E2:
162 49 : return 18;
163 : case INV_W5W7:
164 2688630 : return 24;
165 : case INV_W3W13:
166 14 : return 28;
167 : default:
168 0 : pari_err_BUG("modinv_height_factor");
169 : }
170 0 : return 0;
171 : }
172 :
173 : long
174 2373910 : disc_best_modinv(long D)
175 : {
176 : long ret;
177 2373910 : ret = INV_F; if(modinv_good_disc(ret, D)) return ret;
178 1906282 : ret = INV_W2W3; if(modinv_good_disc(ret, D)) return ret;
179 1906282 : ret = INV_W2W5; if(modinv_good_disc(ret, D)) return ret;
180 1540343 : ret = INV_W2W7; if(modinv_good_disc(ret, D)) return ret;
181 1417500 : ret = INV_W2W13; if(modinv_good_disc(ret, D)) return ret;
182 1043217 : ret = INV_W3W3; if(modinv_good_disc(ret, D)) return ret;
183 809571 : ret = INV_W2W3E2;if(modinv_good_disc(ret, D)) return ret;
184 719789 : ret = INV_W3W5; if(modinv_good_disc(ret, D)) return ret;
185 719719 : ret = INV_W3W7; if(modinv_good_disc(ret, D)) return ret;
186 634613 : ret = INV_W3W13; if(modinv_good_disc(ret, D)) return ret;
187 634613 : ret = INV_W2W5E2;if(modinv_good_disc(ret, D)) return ret;
188 614740 : ret = INV_F3; if(modinv_good_disc(ret, D)) return ret;
189 577066 : ret = INV_W2W7E2;if(modinv_good_disc(ret, D)) return ret;
190 468370 : ret = INV_W5W7; if(modinv_good_disc(ret, D)) return ret;
191 383341 : ret = INV_W3W3E2;if(modinv_good_disc(ret, D)) return ret;
192 383341 : ret = INV_G2; if(modinv_good_disc(ret, D)) return ret;
193 198597 : return INV_J;
194 : }
195 :
196 : INLINE long
197 35448 : modinv_sparse_factor(long inv)
198 : {
199 35448 : switch (inv) {
200 : case INV_G2:
201 : case INV_F8:
202 : case INV_W3W5:
203 : case INV_W2W5E2:
204 : case INV_W3W3E2:
205 3501 : return 3;
206 : case INV_F:
207 736 : return 24;
208 : case INV_F2:
209 : case INV_W2W3:
210 357 : return 12;
211 : case INV_F3:
212 112 : return 8;
213 : case INV_F4:
214 : case INV_W2W3E2:
215 : case INV_W2W5:
216 : case INV_W3W3:
217 1572 : return 6;
218 : case INV_W2W7:
219 1043 : return 4;
220 : case INV_W2W7E2:
221 : case INV_W2W13:
222 : case INV_W3W7:
223 2865 : return 2;
224 : }
225 25262 : return 1;
226 : }
227 :
228 : #define IQ_FILTER_1MOD3 1
229 : #define IQ_FILTER_2MOD3 2
230 : #define IQ_FILTER_1MOD4 4
231 : #define IQ_FILTER_3MOD4 8
232 :
233 : INLINE long
234 11914 : modinv_pfilter(long inv)
235 : {
236 11914 : switch (inv) {
237 : case INV_G2:
238 : case INV_W3W3:
239 : case INV_W3W3E2:
240 : case INV_W3W5:
241 : case INV_W2W5:
242 : case INV_W2W3E2:
243 : case INV_W2W5E2:
244 : case INV_W5W7:
245 : case INV_W3W13:
246 2341 : return IQ_FILTER_1MOD3; /* ensure unique cube roots */
247 : case INV_W2W7:
248 : case INV_F3:
249 529 : return IQ_FILTER_1MOD4; /* ensure at most two 4th/8th roots */
250 : case INV_F:
251 : case INV_F2:
252 : case INV_F4:
253 : case INV_F8:
254 : case INV_W2W3:
255 : /* Ensure unique cube roots and at most two 4th/8th roots */
256 1077 : return IQ_FILTER_1MOD3 | IQ_FILTER_1MOD4;
257 : }
258 7967 : return 0;
259 : }
260 :
261 : int
262 302317 : modinv_good_prime(long inv, long p)
263 : {
264 302317 : switch (inv) {
265 : case INV_G2:
266 : case INV_W2W3E2:
267 : case INV_W3W3:
268 : case INV_W3W3E2:
269 : case INV_W3W5:
270 : case INV_W2W5E2:
271 : case INV_W2W5:
272 27858 : return (p % 3) == 2;
273 : case INV_W2W7:
274 : case INV_F3:
275 9512 : return (p & 3) != 1;
276 : case INV_F2:
277 : case INV_F4:
278 : case INV_F8:
279 : case INV_F:
280 : case INV_W2W3:
281 20036 : return ((p % 3) == 2) && (p & 3) != 1;
282 : }
283 244911 : return 1;
284 : }
285 :
286 : /* Returns true if the prime p does not divide the conductor of D
287 : * (assumes p is prime). */
288 : INLINE int
289 3893743 : prime_to_conductor(long D, long p)
290 : {
291 : long b;
292 3893743 : if (p > 2)
293 2390704 : return (D % (p * p));
294 : /* if 2 divides the conductor of D then D=0 mod 16 (when D_0 is 0
295 : * mod 4) or D=4 mod 16 (when D_0 is 1 mod 4) */
296 1503039 : b = D & 0xF;
297 1503039 : return (b && b != 4);
298 :
299 : }
300 :
301 : INLINE GEN
302 3893743 : red_primeform(long D, long p)
303 : {
304 3893743 : pari_sp av = avma;
305 : GEN P;
306 :
307 3893743 : if ( ! prime_to_conductor(D, p))
308 0 : return NULL;
309 3893743 : P = primeform_u(stoi(D), p);
310 : /* Check that P is primitive */
311 3893743 : if (gequal(gel(P, 1), gel(P, 2)) && dvdis(gel(P, 3), p)) {
312 0 : avma = av;
313 0 : return NULL;
314 : }
315 3893743 : return gerepileupto(av, redimag(P));
316 : }
317 :
318 : /* Computes product of primeforms over primes appearing in the prime
319 : * factorization of n (including multiplicity) */
320 : GEN
321 166005 : qfb_nform(long D, long n)
322 : {
323 166005 : pari_sp av = avma;
324 166005 : GEN N = NULL, fact, ps, es;
325 : long i;
326 :
327 166005 : fact = factoru(n);
328 166005 : ps = gel(fact, 1);
329 166005 : es = gel(fact, 2);
330 497959 : for (i = 1; i < lg(ps); ++i) {
331 : long j, e;
332 331954 : GEN Q = red_primeform(D, ps[i]);
333 331954 : if ( ! Q) {
334 0 : avma = av;
335 0 : return NULL;
336 : }
337 331954 : e = es[i];
338 663964 : for (j = 0; j < e; ++j) {
339 332010 : if ((i-1) || j)
340 166005 : N = qficomp(Q, N);
341 : else
342 166005 : N = Q;
343 : }
344 : }
345 166005 : return gerepileupto(av, N);
346 : }
347 :
348 : INLINE int
349 2004114 : qfb_is_two_torsion(GEN x)
350 : {
351 6012342 : return gequal1(gel(x, 1)) || gequal0(gel(x, 2))
352 4007976 : || gequal(gel(x, 1), gel(x, 2)) || gequal(gel(x, 1), gel(x, 3));
353 : }
354 :
355 : /* Returns true iff the products p1*p2, p1*p2^-1, p1^-1*p2, and
356 : * p1^-1*p2^-1 are all distinct in cl(D) */
357 : INLINE int
358 261990 : qfb_distinct_prods(long D, long p1, long p2)
359 : {
360 261990 : pari_sp av = avma;
361 : GEN P1, P2;
362 : int x;
363 :
364 261990 : P1 = red_primeform(D, p1);
365 261990 : if ( ! P1) {
366 0 : avma = av;
367 0 : return 0;
368 : }
369 261990 : P1 = qfisqr(P1);
370 :
371 261990 : P2 = red_primeform(D, p2);
372 261990 : if ( ! P2) {
373 0 : avma = av;
374 0 : return 0;
375 : }
376 261990 : P2 = qfisqr(P2);
377 :
378 524442 : x = ! (gequal(gel(P1, 1), gel(P2, 1))
379 2941 : && (gequal(gel(P1, 2), gel(P2, 2))
380 1323 : || gequal(gel(P1, 2), gneg(gel(P2, 2)))));
381 261990 : avma = av;
382 261990 : return x;
383 : }
384 :
385 : INLINE int
386 6455552 : modinv_double_eta_good_disc(long D, long inv)
387 : {
388 6455552 : pari_sp av = avma;
389 : GEN P;
390 : long i1, i2, p1, p2, N;
391 :
392 : /* By Corollary 3.1 of Enge-Schertz Constructing elliptic curves
393 : * over finite fields using double eta-quotients, we need p1 != p2
394 : * to both be non-inert and prime to the conductor, and if p1=p2=p
395 : * we want p split and prime to the conductor. We exclude the case
396 : * that p1=p2 divides the conductor, even though this does yield
397 : * class invariants */
398 6455552 : N = modinv_degree(&p1, &p2, inv);
399 6455552 : if ( ! N)
400 0 : return 0;
401 6455552 : i1 = kross(D, p1);
402 6455552 : if (i1 < 0)
403 3551354 : return 0;
404 : /* Exclude ramified case for w_{p,p} */
405 2904198 : if (p1 == p2 && !i1)
406 0 : return 0;
407 2904198 : i2 = kross(D, p2);
408 2904198 : if (i2 < 0)
409 1255097 : return 0;
410 : /* this also verifies that p1 is prime to the conductor */
411 1649101 : P = red_primeform(D, p1);
412 1649101 : if ( ! P
413 1649101 : || gequal1(gel(P, 1)) /* don't allow p1 to be principal */
414 : /* if p1 is unramified, require it to have order > 2 */
415 1648499 : || (i1 && qfb_is_two_torsion(P))) {
416 1008 : avma = av;
417 1008 : return 0;
418 : }
419 1648093 : if (p1 == p2) {
420 : /* if p1=p2 we need p1*p1 to be distinct from its inverse */
421 259385 : int not_2tors = ! qfb_is_two_torsion(qfisqr(P));
422 259385 : avma = av;
423 259385 : return not_2tors;
424 : }
425 : /* this also verifies that p2 is prime to the conductor */
426 1388708 : P = red_primeform(D, p2);
427 1388708 : if ( ! P
428 1388708 : || gequal1(gel(P, 1)) /* don't allow p2 to be principal */
429 : /* if p2 is unramified, require it to have order > 2 */
430 1388631 : || (i2 && qfb_is_two_torsion(P))) {
431 742 : avma = av;
432 742 : return 0;
433 : }
434 1387966 : avma = av;
435 :
436 1387966 : if (i1 > 0 && i2 > 0) {
437 : /* if p1 and p2 are split, we also require p1*p2,
438 : * p1*p2^-1,p1^-1*p2, and p1^-1*p2^-1 to be distinct */
439 261990 : if ( ! qfb_distinct_prods(D, p1, p2))
440 2479 : return 0;
441 : }
442 1385487 : if (!i1 && !i2) {
443 : /* if both p1 and p2 are ramified, make sure their product is not
444 : * principal */
445 165683 : P = qfb_nform(D, N);
446 165683 : if (gequal1(gel(P, 1))) {
447 161 : avma = av;
448 161 : return 0;
449 : }
450 165522 : avma = av;
451 : }
452 1385326 : return 1;
453 : }
454 :
455 : /* Assumes D is a good discriminant for inv, which implies that the
456 : * level is prime to the conductor */
457 : long
458 441 : modinv_ramified(long D, long inv)
459 : {
460 : long p1, p2, N;
461 :
462 441 : N = modinv_degree(&p1, &p2, inv);
463 441 : if (N <= 1)
464 0 : return 0;
465 441 : return !(D % p1) && !(D % p2);
466 : }
467 :
468 : int
469 17707254 : modinv_good_disc(long inv, long D)
470 : {
471 17707254 : switch (inv) {
472 : case INV_J:
473 666074 : return 1;
474 : case INV_G2:
475 472864 : return !!(D % 3);
476 : case INV_F3:
477 622832 : return (-D & 7) == 7;
478 : case INV_F:
479 : case INV_F2:
480 : case INV_F4:
481 : case INV_F8:
482 2550434 : return ((-D & 7) == 7) && (D % 3);
483 : case INV_W3W5:
484 762405 : return (D % 3) && modinv_double_eta_good_disc(D, inv);
485 : case INV_W3W3E2:
486 410424 : return (D % 3) && modinv_double_eta_good_disc(D, inv);
487 : case INV_W3W3:
488 1081248 : return (D & 1) && (D % 3) && modinv_double_eta_good_disc(D, inv);
489 : case INV_W2W3E2:
490 833231 : return (D % 3) && modinv_double_eta_good_disc(D, inv);
491 : case INV_W2W3:
492 1926946 : return ((-D & 7) == 7) && (D % 3) && modinv_double_eta_good_disc(D, inv);
493 : case INV_W2W5:
494 1962611 : return ((-D % 80) != 20) && (D % 3) && modinv_double_eta_good_disc(D, inv);
495 : case INV_W2W5E2:
496 664244 : return (D % 3) && modinv_double_eta_good_disc(D, inv);
497 : case INV_W2W7E2:
498 662102 : return ((-D % 112) != 84) && modinv_double_eta_good_disc(D, inv);
499 : case INV_W2W7:
500 1626195 : return ((-D & 7) == 7) && modinv_double_eta_good_disc(D, inv);
501 : case INV_W2W13:
502 1463413 : return ((-D % 208) != 52) && modinv_double_eta_good_disc(D, inv);
503 : case INV_W3W7:
504 820155 : return (D & 1) && (-D % 21) && modinv_double_eta_good_disc(D, inv);
505 : case INV_W5W7:
506 : /* NB: This is a guess; avs doesn't have an entry */
507 537866 : return (D % 3) && modinv_double_eta_good_disc(D, inv);
508 : case INV_W3W13:
509 : /* NB: This is a guess; avs doesn't have an entry */
510 644210 : return (D & 1) && (D % 3) && modinv_double_eta_good_disc(D, inv);
511 : default:
512 0 : pari_err_BUG("modinv_good_discriminant");
513 : }
514 0 : return 0;
515 : }
516 :
517 : int
518 721 : modinv_is_Weber(long inv)
519 : {
520 721 : return inv == INV_F || inv == INV_F2 || inv == INV_F3 || inv == INV_F4
521 721 : || inv == INV_F8;
522 : }
523 :
524 : int
525 190344 : modinv_is_double_eta(long inv)
526 : {
527 190344 : switch (inv) {
528 : case INV_W2W3:
529 : case INV_W2W3E2:
530 : case INV_W2W5:
531 : case INV_W2W5E2:
532 : case INV_W2W7:
533 : case INV_W2W7E2:
534 : case INV_W2W13:
535 : case INV_W3W3:
536 : case INV_W3W3E2:
537 : case INV_W3W5:
538 : case INV_W3W7:
539 : case INV_W5W7:
540 : case INV_W3W13:
541 28208 : return 1;
542 : }
543 162136 : return 0;
544 : }
545 :
546 : /* END Code from "class_inv.h" */
547 :
548 : INLINE int
549 9976 : safe_abs_sqrt(ulong *r, ulong x, ulong p, ulong pi, ulong s2)
550 : {
551 9976 : if (krouu(x, p) == -1) {
552 4888 : if (p%4 == 1) return 0;
553 4888 : x = Fl_neg(x, p);
554 : }
555 9976 : *r = Fl_sqrt_pre_i(x, s2, p, pi);
556 9976 : return 1;
557 : }
558 :
559 : INLINE int
560 5450 : eighth_root(ulong *r, ulong x, ulong p, ulong pi, ulong s2)
561 : {
562 : ulong s;
563 5450 : if (krouu(x, p) == -1)
564 2647 : return 0;
565 2804 : s = Fl_sqrt_pre_i(x, s2, p, pi);
566 2804 : return safe_abs_sqrt(&s, s, p, pi, s2) && safe_abs_sqrt(r, s, p, pi, s2);
567 : }
568 :
569 : INLINE ulong
570 3063 : modinv_f_from_j(ulong j, ulong p, ulong pi, ulong s2, long only_residue)
571 : {
572 3063 : pari_sp av = avma;
573 : GEN pol, rts;
574 : long i;
575 3063 : ulong g2, f = ULONG_MAX;
576 :
577 : /* f^8 must be a root of X^3 - \gamma_2 X - 16 */
578 3063 : g2 = Fl_sqrtl_pre(j, 3, p, pi);
579 :
580 3063 : pol = mkvecsmall5(0UL, Fl_neg(16 % p, p), Fl_neg(g2, p), 0UL, 1UL);
581 3063 : rts = Flx_roots(pol, p);
582 5887 : for (i = 1; i < lg(rts); ++i) {
583 5887 : if (only_residue) {
584 1175 : if (krouu(rts[i], p) != -1) {
585 616 : f = rts[i];
586 616 : break;
587 : } else
588 559 : continue;
589 4712 : } else if (eighth_root(&f, rts[i], p, pi, s2))
590 2447 : break;
591 : }
592 3063 : if (i == lg(rts))
593 0 : pari_err_BUG("modinv_f_from_j");
594 3063 : avma = av;
595 3063 : return f;
596 : }
597 :
598 : INLINE ulong
599 357 : modinv_f3_from_j(ulong j, ulong p, ulong pi, ulong s2)
600 : {
601 357 : pari_sp av = avma;
602 : GEN pol, rts;
603 : long i;
604 357 : ulong f = ULONG_MAX;
605 :
606 1071 : pol = mkvecsmall5(0UL,
607 1071 : Fl_neg(4096 % p, p), Fl_sub(768 % p, j, p), Fl_neg(48 % p, p), 1UL);
608 357 : rts = Flx_roots(pol, p);
609 738 : for (i = 1; i < lg(rts); ++i) {
610 738 : if (eighth_root(&f, rts[i], p, pi, s2))
611 357 : break;
612 : }
613 357 : if (i == lg(rts))
614 0 : pari_err_BUG("modinv_f3_from_j");
615 357 : avma = av;
616 357 : return f;
617 : }
618 :
619 : /* Return the exponent e for the double-eta "invariant" w such that
620 : * w^e is a class invariant. For example w2w3^12 is a class
621 : * invariant, so double_eta_exponent(INV_W2W3) is 12 and
622 : * double_eta_exponent(INV_W2W3E2) is 6. */
623 : INLINE ulong
624 49197 : double_eta_exponent(long inv)
625 : {
626 49197 : switch (inv) {
627 : case INV_W2W3:
628 2925 : return 12;
629 : case INV_W2W3E2:
630 : case INV_W2W5:
631 : case INV_W3W3:
632 11287 : return 6;
633 : case INV_W2W7:
634 9031 : return 4;
635 : case INV_W3W5:
636 : case INV_W2W5E2:
637 : case INV_W3W3E2:
638 5313 : return 3;
639 : case INV_W2W7E2:
640 : case INV_W2W13:
641 : case INV_W3W7:
642 14295 : return 2;
643 : default:
644 6346 : return 1;
645 : }
646 : }
647 :
648 : INLINE ulong
649 84 : weber_exponent(long inv)
650 : {
651 84 : switch (inv)
652 : {
653 77 : case INV_F: return 24;
654 0 : case INV_F2: return 12;
655 7 : case INV_F3: return 8;
656 0 : case INV_F4: return 6;
657 0 : case INV_F8: return 3;
658 0 : default: return 1;
659 : }
660 : }
661 :
662 : INLINE ulong
663 26730 : double_eta_power(long inv, ulong w, ulong p, ulong pi)
664 : {
665 26730 : return Fl_powu_pre(w, double_eta_exponent(inv), p, pi);
666 : }
667 :
668 : static GEN
669 119 : double_eta_raw_to_Fp(GEN f, GEN p)
670 : {
671 119 : GEN u = FpX_red(RgV_to_RgX(gel(f,1), 0), p);
672 119 : GEN v = FpX_red(RgV_to_RgX(gel(f,2), 0), p);
673 119 : return mkvec3(u, v, gel(f,3));
674 : }
675 :
676 : /* Given a root x of polclass(D, inv) modulo N,
677 : returns a root of polclass(D, 0) modulo N
678 : by plugging x to a modular polynomial.
679 : For double-eta quotients,
680 : this is done by plugging x into the
681 : modular polynomial Phi(w, j) where
682 : w = INV_WpWq.
683 : More information on
684 : Enge, Morain 2013: Generalised Weber Functions. */
685 : GEN
686 665 : Fp_modinv_to_j(GEN x, long inv, GEN p)
687 : {
688 665 : switch(inv)
689 : {
690 : case INV_J:
691 182 : return Fp_red(x, p);
692 : case INV_G2:
693 280 : return Fp_powu(x, 3, p);
694 : case INV_F: case INV_F2: case INV_F3: case INV_F4: case INV_F8:
695 : {
696 84 : GEN xe = Fp_powu(x, weber_exponent(inv), p);
697 84 : return Fp_div(Fp_powu(subiu(xe, 16), 3, p), xe, p);
698 : }
699 : default:
700 119 : if(modinv_is_double_eta(inv))
701 : {
702 119 : GEN xe = Fp_powu(x, double_eta_exponent(inv), p);
703 119 : GEN uvk = double_eta_raw_to_Fp(double_eta_raw(inv), p);
704 119 : GEN J0 = FpX_eval(gel(uvk,1), xe, p);
705 119 : GEN J1 = FpX_eval(gel(uvk,2), xe, p);
706 119 : GEN J2 = Fp_pow(xe, gel(uvk,3), p);
707 119 : GEN phi = mkvec3(J0, J1, J2);
708 119 : return FpX_oneroot(RgX_to_FpX(RgV_to_RgX(phi,1), p),p);
709 : }
710 0 : pari_err_BUG("Fp_modinv_to_j"); return NULL; /* LCOV_EXCL_LIpE */
711 : }
712 : }
713 :
714 : /* Assuming p = 2 (mod 3) and p = 3 (mod 4): if the two 12th roots of
715 : * x (mod p) exist, set *r to one of them and return 1, otherwise
716 : * return 0 (without touching *r). */
717 : INLINE int
718 1056 : twelth_root(ulong *r, ulong x, ulong p, ulong pi, ulong s2)
719 : {
720 : register ulong t;
721 :
722 1056 : t = Fl_sqrtl_pre(x, 3, p, pi);
723 1056 : if (krouu(t, p) == -1)
724 48 : return 0;
725 1008 : t = Fl_sqrt_pre_i(t, s2, p, pi);
726 1008 : return safe_abs_sqrt(r, t, p, pi, s2);
727 : }
728 :
729 : INLINE int
730 5288 : sixth_root(ulong *r, ulong x, ulong p, ulong pi, ulong s2)
731 : {
732 : register ulong t;
733 :
734 5288 : t = Fl_sqrtl_pre(x, 3, p, pi);
735 5288 : if (krouu(t, p) == -1)
736 192 : return 0;
737 5096 : *r = Fl_sqrt_pre_i(t, s2, p, pi);
738 5096 : return 1;
739 : }
740 :
741 : INLINE int
742 3673 : fourth_root(ulong *r, ulong x, ulong p, ulong pi, ulong s2)
743 : {
744 : register ulong s;
745 3673 : if (krouu(x, p) == -1)
746 314 : return 0;
747 3360 : s = Fl_sqrt_pre_i(x, s2, p, pi);
748 3360 : return safe_abs_sqrt(r, s, p, pi, s2);
749 : }
750 :
751 : INLINE int
752 22348 : double_eta_root(long inv, ulong *r, ulong w, ulong p, ulong pi, ulong s2)
753 : {
754 22348 : switch (double_eta_exponent(inv)) {
755 1056 : case 12: return twelth_root(r, w, p, pi, s2);
756 5288 : case 6: return sixth_root(r, w, p, pi, s2);
757 3673 : case 4: return fourth_root(r, w, p, pi, s2);
758 2265 : case 3: *r = Fl_sqrtl_pre(w, 3, p, pi); return 1;
759 7303 : case 2: return krouu(w, p) != -1 && !!(*r = Fl_sqrt_pre_i(w, s2, p, pi));
760 2763 : case 1: *r = w; return 1;
761 : }
762 0 : pari_err_BUG("double_eta_root");
763 0 : return 0;
764 : }
765 :
766 :
767 : /* TODO: Organise things better so that this forward decl is unnecessary. */
768 : static GEN Flx_double_eta_xpoly(GEN f, ulong j, ulong p, ulong pi);
769 : static GEN Flx_double_eta_jpoly(GEN f, ulong x, ulong p, ulong pi);
770 :
771 : /* reduce F = double_eta_raw(inv) mod p */
772 : static GEN
773 26552 : double_eta_raw_to_Fl(GEN f, ulong p)
774 : {
775 26552 : GEN u = ZV_to_Flv(gel(f,1), p);
776 26552 : GEN v = ZV_to_Flv(gel(f,2), p);
777 26551 : return mkvec3(u, v, gel(f,3));
778 : }
779 : /* double_eta_raw(inv) mod p */
780 : static GEN
781 21470 : double_eta_Fl(long inv, ulong p)
782 : {
783 21470 : return double_eta_raw_to_Fl(double_eta_raw(inv), p);
784 : }
785 :
786 : /* Go through the roots of Psi(X,j) until one has an
787 : * double_eta_exponent(inv)-th root, and return that root. F = double_eta_Fl(inv,p) */
788 : INLINE ulong
789 5522 : modinv_double_eta_from_j(GEN F, long inv, ulong j, ulong p, ulong pi, ulong s2)
790 : {
791 5522 : pari_sp av = avma;
792 : long i;
793 5522 : ulong f = ULONG_MAX;
794 5522 : GEN a = Flx_double_eta_xpoly(F, j, p, pi);
795 5522 : a = Flx_roots(a, p);
796 6400 : for (i = 1; i < lg(a); ++i) {
797 6400 : if (double_eta_root(inv, &f, uel(a, i), p, pi, s2))
798 5522 : break;
799 : }
800 5522 : if (i == lg(a))
801 0 : pari_err_BUG("modinv_double_eta_from_j");
802 5522 : avma = av;
803 5522 : return f;
804 : }
805 :
806 : /* TODO: Check whether I can use this to refactor something */
807 : static long
808 10426 : modinv_double_eta_from_2j(
809 : ulong *r, long inv, ulong j1, ulong j2, ulong p, ulong pi, ulong s2)
810 : {
811 10426 : pari_sp av = avma;
812 10426 : GEN f, g, d, F = double_eta_Fl(inv, p);
813 10426 : if (j2 == j1)
814 0 : pari_err_BUG("modinv_double_eta_from_2j");
815 :
816 10426 : f = Flx_double_eta_xpoly(F, j1, p, pi);
817 10426 : g = Flx_double_eta_xpoly(F, j2, p, pi);
818 10426 : d = Flx_gcd(f, g, p);
819 :
820 : /* NB: Morally the next conditional should be written as follows,
821 : * but, I think because of the case when j1 or j2 may not have the
822 : * correct endomorphism ring, we need to use the less strict
823 : * conditional underneath. */
824 : #if 0
825 : if (degpol(d) != 1
826 : || (*r = Flx_oneroot(d, p)) == p
827 : || ! double_eta_root(inv, r, *r, p, pi, s2)) {
828 : pari_err_BUG("modinv_double_eta_from_2j");
829 : }
830 : #endif
831 10426 : if (degpol(d) > 2
832 10426 : || (*r = Flx_oneroot(d, p)) == p
833 10426 : || ! double_eta_root(inv, r, *r, p, pi, s2)) {
834 6 : return 0;
835 : }
836 10420 : avma = av;
837 10420 : return 1;
838 : }
839 :
840 : long
841 10452 : modfn_unambiguous_root(ulong *r, long inv, ulong j0, norm_eqn_t ne, GEN jdb)
842 : {
843 10452 : pari_sp av = avma;
844 10452 : long p1, p2, v = ne->v, p1_depth;
845 10452 : ulong j1, p = ne->p, pi = ne->pi, s2 = ne->s2;
846 : GEN phi;
847 :
848 10452 : (void) modinv_degree(&p1, &p2, inv);
849 10452 : p1_depth = u_lval(v, p1);
850 :
851 10452 : phi = polmodular_db_getp(jdb, p1, p);
852 10452 : if ( ! next_surface_nbr(&j1, phi, p1, p1_depth, j0, NULL, p, pi))
853 0 : pari_err_BUG("modfn_unambiguous_root");
854 10452 : if (p2 == p1) {
855 1306 : if ( ! next_surface_nbr(&j1, phi, p1, p1_depth, j1, &j0, p, pi))
856 0 : pari_err_BUG("modfn_unambiguous_root");
857 : } else {
858 9146 : long p2_depth = u_lval(v, p2);
859 9146 : phi = polmodular_db_getp(jdb, p2, p);
860 9146 : if ( ! next_surface_nbr(&j1, phi, p2, p2_depth, j1, NULL, p, pi))
861 0 : pari_err_BUG("modfn_unambiguous_root");
862 : }
863 10452 : avma = av;
864 10452 : return j1 != j0 && modinv_double_eta_from_2j(r, inv, j0, j1, p, pi, s2);
865 : }
866 :
867 : ulong
868 154530 : modfn_root(ulong j, norm_eqn_t ne, long inv)
869 : {
870 154530 : ulong f, p = ne->p, pi = ne->pi, s2 = ne->s2;
871 154530 : switch (inv) {
872 : case INV_J:
873 145941 : return j;
874 : case INV_G2:
875 5170 : return Fl_sqrtl_pre(j, 3, p, pi);
876 : case INV_F:
877 1698 : return modinv_f_from_j(j, p, pi, s2, 0);
878 : case INV_F2:
879 195 : f = modinv_f_from_j(j, p, pi, s2, 0);
880 196 : return Fl_sqr_pre(f, p, pi);
881 : case INV_F3:
882 357 : return modinv_f3_from_j(j, p, pi, s2);
883 : case INV_F4:
884 553 : f = modinv_f_from_j(j, p, pi, s2, 0);
885 553 : return Fl_sqr_pre(Fl_sqr_pre(f, p, pi), p, pi);
886 : case INV_F8:
887 616 : return modinv_f_from_j(j, p, pi, s2, 1);
888 : }
889 0 : if (modinv_is_double_eta(inv))
890 : {
891 0 : pari_sp av = avma;
892 0 : ulong f = modinv_double_eta_from_j(double_eta_Fl(inv,p), inv, j, p, pi, s2);
893 0 : avma = av; return f;
894 : }
895 0 : pari_err_BUG("modfn_root");
896 0 : return ULONG_MAX;
897 : }
898 :
899 : INLINE ulong
900 1942 : modinv_j_from_f(ulong x, ulong n, ulong p, ulong pi)
901 : {
902 : /* If x satisfies (X^24 - 16)^3 - X^24 * j = 0
903 : * then j = (x^24 - 16)^3 / x^24 */
904 1942 : ulong x24 = Fl_powu_pre(x, 24 / n, p, pi);
905 1942 : return Fl_div(Fl_powu_pre(Fl_sub(x24, 16 % p, p), 3, p, pi), x24, p);
906 : }
907 :
908 : /* TODO: Check whether I can use this to refactor something
909 : * F = double_eta_raw(inv) */
910 : long
911 5082 : modinv_j_from_2double_eta(
912 : GEN F, long inv, ulong *j, ulong x0, ulong x1, ulong p, ulong pi)
913 : {
914 : GEN f, g, d;
915 :
916 5082 : x0 = double_eta_power(inv, x0, p, pi);
917 5082 : x1 = double_eta_power(inv, x1, p, pi);
918 5082 : F = double_eta_raw_to_Fl(F, p);
919 5082 : f = Flx_double_eta_jpoly(F, x0, p, pi);
920 5082 : g = Flx_double_eta_jpoly(F, x1, p, pi);
921 5082 : d = Flx_gcd(f, g, p);
922 5082 : if (degpol(d) > 1)
923 0 : pari_err_BUG("modinv_j_from_2double_eta");
924 5082 : else if (degpol(d) < 1)
925 1869 : return 0;
926 3213 : if (j)
927 0 : *j = Flx_deg1_root(d, p);
928 3213 : return 1;
929 : }
930 :
931 : INLINE ulong
932 49935 : modfn_preimage(ulong x, norm_eqn_t ne, long inv)
933 : {
934 49935 : ulong p = ne->p, pi = ne->pi;
935 49935 : switch (inv) {
936 : case INV_J:
937 44291 : return x;
938 : case INV_G2:
939 3702 : return Fl_powu_pre(x, 3, p, pi);
940 : case INV_F:
941 : /* NB: We could replace these with a single call
942 : *
943 : * modinv_j_from_f(x, inv, p, pi),
944 : *
945 : * but it's probably better to avoid the dependence on the actual
946 : * value of inv and assume it could change. */
947 738 : return modinv_j_from_f(x, 1, p, pi);
948 : case INV_F2:
949 196 : return modinv_j_from_f(x, 2, p, pi);
950 : case INV_F3:
951 168 : return modinv_j_from_f(x, 3, p, pi);
952 : case INV_F4:
953 392 : return modinv_j_from_f(x, 4, p, pi);
954 : case INV_F8:
955 448 : return modinv_j_from_f(x, 8, p, pi);
956 : }
957 : /* NB: This function should never be called if modinv_double_eta(inv) is
958 : * true */
959 0 : pari_err_BUG("modfn_preimage");
960 0 : return ULONG_MAX;
961 : }
962 :
963 :
964 : /**
965 : * SECTION: class group bb_group.
966 : */
967 :
968 : INLINE GEN
969 159350 : mkqfis(long a, long b, long c)
970 : {
971 159350 : retmkqfi(stoi(a), stoi(b), stoi(c));
972 : }
973 :
974 : /**
975 : * SECTION: Fixed-length dot-product-like functions on Fl's with
976 : * precomputed inverse.
977 : */
978 :
979 : /* Computes x0y1 + y0x1 (mod p); assumes p < 2^63. */
980 : INLINE ulong
981 47193966 : Fl_addmul2(
982 : ulong x0, ulong x1, ulong y0, ulong y1,
983 : ulong p, ulong pi)
984 : {
985 47193966 : return Fl_addmulmul_pre(x0,y1,y0,x1,p,pi);
986 : }
987 :
988 :
989 : /* Computes x0y2 + x1y1 + x2y0 (mod p); assumes p < 2^62. */
990 : INLINE ulong
991 8650649 : Fl_addmul3(
992 : ulong x0, ulong x1, ulong x2, ulong y0, ulong y1, ulong y2,
993 : ulong p, ulong pi)
994 : {
995 : ulong l0, l1, h0, h1;
996 : LOCAL_OVERFLOW;
997 : LOCAL_HIREMAINDER;
998 8650649 : l0 = mulll(x0, y2); h0 = hiremainder;
999 8650649 : l1 = mulll(x1, y1); h1 = hiremainder;
1000 8650649 : l1 = addll(l0, l1); h1 = addllx(h0, h1);
1001 8650649 : l0 = mulll(x2, y0); h0 = hiremainder;
1002 8650649 : l1 = addll(l0, l1); h1 = addllx(h0, h1);
1003 8650649 : return remll_pre(h1, l1, p, pi);
1004 : }
1005 :
1006 :
1007 : /* Computes x0y3 + x1y2 + x2y1 + x3y0 (mod p); assumes p < 2^62. */
1008 : INLINE ulong
1009 3964999 : Fl_addmul4(
1010 : ulong x0, ulong x1, ulong x2, ulong x3,
1011 : ulong y0, ulong y1, ulong y2, ulong y3,
1012 : ulong p, ulong pi)
1013 : {
1014 : ulong l0, l1, h0, h1;
1015 : LOCAL_OVERFLOW;
1016 : LOCAL_HIREMAINDER;
1017 3964999 : l0 = mulll(x0, y3); h0 = hiremainder;
1018 3964999 : l1 = mulll(x1, y2); h1 = hiremainder;
1019 3964999 : l1 = addll(l0, l1); h1 = addllx(h0, h1);
1020 3964999 : l0 = mulll(x2, y1); h0 = hiremainder;
1021 3964999 : l1 = addll(l0, l1); h1 = addllx(h0, h1);
1022 3964999 : l0 = mulll(x3, y0); h0 = hiremainder;
1023 3964999 : l1 = addll(l0, l1); h1 = addllx(h0, h1);
1024 3964999 : return remll_pre(h1, l1, p, pi);
1025 : }
1026 :
1027 : /* Computes x0y4 + x1y3 + x2y2 + x3y1 + x4y0 (mod p); assumes p < 2^62. */
1028 : INLINE ulong
1029 19811790 : Fl_addmul5(
1030 : ulong x0, ulong x1, ulong x2, ulong x3, ulong x4,
1031 : ulong y0, ulong y1, ulong y2, ulong y3, ulong y4,
1032 : ulong p, ulong pi)
1033 : {
1034 : ulong l0, l1, h0, h1;
1035 : LOCAL_OVERFLOW;
1036 : LOCAL_HIREMAINDER;
1037 19811790 : l0 = mulll(x0, y4); h0 = hiremainder;
1038 19811790 : l1 = mulll(x1, y3); h1 = hiremainder;
1039 19811790 : l1 = addll(l0, l1); h1 = addllx(h0, h1);
1040 19811790 : l0 = mulll(x2, y2); h0 = hiremainder;
1041 19811790 : l1 = addll(l0, l1); h1 = addllx(h0, h1);
1042 19811790 : l0 = mulll(x3, y1); h0 = hiremainder;
1043 19811790 : l1 = addll(l0, l1); h1 = addllx(h0, h1);
1044 19811790 : l0 = mulll(x4, y0); h0 = hiremainder;
1045 19811790 : l1 = addll(l0, l1); h1 = addllx(h0, h1);
1046 19811790 : return remll_pre(h1, l1, p, pi);
1047 : }
1048 :
1049 : /*
1050 : * A polmodular database for a given class invariant consists of a
1051 : * t_VEC whose L-th entry is 0 or a GEN pointing to Phi_L. This
1052 : * function produces a pair of databases corresponding to the
1053 : * j-invariant and inv.
1054 : */
1055 : GEN
1056 15155 : polmodular_db_init(long inv)
1057 : {
1058 : GEN res, jdb, fdb;
1059 : enum { DEFAULT_MODPOL_DB_LEN = 32 };
1060 :
1061 15155 : res = cgetg_block(2 + 1, t_VEC);
1062 15155 : jdb = zerovec_block(DEFAULT_MODPOL_DB_LEN);
1063 15155 : gel(res, 1) = jdb;
1064 15155 : if (inv != INV_J) {
1065 644 : fdb = zerovec_block(DEFAULT_MODPOL_DB_LEN);
1066 644 : gel(res, 2) = fdb;
1067 : } else {
1068 14511 : gel(res, 2) = gen_0;
1069 : }
1070 15155 : return res;
1071 : }
1072 :
1073 : void
1074 19737 : polmodular_db_add_level(GEN *DB, long L, long inv)
1075 : {
1076 : long max_L;
1077 : GEN db;
1078 :
1079 19737 : if (inv == INV_J) {
1080 16686 : db = gel(*DB, 1);
1081 : } else {
1082 3051 : db = gel(*DB, 2);
1083 3051 : if ( isintzero(db))
1084 0 : pari_err_BUG("polmodular_db_add_level");
1085 : }
1086 :
1087 19737 : max_L = lg(db) - 1;
1088 19737 : if (L > max_L) {
1089 : GEN newdb;
1090 50 : long i, newlen = 2 * L;
1091 :
1092 50 : newdb = cgetg_block(newlen + 1, t_VEC);
1093 1650 : for (i = 1; i <= max_L; ++i)
1094 1600 : gel(newdb, i) = gel(db, i);
1095 3326 : for (i = max_L + 1; i <= newlen; ++i)
1096 3276 : gel(newdb, i) = gen_0;
1097 50 : killblock(db);
1098 : /* NB: Uses the fact that INV_J == 0 */
1099 50 : gel(*DB, 2 - !inv) = db = newdb;
1100 : }
1101 19737 : if ( isintzero(gel(db, L))) {
1102 5828 : pari_sp av = avma;
1103 5828 : gel(db, L) = gclone(polmodular0_ZM(L, inv, NULL, NULL, 0, DB));
1104 5828 : avma = av;
1105 : }
1106 19737 : }
1107 :
1108 :
1109 : void
1110 3581 : polmodular_db_add_levels(GEN *db, long *levels, long k, long inv)
1111 : {
1112 : long i;
1113 7914 : for (i = 0; i < k; ++i)
1114 4333 : polmodular_db_add_level(db, levels[i], inv);
1115 3581 : }
1116 :
1117 : GEN
1118 290663 : polmodular_db_for_inv(GEN db, long inv)
1119 : {
1120 290663 : if (inv == INV_J)
1121 268821 : return gel(db, 1);
1122 21842 : return gel(db, 2);
1123 : }
1124 :
1125 : /* TODO: Also cache modpoly mod p for most recent p (avoid repeated
1126 : * reductions in, for example, polclass.c:oneroot_of_classpoly(). */
1127 : GEN
1128 423039 : polmodular_db_getp(GEN db, long L, ulong p)
1129 : {
1130 423039 : GEN f = gel(db, L);
1131 423039 : if (isintzero(f))
1132 0 : pari_err_BUG("polmodular_db_getp");
1133 423043 : return ZM_to_Flm(f, p);
1134 : }
1135 :
1136 :
1137 : /**
1138 : * SECTION: Table of discriminants to use.
1139 : */
1140 : typedef struct {
1141 : long L; /* modpoly level */
1142 : long D0; /* fundamental discriminant */
1143 : long D1; /* chosen discriminant */
1144 : long L0; /* first generator norm */
1145 : long L1; /* second generator norm */
1146 : long n1; /* order of L0 in cl(D1) */
1147 : long n2; /* order of L0 in cl(D2) where D2 = L^2 D1 */
1148 : long nprimes; /* number of primes needed for D1 */
1149 : long dl1; /* m such that L0^m = L in cl(D1) */
1150 : long dl2_0; /* These two are (m, n) such that L0^m L1^n = form of norm L^2 in D2 */
1151 : long dl2_1; /* This n is always 1 or 0. */
1152 : long cost; /* cost to enumerate subgroup of cl(L^2D): subgroup size is n2 if L1=0, 2*n2 o.w. */
1153 : long bits;
1154 : ulong *primes;
1155 : ulong *traces;
1156 : long inv;
1157 : } modpoly_disc_info;
1158 :
1159 :
1160 : static void
1161 2451 : modpoly_disc_info_clear(modpoly_disc_info *dinfo)
1162 : {
1163 2451 : killblock((GEN) dinfo->primes);
1164 2451 : killblock((GEN) dinfo->traces);
1165 2451 : }
1166 :
1167 : #define MODPOLY_MAX_DCNT 64
1168 :
1169 :
1170 : /* Flag for last parameter of discriminant_with_classno_at_least.
1171 : * Warning: ignoring the sparse factor makes everything slower by
1172 : * something like (sparse factor)^3. */
1173 : #define USE_SPARSE_FACTOR 0
1174 : #define IGNORE_SPARSE_FACTOR 1
1175 :
1176 : static long
1177 : discriminant_with_classno_at_least(
1178 : modpoly_disc_info Ds[MODPOLY_MAX_DCNT], long L, long inv,
1179 : long ignore_sparse);
1180 :
1181 :
1182 : /**
1183 : * SECTION: Hard-coded evaluation functions for modular polynomials of
1184 : * small level.
1185 : */
1186 :
1187 : /*
1188 : * Based on phi2_eval_ff() in Sutherland's classpoly programme.
1189 : * Calculates Phi_2(X, j) (mod p) with 6M+7A (4 reductions, not
1190 : * counting those for Phi_2).
1191 : */
1192 : INLINE GEN
1193 22183108 : Flm_Fl_phi2_evalx(GEN phi2, ulong j, ulong p, ulong pi)
1194 : {
1195 22183108 : GEN res = cgetg(6, t_VECSMALL);
1196 : ulong j2, t1;
1197 :
1198 22174155 : res[1] = 0; /* variable name */
1199 :
1200 22174155 : j2 = Fl_sqr_pre(j, p, pi);
1201 22188860 : t1 = Fl_add(j, coeff(phi2, 3, 1), p);
1202 22176047 : t1 = Fl_addmul2(j, j2, t1, coeff(phi2, 2, 1), p, pi);
1203 22186259 : res[2] = Fl_add(t1, coeff(phi2, 1, 1), p);
1204 :
1205 22173152 : t1 = Fl_addmul2(j, j2, coeff(phi2, 3, 2), coeff(phi2, 2, 2), p, pi);
1206 22187231 : res[3] = Fl_add(t1, coeff(phi2, 2, 1), p);
1207 :
1208 22175214 : t1 = Fl_mul_pre(j, coeff(phi2, 3, 2), p, pi);
1209 22184108 : t1 = Fl_add(t1, coeff(phi2, 3, 1), p);
1210 22173283 : res[4] = Fl_sub(t1, j2, p);
1211 :
1212 22173126 : res[5] = 1;
1213 22173126 : return res;
1214 : }
1215 :
1216 :
1217 : /*
1218 : * Based on phi3_eval_ff() in Sutherland's classpoly programme.
1219 : * Calculates Phi_3(X, j) (mod p) with 13M+13A (6 reductions, not
1220 : * counting those for Phi_3).
1221 : */
1222 : INLINE GEN
1223 2885797 : Flm_Fl_phi3_evalx(GEN phi3, ulong j, ulong p, ulong pi)
1224 : {
1225 2885797 : GEN res = cgetg(7, t_VECSMALL);
1226 : ulong j2, j3, t1;
1227 :
1228 2885235 : res[1] = 0; /* variable name */
1229 :
1230 2885235 : j2 = Fl_sqr_pre(j, p, pi);
1231 2886725 : j3 = Fl_mul_pre(j, j2, p, pi);
1232 :
1233 2887327 : t1 = Fl_add(j, coeff(phi3, 4, 1), p);
1234 8658777 : res[2] = Fl_addmul3(j, j2, j3, t1,
1235 5772518 : coeff(phi3, 3, 1), coeff(phi3, 2, 1), p, pi);
1236 :
1237 5774358 : t1 = Fl_addmul3(j, j2, j3, coeff(phi3, 4, 2),
1238 5774358 : coeff(phi3, 3, 2), coeff(phi3, 2, 2), p, pi);
1239 2887210 : res[3] = Fl_add(t1, coeff(phi3, 2, 1), p);
1240 :
1241 5772640 : t1 = Fl_addmul3(j, j2, j3, coeff(phi3, 4, 3),
1242 5772640 : coeff(phi3, 3, 3), coeff(phi3, 3, 2), p, pi);
1243 2887340 : res[4] = Fl_add(t1, coeff(phi3, 3, 1), p);
1244 :
1245 2886723 : t1 = Fl_addmul2(j, j2, coeff(phi3, 4, 3), coeff(phi3, 4, 2), p, pi);
1246 2887133 : t1 = Fl_add(t1, coeff(phi3, 4, 1), p);
1247 2886316 : res[5] = Fl_sub(t1, j3, p);
1248 :
1249 2886420 : res[6] = 1;
1250 :
1251 2886420 : return res;
1252 : }
1253 :
1254 :
1255 : /*
1256 : * Based on phi5_eval_ff() in Sutherland's classpoly programme.
1257 : * Calculates Phi_5(X, j) (mod p) with 33M+31A (10 reductions, not
1258 : * counting those for Phi_5).
1259 : */
1260 : INLINE GEN
1261 3964093 : Flm_Fl_phi5_evalx(GEN phi5, ulong j, ulong p, ulong pi)
1262 : {
1263 3964093 : GEN res = cgetg(9, t_VECSMALL);
1264 : ulong j2, j3, j4, j5, t1;
1265 :
1266 3963366 : res[1] = 0; /* variable name */
1267 :
1268 3963366 : j2 = Fl_sqr_pre(j, p, pi);
1269 3964915 : j3 = Fl_mul_pre(j, j2, p, pi);
1270 3965475 : j4 = Fl_sqr_pre(j2, p, pi);
1271 3965566 : j5 = Fl_mul_pre(j, j4, p, pi);
1272 :
1273 3965764 : t1 = Fl_add(j, coeff(phi5, 6, 1), p);
1274 15857768 : t1 = Fl_addmul5(j, j2, j3, j4, j5, t1,
1275 7928884 : coeff(phi5, 5, 1), coeff(phi5, 4, 1),
1276 7928884 : coeff(phi5, 3, 1), coeff(phi5, 2, 1),
1277 : p, pi);
1278 3965280 : res[2] = Fl_add(t1, coeff(phi5, 1, 1), p);
1279 :
1280 19821100 : t1 = Fl_addmul5(j, j2, j3, j4, j5,
1281 7928440 : coeff(phi5, 6, 2), coeff(phi5, 5, 2),
1282 11892660 : coeff(phi5, 4, 2), coeff(phi5, 3, 2), coeff(phi5, 2, 2),
1283 : p, pi);
1284 3965288 : res[3] = Fl_add(t1, coeff(phi5, 2, 1), p);
1285 :
1286 19821480 : t1 = Fl_addmul5(j, j2, j3, j4, j5,
1287 7928592 : coeff(phi5, 6, 3), coeff(phi5, 5, 3),
1288 11892888 : coeff(phi5, 4, 3), coeff(phi5, 3, 3), coeff(phi5, 3, 2),
1289 : p, pi);
1290 3965455 : res[4] = Fl_add(t1, coeff(phi5, 3, 1), p);
1291 :
1292 19822565 : t1 = Fl_addmul5(j, j2, j3, j4, j5,
1293 7929026 : coeff(phi5, 6, 4), coeff(phi5, 5, 4),
1294 11893539 : coeff(phi5, 4, 4), coeff(phi5, 4, 3), coeff(phi5, 4, 2),
1295 : p, pi);
1296 3965634 : res[5] = Fl_add(t1, coeff(phi5, 4, 1), p);
1297 :
1298 19823600 : t1 = Fl_addmul5(j, j2, j3, j4, j5,
1299 7929440 : coeff(phi5, 6, 5), coeff(phi5, 5, 5),
1300 11894160 : coeff(phi5, 5, 4), coeff(phi5, 5, 3), coeff(phi5, 5, 2),
1301 : p, pi);
1302 3965703 : res[6] = Fl_add(t1, coeff(phi5, 5, 1), p);
1303 :
1304 15859904 : t1 = Fl_addmul4(j, j2, j3, j4,
1305 7929952 : coeff(phi5, 6, 5), coeff(phi5, 6, 4),
1306 7929952 : coeff(phi5, 6, 3), coeff(phi5, 6, 2),
1307 : p, pi);
1308 3965543 : t1 = Fl_add(t1, coeff(phi5, 6, 1), p);
1309 3964427 : res[7] = Fl_sub(t1, j5, p);
1310 :
1311 3964415 : res[8] = 1;
1312 :
1313 3964415 : return res;
1314 : }
1315 :
1316 : GEN
1317 36030124 : Flm_Fl_polmodular_evalx(GEN phi, long L, ulong j, ulong p, ulong pi)
1318 : {
1319 36030124 : switch (L) {
1320 22186887 : case 2: return Flm_Fl_phi2_evalx(phi, j, p, pi);
1321 2885892 : case 3: return Flm_Fl_phi3_evalx(phi, j, p, pi);
1322 3964377 : case 5: return Flm_Fl_phi5_evalx(phi, j, p, pi);
1323 : default: { /* not GC clean, but gerepileupto-safe */
1324 6992968 : GEN j_powers = Fl_powers_pre(j, L + 1, p, pi);
1325 6987745 : return Flm_Flc_mul_pre_Flx(phi, j_powers, p, pi, 0);
1326 : }
1327 : }
1328 : }
1329 :
1330 : /**
1331 : * SECTION: Velu's formula for the codmain curve in the case of small
1332 : * prime base field.
1333 : */
1334 :
1335 : INLINE ulong
1336 1431162 : Fl_mul4(ulong x, ulong p)
1337 : {
1338 1431162 : return Fl_double(Fl_double(x, p), p);
1339 : }
1340 :
1341 : INLINE ulong
1342 73048 : Fl_mul5(ulong x, ulong p)
1343 : {
1344 73048 : return Fl_add(x, Fl_mul4(x, p), p);
1345 : }
1346 :
1347 : INLINE ulong
1348 715622 : Fl_mul8(ulong x, ulong p)
1349 : {
1350 715622 : return Fl_double(Fl_mul4(x, p), p);
1351 : }
1352 :
1353 : INLINE ulong
1354 642614 : Fl_mul6(ulong x, ulong p)
1355 : {
1356 642614 : return Fl_sub(Fl_mul8(x, p), Fl_double(x, p), p);
1357 : }
1358 :
1359 : INLINE ulong
1360 73048 : Fl_mul7(ulong x, ulong p)
1361 : {
1362 73048 : return Fl_sub(Fl_mul8(x, p), x, p);
1363 : }
1364 :
1365 : /*
1366 : * Given an elliptic curve E = [a4, a6] over F_p and a non-zero point
1367 : * pt on E, return the quotient E' = E/<P> = [a4_img, a6_img].
1368 : */
1369 : static void
1370 73044 : Fle_quotient_from_kernel_generator(
1371 : ulong *a4_img, ulong *a6_img, ulong a4, ulong a6, GEN pt, ulong p, ulong pi)
1372 : {
1373 73044 : pari_sp av = avma;
1374 73044 : ulong t = 0, w = 0;
1375 : GEN Q;
1376 : ulong xQ, yQ, tQ, uQ;
1377 :
1378 73044 : Q = gcopy(pt);
1379 : /* Note that, as L is odd, say L = 2n + 1, we necessarily have
1380 : * [(L - 1)/2]P = [n]P = [n - L]P = -[n + 1]P = -[(L + 1)/2]P. This is
1381 : * what the condition Q[1] != xQ tests, so the loop will execute n times. */
1382 : do {
1383 642575 : xQ = uel(Q, 1);
1384 642575 : yQ = uel(Q, 2);
1385 : /* tQ = 6 xQ^2 + b2 xQ + b4
1386 : * = 6 xQ^2 + 2 a4 (since b2 = 0 and b4 = 2 a4) */
1387 642575 : tQ = Fl_add(Fl_mul6(Fl_sqr_pre(xQ, p, pi), p), Fl_double(a4, p), p);
1388 642542 : uQ = Fl_add(Fl_mul4(Fl_sqr_pre(yQ, p, pi), p),
1389 : Fl_mul_pre(tQ, xQ, p, pi), p);
1390 :
1391 642523 : t = Fl_add(t, tQ, p);
1392 642467 : w = Fl_add(w, uQ, p);
1393 642461 : Q = gerepileupto(av, Fle_add(pt, Q, a4, p));
1394 642578 : } while (uel(Q, 1) != xQ);
1395 :
1396 73048 : avma = av;
1397 : /* a4_img = a4 - 5 * t */
1398 73048 : *a4_img = Fl_sub(a4, Fl_mul5(t, p), p);
1399 : /* a6_img = a6 - b2 * t - 7 * w = a6 - 7 * w
1400 : * (since a1 = a2 = 0 ==> b2 = 0) */
1401 73048 : *a6_img = Fl_sub(a6, Fl_mul7(w, p), p);
1402 73046 : }
1403 :
1404 :
1405 : /**
1406 : * SECTION: Calculation of modular polynomials.
1407 : */
1408 :
1409 : /*
1410 : * Given an elliptic curve [a4, a6] over FF_p, try to find a
1411 : * non-trivial L-torsion point on the curve by considering n times a
1412 : * random point; val controls the maximum L-valuation expected of n
1413 : * times a random point.
1414 : */
1415 : static GEN
1416 106339 : find_L_tors_point(
1417 : ulong *ival,
1418 : ulong a4, ulong a6, ulong p, ulong pi,
1419 : ulong n, ulong L, ulong val)
1420 : {
1421 106339 : pari_sp av = avma;
1422 : ulong i;
1423 : GEN P, Q;
1424 : do {
1425 107353 : Q = random_Flj_pre(a4, a6, p, pi);
1426 107344 : P = Flj_mulu_pre(Q, n, a4, p, pi);
1427 107351 : } while (P[3] == 0);
1428 :
1429 206397 : for (i = 0; i < val; ++i) {
1430 173098 : Q = Flj_mulu_pre(P, L, a4, p, pi);
1431 173102 : if (Q[3] == 0)
1432 73042 : break;
1433 100060 : P = Q;
1434 : }
1435 106341 : if (ival)
1436 99791 : *ival = i;
1437 106341 : return gerepilecopy(av, P);
1438 : }
1439 :
1440 :
1441 : static GEN
1442 66503 : select_curve_with_L_tors_point(
1443 : ulong *a4, ulong *a6,
1444 : ulong L, ulong j, ulong n, ulong card, ulong val,
1445 : norm_eqn_t ne)
1446 : {
1447 66503 : pari_sp av = avma;
1448 : ulong A4, A4t, A6, A6t;
1449 66503 : ulong p = ne->p, pi = ne->pi;
1450 : GEN P;
1451 66503 : if (card % L != 0) {
1452 0 : pari_err_BUG("select_curve_with_L_tors_point: "
1453 : "Cardinality not divisible by L");
1454 : }
1455 :
1456 66503 : Fl_ellj_to_a4a6(j, p, &A4, &A6);
1457 66501 : Fl_elltwist_disc(A4, A6, ne->T, p, &A4t, &A6t);
1458 :
1459 : /* Either E = [a4, a6] or its twist has cardinality divisible by L
1460 : * because of the choice of p and t earlier on. We find out which
1461 : * by attempting to find a point of order L on each. See bot p16 of
1462 : * Sutherland 2012. */
1463 : while (1) {
1464 : ulong i;
1465 99798 : P = find_L_tors_point(&i, A4, A6, p, pi, n, L, val);
1466 99799 : if (i < val)
1467 66503 : break;
1468 33296 : avma = av;
1469 33296 : lswap(A4, A4t);
1470 33296 : lswap(A6, A6t);
1471 33296 : }
1472 :
1473 66503 : *a4 = A4;
1474 66503 : *a6 = A6;
1475 133004 : return gerepilecopy(av, P);
1476 : }
1477 :
1478 :
1479 : /*
1480 : * Return 1 if the L-Sylow subgroup of the curve [a4, a6] (mod p) is
1481 : * cyclic, return 0 if it is not cyclic with "high" probability (I
1482 : * guess around 1/L^3 chance it is still cyclic when we return 0).
1483 : *
1484 : * This code is based on Sutherland's
1485 : * velu.c:velu_verify_Sylow_cyclic() in classpoly-1.0.1.
1486 : */
1487 : INLINE long
1488 37032 : verify_L_sylow_is_cyclic(
1489 : long e, ulong a4, ulong a6, ulong p, ulong pi)
1490 : {
1491 : /* Number of times to try to find a point with maximal order in the
1492 : * L-Sylow subgroup. */
1493 : enum { N_RETRIES = 3 };
1494 37032 : pari_sp av = avma;
1495 37032 : long i, res = 0;
1496 : GEN P;
1497 59975 : for (i = 0; i < N_RETRIES; ++i) {
1498 53432 : P = random_Flj_pre(a4, a6, p, pi);
1499 53426 : P = Flj_mulu_pre(P, e, a4, p, pi);
1500 53433 : if (P[3] != 0) {
1501 30490 : res = 1;
1502 30490 : break;
1503 : }
1504 : }
1505 37033 : avma = av;
1506 37033 : return res;
1507 : }
1508 :
1509 :
1510 : static ulong
1511 66503 : find_noniso_L_isogenous_curve(
1512 : ulong L, ulong n,
1513 : norm_eqn_t ne, long e, ulong val, ulong a4, ulong a6, GEN init_pt, long verify)
1514 : {
1515 : pari_sp ltop, av;
1516 66503 : ulong p = ne->p, pi = ne->pi, j_res = 0;
1517 66503 : GEN pt = init_pt;
1518 66503 : ltop = av = avma;
1519 : while (1) {
1520 : /* c. Use Velu to calculate L-isogenous curve E' = E/<P> */
1521 : ulong a4_img, a6_img;
1522 73048 : ulong z2 = Fl_sqr_pre(pt[3], p, pi);
1523 73045 : pt = mkvecsmall2(Fl_div(pt[1], z2, p),
1524 73049 : Fl_div(pt[2], Fl_mul_pre(z2, pt[3], p, pi), p));
1525 73045 : Fle_quotient_from_kernel_generator(&a4_img, &a6_img,
1526 : a4, a6, pt, p, pi);
1527 :
1528 : /* d. If j(E') = j_res has a different endo ring to j(E), then
1529 : * return j(E'). Otherwise, go to b. */
1530 73045 : if (!verify || verify_L_sylow_is_cyclic(e, a4_img, a6_img, p, pi)) {
1531 66498 : j_res = Fl_ellj_pre(a4_img, a6_img, p, pi);
1532 66502 : break;
1533 : }
1534 :
1535 : /* b. Generate random point P on E of order L */
1536 6545 : avma = av;
1537 6545 : pt = find_L_tors_point(NULL, a4, a6, p, pi, n, L, val);
1538 6545 : }
1539 :
1540 66502 : avma = ltop;
1541 133004 : return j_res;
1542 : }
1543 :
1544 : /*
1545 : * Given a prime L and a j-invariant j (mod p), return the j-invariant
1546 : * of a curve which has a different endomorphism ring to j and is
1547 : * L-isogenous to j.
1548 : */
1549 : INLINE ulong
1550 66500 : compute_L_isogenous_curve(
1551 : ulong L, ulong n, norm_eqn_t ne,
1552 : ulong j, ulong card, ulong val, long verify)
1553 : {
1554 : ulong a4, a6;
1555 : long e;
1556 : GEN pt;
1557 :
1558 66500 : if (ne->p < 5 || j == 0 || j == 1728 % ne->p)
1559 0 : pari_err_BUG("compute_L_isogenous_curve");
1560 66503 : pt = select_curve_with_L_tors_point(&a4, &a6, L, j, n, card, val, ne);
1561 66500 : e = card / L;
1562 66500 : if (e * L != card) pari_err_BUG("compute_L_isogenous_curve");
1563 :
1564 66500 : return find_noniso_L_isogenous_curve(L, n, ne, e, val, a4, a6, pt, verify);
1565 : }
1566 :
1567 : INLINE GEN
1568 30490 : get_Lsqr_cycle(const modpoly_disc_info *dinfo)
1569 : {
1570 30490 : long i, n1 = dinfo->n1, L = dinfo->L;
1571 30490 : GEN cyc = cgetg(L, t_VECSMALL);
1572 30491 : cyc[1] = 0;
1573 266279 : for (i = 2; i <= L / 2; ++i)
1574 235788 : cyc[i] = cyc[i - 1] + n1;
1575 30491 : if ( ! dinfo->L1) {
1576 109784 : for ( ; i < L; ++i)
1577 97764 : cyc[i] = cyc[i - 1] + n1;
1578 : } else {
1579 18471 : cyc[L - 1] = 2 * dinfo->n2 - n1 / 2;
1580 168515 : for (i = L - 2; i > L / 2; --i)
1581 150044 : cyc[i] = cyc[i + 1] - n1;
1582 : }
1583 30491 : return cyc;
1584 : }
1585 :
1586 : INLINE void
1587 424444 : update_Lsqr_cycle(GEN cyc, const modpoly_disc_info *dinfo)
1588 : {
1589 424444 : long i, L = dinfo->L;
1590 13091865 : for (i = 1; i < L; ++i)
1591 12667421 : ++cyc[i];
1592 424444 : if (dinfo->L1 && cyc[L - 1] == 2 * dinfo->n2) {
1593 16857 : long n1 = dinfo->n1;
1594 156221 : for (i = L / 2 + 1; i < L; ++i)
1595 139364 : cyc[i] -= n1;
1596 : }
1597 424444 : }
1598 :
1599 :
1600 : static ulong
1601 30425 : oneroot_of_classpoly(
1602 : GEN hilb, GEN factu, norm_eqn_t ne, GEN jdb)
1603 : {
1604 30425 : pari_sp av = avma;
1605 30425 : ulong j0, p = ne->p, pi = ne->pi;
1606 30425 : long i, nfactors = lg(gel(factu, 1)) - 1;
1607 30425 : GEN hilbp = ZX_to_Flx(hilb, p);
1608 :
1609 : /* TODO: Work out how to use hilb with better invariant */
1610 30346 : j0 = Flx_oneroot_split(hilbp, p);
1611 30488 : if (j0 == p) {
1612 0 : pari_err_BUG("oneroot_of_classpoly: "
1613 : "Didn't find a root of the class polynomial");
1614 : }
1615 31396 : for (i = 1; i <= nfactors; ++i) {
1616 908 : long L = gel(factu, 1)[i];
1617 908 : long val = gel(factu, 2)[i];
1618 908 : GEN phi = polmodular_db_getp(jdb, L, p);
1619 908 : val += z_lval(ne->v, L);
1620 908 : j0 = descend_volcano(phi, j0, p, pi, 0, L, val, val);
1621 908 : avma = av;
1622 : }
1623 30488 : avma = av;
1624 30488 : return j0;
1625 : }
1626 :
1627 : /* TODO: Precompute the classgp_pcp_t structs and link them to dinfo */
1628 : INLINE void
1629 30488 : make_pcp_surface(const modpoly_disc_info *dinfo, classgp_pcp_t G)
1630 : {
1631 30488 : long k = 1, datalen = 3 * k;
1632 :
1633 30488 : memset(G, 0, sizeof *G);
1634 :
1635 30488 : G->_data = cgetg(datalen + 1, t_VECSMALL);
1636 30489 : G->L = G->_data + 1;
1637 30489 : G->n = G->L + k;
1638 30489 : G->o = G->L + k;
1639 :
1640 30489 : G->k = k;
1641 30489 : G->h = G->enum_cnt = dinfo->n1;
1642 30489 : G->L[0] = dinfo->L0;
1643 30489 : G->n[0] = dinfo->n1;
1644 30489 : G->o[0] = dinfo->n1;
1645 30489 : }
1646 :
1647 : INLINE void
1648 30489 : make_pcp_floor(const modpoly_disc_info *dinfo, classgp_pcp_t G)
1649 : {
1650 30489 : long k = dinfo->L1 ? 2 : 1, datalen = 3 * k;
1651 :
1652 30489 : memset(G, 0, sizeof *G);
1653 30489 : G->_data = cgetg(datalen + 1, t_VECSMALL);
1654 30489 : G->L = G->_data + 1;
1655 30489 : G->n = G->L + k;
1656 30489 : G->o = G->L + k;
1657 :
1658 30489 : G->k = k;
1659 30489 : G->h = G->enum_cnt = dinfo->n2 * k;
1660 30489 : G->L[0] = dinfo->L0;
1661 30489 : G->n[0] = dinfo->n2;
1662 30489 : G->o[0] = dinfo->n2;
1663 30489 : if (dinfo->L1) {
1664 18471 : G->L[1] = dinfo->L1;
1665 18471 : G->n[1] = 2;
1666 18471 : G->o[1] = 2;
1667 : }
1668 30489 : }
1669 :
1670 : INLINE GEN
1671 30489 : enum_volcano_surface(
1672 : const modpoly_disc_info *dinfo, norm_eqn_t ne, ulong j0, GEN fdb)
1673 : {
1674 30489 : pari_sp av = avma;
1675 : classgp_pcp_t G;
1676 30489 : make_pcp_surface(dinfo, G);
1677 30490 : return gerepileupto(av, enum_roots(j0, ne, fdb, G));
1678 : }
1679 :
1680 :
1681 : INLINE GEN
1682 30490 : enum_volcano_floor(
1683 : long L, norm_eqn_t ne, ulong j0_pr, GEN fdb,
1684 : const modpoly_disc_info *dinfo)
1685 : {
1686 30490 : pari_sp av = avma;
1687 : /* L^2 D is the discriminant for the order R = Z + L OO. */
1688 30490 : long DR = L * L * ne->D;
1689 30490 : long R_cond = L * ne->u; /* conductor(DR); */
1690 30490 : long w = R_cond * ne->v;
1691 : /* TODO: Calculate these once and for all in polmodular0_ZM(). */
1692 : classgp_pcp_t G;
1693 : norm_eqn_t eqn;
1694 30490 : memcpy(eqn, ne, sizeof *ne);
1695 30490 : eqn->D = DR;
1696 30490 : eqn->u = R_cond;
1697 30490 : eqn->v = w;
1698 30490 : make_pcp_floor(dinfo, G);
1699 30489 : return gerepileupto(av, enum_roots(j0_pr, eqn, fdb, G));
1700 : }
1701 :
1702 : INLINE void
1703 14829 : carray_reverse_inplace(long *arr, long n)
1704 : {
1705 14829 : long lim = n>>1, i;
1706 14829 : --n;
1707 154320 : for (i = 0; i < lim; i++)
1708 139491 : lswap(arr[i], arr[n - i]);
1709 14829 : }
1710 :
1711 : INLINE void
1712 454918 : append_neighbours(GEN rts, GEN surface_js, long njs, long L, long m, long i)
1713 : {
1714 454918 : long r_idx = (((i - 1) + m) % njs) + 1; /* (i + m) % njs */
1715 454918 : long l_idx = smodss((i - 1) - m, njs) + 1; /* (i - m) % njs */
1716 454928 : rts[L] = surface_js[l_idx];
1717 454928 : rts[L + 1] = surface_js[r_idx];
1718 454928 : }
1719 :
1720 : INLINE GEN
1721 32687 : roots_to_coeffs(GEN rts, ulong p, long L)
1722 : {
1723 32687 : long i, k, lrts= lg(rts);
1724 32687 : GEN M = cgetg(L+2+1, t_MAT);
1725 736221 : for (i = 1; i <= L+2; ++i)
1726 703521 : gel(M, i) = cgetg(lrts, t_VECSMALL);
1727 512758 : for (i = 1; i < lrts; ++i) {
1728 480073 : pari_sp av = avma;
1729 480073 : GEN modpol = Flv_roots_to_pol(gel(rts, i), p, 0);
1730 16439138 : for (k = 1; k <= L + 2; ++k)
1731 15959080 : mael(M, k, i) = modpol[k + 1];
1732 480058 : avma = av;
1733 : }
1734 32685 : return M;
1735 : }
1736 :
1737 :
1738 : /* NB: Assumes indices are offset at 0, not at 1 like in GENs;
1739 : * i.e. indices[i] will pick out v[indices[i] + 1] from v. */
1740 : INLINE void
1741 454934 : vecsmall_pick(GEN res, GEN v, GEN indices)
1742 : {
1743 : long i;
1744 13654449 : for (i = 1; i < lg(indices); ++i)
1745 13199515 : res[i] = v[indices[i] + 1];
1746 454934 : }
1747 :
1748 :
1749 : /* The first element of surface_js must lie above the first element of
1750 : * floor_js. Will reverse surface_js if it is not oriented in the
1751 : * same direction as floor_js. */
1752 : INLINE GEN
1753 30491 : root_matrix(
1754 : long L, const modpoly_disc_info *dinfo,
1755 : long njinvs, GEN surface_js, GEN floor_js,
1756 : ulong n, ulong card, ulong val, norm_eqn_t ne)
1757 : {
1758 : pari_sp av;
1759 30491 : long i, m = dinfo->dl1, njs = lg(surface_js) - 1, inv = dinfo->inv, rev;
1760 30491 : GEN rt_mat = zero_Flm_copy(L + 1, njinvs), rts, cyc;
1761 30490 : av = avma;
1762 :
1763 30490 : i = 1;
1764 30490 : cyc = get_Lsqr_cycle(dinfo);
1765 30490 : rts = gel(rt_mat, i);
1766 30490 : vecsmall_pick(rts, floor_js, cyc);
1767 30490 : append_neighbours(rts, surface_js, njs, L, m, i);
1768 :
1769 30491 : i = 2;
1770 30491 : update_Lsqr_cycle(cyc, dinfo);
1771 30491 : rts = gel(rt_mat, i);
1772 30491 : vecsmall_pick(rts, floor_js, cyc);
1773 :
1774 : /* Fix orientation if necessary */
1775 30491 : if (modinv_is_double_eta(inv)) {
1776 : /* TODO: There is potential for refactoring between this,
1777 : * double_eta_initial_js and modfn_preimage. */
1778 5522 : pari_sp av0 = avma;
1779 5522 : ulong p = ne->p, pi = ne->pi, j;
1780 5522 : GEN F = double_eta_Fl(inv, p);
1781 5522 : pari_sp av = avma;
1782 5522 : ulong r1 = double_eta_power(inv, uel(rts, 1), p, pi);
1783 5522 : GEN r, f = Flx_double_eta_jpoly(F, r1, p, pi);
1784 5522 : if ((j = Flx_oneroot(f, p)) == p) pari_err_BUG("root_matrix");
1785 5522 : j = compute_L_isogenous_curve(L, n, ne, j, card, val, 0);
1786 5522 : avma = av;
1787 5522 : r1 = double_eta_power(inv, uel(surface_js, i), p, pi);
1788 5522 : f = Flx_double_eta_jpoly(F, r1, p, pi);
1789 5522 : r = Flx_roots(f, p);
1790 5522 : if (glength(r) != 2) pari_err_BUG("root_matrix");
1791 5522 : rev = (j != uel(r, 1)) && (j != uel(r, 2));
1792 5522 : avma = av0;
1793 : } else {
1794 : ulong j1pr, j1;
1795 24968 : j1pr = modfn_preimage(uel(rts, 1), ne, dinfo->inv);
1796 24968 : j1 = compute_L_isogenous_curve(L, n, ne, j1pr, card, val, 0);
1797 24969 : rev = j1 != modfn_preimage(uel(surface_js, i), ne, dinfo->inv);
1798 : }
1799 30489 : if (rev)
1800 14829 : carray_reverse_inplace(surface_js + 2, njs - 1);
1801 30489 : append_neighbours(rts, surface_js, njs, L, m, i);
1802 :
1803 424432 : for (i = 3; i <= njinvs; ++i) {
1804 393941 : update_Lsqr_cycle(cyc, dinfo);
1805 393953 : rts = gel(rt_mat, i);
1806 393953 : vecsmall_pick(rts, floor_js, cyc);
1807 393945 : append_neighbours(rts, surface_js, njs, L, m, i);
1808 : }
1809 30491 : avma = av;
1810 30491 : return rt_mat;
1811 : }
1812 :
1813 :
1814 : INLINE void
1815 33002 : interpolate_coeffs(GEN phi_modp, ulong p, GEN j_invs, GEN coeff_mat)
1816 : {
1817 33002 : pari_sp av = avma;
1818 : long i;
1819 33002 : GEN pols = Flv_Flm_polint(j_invs, coeff_mat, p, 0);
1820 738698 : for (i = 1; i < lg(pols); ++i) {
1821 705695 : GEN pol = gel(pols, i);
1822 705695 : long k, maxk = lg(pol);
1823 15630198 : for (k = 2; k < maxk; ++k)
1824 14924503 : coeff(phi_modp, k - 1, i) = pol[k];
1825 : }
1826 33003 : avma = av;
1827 33003 : }
1828 :
1829 :
1830 : INLINE long
1831 329084 : Flv_lastnonzero(GEN v)
1832 : {
1833 : long i;
1834 24771334 : for (i = lg(v) - 1; i > 0; --i)
1835 24770527 : if (v[i])
1836 328277 : break;
1837 329084 : return i;
1838 : }
1839 :
1840 :
1841 : /*
1842 : * Assuming the matrix of coefficients in phi corresponds to
1843 : * polynomials phi_k^* satisfying Y^c phi_k^*(Y^s) for c in {0, 1,
1844 : * ..., s} satisfying c + Lk = L + 1 (mod s), change phi so that the
1845 : * coefficients are for the polynomials Y^c phi_k^*(Y^s) (s is the
1846 : * sparsity factor).
1847 : */
1848 : INLINE void
1849 9719 : inflate_polys(GEN phi, long L, long s)
1850 : {
1851 9719 : long k, deg = L + 1;
1852 : long maxr;
1853 9719 : maxr = nbrows(phi);
1854 348520 : for (k = 0; k <= deg; ) {
1855 329081 : long i, c = smodss(L * (1 - k) + 1, s);
1856 : /* TODO: We actually know that the last non-zero element of
1857 : * gel(phi, k) can't be later than index n + 1, where n is about
1858 : * (L + 1)/s. */
1859 329072 : ++k;
1860 5420638 : for (i = Flv_lastnonzero(gel(phi, k)); i > 0; --i) {
1861 5091566 : long r = c + (i - 1) * s + 1;
1862 5091566 : if (r > maxr) {
1863 63175 : coeff(phi, i, k) = 0;
1864 63175 : continue;
1865 : }
1866 5028391 : if (r != i) {
1867 4926381 : coeff(phi, r, k) = coeff(phi, i, k);
1868 4926381 : coeff(phi, i, k) = 0;
1869 : }
1870 : }
1871 : }
1872 9720 : }
1873 :
1874 : INLINE void
1875 39570 : Flv_powu_inplace_pre(GEN v, ulong n, ulong p, ulong pi)
1876 : {
1877 : long i;
1878 326383 : for (i = 1; i < lg(v); ++i)
1879 286812 : v[i] = Fl_powu_pre(v[i], n, p, pi);
1880 39571 : }
1881 :
1882 : INLINE void
1883 9720 : normalise_coeffs(GEN coeffs, GEN js, long L, long s, ulong p, ulong pi)
1884 : {
1885 9720 : pari_sp av = avma;
1886 : long k;
1887 : GEN pows, modinv_js;
1888 :
1889 : /* NB: In fact it would be correct to return the coefficients "as
1890 : * is" when s = 1, but we make that an error anyway since this
1891 : * function should never be called with s = 1. */
1892 9720 : if (s <= 1) pari_err_BUG("normalise_coeffs");
1893 :
1894 : /* pows[i + 1] contains 1 / js[i + 1]^i for i = 0, ..., s - 1. */
1895 9720 : pows = cgetg(s + 1, t_VEC);
1896 9720 : gel(pows, 1) = const_vecsmall(lg(js) - 1, 1);
1897 9720 : modinv_js = Flv_inv_pre(js, p, pi);
1898 9719 : gel(pows, 2) = modinv_js;
1899 37435 : for (k = 3; k <= s; ++k) {
1900 27716 : gel(pows, k) = gcopy(modinv_js);
1901 27715 : Flv_powu_inplace_pre(gel(pows, k), k - 1, p, pi);
1902 : }
1903 :
1904 : /* For each column of coefficients coeffs[k] = [a0 .. an],
1905 : * replace ai by ai / js[i]^c.
1906 : * Said in another way, normalise each row i of coeffs by
1907 : * dividing through by js[i - 1]^c (where c depends on i). */
1908 338790 : for (k = 1; k < lg(coeffs); ++k) {
1909 329070 : long i, c = smodss(L * (1 - (k - 1)) + 1, s);
1910 329011 : GEN col = gel(coeffs, k), C = gel(pows, c + 1);
1911 5786141 : for (i = 1; i < lg(col); ++i)
1912 5457070 : col[i] = Fl_mul_pre(col[i], C[i], p, pi);
1913 : }
1914 9720 : avma = av;
1915 9720 : }
1916 :
1917 : INLINE void
1918 5522 : double_eta_initial_js(
1919 : ulong *x0, ulong *x0pr, ulong j0, ulong j0pr, norm_eqn_t ne,
1920 : long inv, ulong L, ulong n, ulong card, ulong val)
1921 : {
1922 5522 : pari_sp av0 = avma;
1923 5522 : ulong p = ne->p, pi = ne->pi, s2 = ne->s2;
1924 5522 : GEN F = double_eta_Fl(inv, p);
1925 5522 : pari_sp av = avma;
1926 : ulong j1pr, j1, r, t;
1927 : GEN f, g;
1928 :
1929 5522 : *x0pr = modinv_double_eta_from_j(F, inv, j0pr, p, pi, s2);
1930 5522 : t = double_eta_power(inv, *x0pr, p, pi);
1931 5522 : f = Flx_div_by_X_x(Flx_double_eta_jpoly(F, t, p, pi), j0pr, p, &r);
1932 5522 : if (r) pari_err_BUG("double_eta_initial_js");
1933 5522 : j1pr = Flx_deg1_root(f, p);
1934 5522 : avma = av;
1935 :
1936 5522 : j1 = compute_L_isogenous_curve(L, n, ne, j1pr, card, val, 0);
1937 5522 : f = Flx_double_eta_xpoly(F, j0, p, pi);
1938 5522 : g = Flx_double_eta_xpoly(F, j1, p, pi);
1939 : /* x0 is the unique common root of f and g */
1940 5522 : *x0 = Flx_deg1_root(Flx_gcd(f, g, p), p);
1941 5522 : avma = av0;
1942 :
1943 5522 : if ( ! double_eta_root(inv, x0, *x0, p, pi, s2))
1944 0 : pari_err_BUG("double_eta_initial_js");
1945 5522 : }
1946 :
1947 : /*
1948 : * This is Sutherland 2012, Algorithm 2.1, p16.
1949 : */
1950 : static GEN
1951 30448 : polmodular_split_p_Flm(
1952 : ulong L, GEN hilb, GEN factu, norm_eqn_t ne, GEN db,
1953 : const modpoly_disc_info *dinfo)
1954 : {
1955 30448 : pari_sp ltop = avma;
1956 : ulong j0, j0_rt, j0pr, j0pr_rt;
1957 30448 : ulong n, card, val, p = ne->p, pi = ne->pi;
1958 30448 : long s = modinv_sparse_factor(dinfo->inv);
1959 30449 : long nj_selected = ceil((L + 1)/(double)s) + 1;
1960 : GEN surface_js, floor_js, rts;
1961 : GEN phi_modp;
1962 : GEN jdb, fdb;
1963 30449 : long switched_signs = 0;
1964 :
1965 30449 : jdb = polmodular_db_for_inv(db, INV_J);
1966 30436 : fdb = polmodular_db_for_inv(db, dinfo->inv);
1967 :
1968 : /* Precomputation */
1969 30441 : card = p + 1 - ne->t;
1970 30441 : val = u_lvalrem(card, L, &n); /* n = card / L^{v_L(card)} */
1971 :
1972 30448 : j0 = oneroot_of_classpoly(hilb, factu, ne, jdb);
1973 30488 : j0pr = compute_L_isogenous_curve(L, n, ne, j0, card, val, 1);
1974 30490 : if (modinv_is_double_eta(dinfo->inv)) {
1975 5522 : double_eta_initial_js(&j0_rt, &j0pr_rt, j0, j0pr, ne, dinfo->inv,
1976 : L, n, card, val);
1977 : } else {
1978 24968 : j0_rt = modfn_root(j0, ne, dinfo->inv);
1979 24968 : j0pr_rt = modfn_root(j0pr, ne, dinfo->inv);
1980 : }
1981 30489 : surface_js = enum_volcano_surface(dinfo, ne, j0_rt, fdb);
1982 30490 : floor_js = enum_volcano_floor(L, ne, j0pr_rt, fdb, dinfo);
1983 30491 : rts = root_matrix(L, dinfo, nj_selected, surface_js, floor_js,
1984 : n, card, val, ne);
1985 : do {
1986 32688 : pari_sp btop = avma;
1987 : long i;
1988 : GEN coeffs, surf;
1989 :
1990 32688 : coeffs = roots_to_coeffs(rts, p, L);
1991 32682 : surf = vecsmall_shorten(surface_js, nj_selected);
1992 32686 : if (s > 1) {
1993 9720 : normalise_coeffs(coeffs, surf, L, s, p, pi);
1994 9720 : Flv_powu_inplace_pre(surf, s, p, pi);
1995 : }
1996 32686 : phi_modp = zero_Flm_copy(L + 2, L + 2);
1997 32688 : interpolate_coeffs(phi_modp, p, surf, coeffs);
1998 32685 : if (s > 1)
1999 9719 : inflate_polys(phi_modp, L, s);
2000 :
2001 : /* TODO: Calculate just this coefficient of X^L Y^L, so we can do
2002 : * this test, then calculate the other coefficients; at the moment
2003 : * we are sometimes doing all the roots-to-coeffs, normalisation
2004 : * and interpolation work twice. */
2005 32684 : if (ucoeff(phi_modp, L + 1, L + 1) == p - 1)
2006 30487 : break;
2007 :
2008 2197 : if (switched_signs) pari_err_BUG("polmodular_split_p_Flm");
2009 :
2010 2197 : avma = btop;
2011 27394 : for (i = 1; i < lg(rts); ++i) {
2012 25197 : surface_js[i] = Fl_neg(surface_js[i], p);
2013 25197 : coeff(rts, L, i) = Fl_neg(coeff(rts, L, i), p);
2014 25197 : coeff(rts, L + 1, i) = Fl_neg(coeff(rts, L + 1, i), p);
2015 : }
2016 2197 : switched_signs = 1;
2017 2197 : } while (1);
2018 30487 : dbg_printf(4)(" Phi_%lu(X, Y) (mod %lu) = %Ps\n", L, p, phi_modp);
2019 :
2020 30487 : return gerepileupto(ltop, phi_modp);
2021 : }
2022 :
2023 :
2024 : INLINE void
2025 2451 : norm_eqn_init(norm_eqn_t ne, long D, long u)
2026 : {
2027 2451 : memset(ne, 0, sizeof(*ne));
2028 2451 : ne->D = D;
2029 2451 : ne->u = u;
2030 2451 : }
2031 :
2032 :
2033 : INLINE void
2034 30316 : norm_eqn_update(norm_eqn_t ne, ulong t, ulong p, long L)
2035 : {
2036 : long res;
2037 : ulong vL_sqr, vL;
2038 :
2039 30316 : ne->t = t;
2040 30316 : ne->p = p;
2041 30316 : ne->pi = get_Fl_red(p);
2042 30314 : ne->s2 = Fl_2gener_pre(p, ne->pi);
2043 :
2044 30459 : vL_sqr = (4 * p - t * t) / -ne->D;
2045 30459 : res = uissquareall(vL_sqr, &vL);
2046 30427 : if ( ! res || vL % L) pari_err_BUG("norm_eqn_update");
2047 30439 : ne->v = vL;
2048 :
2049 : /* Select twisting parameter. */
2050 60673 : do ne->T = random_Fl(p);
2051 60732 : while (krouu(ne->T, p) != -1);
2052 30461 : }
2053 :
2054 : INLINE void
2055 2407 : Flv_deriv_pre_inplace(GEN v, long deg, ulong p, ulong pi)
2056 : {
2057 2407 : long i, ln = lg(v), d = deg % p;
2058 55817 : for (i = ln - 1; i > 1; --i, --d)
2059 53409 : v[i] = Fl_mul_pre(v[i - 1], d, p, pi);
2060 2408 : v[1] = 0;
2061 2408 : }
2062 :
2063 : /* NB: Deliberately leave a dirty stack, since the result must be
2064 : * gerepileupto'd straight away in any case. */
2065 : INLINE GEN
2066 2408 : eval_modpoly_modp(
2067 : GEN modpoly_modp, GEN j_powers, norm_eqn_t ne, int compute_derivs)
2068 : {
2069 2408 : ulong p = ne->p, pi = ne->pi;
2070 2408 : long L = lg(j_powers) - 3;
2071 2408 : GEN j_pows_p = ZV_to_Flv(j_powers, p);
2072 2408 : GEN tmp = cgetg(2 + 2 * compute_derivs, t_VEC);
2073 : /* We wrap the result in this t_VEC modpoly_modp to trick the
2074 : * ZM_*_CRT() functions into thinking it's a matrix. */
2075 2407 : gel(tmp, 1) = Flm_Flc_mul_pre(modpoly_modp, j_pows_p, p, pi);
2076 2408 : if (compute_derivs) {
2077 1204 : Flv_deriv_pre_inplace(j_pows_p, L + 1, p, pi);
2078 1204 : gel(tmp, 2) = Flm_Flc_mul_pre(modpoly_modp, j_pows_p, p, pi);
2079 1204 : Flv_deriv_pre_inplace(j_pows_p, L + 1, p, pi);
2080 1204 : gel(tmp, 3) = Flm_Flc_mul_pre(modpoly_modp, j_pows_p, p, pi);
2081 : }
2082 2408 : return tmp;
2083 : }
2084 :
2085 : /* Parallel interface */
2086 :
2087 : static void
2088 30321 : vne_to_ne(norm_eqn_t ne, GEN vne)
2089 : {
2090 30321 : ne->D = vne[1];
2091 30321 : ne->u = vne[2];
2092 30321 : }
2093 :
2094 : static GEN
2095 2451 : ne_to_vne(norm_eqn_t ne)
2096 : {
2097 2451 : return mkvecsmall2(ne->D, ne->u);
2098 : }
2099 :
2100 : static void
2101 30329 : vinfo_to_dinfo(modpoly_disc_info *dinfo, GEN vinfo)
2102 : {
2103 30329 : dinfo->L = vinfo[1];
2104 30329 : dinfo->D0 = vinfo[2];
2105 30329 : dinfo->D1 = vinfo[3];
2106 30329 : dinfo->L0 = vinfo[4];
2107 30329 : dinfo->L1 = vinfo[5];
2108 30329 : dinfo->n1 = vinfo[6];
2109 30329 : dinfo->n2 = vinfo[7];
2110 30329 : dinfo->nprimes = vinfo[8];
2111 30329 : dinfo->dl1 = vinfo[9];
2112 30329 : dinfo->dl2_0 = vinfo[10];
2113 30329 : dinfo->dl2_1 = vinfo[11];
2114 30329 : dinfo->inv = vinfo[12];
2115 30329 : }
2116 :
2117 : static GEN
2118 2451 : dinfo_to_vinfo(const modpoly_disc_info *dinfo)
2119 : {
2120 2451 : return mkvecsmalln(12, dinfo->L, dinfo->D0, dinfo->D1, dinfo->L0, dinfo->L1,
2121 : dinfo->n1, dinfo->n2, dinfo->nprimes, dinfo->dl1, dinfo->dl2_0,
2122 : dinfo->dl2_1, dinfo->inv);
2123 : }
2124 :
2125 : GEN
2126 30330 : polmodular_worker(ulong p, ulong t,
2127 : ulong L, GEN hilb, GEN factu, GEN vne, GEN vinfo,
2128 : long compute_derivs, GEN j_powers, GEN fdb)
2129 : {
2130 30330 : pari_sp av = avma;
2131 : GEN modpoly_modp;
2132 : modpoly_disc_info dinfo;
2133 : norm_eqn_t ne;
2134 30330 : vinfo_to_dinfo(&dinfo, vinfo);
2135 30396 : vne_to_ne(ne, vne);
2136 30323 : norm_eqn_update(ne, t, p, L);
2137 30460 : modpoly_modp = polmodular_split_p_Flm(L, hilb, factu, ne, fdb, &dinfo);
2138 30486 : if (!isintzero(j_powers)) {
2139 2408 : modpoly_modp = eval_modpoly_modp(modpoly_modp, j_powers, ne, compute_derivs);
2140 2408 : modpoly_modp = gerepileupto(av, modpoly_modp);
2141 : }
2142 30487 : return modpoly_modp;
2143 : }
2144 :
2145 : static GEN
2146 17941 : sympol_to_ZM(GEN phi, long L)
2147 : {
2148 17941 : pari_sp av = avma;
2149 17941 : GEN res = zeromatcopy(L + 2, L + 2);
2150 17941 : long i, j, c = 1;
2151 77833 : for (i = 1; i <= L + 1; ++i) {
2152 197099 : for (j = 1; j <= i; ++j, ++c)
2153 137207 : gcoeff(res, i, j) = gcoeff(res, j, i) = gel(phi, c);
2154 : }
2155 17941 : gcoeff(res, L + 2, 1) = gcoeff(res, 1, L + 2) = gen_1;
2156 17941 : return gerepilecopy(av, res);
2157 : }
2158 :
2159 : static GEN
2160 : polmodular_small_ZM(long L, long inv, GEN *db);
2161 :
2162 : INLINE long
2163 20619 : modinv_max_internal_level(long inv)
2164 : {
2165 20619 : switch (inv) {
2166 : case INV_J:
2167 18222 : return 5;
2168 : case INV_G2:
2169 203 : return 2;
2170 : case INV_F:
2171 : case INV_F2:
2172 : case INV_F4:
2173 : case INV_F8:
2174 485 : return 5;
2175 : case INV_W2W5:
2176 : case INV_W2W5E2:
2177 238 : return 7;
2178 : case INV_W2W3:
2179 : case INV_W2W3E2:
2180 : case INV_W3W3:
2181 : case INV_W3W7:
2182 476 : return 5;
2183 : case INV_W3W3E2:
2184 63 : return 2;
2185 : case INV_F3:
2186 : case INV_W2W7:
2187 : case INV_W2W7E2:
2188 : case INV_W2W13:
2189 638 : return 3;
2190 : case INV_W3W5:
2191 : case INV_W5W7:
2192 : case INV_W3W13:
2193 294 : return 2;
2194 : }
2195 0 : pari_err_BUG("modinv_max_internal_level");
2196 0 : return LONG_MAX;
2197 : }
2198 :
2199 : GEN
2200 20514 : polmodular0_ZM(
2201 : long L, long inv, GEN J, GEN Q, int compute_derivs, GEN *db)
2202 : {
2203 20514 : pari_sp ltop = avma;
2204 20514 : long k, d, Dcnt, nprimes = 0;
2205 : GEN modpoly, plist;
2206 : modpoly_disc_info Ds[MODPOLY_MAX_DCNT];
2207 :
2208 20514 : long lvl = modinv_level(inv);
2209 20514 : if (cgcd(L, lvl) != 1) {
2210 7 : pari_err_DOMAIN("polmodular0_ZM", "invariant",
2211 : "incompatible with", stoi(L), stoi(lvl));
2212 : }
2213 :
2214 20507 : dbg_printf(1)("Calculating modular polynomial of level %lu for invariant %d\n", L, inv);
2215 20507 : if (L <= modinv_max_internal_level(inv))
2216 18074 : return polmodular_small_ZM(L, inv, db);
2217 :
2218 2433 : Dcnt = discriminant_with_classno_at_least(Ds, L, inv, USE_SPARSE_FACTOR);
2219 2433 : for (d = 0; d < Dcnt; ++d) nprimes += Ds[d].nprimes;
2220 2433 : modpoly = cgetg(nprimes+1, t_VEC);
2221 2433 : plist = cgetg(nprimes+1, t_VECSMALL);
2222 2433 : k = 1;
2223 :
2224 4884 : for (d = 0; d < Dcnt; ++d) {
2225 : long D, DK, i;
2226 : ulong cond;
2227 : GEN j_powers, factu, hilb;
2228 : norm_eqn_t ne;
2229 2451 : modpoly_disc_info *dinfo = &Ds[d];
2230 : GEN worker;
2231 : struct pari_mt pt;
2232 2451 : long pending = 0;
2233 :
2234 2451 : polmodular_db_add_level(db, dinfo->L0, inv);
2235 2451 : if (dinfo->L1)
2236 1155 : polmodular_db_add_level(db, dinfo->L1, inv);
2237 :
2238 2451 : D = dinfo->D1;
2239 2451 : DK = dinfo->D0;
2240 2451 : cond = sqrt((double)(D / DK));
2241 2451 : factu = factoru(cond);
2242 2451 : dbg_printf(1)("Selected discriminant D = %ld = %ld^2 * %ld.\n",
2243 : D, cond, DK);
2244 :
2245 2451 : hilb = polclass0(DK, INV_J, 0, db);
2246 2451 : norm_eqn_init(ne, D, cond);
2247 :
2248 2451 : if (cond > 1)
2249 31 : polmodular_db_add_levels(db, zv_to_longptr(gel(factu, 1)), glength(gel(factu, 1)), INV_J);
2250 :
2251 2451 : dbg_printf(1)("D = %ld, L0 = %lu, L1 = %lu, ", dinfo->D1, dinfo->L0, dinfo->L1);
2252 2451 : dbg_printf(1)("n1 = %lu, n2 = %lu, dl1 = %lu, dl2_0 = %lu, dl2_1 = %lu\n",
2253 : dinfo->n1, dinfo->n2, dinfo->dl1, dinfo->dl2_0, dinfo->dl2_1);
2254 2451 : dbg_printf(0)("Calculating modular polynomial of level %lu:", L);
2255 :
2256 2451 : j_powers = gen_0;
2257 2451 : if (J) {
2258 56 : compute_derivs = !!compute_derivs;
2259 56 : j_powers = Fp_powers(J, L + 1, Q);
2260 : }
2261 2451 : worker = strtoclosure("_polmodular_worker", 8, utoi(L), hilb, factu, ne_to_vne(ne),
2262 : dinfo_to_vinfo(dinfo),
2263 : stoi(compute_derivs), j_powers, *db);
2264 2451 : mt_queue_start_lim(&pt, worker, dinfo->nprimes);
2265 36069 : for (i = 0; i < dinfo->nprimes || pending; ++i)
2266 : {
2267 : GEN done;
2268 : long workid;
2269 33618 : ulong p = dinfo->primes[i];
2270 33618 : ulong t = dinfo->traces[i];
2271 33618 : mt_queue_submit(&pt, p, i < dinfo-> nprimes? mkvec2(utoi(p), utoi(t)): NULL);
2272 33618 : done = mt_queue_get(&pt, &workid, &pending);
2273 33618 : if (done)
2274 : {
2275 30491 : gel(modpoly, k) = done;
2276 30491 : plist[k] = workid;
2277 30491 : k++;
2278 30491 : dbg_printf(0)(" %ld%%", k*100/nprimes);
2279 : }
2280 : }
2281 2451 : dbg_printf(0)("\n");
2282 2451 : mt_queue_end(&pt);
2283 2451 : modpoly_disc_info_clear(dinfo);
2284 : }
2285 2433 : modpoly = nmV_chinese_center(modpoly, plist, NULL);
2286 2433 : if (J)
2287 56 : return gerepileupto(ltop, FpM_red(modpoly, Q));
2288 2377 : return gerepileupto(ltop, modpoly);
2289 : }
2290 :
2291 :
2292 : GEN
2293 14511 : polmodular_ZM(long L, long inv)
2294 : {
2295 : GEN db, Phi;
2296 :
2297 14511 : if (L < 2)
2298 7 : pari_err_DOMAIN("polmodular_ZM", "L", "<", gen_2, stoi(L));
2299 :
2300 : /* TODO: Handle non-prime L. This is Algorithm 1.1 and Corollary
2301 : * 3.4 in Sutherland, "Class polynomials for nonholomorphic modular
2302 : * functions". */
2303 14504 : if ( ! uisprime(L))
2304 7 : pari_err_IMPL("composite level");
2305 :
2306 14497 : db = polmodular_db_init(inv);
2307 14497 : Phi = polmodular0_ZM(L, inv, NULL, NULL, 0, &db);
2308 14490 : gunclone_deep(db);
2309 :
2310 14490 : return Phi;
2311 : }
2312 :
2313 :
2314 : GEN
2315 14434 : polmodular_ZXX(long L, long inv, long vx, long vy)
2316 : {
2317 14434 : pari_sp av = avma;
2318 14434 : GEN phi = polmodular_ZM(L, inv);
2319 :
2320 14413 : if (vx < 0) vx = 0;
2321 14413 : if (vy < 0) vy = 1;
2322 14413 : if (varncmp(vx, vy) >= 0)
2323 14 : pari_err_PRIORITY("polmodular_ZXX", pol_x(vx), "<=", vy);
2324 :
2325 14399 : return gerepilecopy(av, RgM_to_RgXX(phi, vx, vy));
2326 : }
2327 :
2328 : INLINE GEN
2329 56 : FpV_deriv(GEN v, long deg, GEN P)
2330 : {
2331 56 : long i, ln = lg(v);
2332 56 : GEN dv = cgetg(ln, t_VEC);
2333 392 : for (i = ln - 1; i > 1; --i, --deg)
2334 336 : gel(dv, i) = Fp_mulu(gel(v, i - 1), deg, P);
2335 56 : gel(dv, 1) = gen_0;
2336 56 : return dv;
2337 : }
2338 :
2339 : GEN
2340 112 : Fp_polmodular_evalx(
2341 : long L, long inv, GEN J, GEN P, long v, int compute_derivs)
2342 : {
2343 112 : pari_sp av = avma;
2344 : GEN db, phi;
2345 :
2346 112 : if (L <= modinv_max_internal_level(inv)) {
2347 : GEN tmp;
2348 56 : GEN phi = RgM_to_FpM(polmodular_ZM(L, inv), P);
2349 56 : GEN j_powers = Fp_powers(J, L + 1, P);
2350 56 : GEN modpol = RgV_to_RgX(FpM_FpC_mul(phi, j_powers, P), v);
2351 56 : if (compute_derivs) {
2352 28 : tmp = cgetg(4, t_VEC);
2353 28 : gel(tmp, 1) = modpol;
2354 28 : j_powers = FpV_deriv(j_powers, L + 1, P);
2355 28 : gel(tmp, 2) = RgV_to_RgX(FpM_FpC_mul(phi, j_powers, P), v);
2356 28 : j_powers = FpV_deriv(j_powers, L + 1, P);
2357 28 : gel(tmp, 3) = RgV_to_RgX(FpM_FpC_mul(phi, j_powers, P), v);
2358 : } else {
2359 28 : tmp = modpol;
2360 : }
2361 56 : return gerepilecopy(av, tmp);
2362 : }
2363 :
2364 56 : db = polmodular_db_init(inv);
2365 56 : phi = polmodular0_ZM(L, inv, J, P, compute_derivs, &db);
2366 56 : phi = RgM_to_RgXV(phi, v);
2367 56 : gunclone_deep(db);
2368 56 : return gerepilecopy(av, compute_derivs ? phi : gel(phi, 1));
2369 : }
2370 :
2371 : GEN
2372 609 : polmodular(long L, long inv, GEN x, long v, long compute_derivs)
2373 : {
2374 609 : pari_sp av = avma;
2375 : long tx;
2376 609 : GEN J = NULL, P = NULL, res = NULL, one = NULL;
2377 :
2378 609 : check_modinv(inv);
2379 602 : if ( ! x || gequalX(x)) {
2380 476 : long xv = 0;
2381 476 : if (x)
2382 42 : xv = varn(x);
2383 476 : if (compute_derivs)
2384 7 : pari_err_FLAG("polmodular");
2385 469 : return polmodular_ZXX(L, inv, xv, v);
2386 : }
2387 :
2388 126 : tx = typ(x);
2389 126 : if (tx == t_INTMOD) {
2390 56 : J = gel(x, 2);
2391 56 : P = gel(x, 1);
2392 56 : one = mkintmod(gen_1, P);
2393 70 : } else if (tx == t_FFELT) {
2394 63 : J = FF_to_FpXQ_i(x);
2395 63 : if (degpol(J) > 0)
2396 7 : pari_err_DOMAIN("polmodular", "x", "not in prime subfield ", gen_0, x);
2397 56 : J = constant_coeff(J);
2398 56 : P = FF_p_i(x);
2399 56 : one = p_to_FF(P, 0);
2400 : } else {
2401 7 : pari_err_TYPE("polmodular", x);
2402 : }
2403 :
2404 112 : if (v < 0) v = 1;
2405 112 : res = Fp_polmodular_evalx(L, inv, J, P, v, compute_derivs);
2406 112 : res = gmul(res, one);
2407 112 : return gerepileupto(av, res);
2408 : }
2409 :
2410 : /**
2411 : * SECTION: Modular polynomials of level <= MAX_INTERNAL_MODPOLY_LEVEL.
2412 : */
2413 :
2414 : /*
2415 : * These functions return a vector of unique coefficients of classical
2416 : * modular polynomials \Phi_L(X, Y) of small level L. The number of
2417 : * such coefficients is (L + 1)(L + 2)/2 since \Phi is symmetric.
2418 : * (Note that we omit the (common) coefficient of X^{L + 1} and Y^{L +
2419 : * 1} since it is always 1.) See sympol_to_ZM() for how to interpret
2420 : * the coefficients, and use that function to get the corresponding
2421 : * full (desymmetrised) matrix of coefficients.
2422 : */
2423 :
2424 : /*
2425 : * Phi2, the modular polynomial of level 2:
2426 : *
2427 : * X^3
2428 : * + X^2 * (-Y^2 + 1488*Y - 162000)
2429 : * + X * (1488*Y^2 + 40773375*Y + 8748000000)
2430 : * + Y^3 - 162000*Y^2 + 8748000000*Y - 157464000000000
2431 : *
2432 : * [[3, 0, 1],
2433 : * [2, 2, -1],
2434 : * [2, 1, 1488],
2435 : * [2, 0, -162000],
2436 : * [1, 1, 40773375],
2437 : * [1, 0, 8748000000],
2438 : * [0, 0, -157464000000000]],
2439 : */
2440 :
2441 : static GEN
2442 14952 : phi2_ZV(void)
2443 : {
2444 14952 : GEN phi2 = cgetg(7, t_VEC);
2445 14952 : gel(phi2, 1) = uu32toi(36662, 1908994048);
2446 14952 : setsigne(gel(phi2, 1), -1);
2447 14952 : gel(phi2, 2) = uu32toi(2, 158065408);
2448 14952 : gel(phi2, 3) = stoi(40773375);
2449 14952 : gel(phi2, 4) = stoi(-162000);
2450 14952 : gel(phi2, 5) = stoi(1488);
2451 14952 : gel(phi2, 6) = gen_m1;
2452 14952 : return phi2;
2453 : }
2454 :
2455 :
2456 : /*
2457 : * L = 3
2458 : *
2459 : * [4, 0, 1],
2460 : * [3, 3, -1],
2461 : * [3, 2, 2232],
2462 : * [3, 1, -1069956],
2463 : * [3, 0, 36864000],
2464 : * [2, 2, 2587918086],
2465 : * [2, 1, 8900222976000],
2466 : * [2, 0, 452984832000000],
2467 : * [1, 1, -770845966336000000],
2468 : * [1, 0, 1855425871872000000000]
2469 : * [0, 0, 0]
2470 : *
2471 : * X^4
2472 : * + X^3 (-Y^3 + 2232*Y^2 - 1069956*Y + 36864000)
2473 : * + X^2 (2232*Y^3 + 2587918086*Y^2 + 8900222976000*Y + 452984832000000)
2474 : * + X (-1069956*Y^3 + 8900222976000*Y^2 - 770845966336000000*Y + 1855425871872000000000)
2475 : * + Y^4 + 36864000*Y^3 + 452984832000000*Y^2 + 1855425871872000000000*Y
2476 : *
2477 : * 1855425871872000000000 == 2^32 * (100 * 2^32 + 2503270400)
2478 : */
2479 : static GEN
2480 1057 : phi3_ZV(void)
2481 : {
2482 1057 : GEN phi3 = cgetg(11, t_VEC);
2483 1057 : pari_sp av = avma;
2484 1057 : gel(phi3, 1) = gen_0;
2485 1057 : gel(phi3, 2) = gerepileupto(av, shifti(uu32toi(100, 2503270400UL), 32));
2486 1057 : gel(phi3, 3) = uu32toi(179476562, 2147483648UL);
2487 1057 : setsigne(gel(phi3, 3), -1);
2488 1057 : gel(phi3, 4) = uu32toi(105468, 3221225472UL);
2489 1057 : gel(phi3, 5) = uu32toi(2072, 1050738688);
2490 1057 : gel(phi3, 6) = utoi(2587918086UL);
2491 1057 : gel(phi3, 7) = stoi(36864000);
2492 1057 : gel(phi3, 8) = stoi(-1069956);
2493 1057 : gel(phi3, 9) = stoi(2232);
2494 1057 : gel(phi3, 10) = gen_m1;
2495 1057 : return phi3;
2496 : }
2497 :
2498 :
2499 : static GEN
2500 1071 : phi5_ZV(void)
2501 : {
2502 1071 : GEN phi5 = cgetg(22, t_VEC);
2503 1071 : gel(phi5, 1) = mkintn(5, 0x18c2cc9cUL, 0x484382b2UL, 0xdc000000UL, 0x0UL, 0x0UL);
2504 1071 : gel(phi5, 2) = mkintn(5, 0x2638fUL, 0x2ff02690UL, 0x68026000UL, 0x0UL, 0x0UL);
2505 1071 : gel(phi5, 3) = mkintn(5, 0x308UL, 0xac9d9a4UL, 0xe0fdab12UL, 0xc0000000UL, 0x0UL);
2506 1071 : setsigne(gel(phi5, 3), -1);
2507 1071 : gel(phi5, 4) = mkintn(5, 0x13UL, 0xaae09f9dUL, 0x1b5ef872UL, 0x30000000UL, 0x0UL);
2508 1071 : gel(phi5, 5) = mkintn(4, 0x1b802fa9UL, 0x77ba0653UL, 0xd2f78000UL, 0x0UL);
2509 1071 : gel(phi5, 6) = mkintn(4, 0xfbfdUL, 0x278e4756UL, 0xdf08a7c4UL, 0x40000000UL);
2510 1071 : gel(phi5, 7) = mkintn(4, 0x35f922UL, 0x62ccea6fUL, 0x153d0000UL, 0x0UL);
2511 1071 : gel(phi5, 8) = mkintn(4, 0x97dUL, 0x29203fafUL, 0xc3036909UL, 0x80000000UL);
2512 1071 : setsigne(gel(phi5, 8), -1);
2513 1071 : gel(phi5, 9) = mkintn(3, 0x56e9e892UL, 0xd7781867UL, 0xf2ea0000UL);
2514 1071 : gel(phi5, 10) = mkintn(3, 0x5d6dUL, 0xe0a58f4eUL, 0x9ee68c14UL);
2515 1071 : setsigne(gel(phi5, 10), -1);
2516 1071 : gel(phi5, 11) = mkintn(3, 0x1100dUL, 0x85cea769UL, 0x40000000UL);
2517 1071 : gel(phi5, 12) = mkintn(3, 0x1b38UL, 0x43cf461fUL, 0x3a900000UL);
2518 1071 : gel(phi5, 13) = mkintn(3, 0x14UL, 0xc45a616eUL, 0x4801680fUL);
2519 1071 : gel(phi5, 14) = uu32toi(0x17f4350UL, 0x493ca3e0UL);
2520 1071 : gel(phi5, 15) = uu32toi(0x183UL, 0xe54ce1f8UL);
2521 1071 : gel(phi5, 16) = uu32toi(0x1c9UL, 0x18860000UL);
2522 1071 : gel(phi5, 17) = uu32toi(0x39UL, 0x6f7a2206UL);
2523 1071 : setsigne(gel(phi5, 17), -1);
2524 1071 : gel(phi5, 18) = stoi(2028551200);
2525 1071 : gel(phi5, 19) = stoi(-4550940);
2526 1071 : gel(phi5, 20) = stoi(3720);
2527 1071 : gel(phi5, 21) = gen_m1;
2528 1071 : return phi5;
2529 : }
2530 :
2531 : static GEN
2532 189 : phi5_f_ZV(void)
2533 : {
2534 189 : GEN phi = zerovec(21);
2535 189 : gel(phi, 3) = stoi(4);
2536 189 : gel(phi, 21) = gen_m1;
2537 189 : return phi;
2538 : }
2539 :
2540 : static GEN
2541 21 : phi3_f3_ZV(void)
2542 : {
2543 21 : GEN phi = zerovec(10);
2544 21 : gel(phi, 3) = stoi(8);
2545 21 : gel(phi, 10) = gen_m1;
2546 21 : return phi;
2547 : }
2548 :
2549 : static GEN
2550 63 : phi2_g2_ZV(void)
2551 : {
2552 63 : GEN phi = zerovec(6);
2553 63 : gel(phi, 1) = stoi(-54000);
2554 63 : gel(phi, 3) = stoi(495);
2555 63 : gel(phi, 6) = gen_m1;
2556 63 : return phi;
2557 : }
2558 :
2559 : static GEN
2560 70 : phi5_w2w3_ZV(void)
2561 : {
2562 70 : GEN phi = zerovec(21);
2563 70 : gel(phi, 3) = gen_m1;
2564 70 : gel(phi, 10) = stoi(5);
2565 70 : gel(phi, 21) = gen_m1;
2566 70 : return phi;
2567 : }
2568 :
2569 : static GEN
2570 98 : phi7_w2w5_ZV(void)
2571 : {
2572 98 : GEN phi = zerovec(36);
2573 98 : gel(phi, 3) = gen_m1;
2574 98 : gel(phi, 15) = stoi(56);
2575 98 : gel(phi, 19) = stoi(42);
2576 98 : gel(phi, 24) = stoi(21);
2577 98 : gel(phi, 30) = stoi(7);
2578 98 : gel(phi, 36) = gen_m1;
2579 98 : return phi;
2580 : }
2581 :
2582 : static GEN
2583 56 : phi5_w3w3_ZV(void)
2584 : {
2585 56 : GEN phi = zerovec(21);
2586 56 : gel(phi, 3) = stoi(9);
2587 56 : gel(phi, 6) = stoi(-15);
2588 56 : gel(phi, 15) = stoi(5);
2589 56 : gel(phi, 21) = gen_m1;
2590 56 : return phi;
2591 : }
2592 :
2593 : static GEN
2594 154 : phi3_w2w7_ZV(void)
2595 : {
2596 154 : GEN phi = zerovec(10);
2597 154 : gel(phi, 3) = gen_m1;
2598 154 : gel(phi, 6) = stoi(3);
2599 154 : gel(phi, 10) = gen_m1;
2600 154 : return phi;
2601 : }
2602 :
2603 : static GEN
2604 35 : phi2_w3w5_ZV(void)
2605 : {
2606 35 : GEN phi = zerovec(6);
2607 35 : gel(phi, 3) = gen_1;
2608 35 : gel(phi, 6) = gen_m1;
2609 35 : return phi;
2610 : }
2611 :
2612 : static GEN
2613 49 : phi5_w3w7_ZV(void)
2614 : {
2615 49 : GEN phi = zerovec(21);
2616 49 : gel(phi, 3) = gen_m1;
2617 49 : gel(phi, 6) = stoi(10);
2618 49 : gel(phi, 8) = stoi(5);
2619 49 : gel(phi, 10) = stoi(35);
2620 49 : gel(phi, 13) = stoi(20);
2621 49 : gel(phi, 15) = stoi(10);
2622 49 : gel(phi, 17) = stoi(5);
2623 49 : gel(phi, 19) = stoi(5);
2624 49 : gel(phi, 21) = gen_m1;
2625 49 : return phi;
2626 : }
2627 :
2628 : static GEN
2629 42 : phi3_w2w13_ZV(void)
2630 : {
2631 42 : GEN phi = zerovec(10);
2632 42 : gel(phi, 3) = gen_m1;
2633 42 : gel(phi, 6) = stoi(3);
2634 42 : gel(phi, 8) = stoi(3);
2635 42 : gel(phi, 10) = gen_m1;
2636 42 : return phi;
2637 : }
2638 :
2639 : static GEN
2640 21 : phi2_w3w3e2_ZV(void)
2641 : {
2642 21 : GEN phi = zerovec(6);
2643 21 : gel(phi, 3) = stoi(3);
2644 21 : gel(phi, 6) = gen_m1;
2645 21 : return phi;
2646 : }
2647 :
2648 : static GEN
2649 49 : phi2_w5w7_ZV(void)
2650 : {
2651 49 : GEN phi = zerovec(6);
2652 49 : gel(phi, 3) = gen_1;
2653 49 : gel(phi, 5) = gen_2;
2654 49 : gel(phi, 6) = gen_m1;
2655 49 : return phi;
2656 : }
2657 :
2658 : static GEN
2659 14 : phi2_w3w13_ZV(void)
2660 : {
2661 14 : GEN phi = zerovec(6);
2662 14 : gel(phi, 3) = gen_m1;
2663 14 : gel(phi, 5) = gen_2;
2664 14 : gel(phi, 6) = gen_m1;
2665 14 : return phi;
2666 : }
2667 :
2668 : INLINE long
2669 133 : modinv_parent(long inv)
2670 : {
2671 133 : switch (inv) {
2672 : case INV_F2:
2673 : case INV_F4:
2674 : case INV_F8:
2675 42 : return INV_F;
2676 : case INV_W2W3E2:
2677 21 : return INV_W2W3;
2678 : case INV_W2W5E2:
2679 21 : return INV_W2W5;
2680 : case INV_W2W7E2:
2681 49 : return INV_W2W7;
2682 : case INV_W3W3E2:
2683 0 : return INV_W3W3;
2684 : default:
2685 0 : pari_err_BUG("modinv_parent");
2686 : }
2687 0 : return -1;
2688 : }
2689 :
2690 : /* TODO: Think of a better name than "parent power"; sheesh. */
2691 : INLINE long
2692 133 : modinv_parent_power(long inv)
2693 : {
2694 133 : switch (inv) {
2695 : case INV_F4:
2696 14 : return 4;
2697 : case INV_F8:
2698 14 : return 8;
2699 : case INV_F2:
2700 : case INV_W2W3E2:
2701 : case INV_W2W5E2:
2702 : case INV_W2W7E2:
2703 : case INV_W3W3E2:
2704 105 : return 2;
2705 : default:
2706 0 : pari_err_BUG("modinv_parent_power");
2707 : }
2708 0 : return -1;
2709 : }
2710 :
2711 : static GEN
2712 133 : polmodular0_powerup_ZM(long L, long inv, GEN *db)
2713 : {
2714 133 : pari_sp ltop = avma, av;
2715 : long s, D, nprimes, N;
2716 : GEN mp, pol, P, H;
2717 :
2718 133 : long parent = modinv_parent(inv);
2719 133 : long e = modinv_parent_power(inv);
2720 :
2721 : modpoly_disc_info Ds[MODPOLY_MAX_DCNT];
2722 : /* FIXME: We throw away the table of fundamental discriminants here. */
2723 133 : long nDs = discriminant_with_classno_at_least(Ds, L, inv, IGNORE_SPARSE_FACTOR);
2724 133 : if (nDs != 1) pari_err_BUG("polmodular0_powerup_ZM");
2725 133 : D = Ds[0].D1;
2726 133 : nprimes = Ds[0].nprimes + 1;
2727 133 : mp = polmodular0_ZM(L, parent, NULL, NULL, 0, db);
2728 133 : H = polclass0(D, parent, 0, db);
2729 :
2730 133 : N = L + 2;
2731 133 : if (degpol(H) < N) pari_err_BUG("polmodular0_powerup_ZM");
2732 :
2733 133 : av = avma;
2734 133 : pol = ZM_init_CRT(zero_Flm_copy(N, L + 2), 1);
2735 133 : P = gen_1;
2736 448 : for (s = 1; s < nprimes; ++s) {
2737 : pari_sp av1, av2;
2738 315 : ulong p = Ds[0].primes[s-1], pi = get_Fl_red(p);
2739 : long i;
2740 : GEN Hrts, js, Hp, Phip, coeff_mat, phi_modp;
2741 :
2742 315 : phi_modp = zero_Flm_copy(N, L + 2);
2743 315 : av1 = avma;
2744 315 : Hp = ZX_to_Flx(H, p);
2745 315 : Hrts = Flx_roots(Hp, p);
2746 315 : if (glength(Hrts) < N) pari_err_BUG("polmodular0_powerup_ZM");
2747 315 : js = cgetg(N + 1, t_VECSMALL);
2748 2450 : for (i = 1; i <= N; ++i)
2749 2135 : uel(js, i) = Fl_powu_pre(uel(Hrts, i), e, p, pi);
2750 :
2751 315 : Phip = ZM_to_Flm(mp, p);
2752 315 : coeff_mat = zero_Flm_copy(N, L + 2);
2753 315 : av2 = avma;
2754 2450 : for (i = 1; i <= N; ++i) {
2755 : long k;
2756 : GEN phi_at_ji, mprts;
2757 :
2758 2135 : phi_at_ji = Flm_Fl_polmodular_evalx(Phip, L, uel(Hrts, i), p, pi);
2759 2135 : mprts = Flx_roots(phi_at_ji, p);
2760 2135 : if (lg(mprts) != L + 2) pari_err_BUG("polmodular0_powerup_ZM");
2761 :
2762 2135 : Flv_powu_inplace_pre(mprts, e, p, pi);
2763 2135 : phi_at_ji = Flv_roots_to_pol(mprts, p, 0);
2764 :
2765 17234 : for (k = 1; k <= L + 2; ++k)
2766 15099 : ucoeff(coeff_mat, i, k) = uel(phi_at_ji, k + 1);
2767 2135 : avma = av2;
2768 : }
2769 :
2770 315 : interpolate_coeffs(phi_modp, p, js, coeff_mat);
2771 315 : avma = av1;
2772 :
2773 315 : (void) ZM_incremental_CRT(&pol, phi_modp, &P, p);
2774 315 : if (gc_needed(av, 2))
2775 0 : gerepileall(av, 2, &pol, &P);
2776 : }
2777 :
2778 133 : return gerepileupto(ltop, pol);
2779 : }
2780 :
2781 : /* Returns the modular polynomial with the smallest level for the
2782 : * given invariant, except if inv is INV_J, in which case return the
2783 : * modular polynomial of level L in {2,3,5}. NULL is returned if the
2784 : * modular polynomial can be calculated using polmodular0_powerup_ZM. */
2785 : INLINE GEN
2786 18074 : internal_db(long L, long inv)
2787 : {
2788 18074 : switch (inv) {
2789 : case INV_J: {
2790 17080 : switch (L) {
2791 14952 : case 2: return phi2_ZV();
2792 1057 : case 3: return phi3_ZV();
2793 1071 : case 5: return phi5_ZV();
2794 0 : default: break;
2795 : }
2796 : }
2797 189 : case INV_F: return phi5_f_ZV();
2798 14 : case INV_F2: return NULL;
2799 21 : case INV_F3: return phi3_f3_ZV();
2800 14 : case INV_F4: return NULL;
2801 63 : case INV_G2: return phi2_g2_ZV();
2802 70 : case INV_W2W3: return phi5_w2w3_ZV();
2803 14 : case INV_F8: return NULL;
2804 56 : case INV_W3W3: return phi5_w3w3_ZV();
2805 98 : case INV_W2W5: return phi7_w2w5_ZV();
2806 154 : case INV_W2W7: return phi3_w2w7_ZV();
2807 35 : case INV_W3W5: return phi2_w3w5_ZV();
2808 49 : case INV_W3W7: return phi5_w3w7_ZV();
2809 21 : case INV_W2W3E2: return NULL;
2810 21 : case INV_W2W5E2: return NULL;
2811 42 : case INV_W2W13: return phi3_w2w13_ZV();
2812 49 : case INV_W2W7E2: return NULL;
2813 21 : case INV_W3W3E2: return phi2_w3w3e2_ZV();
2814 49 : case INV_W5W7: return phi2_w5w7_ZV();
2815 14 : case INV_W3W13: return phi2_w3w13_ZV();
2816 : }
2817 0 : pari_err_BUG("internal_db");
2818 0 : return NULL;
2819 : }
2820 :
2821 : /* NB: Should only be called if L <= modinv_max_internal_level(inv) */
2822 : static GEN
2823 18074 : polmodular_small_ZM(long L, long inv, GEN *db)
2824 : {
2825 18074 : GEN f = internal_db(L, inv);
2826 18074 : if ( ! f)
2827 133 : return polmodular0_powerup_ZM(L, inv, db);
2828 17941 : return sympol_to_ZM(f, L);
2829 : }
2830 :
2831 :
2832 : /* Each function phi_w?w?_j() returns a vector V containing two
2833 : * vectors u and v, and a scalar k, which together represent the
2834 : * bivariate polnomial
2835 : *
2836 : * phi(X, Y) = \sum_i u[i] X^i + Y \sum_i v[i] X^i + Y^2 X^k
2837 : */
2838 : static GEN
2839 1449 : phi_w2w3_j(void)
2840 : {
2841 : GEN phi, phi0, phi1;
2842 1449 : phi = cgetg(4, t_VEC);
2843 :
2844 1449 : phi0 = cgetg(14, t_VEC);
2845 1449 : gel(phi0, 1) = gen_1;
2846 1449 : gel(phi0, 2) = utoi(0x3cUL); setsigne(gel(phi0, 2), -1);
2847 1449 : gel(phi0, 3) = utoi(0x702UL);
2848 1449 : gel(phi0, 4) = utoi(0x797cUL); setsigne(gel(phi0, 4), -1);
2849 1449 : gel(phi0, 5) = utoi(0x5046fUL);
2850 1449 : gel(phi0, 6) = utoi(0x1be0b8UL); setsigne(gel(phi0, 6), -1);
2851 1449 : gel(phi0, 7) = utoi(0x28ef9cUL);
2852 1449 : gel(phi0, 8) = utoi(0x15e2968UL);
2853 1449 : gel(phi0, 9) = utoi(0x1b8136fUL);
2854 1449 : gel(phi0, 10) = utoi(0xa67674UL);
2855 1449 : gel(phi0, 11) = utoi(0x23982UL);
2856 1449 : gel(phi0, 12) = utoi(0x294UL);
2857 1449 : gel(phi0, 13) = gen_1;
2858 :
2859 1449 : phi1 = cgetg(13, t_VEC);
2860 1449 : gel(phi1, 1) = gen_0;
2861 1449 : gel(phi1, 2) = gen_0;
2862 1449 : gel(phi1, 3) = gen_m1;
2863 1449 : gel(phi1, 4) = utoi(0x23UL);
2864 1449 : gel(phi1, 5) = utoi(0xaeUL); setsigne(gel(phi1, 5), -1);
2865 1449 : gel(phi1, 6) = utoi(0x5b8UL); setsigne(gel(phi1, 6), -1);
2866 1449 : gel(phi1, 7) = utoi(0x12d7UL);
2867 1449 : gel(phi1, 8) = utoi(0x7c86UL); setsigne(gel(phi1, 8), -1);
2868 1449 : gel(phi1, 9) = utoi(0x37c8UL);
2869 1449 : gel(phi1, 10) = utoi(0x69cUL); setsigne(gel(phi1, 10), -1);
2870 1449 : gel(phi1, 11) = utoi(0x48UL);
2871 1449 : gel(phi1, 12) = gen_m1;
2872 :
2873 1449 : gel(phi, 1) = phi0;
2874 1449 : gel(phi, 2) = phi1;
2875 1449 : gel(phi, 3) = stoi(5);
2876 :
2877 1449 : return phi;
2878 : }
2879 :
2880 : static GEN
2881 2495 : phi_w3w3_j(void)
2882 : {
2883 : GEN phi, phi0, phi1;
2884 2495 : phi = cgetg(4, t_VEC);
2885 :
2886 2495 : phi0 = cgetg(14, t_VEC);
2887 2495 : gel(phi0, 1) = utoi(0x2d9UL);
2888 2495 : gel(phi0, 2) = utoi(0x4fbcUL);
2889 2495 : gel(phi0, 3) = utoi(0x5828aUL);
2890 2495 : gel(phi0, 4) = utoi(0x3a7a3cUL);
2891 2495 : gel(phi0, 5) = utoi(0x1bd8edfUL);
2892 2495 : gel(phi0, 6) = utoi(0x8348838UL);
2893 2494 : gel(phi0, 7) = utoi(0x1983f8acUL);
2894 2495 : gel(phi0, 8) = utoi(0x14e4e098UL);
2895 2495 : gel(phi0, 9) = utoi(0x69ed1a7UL);
2896 2495 : gel(phi0, 10) = utoi(0xc3828cUL);
2897 2495 : gel(phi0, 11) = utoi(0x2696aUL);
2898 2495 : gel(phi0, 12) = utoi(0x2acUL);
2899 2495 : gel(phi0, 13) = gen_1;
2900 :
2901 2495 : phi1 = cgetg(13, t_VEC);
2902 2495 : gel(phi1, 1) = gen_0;
2903 2495 : gel(phi1, 2) = utoi(0x1bUL); setsigne(gel(phi1, 2), -1);
2904 2495 : gel(phi1, 3) = utoi(0x5d6UL); setsigne(gel(phi1, 3), -1);
2905 2495 : gel(phi1, 4) = utoi(0x1c7bUL); setsigne(gel(phi1, 4), -1);
2906 2495 : gel(phi1, 5) = utoi(0x7980UL);
2907 2495 : gel(phi1, 6) = utoi(0x12168UL);
2908 2495 : gel(phi1, 7) = utoi(0x3528UL); setsigne(gel(phi1, 7), -1);
2909 2495 : gel(phi1, 8) = utoi(0x6174UL); setsigne(gel(phi1, 8), -1);
2910 2495 : gel(phi1, 9) = utoi(0x2208UL);
2911 2495 : gel(phi1, 10) = utoi(0x41dUL); setsigne(gel(phi1, 10), -1);
2912 2495 : gel(phi1, 11) = utoi(0x36UL);
2913 2495 : gel(phi1, 12) = gen_m1;
2914 :
2915 2495 : gel(phi, 1) = phi0;
2916 2495 : gel(phi, 2) = phi1;
2917 2495 : gel(phi, 3) = utoi(2);
2918 :
2919 2495 : return phi;
2920 : }
2921 :
2922 : static GEN
2923 3500 : phi_w2w5_j(void)
2924 : {
2925 : GEN phi, phi0, phi1;
2926 3500 : phi = cgetg(4, t_VEC);
2927 :
2928 3500 : phi0 = cgetg(20, t_VEC);
2929 3500 : gel(phi0, 1) = gen_1;
2930 3500 : gel(phi0, 2) = utoi(0x2aUL); setsigne(gel(phi0, 2), -1);
2931 3500 : gel(phi0, 3) = utoi(0x549UL);
2932 3500 : gel(phi0, 4) = utoi(0x6530UL); setsigne(gel(phi0, 4), -1);
2933 3500 : gel(phi0, 5) = utoi(0x60504UL);
2934 3500 : gel(phi0, 6) = utoi(0x3cbbc8UL); setsigne(gel(phi0, 6), -1);
2935 3500 : gel(phi0, 7) = utoi(0x1d1ee74UL);
2936 3500 : gel(phi0, 8) = utoi(0x7ef9ab0UL); setsigne(gel(phi0, 8), -1);
2937 3500 : gel(phi0, 9) = utoi(0x12b888beUL);
2938 3500 : gel(phi0, 10) = utoi(0x15fa174cUL); setsigne(gel(phi0, 10), -1);
2939 3500 : gel(phi0, 11) = utoi(0x615d9feUL);
2940 3500 : gel(phi0, 12) = utoi(0xbeca070UL);
2941 3500 : gel(phi0, 13) = utoi(0x88de74cUL); setsigne(gel(phi0, 13), -1);
2942 3500 : gel(phi0, 14) = utoi(0x2b3a268UL); setsigne(gel(phi0, 14), -1);
2943 3500 : gel(phi0, 15) = utoi(0x24b3244UL);
2944 3500 : gel(phi0, 16) = utoi(0xb56270UL);
2945 3500 : gel(phi0, 17) = utoi(0x25989UL);
2946 3500 : gel(phi0, 18) = utoi(0x2a6UL);
2947 3500 : gel(phi0, 19) = gen_1;
2948 :
2949 3500 : phi1 = cgetg(19, t_VEC);
2950 3500 : gel(phi1, 1) = gen_0;
2951 3500 : gel(phi1, 2) = gen_0;
2952 3500 : gel(phi1, 3) = gen_m1;
2953 3500 : gel(phi1, 4) = utoi(0x1eUL);
2954 3500 : gel(phi1, 5) = utoi(0xffUL); setsigne(gel(phi1, 5), -1);
2955 3500 : gel(phi1, 6) = utoi(0x243UL);
2956 3500 : gel(phi1, 7) = utoi(0xf3UL); setsigne(gel(phi1, 7), -1);
2957 3500 : gel(phi1, 8) = utoi(0x5c4UL); setsigne(gel(phi1, 8), -1);
2958 3500 : gel(phi1, 9) = utoi(0x107bUL);
2959 3500 : gel(phi1, 10) = utoi(0x11b2fUL); setsigne(gel(phi1, 10), -1);
2960 3500 : gel(phi1, 11) = utoi(0x48fa8UL);
2961 3500 : gel(phi1, 12) = utoi(0x6ff7cUL); setsigne(gel(phi1, 12), -1);
2962 3500 : gel(phi1, 13) = utoi(0x4bf48UL);
2963 3500 : gel(phi1, 14) = utoi(0x187efUL); setsigne(gel(phi1, 14), -1);
2964 3500 : gel(phi1, 15) = utoi(0x404cUL);
2965 3500 : gel(phi1, 16) = utoi(0x582UL); setsigne(gel(phi1, 16), -1);
2966 3500 : gel(phi1, 17) = utoi(0x3cUL);
2967 3500 : gel(phi1, 18) = gen_m1;
2968 :
2969 3500 : gel(phi, 1) = phi0;
2970 3500 : gel(phi, 2) = phi1;
2971 3500 : gel(phi, 3) = utoi(7);
2972 :
2973 3500 : return phi;
2974 : }
2975 :
2976 : static GEN
2977 5665 : phi_w2w7_j(void)
2978 : {
2979 : GEN phi, phi0, phi1;
2980 5665 : phi = cgetg(4, t_VEC);
2981 :
2982 5665 : phi0 = cgetg(26, t_VEC);
2983 5665 : gel(phi0, 1) = gen_1;
2984 5665 : gel(phi0, 2) = utoi(0x24UL); setsigne(gel(phi0, 2), -1);
2985 5665 : gel(phi0, 3) = utoi(0x4ceUL);
2986 5665 : gel(phi0, 4) = utoi(0x5d60UL); setsigne(gel(phi0, 4), -1);
2987 5665 : gel(phi0, 5) = utoi(0x62b05UL);
2988 5665 : gel(phi0, 6) = utoi(0x47be78UL); setsigne(gel(phi0, 6), -1);
2989 5665 : gel(phi0, 7) = utoi(0x2a3880aUL);
2990 5665 : gel(phi0, 8) = utoi(0x114bccf4UL); setsigne(gel(phi0, 8), -1);
2991 5665 : gel(phi0, 9) = utoi(0x4b95e79aUL);
2992 5665 : gel(phi0, 10) = utoi(0xe2cfee1cUL); setsigne(gel(phi0, 10), -1);
2993 5665 : gel(phi0, 11) = uu32toi(0x1UL, 0xe43d1126UL);
2994 5665 : gel(phi0, 12) = uu32toi(0x2UL, 0xf04dc6f8UL); setsigne(gel(phi0, 12), -1);
2995 5665 : gel(phi0, 13) = uu32toi(0x3UL, 0x5384987dUL);
2996 5665 : gel(phi0, 14) = uu32toi(0x2UL, 0xa5ccbe18UL); setsigne(gel(phi0, 14), -1);
2997 5665 : gel(phi0, 15) = uu32toi(0x1UL, 0x4c52c8a6UL);
2998 5665 : gel(phi0, 16) = utoi(0x2643fdecUL); setsigne(gel(phi0, 16), -1);
2999 5665 : gel(phi0, 17) = utoi(0x49f5ab66UL); setsigne(gel(phi0, 17), -1);
3000 5665 : gel(phi0, 18) = utoi(0x33074d3cUL);
3001 5665 : gel(phi0, 19) = utoi(0x6a3e376UL); setsigne(gel(phi0, 19), -1);
3002 5665 : gel(phi0, 20) = utoi(0x675aa58UL); setsigne(gel(phi0, 20), -1);
3003 5665 : gel(phi0, 21) = utoi(0x2674005UL);
3004 5665 : gel(phi0, 22) = utoi(0xba5be0UL);
3005 5665 : gel(phi0, 23) = utoi(0x2644eUL);
3006 5665 : gel(phi0, 24) = utoi(0x2acUL);
3007 5665 : gel(phi0, 25) = gen_1;
3008 :
3009 5665 : phi1 = cgetg(25, t_VEC);
3010 5665 : gel(phi1, 1) = gen_0;
3011 5665 : gel(phi1, 2) = gen_0;
3012 5665 : gel(phi1, 3) = gen_m1;
3013 5665 : gel(phi1, 4) = utoi(0x1cUL);
3014 5665 : gel(phi1, 5) = utoi(0x10aUL); setsigne(gel(phi1, 5), -1);
3015 5665 : gel(phi1, 6) = utoi(0x3f0UL);
3016 5665 : gel(phi1, 7) = utoi(0x5d3UL); setsigne(gel(phi1, 7), -1);
3017 5665 : gel(phi1, 8) = utoi(0x3efUL);
3018 5665 : gel(phi1, 9) = utoi(0x102UL); setsigne(gel(phi1, 9), -1);
3019 5665 : gel(phi1, 10) = utoi(0x5c8UL); setsigne(gel(phi1, 10), -1);
3020 5665 : gel(phi1, 11) = utoi(0x102fUL);
3021 5665 : gel(phi1, 12) = utoi(0x13f8aUL); setsigne(gel(phi1, 12), -1);
3022 5665 : gel(phi1, 13) = utoi(0x86538UL);
3023 5665 : gel(phi1, 14) = utoi(0x1bbd10UL); setsigne(gel(phi1, 14), -1);
3024 5665 : gel(phi1, 15) = utoi(0x3614e8UL);
3025 5665 : gel(phi1, 16) = utoi(0x42f793UL); setsigne(gel(phi1, 16), -1);
3026 5665 : gel(phi1, 17) = utoi(0x364698UL);
3027 5665 : gel(phi1, 18) = utoi(0x1c7a10UL); setsigne(gel(phi1, 18), -1);
3028 5665 : gel(phi1, 19) = utoi(0x97cc8UL);
3029 5665 : gel(phi1, 20) = utoi(0x1fc8aUL); setsigne(gel(phi1, 20), -1);
3030 5665 : gel(phi1, 21) = utoi(0x4210UL);
3031 5665 : gel(phi1, 22) = utoi(0x524UL); setsigne(gel(phi1, 22), -1);
3032 5665 : gel(phi1, 23) = utoi(0x38UL);
3033 5665 : gel(phi1, 24) = gen_m1;
3034 :
3035 5665 : gel(phi, 1) = phi0;
3036 5665 : gel(phi, 2) = phi1;
3037 5665 : gel(phi, 3) = utoi(9);
3038 :
3039 5665 : return phi;
3040 : }
3041 :
3042 : static GEN
3043 1876 : phi_w2w13_j(void)
3044 : {
3045 : GEN phi, phi0, phi1;
3046 1876 : phi = cgetg(4, t_VEC);
3047 :
3048 1876 : phi0 = cgetg(44, t_VEC);
3049 1876 : gel(phi0, 1) = gen_1;
3050 1876 : gel(phi0, 2) = utoi(0x1eUL); setsigne(gel(phi0, 2), -1);
3051 1876 : gel(phi0, 3) = utoi(0x45fUL);
3052 1876 : gel(phi0, 4) = utoi(0x5590UL); setsigne(gel(phi0, 4), -1);
3053 1876 : gel(phi0, 5) = utoi(0x64407UL);
3054 1876 : gel(phi0, 6) = utoi(0x53a792UL); setsigne(gel(phi0, 6), -1);
3055 1876 : gel(phi0, 7) = utoi(0x3b21af3UL);
3056 1876 : gel(phi0, 8) = utoi(0x20d056d0UL); setsigne(gel(phi0, 8), -1);
3057 1876 : gel(phi0, 9) = utoi(0xe02db4a6UL);
3058 1876 : gel(phi0, 10) = uu32toi(0x4UL, 0xb23400b0UL); setsigne(gel(phi0, 10), -1);
3059 1876 : gel(phi0, 11) = uu32toi(0x14UL, 0x57fbb906UL);
3060 1876 : gel(phi0, 12) = uu32toi(0x49UL, 0xcf80c00UL); setsigne(gel(phi0, 12), -1);
3061 1876 : gel(phi0, 13) = uu32toi(0xdeUL, 0x84ff421UL);
3062 1876 : gel(phi0, 14) = uu32toi(0x244UL, 0xc500c156UL); setsigne(gel(phi0, 14), -1);
3063 1876 : gel(phi0, 15) = uu32toi(0x52cUL, 0x79162979UL);
3064 1876 : gel(phi0, 16) = uu32toi(0xa64UL, 0x8edc5650UL); setsigne(gel(phi0, 16), -1);
3065 1876 : gel(phi0, 17) = uu32toi(0x1289UL, 0x4225bb41UL);
3066 1876 : gel(phi0, 18) = uu32toi(0x1d89UL, 0x2a15229aUL); setsigne(gel(phi0, 18), -1);
3067 1876 : gel(phi0, 19) = uu32toi(0x2a3eUL, 0x4539f1ebUL);
3068 1876 : gel(phi0, 20) = uu32toi(0x366aUL, 0xa5ea1130UL); setsigne(gel(phi0, 20), -1);
3069 1876 : gel(phi0, 21) = uu32toi(0x3f47UL, 0xa19fecb4UL);
3070 1876 : gel(phi0, 22) = uu32toi(0x4282UL, 0x91a3c4a0UL); setsigne(gel(phi0, 22), -1);
3071 1876 : gel(phi0, 23) = uu32toi(0x3f30UL, 0xbaa305b4UL);
3072 1876 : gel(phi0, 24) = uu32toi(0x3635UL, 0xd11c2530UL); setsigne(gel(phi0, 24), -1);
3073 1876 : gel(phi0, 25) = uu32toi(0x29e2UL, 0x89df27ebUL);
3074 1876 : gel(phi0, 26) = uu32toi(0x1d03UL, 0x6509d48aUL); setsigne(gel(phi0, 26), -1);
3075 1876 : gel(phi0, 27) = uu32toi(0x11e2UL, 0x272cc601UL);
3076 1876 : gel(phi0, 28) = uu32toi(0x9b0UL, 0xacd58ff0UL); setsigne(gel(phi0, 28), -1);
3077 1876 : gel(phi0, 29) = uu32toi(0x485UL, 0x608d7db9UL);
3078 1876 : gel(phi0, 30) = uu32toi(0x1bfUL, 0xa941546UL); setsigne(gel(phi0, 30), -1);
3079 1876 : gel(phi0, 31) = uu32toi(0x82UL, 0x56e48b21UL);
3080 1876 : gel(phi0, 32) = uu32toi(0x13UL, 0xc36b2340UL); setsigne(gel(phi0, 32), -1);
3081 1876 : gel(phi0, 33) = uu32toi(0x5UL, 0x6637257aUL); setsigne(gel(phi0, 33), -1);
3082 1876 : gel(phi0, 34) = uu32toi(0x5UL, 0x40f70bd0UL);
3083 1876 : gel(phi0, 35) = uu32toi(0x1UL, 0xf70842daUL); setsigne(gel(phi0, 35), -1);
3084 1876 : gel(phi0, 36) = utoi(0x53eea5f0UL);
3085 1876 : gel(phi0, 37) = utoi(0xda17bf3UL);
3086 1876 : gel(phi0, 38) = utoi(0xaf246c2UL); setsigne(gel(phi0, 38), -1);
3087 1876 : gel(phi0, 39) = utoi(0x278f847UL);
3088 1876 : gel(phi0, 40) = utoi(0xbf5550UL);
3089 1876 : gel(phi0, 41) = utoi(0x26f1fUL);
3090 1876 : gel(phi0, 42) = utoi(0x2b2UL);
3091 1876 : gel(phi0, 43) = gen_1;
3092 :
3093 1876 : phi1 = cgetg(43, t_VEC);
3094 1876 : gel(phi1, 1) = gen_0;
3095 1876 : gel(phi1, 2) = gen_0;
3096 1876 : gel(phi1, 3) = gen_m1;
3097 1876 : gel(phi1, 4) = utoi(0x1aUL);
3098 1876 : gel(phi1, 5) = utoi(0x111UL); setsigne(gel(phi1, 5), -1);
3099 1876 : gel(phi1, 6) = utoi(0x5e4UL);
3100 1876 : gel(phi1, 7) = utoi(0x1318UL); setsigne(gel(phi1, 7), -1);
3101 1876 : gel(phi1, 8) = utoi(0x2804UL);
3102 1876 : gel(phi1, 9) = utoi(0x3cd6UL); setsigne(gel(phi1, 9), -1);
3103 1876 : gel(phi1, 10) = utoi(0x467cUL);
3104 1876 : gel(phi1, 11) = utoi(0x3cd6UL); setsigne(gel(phi1, 11), -1);
3105 1876 : gel(phi1, 12) = utoi(0x2804UL);
3106 1876 : gel(phi1, 13) = utoi(0x1318UL); setsigne(gel(phi1, 13), -1);
3107 1876 : gel(phi1, 14) = utoi(0x5e3UL);
3108 1876 : gel(phi1, 15) = utoi(0x10dUL); setsigne(gel(phi1, 15), -1);
3109 1876 : gel(phi1, 16) = utoi(0x5ccUL); setsigne(gel(phi1, 16), -1);
3110 1876 : gel(phi1, 17) = utoi(0x100bUL);
3111 1876 : gel(phi1, 18) = utoi(0x160e1UL); setsigne(gel(phi1, 18), -1);
3112 1876 : gel(phi1, 19) = utoi(0xd2cb0UL);
3113 1876 : gel(phi1, 20) = utoi(0x4c85fcUL); setsigne(gel(phi1, 20), -1);
3114 1876 : gel(phi1, 21) = utoi(0x137cb98UL);
3115 1876 : gel(phi1, 22) = utoi(0x3c75568UL); setsigne(gel(phi1, 22), -1);
3116 1876 : gel(phi1, 23) = utoi(0x95c69c8UL);
3117 1876 : gel(phi1, 24) = utoi(0x131557bcUL); setsigne(gel(phi1, 24), -1);
3118 1876 : gel(phi1, 25) = utoi(0x20aacfd0UL);
3119 1876 : gel(phi1, 26) = utoi(0x2f9164e6UL); setsigne(gel(phi1, 26), -1);
3120 1876 : gel(phi1, 27) = utoi(0x3b6a5e40UL);
3121 1876 : gel(phi1, 28) = utoi(0x3ff54344UL); setsigne(gel(phi1, 28), -1);
3122 1876 : gel(phi1, 29) = utoi(0x3b6a9140UL);
3123 1876 : gel(phi1, 30) = utoi(0x2f927fa6UL); setsigne(gel(phi1, 30), -1);
3124 1876 : gel(phi1, 31) = utoi(0x20ae6450UL);
3125 1876 : gel(phi1, 32) = utoi(0x131cd87cUL); setsigne(gel(phi1, 32), -1);
3126 1876 : gel(phi1, 33) = utoi(0x967d1e8UL);
3127 1876 : gel(phi1, 34) = utoi(0x3d48ca8UL); setsigne(gel(phi1, 34), -1);
3128 1876 : gel(phi1, 35) = utoi(0x14333b8UL);
3129 1876 : gel(phi1, 36) = utoi(0x5406bcUL); setsigne(gel(phi1, 36), -1);
3130 1876 : gel(phi1, 37) = utoi(0x10c130UL);
3131 1876 : gel(phi1, 38) = utoi(0x27ba1UL); setsigne(gel(phi1, 38), -1);
3132 1876 : gel(phi1, 39) = utoi(0x433cUL);
3133 1876 : gel(phi1, 40) = utoi(0x4c6UL); setsigne(gel(phi1, 40), -1);
3134 1876 : gel(phi1, 41) = utoi(0x34UL);
3135 1876 : gel(phi1, 42) = gen_m1;
3136 :
3137 1876 : gel(phi, 1) = phi0;
3138 1876 : gel(phi, 2) = phi1;
3139 1876 : gel(phi, 3) = utoi(15);
3140 :
3141 1876 : return phi;
3142 : }
3143 :
3144 : static GEN
3145 1107 : phi_w3w5_j(void)
3146 : {
3147 : GEN phi, phi0, phi1;
3148 1107 : phi = cgetg(4, t_VEC);
3149 :
3150 1107 : phi0 = cgetg(26, t_VEC);
3151 1107 : gel(phi0, 1) = gen_1;
3152 1107 : gel(phi0, 2) = utoi(0x18UL);
3153 1107 : gel(phi0, 3) = utoi(0xb4UL);
3154 1107 : gel(phi0, 4) = utoi(0x178UL); setsigne(gel(phi0, 4), -1);
3155 1107 : gel(phi0, 5) = utoi(0x2d7eUL); setsigne(gel(phi0, 5), -1);
3156 1107 : gel(phi0, 6) = utoi(0x89b8UL); setsigne(gel(phi0, 6), -1);
3157 1107 : gel(phi0, 7) = utoi(0x35c24UL);
3158 1107 : gel(phi0, 8) = utoi(0x128a18UL);
3159 1107 : gel(phi0, 9) = utoi(0x12a911UL); setsigne(gel(phi0, 9), -1);
3160 1107 : gel(phi0, 10) = utoi(0xcc0190UL); setsigne(gel(phi0, 10), -1);
3161 1107 : gel(phi0, 11) = utoi(0x94368UL);
3162 1107 : gel(phi0, 12) = utoi(0x1439d0UL);
3163 1107 : gel(phi0, 13) = utoi(0x96f931cUL);
3164 1107 : gel(phi0, 14) = utoi(0x1f59ff0UL); setsigne(gel(phi0, 14), -1);
3165 1107 : gel(phi0, 15) = utoi(0x20e7e8UL);
3166 1107 : gel(phi0, 16) = utoi(0x25fdf150UL); setsigne(gel(phi0, 16), -1);
3167 1107 : gel(phi0, 17) = utoi(0x7091511UL); setsigne(gel(phi0, 17), -1);
3168 1107 : gel(phi0, 18) = utoi(0x1ef52f8UL);
3169 1107 : gel(phi0, 19) = utoi(0x341f2de4UL);
3170 1107 : gel(phi0, 20) = utoi(0x25d72c28UL);
3171 1107 : gel(phi0, 21) = utoi(0x95d2082UL);
3172 1107 : gel(phi0, 22) = utoi(0xd2d828UL);
3173 1107 : gel(phi0, 23) = utoi(0x281f4UL);
3174 1107 : gel(phi0, 24) = utoi(0x2b8UL);
3175 1107 : gel(phi0, 25) = gen_1;
3176 :
3177 1107 : phi1 = cgetg(25, t_VEC);
3178 1107 : gel(phi1, 1) = gen_0;
3179 1107 : gel(phi1, 2) = gen_0;
3180 1107 : gel(phi1, 3) = gen_0;
3181 1107 : gel(phi1, 4) = gen_1;
3182 1107 : gel(phi1, 5) = utoi(0xfUL);
3183 1107 : gel(phi1, 6) = utoi(0x2eUL);
3184 1107 : gel(phi1, 7) = utoi(0x1fUL); setsigne(gel(phi1, 7), -1);
3185 1107 : gel(phi1, 8) = utoi(0x2dUL); setsigne(gel(phi1, 8), -1);
3186 1107 : gel(phi1, 9) = utoi(0x5caUL); setsigne(gel(phi1, 9), -1);
3187 1107 : gel(phi1, 10) = utoi(0x358UL); setsigne(gel(phi1, 10), -1);
3188 1107 : gel(phi1, 11) = utoi(0x2f1cUL);
3189 1107 : gel(phi1, 12) = utoi(0xd8eaUL);
3190 1107 : gel(phi1, 13) = utoi(0x38c70UL); setsigne(gel(phi1, 13), -1);
3191 1107 : gel(phi1, 14) = utoi(0x1a964UL); setsigne(gel(phi1, 14), -1);
3192 1107 : gel(phi1, 15) = utoi(0x93512UL);
3193 1107 : gel(phi1, 16) = utoi(0x58f2UL); setsigne(gel(phi1, 16), -1);
3194 1107 : gel(phi1, 17) = utoi(0x5af1eUL); setsigne(gel(phi1, 17), -1);
3195 1107 : gel(phi1, 18) = utoi(0x1afb8UL);
3196 1107 : gel(phi1, 19) = utoi(0xc084UL);
3197 1107 : gel(phi1, 20) = utoi(0x7fcbUL); setsigne(gel(phi1, 20), -1);
3198 1107 : gel(phi1, 21) = utoi(0x1c89UL);
3199 1107 : gel(phi1, 22) = utoi(0x32aUL); setsigne(gel(phi1, 22), -1);
3200 1107 : gel(phi1, 23) = utoi(0x2dUL);
3201 1107 : gel(phi1, 24) = gen_m1;
3202 :
3203 1107 : gel(phi, 1) = phi0;
3204 1107 : gel(phi, 2) = phi1;
3205 1107 : gel(phi, 3) = utoi(8);
3206 :
3207 1107 : return phi;
3208 : }
3209 :
3210 : static GEN
3211 2944 : phi_w3w7_j(void)
3212 : {
3213 : GEN phi, phi0, phi1;
3214 2944 : phi = cgetg(4, t_VEC);
3215 :
3216 2944 : phi0 = cgetg(34, t_VEC);
3217 2944 : gel(phi0, 1) = gen_1;
3218 2944 : gel(phi0, 2) = utoi(0x14UL); setsigne(gel(phi0, 2), -1);
3219 2944 : gel(phi0, 3) = utoi(0x82UL);
3220 2944 : gel(phi0, 4) = utoi(0x1f8UL);
3221 2944 : gel(phi0, 5) = utoi(0x2a45UL); setsigne(gel(phi0, 5), -1);
3222 2944 : gel(phi0, 6) = utoi(0x9300UL);
3223 2944 : gel(phi0, 7) = utoi(0x32abeUL);
3224 2944 : gel(phi0, 8) = utoi(0x19c91cUL); setsigne(gel(phi0, 8), -1);
3225 2944 : gel(phi0, 9) = utoi(0xc1ba9UL);
3226 2944 : gel(phi0, 10) = utoi(0x1788f68UL);
3227 2944 : gel(phi0, 11) = utoi(0x2b1989cUL); setsigne(gel(phi0, 11), -1);
3228 2944 : gel(phi0, 12) = utoi(0x7a92408UL); setsigne(gel(phi0, 12), -1);
3229 2944 : gel(phi0, 13) = utoi(0x1238d56eUL);
3230 2944 : gel(phi0, 14) = utoi(0x13dd66a0UL);
3231 2944 : gel(phi0, 15) = utoi(0x2dbedca8UL); setsigne(gel(phi0, 15), -1);
3232 2944 : gel(phi0, 16) = utoi(0x34282eb8UL); setsigne(gel(phi0, 16), -1);
3233 2944 : gel(phi0, 17) = utoi(0x2c2a54d2UL);
3234 2944 : gel(phi0, 18) = utoi(0x98db81a8UL);
3235 2944 : gel(phi0, 19) = utoi(0x4088be8UL); setsigne(gel(phi0, 19), -1);
3236 2944 : gel(phi0, 20) = utoi(0xe424a220UL); setsigne(gel(phi0, 20), -1);
3237 2944 : gel(phi0, 21) = utoi(0x67bbb232UL); setsigne(gel(phi0, 21), -1);
3238 2944 : gel(phi0, 22) = utoi(0x7dd8bb98UL);
3239 2944 : gel(phi0, 23) = uu32toi(0x1UL, 0xcaff744UL);
3240 2944 : gel(phi0, 24) = utoi(0x1d46a378UL); setsigne(gel(phi0, 24), -1);
3241 2944 : gel(phi0, 25) = utoi(0x82fa50f7UL); setsigne(gel(phi0, 25), -1);
3242 2944 : gel(phi0, 26) = utoi(0x700ef38cUL); setsigne(gel(phi0, 26), -1);
3243 2944 : gel(phi0, 27) = utoi(0x20aa202eUL);
3244 2944 : gel(phi0, 28) = utoi(0x299b3440UL);
3245 2944 : gel(phi0, 29) = utoi(0xa476c4bUL);
3246 2944 : gel(phi0, 30) = utoi(0xd80558UL);
3247 2944 : gel(phi0, 31) = utoi(0x28a32UL);
3248 2944 : gel(phi0, 32) = utoi(0x2bcUL);
3249 2944 : gel(phi0, 33) = gen_1;
3250 :
3251 2944 : phi1 = cgetg(33, t_VEC);
3252 2944 : gel(phi1, 1) = gen_0;
3253 2944 : gel(phi1, 2) = gen_0;
3254 2944 : gel(phi1, 3) = gen_0;
3255 2944 : gel(phi1, 4) = gen_m1;
3256 2944 : gel(phi1, 5) = utoi(0xeUL);
3257 2944 : gel(phi1, 6) = utoi(0x31UL); setsigne(gel(phi1, 6), -1);
3258 2944 : gel(phi1, 7) = utoi(0xeUL); setsigne(gel(phi1, 7), -1);
3259 2944 : gel(phi1, 8) = utoi(0x99UL);
3260 2944 : gel(phi1, 9) = utoi(0x8UL); setsigne(gel(phi1, 9), -1);
3261 2944 : gel(phi1, 10) = utoi(0x2eUL); setsigne(gel(phi1, 10), -1);
3262 2944 : gel(phi1, 11) = utoi(0x5ccUL); setsigne(gel(phi1, 11), -1);
3263 2944 : gel(phi1, 12) = utoi(0x308UL);
3264 2944 : gel(phi1, 13) = utoi(0x2904UL);
3265 2944 : gel(phi1, 14) = utoi(0x15700UL); setsigne(gel(phi1, 14), -1);
3266 2944 : gel(phi1, 15) = utoi(0x2b9ecUL); setsigne(gel(phi1, 15), -1);
3267 2944 : gel(phi1, 16) = utoi(0xf0966UL);
3268 2944 : gel(phi1, 17) = utoi(0xb3cc8UL);
3269 2944 : gel(phi1, 18) = utoi(0x38241cUL); setsigne(gel(phi1, 18), -1);
3270 2944 : gel(phi1, 19) = utoi(0x8604cUL); setsigne(gel(phi1, 19), -1);
3271 2944 : gel(phi1, 20) = utoi(0x578a64UL);
3272 2944 : gel(phi1, 21) = utoi(0x11a798UL); setsigne(gel(phi1, 21), -1);
3273 2944 : gel(phi1, 22) = utoi(0x39c85eUL); setsigne(gel(phi1, 22), -1);
3274 2944 : gel(phi1, 23) = utoi(0x1a5084UL);
3275 2944 : gel(phi1, 24) = utoi(0xcdeb4UL);
3276 2944 : gel(phi1, 25) = utoi(0xb0364UL); setsigne(gel(phi1, 25), -1);
3277 2944 : gel(phi1, 26) = utoi(0x129d4UL);
3278 2944 : gel(phi1, 27) = utoi(0x126fcUL);
3279 2944 : gel(phi1, 28) = utoi(0x8649UL); setsigne(gel(phi1, 28), -1);
3280 2944 : gel(phi1, 29) = utoi(0x1aa2UL);
3281 2944 : gel(phi1, 30) = utoi(0x2dfUL); setsigne(gel(phi1, 30), -1);
3282 2944 : gel(phi1, 31) = utoi(0x2aUL);
3283 2944 : gel(phi1, 32) = gen_m1;
3284 :
3285 2944 : gel(phi, 1) = phi0;
3286 2944 : gel(phi, 2) = phi1;
3287 2944 : gel(phi, 3) = utoi(10);
3288 :
3289 2944 : return phi;
3290 : }
3291 :
3292 : static GEN
3293 210 : phi_w3w13_j(void)
3294 : {
3295 : GEN phi, phi0, phi1;
3296 210 : phi = cgetg(4, t_VEC);
3297 :
3298 210 : phi0 = cgetg(58, t_VEC);
3299 210 : gel(phi0, 1) = gen_1;
3300 210 : gel(phi0, 2) = utoi(0x10UL); setsigne(gel(phi0, 2), -1);
3301 210 : gel(phi0, 3) = utoi(0x58UL);
3302 210 : gel(phi0, 4) = utoi(0x258UL);
3303 210 : gel(phi0, 5) = utoi(0x270cUL); setsigne(gel(phi0, 5), -1);
3304 210 : gel(phi0, 6) = utoi(0x9c00UL);
3305 210 : gel(phi0, 7) = utoi(0x2b40cUL);
3306 210 : gel(phi0, 8) = utoi(0x20e250UL); setsigne(gel(phi0, 8), -1);
3307 210 : gel(phi0, 9) = utoi(0x4f46baUL);
3308 210 : gel(phi0, 10) = utoi(0x1869448UL);
3309 210 : gel(phi0, 11) = utoi(0xa49ab68UL); setsigne(gel(phi0, 11), -1);
3310 210 : gel(phi0, 12) = utoi(0x96c7630UL);
3311 210 : gel(phi0, 13) = utoi(0x4f7e0af6UL);
3312 210 : gel(phi0, 14) = utoi(0xea093590UL); setsigne(gel(phi0, 14), -1);
3313 210 : gel(phi0, 15) = utoi(0x6735bc50UL); setsigne(gel(phi0, 15), -1);
3314 210 : gel(phi0, 16) = uu32toi(0x5UL, 0x971a2e08UL);
3315 210 : gel(phi0, 17) = uu32toi(0x6UL, 0x29c9d965UL); setsigne(gel(phi0, 17), -1);
3316 210 : gel(phi0, 18) = uu32toi(0xdUL, 0xeb9aa360UL); setsigne(gel(phi0, 18), -1);
3317 210 : gel(phi0, 19) = uu32toi(0x26UL, 0xe9c0584UL);
3318 210 : gel(phi0, 20) = uu32toi(0x1UL, 0xb0cadce8UL); setsigne(gel(phi0, 20), -1);
3319 210 : gel(phi0, 21) = uu32toi(0x62UL, 0x73586014UL); setsigne(gel(phi0, 21), -1);
3320 210 : gel(phi0, 22) = uu32toi(0x66UL, 0xaf672e38UL);
3321 210 : gel(phi0, 23) = uu32toi(0x6bUL, 0x93c28cdcUL);
3322 210 : gel(phi0, 24) = uu32toi(0x11eUL, 0x4f633080UL); setsigne(gel(phi0, 24), -1);
3323 210 : gel(phi0, 25) = uu32toi(0x3cUL, 0xcc42461bUL);
3324 210 : gel(phi0, 26) = uu32toi(0x17bUL, 0xdec0a78UL);
3325 210 : gel(phi0, 27) = uu32toi(0x166UL, 0x910d8bd0UL); setsigne(gel(phi0, 27), -1);
3326 210 : gel(phi0, 28) = uu32toi(0xd4UL, 0x47873030UL); setsigne(gel(phi0, 28), -1);
3327 210 : gel(phi0, 29) = uu32toi(0x204UL, 0x811828baUL);
3328 210 : gel(phi0, 30) = uu32toi(0x50UL, 0x5d713960UL); setsigne(gel(phi0, 30), -1);
3329 210 : gel(phi0, 31) = uu32toi(0x198UL, 0xa27e42b0UL); setsigne(gel(phi0, 31), -1);
3330 210 : gel(phi0, 32) = uu32toi(0xe1UL, 0x25685138UL);
3331 210 : gel(phi0, 33) = uu32toi(0xe3UL, 0xaa5774bbUL);
3332 210 : gel(phi0, 34) = uu32toi(0xcfUL, 0x392a9a00UL); setsigne(gel(phi0, 34), -1);
3333 210 : gel(phi0, 35) = uu32toi(0x81UL, 0xfb334d04UL); setsigne(gel(phi0, 35), -1);
3334 210 : gel(phi0, 36) = uu32toi(0xabUL, 0x59594a68UL);
3335 210 : gel(phi0, 37) = uu32toi(0x42UL, 0x356993acUL);
3336 210 : gel(phi0, 38) = uu32toi(0x86UL, 0x307ba678UL); setsigne(gel(phi0, 38), -1);
3337 210 : gel(phi0, 39) = uu32toi(0xbUL, 0x7a9e59dcUL); setsigne(gel(phi0, 39), -1);
3338 210 : gel(phi0, 40) = uu32toi(0x4cUL, 0x27935f20UL);
3339 210 : gel(phi0, 41) = uu32toi(0x2UL, 0xe0ac9045UL); setsigne(gel(phi0, 41), -1);
3340 210 : gel(phi0, 42) = uu32toi(0x24UL, 0x14495758UL); setsigne(gel(phi0, 42), -1);
3341 210 : gel(phi0, 43) = utoi(0x20973410UL);
3342 210 : gel(phi0, 44) = uu32toi(0x13UL, 0x99ff4e00UL);
3343 210 : gel(phi0, 45) = uu32toi(0x1UL, 0xa710d34aUL); setsigne(gel(phi0, 45), -1);
3344 210 : gel(phi0, 46) = uu32toi(0x7UL, 0xfe5405c0UL); setsigne(gel(phi0, 46), -1);
3345 210 : gel(phi0, 47) = uu32toi(0x1UL, 0xcdee0f8UL);
3346 210 : gel(phi0, 48) = uu32toi(0x2UL, 0x660c92a8UL);
3347 210 : gel(phi0, 49) = utoi(0x3f13a35aUL);
3348 210 : gel(phi0, 50) = utoi(0xe4eb4ba0UL); setsigne(gel(phi0, 50), -1);
3349 210 : gel(phi0, 51) = utoi(0x6420f4UL); setsigne(gel(phi0, 51), -1);
3350 210 : gel(phi0, 52) = utoi(0x2c624370UL);
3351 210 : gel(phi0, 53) = utoi(0xb31b814UL);
3352 210 : gel(phi0, 54) = utoi(0xdd3ad8UL);
3353 210 : gel(phi0, 55) = utoi(0x29278UL);
3354 210 : gel(phi0, 56) = utoi(0x2c0UL);
3355 210 : gel(phi0, 57) = gen_1;
3356 :
3357 210 : phi1 = cgetg(57, t_VEC);
3358 210 : gel(phi1, 1) = gen_0;
3359 210 : gel(phi1, 2) = gen_0;
3360 210 : gel(phi1, 3) = gen_0;
3361 210 : gel(phi1, 4) = gen_m1;
3362 210 : gel(phi1, 5) = utoi(0xdUL);
3363 210 : gel(phi1, 6) = utoi(0x34UL); setsigne(gel(phi1, 6), -1);
3364 210 : gel(phi1, 7) = utoi(0x1aUL);
3365 210 : gel(phi1, 8) = utoi(0xf7UL);
3366 210 : gel(phi1, 9) = utoi(0x16cUL); setsigne(gel(phi1, 9), -1);
3367 210 : gel(phi1, 10) = utoi(0xddUL); setsigne(gel(phi1, 10), -1);
3368 210 : gel(phi1, 11) = utoi(0x28aUL);
3369 210 : gel(phi1, 12) = utoi(0xddUL); setsigne(gel(phi1, 12), -1);
3370 210 : gel(phi1, 13) = utoi(0x16cUL); setsigne(gel(phi1, 13), -1);
3371 210 : gel(phi1, 14) = utoi(0xf6UL);
3372 210 : gel(phi1, 15) = utoi(0x1dUL);
3373 210 : gel(phi1, 16) = utoi(0x31UL); setsigne(gel(phi1, 16), -1);
3374 210 : gel(phi1, 17) = utoi(0x5ceUL); setsigne(gel(phi1, 17), -1);
3375 210 : gel(phi1, 18) = utoi(0x2e4UL);
3376 210 : gel(phi1, 19) = utoi(0x252cUL);
3377 210 : gel(phi1, 20) = utoi(0x1b34cUL); setsigne(gel(phi1, 20), -1);
3378 210 : gel(phi1, 21) = utoi(0xaf80UL);
3379 210 : gel(phi1, 22) = utoi(0x1cc5f9UL);
3380 210 : gel(phi1, 23) = utoi(0x3e1aa5UL); setsigne(gel(phi1, 23), -1);
3381 210 : gel(phi1, 24) = utoi(0x86d17aUL); setsigne(gel(phi1, 24), -1);
3382 210 : gel(phi1, 25) = utoi(0x2427264UL);
3383 210 : gel(phi1, 26) = utoi(0x691c1fUL); setsigne(gel(phi1, 26), -1);
3384 210 : gel(phi1, 27) = utoi(0x862ad4eUL); setsigne(gel(phi1, 27), -1);
3385 210 : gel(phi1, 28) = utoi(0xab21e1fUL);
3386 210 : gel(phi1, 29) = utoi(0xbc19ddcUL);
3387 210 : gel(phi1, 30) = utoi(0x24331db8UL); setsigne(gel(phi1, 30), -1);
3388 210 : gel(phi1, 31) = utoi(0x972c105UL);
3389 210 : gel(phi1, 32) = utoi(0x363d7107UL);
3390 210 : gel(phi1, 33) = utoi(0x39696450UL); setsigne(gel(phi1, 33), -1);
3391 210 : gel(phi1, 34) = utoi(0x1bce7c48UL); setsigne(gel(phi1, 34), -1);
3392 210 : gel(phi1, 35) = utoi(0x552ecba0UL);
3393 210 : gel(phi1, 36) = utoi(0x1c7771b8UL); setsigne(gel(phi1, 36), -1);
3394 210 : gel(phi1, 37) = utoi(0x393029b8UL); setsigne(gel(phi1, 37), -1);
3395 210 : gel(phi1, 38) = utoi(0x3755be97UL);
3396 210 : gel(phi1, 39) = utoi(0x83402a9UL);
3397 210 : gel(phi1, 40) = utoi(0x24d5be62UL); setsigne(gel(phi1, 40), -1);
3398 210 : gel(phi1, 41) = utoi(0xdb6d90aUL);
3399 210 : gel(phi1, 42) = utoi(0xa0ef177UL);
3400 210 : gel(phi1, 43) = utoi(0x99ff162UL); setsigne(gel(phi1, 43), -1);
3401 210 : gel(phi1, 44) = utoi(0xb09e27UL);
3402 210 : gel(phi1, 45) = utoi(0x26a7adcUL);
3403 210 : gel(phi1, 46) = utoi(0x116e2fcUL); setsigne(gel(phi1, 46), -1);
3404 210 : gel(phi1, 47) = utoi(0x1383b5UL); setsigne(gel(phi1, 47), -1);
3405 210 : gel(phi1, 48) = utoi(0x35a9e7UL);
3406 210 : gel(phi1, 49) = utoi(0x1082a0UL); setsigne(gel(phi1, 49), -1);
3407 210 : gel(phi1, 50) = utoi(0x4696UL); setsigne(gel(phi1, 50), -1);
3408 210 : gel(phi1, 51) = utoi(0x19f98UL);
3409 210 : gel(phi1, 52) = utoi(0x8bb3UL); setsigne(gel(phi1, 52), -1);
3410 210 : gel(phi1, 53) = utoi(0x18bbUL);
3411 210 : gel(phi1, 54) = utoi(0x297UL); setsigne(gel(phi1, 54), -1);
3412 210 : gel(phi1, 55) = utoi(0x27UL);
3413 210 : gel(phi1, 56) = gen_m1;
3414 :
3415 210 : gel(phi, 1) = phi0;
3416 210 : gel(phi, 2) = phi1;
3417 210 : gel(phi, 3) = utoi(16);
3418 :
3419 210 : return phi;
3420 : }
3421 :
3422 : static GEN
3423 2567 : phi_w5w7_j(void)
3424 : {
3425 : GEN phi, phi0, phi1;
3426 2567 : phi = cgetg(4, t_VEC);
3427 :
3428 2567 : phi0 = cgetg(50, t_VEC);
3429 2567 : gel(phi0, 1) = gen_1;
3430 2567 : gel(phi0, 2) = utoi(0xcUL);
3431 2567 : gel(phi0, 3) = utoi(0x2aUL);
3432 2567 : gel(phi0, 4) = utoi(0x10UL);
3433 2567 : gel(phi0, 5) = utoi(0x69UL); setsigne(gel(phi0, 5), -1);
3434 2567 : gel(phi0, 6) = utoi(0x318UL); setsigne(gel(phi0, 6), -1);
3435 2567 : gel(phi0, 7) = utoi(0x148aUL); setsigne(gel(phi0, 7), -1);
3436 2567 : gel(phi0, 8) = utoi(0x17c4UL); setsigne(gel(phi0, 8), -1);
3437 2567 : gel(phi0, 9) = utoi(0x1a73UL);
3438 2567 : gel(phi0, 10) = gen_0;
3439 2567 : gel(phi0, 11) = utoi(0x338a0UL);
3440 2567 : gel(phi0, 12) = utoi(0x61698UL);
3441 2567 : gel(phi0, 13) = utoi(0x96e8UL); setsigne(gel(phi0, 13), -1);
3442 2567 : gel(phi0, 14) = utoi(0x140910UL);
3443 2567 : gel(phi0, 15) = utoi(0x45f6b4UL); setsigne(gel(phi0, 15), -1);
3444 2567 : gel(phi0, 16) = utoi(0x309f50UL); setsigne(gel(phi0, 16), -1);
3445 2567 : gel(phi0, 17) = utoi(0xef9f8bUL); setsigne(gel(phi0, 17), -1);
3446 2567 : gel(phi0, 18) = utoi(0x283167cUL); setsigne(gel(phi0, 18), -1);
3447 2567 : gel(phi0, 19) = utoi(0x625e20aUL);
3448 2567 : gel(phi0, 20) = utoi(0x16186350UL); setsigne(gel(phi0, 20), -1);
3449 2567 : gel(phi0, 21) = utoi(0x46861281UL);
3450 2567 : gel(phi0, 22) = utoi(0x754b96a0UL); setsigne(gel(phi0, 22), -1);
3451 2567 : gel(phi0, 23) = uu32toi(0x1UL, 0x421ca02aUL);
3452 2567 : gel(phi0, 24) = uu32toi(0x2UL, 0xdb76a5cUL); setsigne(gel(phi0, 24), -1);
3453 2567 : gel(phi0, 25) = uu32toi(0x4UL, 0xf6afd8eUL);
3454 2567 : gel(phi0, 26) = uu32toi(0x6UL, 0xaafd3cb4UL); setsigne(gel(phi0, 26), -1);
3455 2567 : gel(phi0, 27) = uu32toi(0x8UL, 0xda2539caUL);
3456 2567 : gel(phi0, 28) = uu32toi(0xfUL, 0x84343790UL); setsigne(gel(phi0, 28), -1);
3457 2567 : gel(phi0, 29) = uu32toi(0xfUL, 0x914ff421UL);
3458 2567 : gel(phi0, 30) = uu32toi(0x19UL, 0x3c123950UL); setsigne(gel(phi0, 30), -1);
3459 2567 : gel(phi0, 31) = uu32toi(0x15UL, 0x381f722aUL);
3460 2567 : gel(phi0, 32) = uu32toi(0x15UL, 0xe01c0c24UL); setsigne(gel(phi0, 32), -1);
3461 2567 : gel(phi0, 33) = uu32toi(0x19UL, 0x3360b375UL);
3462 2567 : gel(phi0, 34) = utoi(0x59fda9c0UL); setsigne(gel(phi0, 34), -1);
3463 2567 : gel(phi0, 35) = uu32toi(0x20UL, 0xff55024cUL);
3464 2567 : gel(phi0, 36) = uu32toi(0x16UL, 0xcc600800UL);
3465 2567 : gel(phi0, 37) = uu32toi(0x24UL, 0x1879c898UL);
3466 2567 : gel(phi0, 38) = uu32toi(0x1cUL, 0x37f97498UL);
3467 2567 : gel(phi0, 39) = uu32toi(0x19UL, 0x39ec4b60UL);
3468 2567 : gel(phi0, 40) = uu32toi(0x10UL, 0x52c660d0UL);
3469 2567 : gel(phi0, 41) = uu32toi(0x9UL, 0xcab00333UL);
3470 2567 : gel(phi0, 42) = uu32toi(0x4UL, 0x7fe69be4UL);
3471 2567 : gel(phi0, 43) = uu32toi(0x1UL, 0xa0c6f116UL);
3472 2567 : gel(phi0, 44) = utoi(0x69244638UL);
3473 2567 : gel(phi0, 45) = utoi(0xed560f7UL);
3474 2567 : gel(phi0, 46) = utoi(0xe7b660UL);
3475 2567 : gel(phi0, 47) = utoi(0x29d8aUL);
3476 2567 : gel(phi0, 48) = utoi(0x2c4UL);
3477 2567 : gel(phi0, 49) = gen_1;
3478 :
3479 2567 : phi1 = cgetg(49, t_VEC);
3480 2567 : gel(phi1, 1) = gen_0;
3481 2567 : gel(phi1, 2) = gen_0;
3482 2567 : gel(phi1, 3) = gen_0;
3483 2567 : gel(phi1, 4) = gen_0;
3484 2567 : gel(phi1, 5) = gen_0;
3485 2567 : gel(phi1, 6) = gen_1;
3486 2567 : gel(phi1, 7) = utoi(0x7UL);
3487 2567 : gel(phi1, 8) = utoi(0x8UL);
3488 2567 : gel(phi1, 9) = utoi(0x9UL); setsigne(gel(phi1, 9), -1);
3489 2567 : gel(phi1, 10) = gen_0;
3490 2567 : gel(phi1, 11) = utoi(0x13UL); setsigne(gel(phi1, 11), -1);
3491 2567 : gel(phi1, 12) = utoi(0x7UL); setsigne(gel(phi1, 12), -1);
3492 2567 : gel(phi1, 13) = utoi(0x5ceUL); setsigne(gel(phi1, 13), -1);
3493 2567 : gel(phi1, 14) = utoi(0xb0UL); setsigne(gel(phi1, 14), -1);
3494 2567 : gel(phi1, 15) = utoi(0x460UL);
3495 2567 : gel(phi1, 16) = utoi(0x194bUL); setsigne(gel(phi1, 16), -1);
3496 2567 : gel(phi1, 17) = utoi(0x87c3UL);
3497 2567 : gel(phi1, 18) = utoi(0x3cdeUL);
3498 2567 : gel(phi1, 19) = utoi(0xd683UL); setsigne(gel(phi1, 19), -1);
3499 2567 : gel(phi1, 20) = utoi(0x6099bUL);
3500 2567 : gel(phi1, 21) = utoi(0x111ea8UL); setsigne(gel(phi1, 21), -1);
3501 2567 : gel(phi1, 22) = utoi(0xfa113UL);
3502 2567 : gel(phi1, 23) = utoi(0x1a6561UL); setsigne(gel(phi1, 23), -1);
3503 2567 : gel(phi1, 24) = utoi(0x1e997UL); setsigne(gel(phi1, 24), -1);
3504 2567 : gel(phi1, 25) = utoi(0x214e54UL);
3505 2567 : gel(phi1, 26) = utoi(0x29c3f4UL); setsigne(gel(phi1, 26), -1);
3506 2567 : gel(phi1, 27) = utoi(0x67e102UL);
3507 2567 : gel(phi1, 28) = utoi(0x227eaaUL); setsigne(gel(phi1, 28), -1);
3508 2567 : gel(phi1, 29) = utoi(0x191d10UL);
3509 2567 : gel(phi1, 30) = utoi(0x1a9cd5UL);
3510 2567 : gel(phi1, 31) = utoi(0x58386fUL); setsigne(gel(phi1, 31), -1);
3511 2567 : gel(phi1, 32) = utoi(0x2e49f6UL);
3512 2567 : gel(phi1, 33) = utoi(0x31194bUL); setsigne(gel(phi1, 33), -1);
3513 2567 : gel(phi1, 34) = utoi(0x9e07aUL);
3514 2567 : gel(phi1, 35) = utoi(0x260d59UL);
3515 2567 : gel(phi1, 36) = utoi(0x189921UL); setsigne(gel(phi1, 36), -1);
3516 2567 : gel(phi1, 37) = utoi(0xeca4aUL);
3517 2567 : gel(phi1, 38) = utoi(0xa3d9cUL); setsigne(gel(phi1, 38), -1);
3518 2567 : gel(phi1, 39) = utoi(0x426daUL); setsigne(gel(phi1, 39), -1);
3519 2567 : gel(phi1, 40) = utoi(0x91875UL);
3520 2567 : gel(phi1, 41) = utoi(0x3b55bUL); setsigne(gel(phi1, 41), -1);
3521 2567 : gel(phi1, 42) = utoi(0x56f4UL); setsigne(gel(phi1, 42), -1);
3522 2567 : gel(phi1, 43) = utoi(0xcd1bUL);
3523 2567 : gel(phi1, 44) = utoi(0x5159UL); setsigne(gel(phi1, 44), -1);
3524 2567 : gel(phi1, 45) = utoi(0x10f4UL);
3525 2567 : gel(phi1, 46) = utoi(0x20dUL); setsigne(gel(phi1, 46), -1);
3526 2567 : gel(phi1, 47) = utoi(0x23UL);
3527 2567 : gel(phi1, 48) = gen_m1;
3528 :
3529 2567 : gel(phi, 1) = phi0;
3530 2567 : gel(phi, 2) = phi1;
3531 2567 : gel(phi, 3) = utoi(12);
3532 :
3533 2567 : return phi;
3534 : }
3535 :
3536 : GEN
3537 21813 : double_eta_raw(long inv)
3538 : {
3539 21813 : switch (inv) {
3540 : case INV_W2W3:
3541 : case INV_W2W3E2:
3542 1449 : return phi_w2w3_j();
3543 : case INV_W3W3:
3544 : case INV_W3W3E2:
3545 2495 : return phi_w3w3_j();
3546 : case INV_W2W5:
3547 : case INV_W2W5E2:
3548 3500 : return phi_w2w5_j();
3549 : case INV_W2W7:
3550 : case INV_W2W7E2:
3551 5665 : return phi_w2w7_j();
3552 : case INV_W3W5:
3553 1107 : return phi_w3w5_j();
3554 : case INV_W3W7:
3555 2944 : return phi_w3w7_j();
3556 : case INV_W2W13:
3557 1876 : return phi_w2w13_j();
3558 : case INV_W3W13:
3559 210 : return phi_w3w13_j();
3560 : case INV_W5W7:
3561 2567 : return phi_w5w7_j();
3562 : default:
3563 0 : pari_err_BUG("double_eta_raw");
3564 : }
3565 0 : return NULL;
3566 : }
3567 :
3568 : /* F = double_eta_Fl(inv, p) */
3569 : static GEN
3570 37417 : Flx_double_eta_xpoly(GEN F, ulong j, ulong p, ulong pi)
3571 : {
3572 37417 : GEN u = gel(F,1), v = gel(F,2), w;
3573 37417 : long i, k = itos(gel(F,3)), lu = lg(u), lv = lg(v), lw = lu + 1;
3574 :
3575 37417 : w = cgetg(lw, t_VECSMALL); /* lu >= max(lv,k) */
3576 37416 : w[1] = 0; /* variable number */
3577 1028535 : for (i = 1; i < lv; i++)
3578 991117 : uel(w, i + 1) = Fl_add(uel(u, i), Fl_mul_pre(j, uel(v, i), p, pi), p);
3579 74836 : for ( ; i < lu; i++)
3580 37418 : uel(w, i + 1) = uel(u, i);
3581 37418 : uel(w, k + 2) = Fl_add(uel(w, k + 2), Fl_sqr_pre(j, p, pi), p);
3582 37418 : return Flx_renormalize(w, lw);
3583 : }
3584 :
3585 : /* F = double_eta_Fl(inv, p) */
3586 : static GEN
3587 26730 : Flx_double_eta_jpoly(GEN F, ulong x, ulong p, ulong pi)
3588 : {
3589 26730 : pari_sp av = avma;
3590 26730 : GEN u = gel(F,1), v = gel(F,2), xs;
3591 26730 : long k = itos(gel(F,3));
3592 : ulong a, b, c;
3593 :
3594 : /* u is always longest and the length is bigger than k */
3595 26730 : xs = Fl_powers_pre(x, lg(u) - 1, p, pi);
3596 26730 : c = Flv_dotproduct_pre(u, xs, p, pi);
3597 26729 : b = Flv_dotproduct_pre(v, xs, p, pi);
3598 26730 : a = uel(xs, k + 1);
3599 26730 : avma = av;
3600 26730 : return mkvecsmall4(0, c, b, a);
3601 : }
3602 :
3603 : /**
3604 : * SECTION: Select discriminant for given modpoly level.
3605 : */
3606 :
3607 : /* require an L1, useful for multi-threading */
3608 : #define MODPOLY_USE_L1 1
3609 : /* no bound on L1 other than the fixed bound MAX_L1 - needed to
3610 : * handle small L for certain invariants (but not for j) */
3611 : #define MODPOLY_NO_MAX_L1 2
3612 : /* don't use any auxilliary primes - needed to handle small L for
3613 : * certain invariants (but not for j) */
3614 : #define MODPOLY_NO_AUX_L 4
3615 : #define MODPOLY_IGNORE_SPARSE_FACTOR 8
3616 :
3617 : INLINE double
3618 2566 : modpoly_height_bound(long L, long inv)
3619 : {
3620 : double nbits, nbits2;
3621 : double c;
3622 : long hf;
3623 :
3624 : /* proven bound (in bits), derived from: 6l*log(l)+16*l+13*sqrt(l)*log(l) */
3625 2566 : nbits = 6.0*L*log2(L)+16/LOG2*L+8.0*sqrt((double)L)*log2(L);
3626 : /* alternative proven bound (in bits), derived from: 6l*log(l)+17*l */
3627 2566 : nbits2 = 6.0*L*log2(L)+17/LOG2*L;
3628 2566 : if ( nbits2 < nbits ) nbits = nbits2;
3629 2566 : hf = modinv_height_factor(inv);
3630 2566 : if (hf > 1) {
3631 : /*
3632 : *** Important *** when dividing by the height factor, we only want to reduce terms related
3633 : to the bound on j (the roots of Phi_l(X,y)), not terms arising from binomial coefficients.
3634 : These arise in lemmas 2 and 3 of the height bound paper, terms of (log 2)*L and 2.085*(L+1)
3635 : which we convert here to binary logs.
3636 : */
3637 : /* this is a massive overestimate, if you care about speed,
3638 : * determine a good height bound empirically as done for INV_F
3639 : * below */
3640 1536 : nbits2 = nbits - 4.01*L -3.0;
3641 1536 : nbits = nbits2/hf + 4.01*L + 3.0;
3642 : }
3643 :
3644 2566 : if (inv == INV_F) {
3645 184 : if (L < 30) c = 45;
3646 35 : else if (L < 100) c = 36;
3647 21 : else if (L < 300) c = 32;
3648 7 : else if (L < 600) c = 26;
3649 0 : else if (L < 1200) c = 24;
3650 0 : else if (L < 2400) c = 22;
3651 0 : else c = 20;
3652 184 : nbits = (6.0*L*log2(L) + c*L)/hf;
3653 : }
3654 :
3655 2566 : return nbits;
3656 : }
3657 :
3658 : /* small enough so we can write the factorization of a smooth in a BIL
3659 : * bit integer */
3660 : #define SMOOTH_PRIMES ((BITS_IN_LONG >> 1) - 1)
3661 :
3662 :
3663 : #define MAX_ATKIN 255
3664 :
3665 : /* Must have primes at least up to MAX_ATKIN */
3666 : static const long PRIMES[] = {
3667 : 0, 2, 3, 5, 7, 11, 13, 17, 19, 23,
3668 : 29, 31, 37, 41, 43, 47, 53, 59, 61, 67,
3669 : 71, 73, 79, 83, 89, 97, 101, 103, 107, 109,
3670 : 113, 127, 131, 137, 139, 149, 151, 157, 163, 167,
3671 : 173, 179, 181, 191, 193, 197, 199, 211, 223, 227,
3672 : 229, 233, 239, 241, 251, 257, 263, 269, 271, 277
3673 : };
3674 :
3675 : #define MAX_L1 255
3676 :
3677 : typedef struct D_entry_struct {
3678 : ulong m;
3679 : long D, h;
3680 : } D_entry;
3681 :
3682 : /* Returns a form that generates the classes of norm p^2 in cl(p^2D)
3683 : * (i.e. one with order p-1), where p is an odd prime that splits in D
3684 : * and does not divide its conductor (but this is not verified) */
3685 : INLINE GEN
3686 57879 : qform_primeform2(long p, long D)
3687 : {
3688 57879 : pari_sp ltop = avma, av;
3689 : GEN M;
3690 : register long k;
3691 :
3692 57879 : M = factor_pn_1(stoi(p), 1);
3693 57879 : av = avma;
3694 117630 : for (k = D & 1; k <= p; k += 2)
3695 : {
3696 : GEN Q;
3697 117630 : long ord, a, b, c = (k * k - D) / 4;
3698 :
3699 117630 : if (!(c % p)) continue;
3700 101471 : a = p * p;
3701 101471 : b = k * p;
3702 101471 : Q = redimag(mkqfis(a, b, c));
3703 : /* TODO: How do we know that Q has order dividing p - 1? If we
3704 : * don't, then the call to gen_order should be replaced with a
3705 : * call to something with fastorder semantics (i.e. return 0 if
3706 : * ord(Q) \ndiv M). */
3707 101471 : ord = itos(qfi_order(Q, M));
3708 101471 : if (ord == p - 1) {
3709 : /* TODO: This check that gen_order returned the correct result
3710 : * should be removed when gen_order is replaced with something
3711 : * with fastorder semantics. */
3712 57879 : GEN tst = gpowgs(Q, p - 1);
3713 57879 : if (qfb_equal1(tst)) { avma = ltop; return mkqfis(a, b, c); }
3714 0 : break;
3715 : }
3716 43592 : avma = av;
3717 : }
3718 0 : avma = ltop; return NULL;
3719 : }
3720 :
3721 : #define divides(a,b) (!((b) % (a)))
3722 :
3723 : /* This gets around the fact that gen_PH_log returns garbage if g
3724 : * doesn't have order n. */
3725 : INLINE GEN
3726 596537 : discrete_log(GEN a, GEN g, long n)
3727 : {
3728 596537 : return qfi_Shanks(a, g, n);
3729 : }
3730 :
3731 :
3732 : /*
3733 : * Output x such that [L0]^x = [L] in cl(D) where n = #cl(D). Return
3734 : * value indicates whether x was found or not.
3735 : */
3736 : INLINE long
3737 159477 : primeform_discrete_log(long *x, long L0, long L, long n, long D)
3738 : {
3739 : GEN X, Q, R, DD;
3740 159477 : pari_sp av = avma;
3741 159477 : long res = 0;
3742 :
3743 159477 : DD = stoi(D);
3744 159477 : Q = redimag(primeform_u(DD, L0));
3745 159477 : R = redimag(primeform_u(DD, L));
3746 159477 : X = discrete_log(R, Q, n);
3747 159477 : if (X) {
3748 60933 : *x = itos(X);
3749 60933 : res = 1;
3750 : }
3751 159477 : avma = av;
3752 159477 : return res;
3753 : }
3754 :
3755 :
3756 : /*
3757 : * Return the norm of a class group generator appropriate for a
3758 : * discriminant that will be used to calculate the modular polynomial
3759 : * of level L and invariant inv. Don't consider norms less than
3760 : * initial_L0.
3761 : */
3762 : static long
3763 2566 : select_L0(long L, long inv, long initial_L0)
3764 : {
3765 : long modinv_N;
3766 : long L0;
3767 :
3768 2566 : modinv_N = modinv_level(inv);
3769 2566 : if (divides(L, modinv_N))
3770 0 : pari_err_BUG("select_L0");
3771 :
3772 : /* TODO: Clean up these anomolous L0 choices */
3773 :
3774 : /* I've no idea why the discriminant-finding code fails with L0=5
3775 : * when L=19 and L=29, nor why L0=7 and L0=11 don't work for L=19
3776 : * either, nor why this happens for the otherwise unrelated
3777 : * invariants Weber-f and (2,3) double-eta. */
3778 2566 : if (inv == INV_W3W3E2 && L == 5)
3779 14 : return 2;
3780 :
3781 2552 : if (inv == INV_F || inv == INV_F2 || inv == INV_F4 || inv == INV_F8
3782 2256 : || inv == INV_W2W3 || inv == INV_W2W3E2
3783 2179 : || inv == INV_W3W3 /* || inv == INV_W3W3E2 */) {
3784 436 : if (L == 19)
3785 50 : return 13;
3786 386 : else if (L == 29 || L == 5)
3787 70 : return 7;
3788 316 : return 5;
3789 : }
3790 2116 : if ((inv == INV_W2W5 || inv == INV_W2W5E2)
3791 140 : && (L == 7 || L == 19))
3792 35 : return 13;
3793 2081 : if ((inv == INV_W2W7 || inv == INV_W2W7E2)
3794 323 : && L == 11)
3795 28 : return 13;
3796 2053 : if (inv == INV_W3W5) {
3797 63 : if (L == 7)
3798 7 : return 13;
3799 56 : else if (L == 17)
3800 0 : return 7;
3801 : }
3802 2046 : if (inv == INV_W3W7) {
3803 161 : if (L == 29 || L == 101)
3804 28 : return 11;
3805 133 : if (L == 11 || L == 19)
3806 35 : return 13;
3807 : }
3808 1983 : if (inv == INV_W5W7 && L == 17)
3809 21 : return 3;
3810 :
3811 : /* L0 is the smallest small prime different from L that doesn't divide modinv_N. */
3812 4877 : for (L0 = unextprime(initial_L0 + 1);
3813 2866 : L0 == L || divides(L0, modinv_N);
3814 953 : L0 = unextprime(L0 + 1))
3815 : ;
3816 1962 : return L0;
3817 : }
3818 :
3819 :
3820 : /*
3821 : * Return the order of [L]^n in cl(D), where #cl(D) = ord.
3822 : */
3823 : INLINE long
3824 851486 : primeform_exp_order(long L, long n, long D, long ord)
3825 : {
3826 851486 : pari_sp av = avma;
3827 851486 : GEN Q = gpowgs(primeform_u(stoi(D), L), n);
3828 851486 : long m = itos(qfi_order(Q, Z_factor(stoi(ord))));
3829 851486 : avma = av;
3830 851486 : return m;
3831 : }
3832 :
3833 : INLINE long
3834 87013 : eql_elt(GEN P, GEN Q, long i)
3835 : {
3836 87013 : return gequal(gel(P, i), gel(Q, i));
3837 : }
3838 :
3839 : /* If an ideal of norm modinv_deg is equivalent to an ideal of
3840 : * norm L0, we have an orientation ambiguity that we need to
3841 : * avoid. Note that we need to check all the possibilities (up
3842 : * to 8), but we can cheaply check inverses (so at most 2) */
3843 : static long
3844 32392 : orientation_ambiguity(long D1, long L0, long modinv_p1, long modinv_p2, long modinv_N)
3845 : {
3846 32392 : pari_sp av = avma;
3847 32392 : long ambiguity = 0;
3848 32392 : GEN D = stoi(D1);
3849 32392 : GEN Q1 = primeform_u(D, modinv_p1), Q2 = NULL;
3850 :
3851 32392 : if (modinv_p2 > 1) {
3852 32392 : if (modinv_p1 != modinv_p2) {
3853 27759 : GEN P2 = primeform_u(D, modinv_p2);
3854 27759 : GEN Q = gsqr(P2);
3855 27759 : GEN R = gsqr(Q1);
3856 : /* check that p1^2 != p2^{+/-2}, since this leads to
3857 : * ambiguities when converting j's to f's */
3858 27759 : if (eql_elt(Q, R, 1)
3859 84 : && (eql_elt(Q, R, 2) || gequal(gel(Q, 2), gneg(gel(R, 2))))) {
3860 0 : dbg_printf(3)("Bad D=%ld, a^2=b^2 problem between modinv_p1=%ld and modinv_p2=%ld\n",
3861 : D1, modinv_p1, modinv_p2);
3862 0 : ambiguity = 1;
3863 : } else {
3864 : /* generate both p1*p2 and p1*p2^{-1} */
3865 27759 : Q2 = gmul(Q1, P2);
3866 27759 : P2 = ginv(P2);
3867 27759 : Q1 = gmul(Q1, P2);
3868 : }
3869 : } else {
3870 4633 : Q1 = gsqr(Q1);
3871 : }
3872 : }
3873 32392 : if ( ! ambiguity) {
3874 32392 : GEN P = gsqr(primeform_u(D, L0));
3875 32392 : if (eql_elt(P, Q1, 1)
3876 31397 : || (modinv_p2 && modinv_p1 != modinv_p2 && eql_elt(P, Q2, 1))) {
3877 1424 : dbg_printf(3)("Bad D=%ld, a=b^{+/-2} problem between modinv_N=%ld and L0=%ld\n",
3878 : D1, modinv_N, L0);
3879 1424 : ambiguity = 1;
3880 : }
3881 : }
3882 32392 : avma = av;
3883 32392 : return ambiguity;
3884 : }
3885 :
3886 : static long
3887 616175 : check_generators(
3888 : long *n1_, long *m_,
3889 : long D, long h, long n, long subgrp_sz, long L0, long L1)
3890 : {
3891 616175 : long n1, m = primeform_exp_order(L0, n, D, h);
3892 616175 : if (m_)
3893 251077 : *m_ = m;
3894 616175 : n1 = n * m;
3895 616175 : if ( ! n1)
3896 0 : pari_err_BUG("check_generators");
3897 616175 : *n1_ = n1;
3898 616175 : if (n1 < subgrp_sz/2 || ( ! L1 && n1 < subgrp_sz)) {
3899 23043 : dbg_printf(3)("Bad D1=%ld with n1=%ld, h1=%ld, L1=%ld: "
3900 : "L0 and L1 don't span subgroup of size d in cl(D1)\n",
3901 : D, n, h, L1);
3902 23043 : return 0;
3903 : }
3904 593132 : if (n1 < subgrp_sz && ! (n1 & 1)) {
3905 : int res;
3906 : /* check whether L1 is generated by L0, using the fact that it has order 2 */
3907 13178 : pari_sp av = avma;
3908 13178 : GEN D1 = stoi(D);
3909 13178 : GEN Q = gpowgs(primeform_u(D1, L0), n1 / 2);
3910 13178 : res = gequal(Q, redimag(primeform_u(D1, L1)));
3911 13178 : avma = av;
3912 13178 : if (res) {
3913 3878 : dbg_printf(3)("Bad D1=%ld, with n1=%ld, h1=%ld, L1=%ld: "
3914 : "L1 generated by L0 in cl(D1)\n", D, n, h, L1);
3915 3878 : return 0;
3916 : }
3917 : }
3918 589254 : return 1;
3919 : }
3920 :
3921 : /*
3922 : * Calculate solutions (p, t) to the norm equation
3923 : *
3924 : * 4 p = t^2 - v^2 L^2 D (*)
3925 : *
3926 : * corresponding to the descriminant described by Dinfo.
3927 : *
3928 : * INPUT:
3929 : * - max: length of primes and traces
3930 : * - xprimes: p to exclude from primes (if they arise)
3931 : * - xcnt: length of xprimes
3932 : * - minbits: sum of log2(p) must be larger than this
3933 : * - Dinfo: discriminant, invariant and L for which we seek solutions
3934 : * to (*)
3935 : *
3936 : * OUTPUT:
3937 : * - primes: array of p in (*)
3938 : * - traces: array of t in (*)
3939 : * - totbits: sum of log2(p) for p in primes.
3940 : *
3941 : * RETURN:
3942 : * - the number of primes and traces found (these are always the same).
3943 : *
3944 : * NOTE: Any of primes, traces and totbits can be zero, in which case
3945 : * these values are discarded after calculation (however it is not
3946 : * permitted for traces to be zero if primes is non-zero). xprimes
3947 : * can be zero, in which case it is treated as empty.
3948 : */
3949 : static long
3950 9346 : modpoly_pickD_primes(
3951 : ulong *primes, ulong *traces, long max, ulong *xprimes, long xcnt,
3952 : long *totbits, long minbits, modpoly_disc_info *Dinfo)
3953 : {
3954 : double bits;
3955 : long D, m, n, vcnt, pfilter, one_prime, inv;
3956 : ulong maxp;
3957 : register ulong a1, a2, v, t, p, a1_start, a1_delta, L0, L1, L, absD;
3958 : /* FF_BITS = BITS_IN_LONG - NAIL_BITS (latter is 2 or 7) */
3959 9346 : ulong FF_BITS = BITS_IN_LONG - 2;
3960 :
3961 9346 : D = Dinfo->D1;
3962 9346 : L0 = Dinfo->L0;
3963 9346 : L1 = Dinfo->L1;
3964 9346 : L = Dinfo->L;
3965 9346 : inv = Dinfo->inv;
3966 :
3967 9346 : absD = -D;
3968 :
3969 : /* make sure pfilter and D don't preclude the possibility of p=(t^2-v^2D)/4 being prime */
3970 9346 : pfilter = modinv_pfilter(inv);
3971 9346 : if ((pfilter & IQ_FILTER_1MOD3) && ! (D % 3))
3972 35 : return 0;
3973 9311 : if ((pfilter & IQ_FILTER_1MOD4) && ! (D & 0xF))
3974 0 : return 0;
3975 :
3976 : /* Naively estimate the number of primes satisfying 4p=t^2-L^2D
3977 : * with t=2 mod L and pfilter using the PNT this is roughly #{t:
3978 : * t^2 < max p and t=2 mod L} / pi(max p) * filter_density, where
3979 : * filter_density is 1, 2, or 4 depending on pfilter. If this
3980 : * quantity is already more than twice the number of bits we need,
3981 : * assume that, barring some obstruction, we should have no problem
3982 : * getting enough primes. In this case we just verify we can get
3983 : * one prime (which should always be true, assuming we chose D
3984 : * properly). */
3985 9311 : one_prime = 0;
3986 9311 : if (totbits)
3987 9311 : *totbits = 0;
3988 9311 : if (max <= 1 && ! one_prime) {
3989 6725 : p = ((pfilter & IQ_FILTER_1MOD3) ? 2 : 1) * ((pfilter & IQ_FILTER_1MOD4) ? 2 : 1);
3990 13450 : one_prime = (1UL << ((FF_BITS+1)/2)) * (log2(L*L*(-D))-1)
3991 6725 : > p*L*minbits*FF_BITS*LOG2;
3992 6725 : if (one_prime && totbits)
3993 6629 : *totbits = minbits+1; /* lie */
3994 : }
3995 :
3996 9311 : m = n = 0;
3997 9311 : bits = 0.0;
3998 9311 : maxp = 0;
3999 23324 : for (v = 1; v < 100 && bits < minbits; v++) {
4000 : /* Don't allow v dividing the conductor. */
4001 20688 : if (ugcd(absD, v) != 1)
4002 98 : continue;
4003 : /* Avoid v dividing the level. */
4004 20590 : if (v > 2 && modinv_is_double_eta(inv) && ugcd(modinv_level(inv), v) != 1)
4005 966 : continue;
4006 : /* can't get odd p with D=1 mod 8 unless v is even */
4007 19624 : if ((v & 1) && (D & 7) == 1)
4008 9944 : continue;
4009 : /* don't use v divisible by 4 for L0=2 (we could remove this restriction, for a cost) */
4010 9680 : if (L0 == 2 && !(v & 3))
4011 111 : continue;
4012 : /* can't get p=3mod4 if v^2D is 0 mod 16 */
4013 9569 : if ((pfilter & IQ_FILTER_1MOD4) && !((v*v*D) & 0xF))
4014 83 : continue;
4015 9486 : if ((pfilter & IQ_FILTER_1MOD3) && !(v%3) )
4016 58 : continue;
4017 : /* avoid L0-volcanos with non-zero height */
4018 9428 : if (L0 != 2 && ! (v % L0))
4019 21 : continue;
4020 : /* ditto for L1 */
4021 9407 : if (L1 && !(v % L1))
4022 0 : continue;
4023 9407 : vcnt = 0;
4024 9407 : if ((v*v*absD)/4 > (1L<<FF_BITS)/(L*L))
4025 46 : break;
4026 9361 : if (((v*D) & 1)) {
4027 0 : a1_start = 1;
4028 0 : a1_delta = 2;
4029 : } else {
4030 9361 : a1_start = (((v*v*D) & 7) ? 2 : 0 );
4031 9361 : a1_delta = 4;
4032 : }
4033 443560 : for (a1 = a1_start; bits < minbits; a1 += a1_delta) {
4034 440938 : a2 = (a1*a1 + v*v*absD) >> 2;
4035 440938 : if ( !(a2 % L))
4036 67144 : continue;
4037 373794 : t = a1*L + 2;
4038 373794 : p = a2*L*L + t - 1;
4039 373794 : if ( !(p & 1))
4040 0 : pari_err_BUG("modpoly_pickD_primes");
4041 : /* double check calculation just in case of overflow or other weirdness */
4042 373794 : if (t*t + v*v*L*L*absD != 4*p)
4043 0 : pari_err_BUG("modpoly_pickD_primes");
4044 373794 : if (p > (1UL<<FF_BITS))
4045 110 : break;
4046 373684 : if (xprimes) {
4047 590886 : while (m < xcnt && xprimes[m] < p)
4048 382 : m++;
4049 295252 : if (m < xcnt && p == xprimes[m]) {
4050 0 : dbg_printf(1)("skipping duplicate prime %ld\n", p);
4051 0 : continue;
4052 : }
4053 : }
4054 373684 : if ( ! uisprime(p))
4055 331598 : continue;
4056 42086 : if ( ! modinv_good_prime(inv, p))
4057 1946 : continue;
4058 40140 : if (primes) {
4059 : /* ulong vLp, vLm; */
4060 30830 : if (n >= max)
4061 0 : goto done;
4062 : /* TODO: Implement test to filter primes that lead to
4063 : * L-valuation != 2 */
4064 30830 : primes[n] = p;
4065 30830 : traces[n] = t;
4066 : }
4067 40140 : n++;
4068 40140 : vcnt++;
4069 40140 : bits += log2(p);
4070 40140 : if (p > maxp)
4071 39271 : maxp = p;
4072 40140 : if (one_prime)
4073 6629 : goto done;
4074 : }
4075 2732 : if (vcnt)
4076 2729 : dbg_printf(3)("%ld primes with v=%ld, maxp=%ld (%.2f bits)\n",
4077 : vcnt, v, maxp, log2(maxp));
4078 : }
4079 : done:
4080 9311 : if ( ! n) {
4081 9 : dbg_printf(3)("check_primes failed completely for D=%ld\n", D);
4082 9 : return 0;
4083 : }
4084 9302 : dbg_printf(3)("D=%ld: Found %ld primes totalling %0.2f of %ld bits\n",
4085 : D, n, bits, minbits);
4086 9302 : if (totbits && ! *totbits)
4087 2673 : *totbits = (long)bits;
4088 9302 : return n;
4089 : }
4090 :
4091 : #define MAX_VOLCANO_FLOOR_SIZE 100000000
4092 :
4093 : static long
4094 2568 : calc_primes_for_discriminants(modpoly_disc_info Ds[], long Dcnt, long L, long minbits)
4095 : {
4096 2568 : pari_sp av = avma;
4097 : long i, j, k, m, n, D1, pcnt, totbits;
4098 : ulong *primes, *Dprimes, *Dtraces;
4099 :
4100 : /* D1 is the discriminant with smallest absolute value among those we found */
4101 2568 : D1 = Ds[0].D1;
4102 6716 : for (i = 1; i < Dcnt; i++) {
4103 4148 : if (Ds[i].D1 > D1)
4104 336 : D1 = Ds[i].D1;
4105 : }
4106 :
4107 : /* n is an upper bound on the number of primes we might get. */
4108 2568 : n = ceil(minbits / (log2(L * L * (-D1)) - 2)) + 1;
4109 2568 : primes = (ulong *) stack_malloc(n * sizeof(*primes));
4110 2568 : Dprimes = (ulong *) stack_malloc(n * sizeof(*Dprimes));
4111 2568 : Dtraces = (ulong *) stack_malloc(n * sizeof(*Dtraces));
4112 2586 : for (i = 0, totbits = 0, pcnt = 0; i < Dcnt && totbits < minbits; i++) {
4113 7758 : Ds[i].nprimes = modpoly_pickD_primes(Dprimes, Dtraces, n, primes, pcnt,
4114 5172 : &Ds[i].bits, minbits - totbits, Ds + i);
4115 2586 : Ds[i].primes = (ulong *) newblock(Ds[i].nprimes);
4116 2586 : memcpy(Ds[i].primes, Dprimes, Ds[i].nprimes * sizeof(*Dprimes));
4117 :
4118 2586 : Ds[i].traces = (ulong *) newblock(Ds[i].nprimes);
4119 2586 : memcpy(Ds[i].traces, Dtraces, Ds[i].nprimes * sizeof(*Dtraces));
4120 :
4121 2586 : totbits += Ds[i].bits;
4122 2586 : pcnt += Ds[i].nprimes;
4123 :
4124 2586 : if (totbits >= minbits || i == Dcnt - 1) {
4125 2568 : Dcnt = i + 1;
4126 2568 : break;
4127 : }
4128 : /* merge lists */
4129 552 : for (j = pcnt - Ds[i].nprimes - 1, k = Ds[i].nprimes - 1, m = pcnt - 1; m >= 0; m--) {
4130 534 : if (k >= 0) {
4131 509 : if (j >= 0 && primes[j] > Dprimes[k])
4132 301 : primes[m] = primes[j--];
4133 : else
4134 208 : primes[m] = Dprimes[k--];
4135 : } else {
4136 25 : primes[m] = primes[j--];
4137 : }
4138 : }
4139 : }
4140 2568 : if (totbits < minbits) {
4141 2 : dbg_printf(1)("Only obtained %ld of %ld bits using %ld discriminants\n",
4142 : totbits, minbits, Dcnt);
4143 2 : Dcnt = 0;
4144 : }
4145 2568 : avma = av;
4146 2568 : return Dcnt;
4147 : }
4148 :
4149 : /*
4150 : * Select discriminant(s) to use when calculating the modular
4151 : * polynomial of level L and invariant inv.
4152 : *
4153 : * INPUT:
4154 : * - L: level of modular polynomial (must be odd)
4155 : * - inv: invariant of modular polynomial
4156 : * - L0: result of select_L0(L, inv)
4157 : * - minbits: height of modular polynomial
4158 : * - flags: see below
4159 : * - tab: result of scanD0(L0)
4160 : * - tablen: length of tab
4161 : *
4162 : * OUTPUT:
4163 : * - Ds: the selected discriminant(s)
4164 : *
4165 : * RETURN:
4166 : * - the number of Ds found
4167 : *
4168 : * The flags parameter is constructed by ORing zero or more of the
4169 : * following values:
4170 : * - MODPOLY_USE_L1: force use of second class group generator
4171 : * - MODPOLY_NO_AUX_L: don't use auxillary class group elements
4172 : * - MODPOLY_IGNORE_SPARSE_FACTOR: obtain D for which h(D) > L + 1
4173 : * rather than h(D) > (L + 1)/s
4174 : */
4175 : static long
4176 2568 : modpoly_pickD(
4177 : modpoly_disc_info Ds[MODPOLY_MAX_DCNT], long L, long inv, long L0, long max_L1, long minbits, long flags,
4178 : D_entry *tab, long tablen)
4179 : {
4180 2568 : pari_sp ltop = avma, btop;
4181 : D_entry D0_entry;
4182 : modpoly_disc_info Dinfo;
4183 : pari_timer T;
4184 : long p, q, L1;
4185 : long tmp, modinv_p1, modinv_p2;
4186 : long modinv_N, modinv_deg, deg, use_L1, twofactor, pfilter, Dcnt;
4187 2568 : long how_many_D0s = 0;
4188 : double D0_bits, p_bits, L_bits;
4189 : long D0, D1, D2, m, h0, h1, h2, n0, n1, n2, dl1, dl2[2], d, cost, enum_cost, H_cost, best_cost, totbits;
4190 : register long i, j, k;
4191 :
4192 2568 : if ( !(L & 1))
4193 0 : pari_err_BUG("modpoly_pickD");
4194 :
4195 2568 : timer_start(&T);
4196 2568 : d = (flags & MODPOLY_IGNORE_SPARSE_FACTOR) ? 1 : modinv_sparse_factor(inv);
4197 : /*d = ui_ceil_ratio(L + 1, d) + 1; */
4198 2568 : tmp = (L + 1) / d;
4199 2568 : d = ((tmp * d < (L + 1)) ? tmp + 1 : tmp);
4200 2568 : d += 1;
4201 2568 : modinv_N = modinv_level(inv);
4202 2568 : modinv_deg = modinv_degree(&modinv_p1, &modinv_p2, inv);
4203 2568 : pfilter = modinv_pfilter(inv);
4204 :
4205 : /* Now set level to 0 unless we will need to compute N-isogenies */
4206 2568 : dbg_printf(1)("Using L0=%ld for L=%ld, d=%ld, modinv_N=%ld, modinv_deg=%ld\n",
4207 : L0, L, d, modinv_N, modinv_deg);
4208 :
4209 : /* We use L1 if (L0|L) == 1 or if we are forced to by flags. */
4210 2568 : use_L1 = (kross(L0,L) > 0 || (flags & MODPOLY_USE_L1));
4211 :
4212 2568 : Dcnt = 0;
4213 2568 : best_cost = 0;
4214 2568 : L_bits = log2(L);
4215 2568 : dbg_printf(3)("use_L1=%ld\n", use_L1);
4216 2568 : dbg_printf(3)("minbits = %ld\n", minbits);
4217 :
4218 2568 : totbits = 0;
4219 :
4220 : /* Iterate over the fundamental discriminants for L0 */
4221 1579091 : while (how_many_D0s < tablen) {
4222 1573955 : D0_entry = tab[how_many_D0s];
4223 1573955 : how_many_D0s++;
4224 : /* extract D0 from D0_entry */
4225 1573955 : D0 = D0_entry.D;
4226 :
4227 1573955 : if ( ! modinv_good_disc(inv, D0))
4228 600962 : continue;
4229 : /* extract class poly degree from D0_entry */
4230 972993 : deg = D0_entry.h;
4231 :
4232 : /* compute class number */
4233 972993 : h0 = ((D0_entry.m & 2) ? 2*deg : deg);
4234 972993 : dbg_printf(3)("D0=%ld\n", D0);
4235 :
4236 : /* don't allow either modinv_p1 or modinv_p2 to ramify */
4237 972993 : if (kross(D0, L) < 1
4238 440826 : || (modinv_p1 > 1 && kross(D0, modinv_p1) < 1)
4239 435560 : || (modinv_p2 > 1 && kross(D0, modinv_p2) < 1)) {
4240 547627 : dbg_printf(3)("Bad D0=%ld due to non-split L or ramified level\n", D0);
4241 547627 : continue;
4242 : }
4243 :
4244 : /* (D0_entry.m & 1) is 1 if ord(L0) < h0, is 0 if ord(L0) == h0.
4245 : * Hence (D0_entry.m & 1) + 1 is 2 if ord(L0) < h0 (hence ==
4246 : * h0/2) and is 1 if ord(L0) == h0. Hence n0 = ord(L0). */
4247 425366 : n0 = h0/((D0_entry.m & 1) + 1);
4248 :
4249 : /* Look for L1: for each smooth prime p */
4250 10332729 : for (i = 1 ; i <= SMOOTH_PRIMES; i++) {
4251 9992030 : if (PRIMES[i] <= L0)
4252 608097 : continue;
4253 : /* If 1 + (D0 | p) == 1, i.e. if (D0 | p) == 0, i.e. if p | D0. */
4254 9383933 : if (((D0_entry.m >> (2*i)) & 3) == 1) {
4255 : /* set L1 = p if it's not L0, it's less than the maximum,
4256 : * it doesn't divide modinv_N, and (L1 | L) == -1. */
4257 : /* XXX: Why do we want (L1 | L) == -1? Presumably so (L^2 v^2 D0 | L1) == -1? */
4258 303198 : L1 = PRIMES[i];
4259 303198 : if (L1 <= max_L1 && (modinv_N % L1) && kross(L1, L) < 0)
4260 84667 : break;
4261 : }
4262 : }
4263 : /* Didn't find suitable L1... */
4264 425366 : if (i > SMOOTH_PRIMES) {
4265 340699 : if (n0 < h0 || use_L1) {
4266 : /* ...though we needed one. */
4267 190055 : dbg_printf(3)("Bad D0=%ld because there is no good L1\n", D0);
4268 190055 : continue;
4269 : }
4270 150644 : L1 = 0;
4271 : }
4272 235311 : dbg_printf(3)("Good D0=%ld with L1=%ld, n0=%ld, h0=%ld, d=%ld\n",
4273 : D0, L1, n0, h0, d);
4274 :
4275 : /* We're finished if we have sufficiently many discriminants that satisfy
4276 : * the cost requirement */
4277 235311 : if (totbits > minbits && best_cost && h0*(L-1) > 3*best_cost)
4278 0 : break;
4279 :
4280 235311 : D0_bits = log2(-D0);
4281 : /* If L^2 D0 is too big to fit in a BIL bit integer, skip D0. */
4282 235311 : if (D0_bits + 2 * L_bits > (BITS_IN_LONG - 1))
4283 0 : continue;
4284 :
4285 : /* m is the order of L0^n0 in L^2 D0? */
4286 235311 : m = primeform_exp_order(L0, n0, L * L * D0, n0 * (L - 1));
4287 235311 : if ( m < (L - 1)/2 ) {
4288 67251 : dbg_printf(3)("Bad D0=%ld because %ld is less than (L-1)/2=%ld\n",
4289 0 : D0, m, (L - 1)/2);
4290 67251 : continue;
4291 : }
4292 : /* Heuristic. Doesn't end up contributing much. */
4293 168060 : H_cost = 2 * deg * deg;
4294 :
4295 : /* 7 = 2^3 - 1 = 0b111, so D0 & 7 == D0 (mod 8).
4296 : * 0xc = 0b1100, so D0_entry.m & 0xc == 1 + (D0 | 2).
4297 : * Altogether, we get:
4298 : * if D0 = 5 (mod 8), then
4299 : * if (D0 | 2) == -1, twofactor = 3
4300 : * otherwise (i.e. (D0 | 2) == 0 or 1), twofactor = 1
4301 : * otherwise
4302 : * twofactor is 0 */
4303 168060 : if ((D0 & 7) == 5)
4304 4917 : twofactor = ((D0_entry.m & 0xc) ? 1 : 3);
4305 : else
4306 163143 : twofactor = 0;
4307 :
4308 168060 : btop = avma;
4309 : /* For each small prime... */
4310 586663 : for (i = 0; i <= SMOOTH_PRIMES; i++) {
4311 586558 : avma = btop;
4312 : /* i = 0 corresponds to 1, which we do not want to skip! (i.e. DK = D) */
4313 586558 : if (i) {
4314 418498 : if (modinv_odd_conductor(inv) && i == 1)
4315 8866 : continue;
4316 409632 : p = PRIMES[i];
4317 : /* Don't allow large factors in the conductor. */
4318 409632 : if (p > max_L1)
4319 75901 : break;
4320 333731 : if (p == L0 || p == L1 || p == L || p == modinv_p1 || p == modinv_p2)
4321 115377 : continue;
4322 218354 : p_bits = log2(p);
4323 : /* h1 is the class number of D1 = q^2 D0, where q = p^j (j defined in the loop below) */
4324 218354 : h1 = h0 * (p - ((D0_entry.m >> (2*i)) & 0x3) + 1);
4325 : /* q is the smallest power of p such that h1 >= d ~ "L + 1". */
4326 218354 : for (j = 1, q = p; h1 < d; j++, q *= p, h1 *= p)
4327 : ;
4328 218354 : D1 = q * q * D0;
4329 : /* can't have D1 = 0 mod 16 and hope to get any primes congruent to 3 mod 4 */
4330 218354 : if ((pfilter & IQ_FILTER_1MOD4) && !(D1 & 0xF))
4331 70 : continue;
4332 : } else {
4333 : /* i = 0, corresponds to "p = 1". */
4334 168060 : h1 = h0;
4335 168060 : D1 = D0;
4336 168060 : p = q = j = 1;
4337 168060 : p_bits = 0;
4338 : }
4339 : /* include a factor of 4 if D1 is 5 mod 8 */
4340 : /* XXX: No idea why he does this. */
4341 386344 : if (twofactor && (q & 1)) {
4342 11869 : if (modinv_odd_conductor(inv))
4343 11351 : continue;
4344 518 : D1 *= 4;
4345 518 : h1 *= twofactor;
4346 : }
4347 : /* heuristic early abort -- we may miss some good D1's, but this saves a *lot* of time */
4348 374993 : if (totbits > minbits && best_cost && h1*(L-1) > 2.2*best_cost)
4349 8261 : continue;
4350 :
4351 : /* log2(D0 * (p^j)^2 * L^2 * twofactor) > (BIL - 1) -- i.e. params too big. */
4352 366732 : if (D0_bits + 2*j*p_bits + 2*L_bits + (twofactor && (q & 1) ? 2.0 : 0.0) > (BITS_IN_LONG - 1))
4353 1634 : continue;
4354 :
4355 365098 : if ( ! check_generators(&n1, NULL, D1, h1, n0, d, L0, L1))
4356 14053 : continue;
4357 :
4358 351045 : if (n1 < h1) {
4359 156938 : if ( ! primeform_discrete_log(&dl1, L0, L, n1, D1))
4360 98544 : continue;
4361 : } else {
4362 194107 : dl1 = -1; /* fill it in later, no point in wasting the effort right now */
4363 : }
4364 252501 : dbg_printf(3)("Good D0=%ld, D1=%ld with q=%ld, L1=%ld, n1=%ld, h1=%ld\n",
4365 : D0, D1, q, L1, n1, h1);
4366 252501 : if (modinv_deg && orientation_ambiguity(D1, L0, modinv_p1, modinv_p2, modinv_N))
4367 1424 : continue;
4368 :
4369 251077 : D2 = L * L * D1;
4370 251077 : h2 = h1 * (L-1);
4371 : /* m is the order of L0^n1 in cl(D2) */
4372 251077 : if ( ! check_generators(&n2, &m, D2, h2, n1, d*(L - 1), L0, L1))
4373 12868 : continue;
4374 :
4375 : /* This restriction on m is not strictly necessary, but simplifies life later. */
4376 238209 : if (m < (L-1)/2 || (!L1 && m < L-1) ) {
4377 116682 : dbg_printf(3)("Bad D2=%ld for D1=%ld, D0=%ld, with n2=%ld, h2=%ld, L1=%ld, "
4378 : "order of L0^n1 in cl(D2) is too small\n", D2, D1, D0, n2, h2, L1);
4379 116682 : continue;
4380 : }
4381 121527 : dl2[0] = n1;
4382 121527 : dl2[1] = 0;
4383 121527 : if (m < L - 1) {
4384 57879 : GEN Q1 = qform_primeform2(L, D1), Q2, X;
4385 57879 : if ( ! Q1)
4386 0 : pari_err_BUG("modpoly_pickD");
4387 57879 : Q2 = primeform_u(stoi(D2), L1);
4388 57879 : Q2 = gmul(Q1, Q2); /* we know this element has order L-1 */
4389 57879 : Q1 = primeform_u(stoi(D2), L0);
4390 57879 : k = ((n2 & 1) ? 2*n2 : n2)/(L-1);
4391 57879 : Q1 = gpowgs(Q1, k);
4392 57879 : X = discrete_log(Q2, Q1, L - 1);
4393 57879 : if ( ! X) {
4394 8202 : dbg_printf(3)("Bad D2=%ld for D1=%ld, D0=%ld, with n2=%ld, h2=%ld, L1=%ld, "
4395 : "form of norm L^2 not generated by L0 and L1\n",
4396 : D2, D1, D0, n2, h2, L1);
4397 8202 : continue;
4398 : }
4399 49677 : dl2[0] = itos(X) * k;
4400 49677 : dl2[1] = 1;
4401 : }
4402 113325 : if ( ! (m < L-1 || n2 < d*(L-1)) && n1 >= d && ! use_L1)
4403 63296 : L1 = 0; /* we don't need L1 */
4404 :
4405 113325 : if ( ! L1 && use_L1) {
4406 0 : dbg_printf(3)("not using D2=%ld for D1=%ld, D0=%ld, with n2=%ld, h2=%ld, L1=%ld, "
4407 : "because we don't need L1 but must use it\n",
4408 : D2, D1, D0, n2, h2, L1);
4409 0 : continue;
4410 : }
4411 : /* don't allow zero dl2[1] with L1 for the moment, since
4412 : * modpoly doesn't handle it - we may change this in the future */
4413 113325 : if (L1 && ! dl2[1])
4414 352 : continue;
4415 112973 : dbg_printf(3)("Good D0=%ld, D1=%ld, D2=%ld with s=%ld^%ld, L1=%ld, dl2=%ld, n2=%ld, h2=%ld\n",
4416 : D0, D1, D2, p, j, L1, dl2[0], n2, h2);
4417 :
4418 : /* This is estimate is heuristic and fiddling with the
4419 : * parameters 5 and 0.25 can change things quite a bit. */
4420 112973 : enum_cost = n2 * (5 * L0 * L0 + 0.25 * L1 * L1);
4421 112973 : cost = enum_cost + H_cost;
4422 112973 : if (best_cost && cost > 2.2*best_cost)
4423 85721 : break;
4424 27252 : if (best_cost && cost >= 0.99*best_cost)
4425 20492 : continue;
4426 :
4427 6760 : Dinfo.L = L;
4428 6760 : Dinfo.D0 = D0;
4429 6760 : Dinfo.D1 = D1;
4430 6760 : Dinfo.L0 = L0;
4431 6760 : Dinfo.L1 = L1;
4432 6760 : Dinfo.n1 = n1;
4433 6760 : Dinfo.n2 = n2;
4434 6760 : Dinfo.dl1 = dl1;
4435 6760 : Dinfo.dl2_0 = dl2[0];
4436 6760 : Dinfo.dl2_1 = dl2[1];
4437 6760 : Dinfo.cost = cost;
4438 6760 : Dinfo.inv = inv;
4439 :
4440 6760 : if ( ! modpoly_pickD_primes (NULL, NULL, 0, NULL, 0, &Dinfo.bits, minbits, &Dinfo))
4441 44 : continue;
4442 6716 : dbg_printf(2)("Best D2=%ld, D1=%ld, D0=%ld with s=%ld^%ld, L1=%ld, "
4443 : "n1=%ld, n2=%ld, cost ratio %.2f, bits=%ld\n",
4444 : D2, D1, D0, p, j, L1, n1, n2,
4445 0 : (double)cost/(d*(L-1)), Dinfo.bits);
4446 : /* Insert Dinfo into the Ds array. Ds is sorted by ascending cost. */
4447 33870 : for (j = 0; j < Dcnt; j++)
4448 31293 : if (Dinfo.cost < Ds[j].cost)
4449 4139 : break;
4450 6716 : if (n2 > MAX_VOLCANO_FLOOR_SIZE && n2*(L1 ? 2 : 1) > 1.2* (d*(L-1)) ) {
4451 0 : dbg_printf(3)("Not using D1=%ld, D2=%ld for space reasons\n", D1, D2);
4452 0 : continue;
4453 : }
4454 6716 : if (j == Dcnt && Dcnt == MODPOLY_MAX_DCNT)
4455 0 : continue;
4456 6716 : totbits += Dinfo.bits;
4457 6716 : if (Dcnt == MODPOLY_MAX_DCNT)
4458 0 : totbits -= Ds[Dcnt-1].bits;
4459 6716 : if (n2 > MAX_VOLCANO_FLOOR_SIZE)
4460 0 : dbg_printf(3)("totbits=%ld, minbits=%ld\n", totbits, minbits);
4461 6716 : if (Dcnt < MODPOLY_MAX_DCNT)
4462 6716 : Dcnt++;
4463 14409 : for (k = Dcnt - 1; k > j; k--)
4464 7693 : Ds[k] = Ds[k - 1];
4465 6716 : Ds[k] = Dinfo;
4466 6716 : if (totbits > minbits)
4467 6700 : best_cost = Ds[Dcnt-1].cost;
4468 : else
4469 16 : best_cost = 0;
4470 : /* if we were able to use D1 with s = 1, there is no point in
4471 : * using any larger D1 for the same D0 */
4472 6716 : if ( ! i)
4473 6333 : break;
4474 : } /* END FOR over small primes */
4475 : } /* END WHILE over D0's */
4476 2568 : dbg_printf(2)(" checked %ld of %ld fundamental discriminants to find suitable "
4477 : "discriminant (Dcnt = %ld)\n", how_many_D0s, tablen, Dcnt);
4478 2568 : if ( ! Dcnt) {
4479 0 : dbg_printf(1)("failed completely for L=%ld\n", L);
4480 0 : return 0;
4481 : }
4482 :
4483 2568 : Dcnt = calc_primes_for_discriminants(Ds, Dcnt, L, minbits);
4484 :
4485 : /* fill in any missing dl1's */
4486 5152 : for (i = 0 ; i < Dcnt; i++) {
4487 2584 : if (Ds[i].dl1 < 0) {
4488 : long dl;
4489 2539 : if ( ! primeform_discrete_log(&dl, L0, L, Ds[i].n1, Ds[i].D1))
4490 : {
4491 0 : pari_err_BUG("modpoly_pickD");
4492 : return -1; /* LCOV_EXCL_LINE */
4493 : }
4494 2539 : Ds[i].dl1 = dl;
4495 : }
4496 : }
4497 2568 : if (DEBUGLEVEL > 1+3) {
4498 0 : err_printf("Selected %ld discriminants using %ld msecs\n", Dcnt, timer_delay(&T));
4499 0 : for (i = 0 ; i < Dcnt ; i++) {
4500 : /* TODO: Reuse the calculation from the D_entry */
4501 0 : GEN H = classno(stoi(Ds[i].D0));
4502 0 : long h0 = itos(H);
4503 0 : err_printf (" D0=%ld, h(D0)=%ld, D=%ld, L0=%ld, L1=%ld, "
4504 : "cost ratio=%.2f, enum ratio=%.2f,",
4505 0 : Ds[i].D0, h0, Ds[i].D1, Ds[i].L0, Ds[i].L1,
4506 0 : (double)Ds[i].cost/(d*(L-1)),
4507 0 : (double)(Ds[i].n2*(Ds[i].L1 ? 2 : 1))/(d*(L-1)));
4508 0 : err_printf (" %ld primes, %ld bits\n", Ds[i].nprimes, Ds[i].bits);
4509 : }
4510 : }
4511 2568 : avma = ltop;
4512 2568 : return Dcnt;
4513 : }
4514 :
4515 : static int
4516 12283032 : _qsort_cmp(const void *a, const void *b)
4517 : {
4518 : D_entry *x, *y;
4519 : long u, v;
4520 :
4521 12283032 : x = (D_entry *)a;
4522 12283032 : y = (D_entry *)b;
4523 : /* u and v are the class numbers of x and y */
4524 12283032 : u = x->h * (!!(x->m & 2) + 1);
4525 12283032 : v = y->h * (!!(y->m & 2) + 1);
4526 : /* Sort by class number */
4527 12283032 : if (u < v)
4528 3734920 : return -1;
4529 8548112 : if (u > v)
4530 5985763 : return 1;
4531 : /* Sort by discriminant (which is < 0, hence the sign reversal) */
4532 2562349 : if (x->D > y->D)
4533 2562349 : return -1;
4534 0 : if (x->D < y->D)
4535 0 : return 1;
4536 0 : return 0;
4537 : }
4538 :
4539 : /*
4540 : * Build a table containing fundamental discriminants less than maxd
4541 : * whose class groups
4542 : * - are cyclic generated by an element of norm L0
4543 : * - have class number at most maxh
4544 : * The table is ordered using _qsort_cmp above, which ranks the
4545 : * discriminants by class number, then by absolute discriminant.
4546 : *
4547 : * INPUT:
4548 : * - maxd: largest allowed discriminant
4549 : * - maxh: largest allowed class number
4550 : * - L0: norm of class group generator
4551 : *
4552 : * OUTPUT:
4553 : * - tablelen: length of return value
4554 : *
4555 : * RETURN:
4556 : * - array of {discriminant D, h(D), kronecker symbols for small p} where D
4557 : * satisfies the properties above
4558 : */
4559 : static D_entry *
4560 2568 : scanD0(long *tablelen, long *minD, long maxD, long maxh, long L0)
4561 : {
4562 : GEN fact, DD, H, ordL, frm;
4563 : pari_sp av;
4564 : long *q, *e;
4565 : D_entry *tab;
4566 : ulong m, x;
4567 : long h, d, D, n;
4568 : long L1, cnt, i, j, k;
4569 :
4570 2568 : if (maxD < 0)
4571 0 : maxD = -maxD;
4572 :
4573 : /* NB: As seen in the loop below, the real class number of D can be */
4574 : /* 2*maxh if cl(D) is cyclic. */
4575 2568 : if (maxh < 0)
4576 0 : pari_err_BUG("scanD0");
4577 :
4578 : /* Not checked, but L0 should be 2, 3, 5 or 7. */
4579 :
4580 2568 : tab = (D_entry *) stack_malloc((maxD/4)*sizeof(*tab)); /* Overestimate */
4581 2568 : cnt = 0;
4582 2568 : av = avma;
4583 :
4584 : /* d = 7, 11, 15, 19, 23, ... */
4585 6420000 : for (d = *minD, cnt = 0; d <= maxD; d += 4) {
4586 6417432 : D = -d;
4587 : /* Check to see if (D | L0) = 1 */
4588 6417432 : if (kross(D, L0) < 1)
4589 3525594 : continue;
4590 :
4591 : /* [q, e] is the factorisation of d. */
4592 2891838 : fact = factoru(d);
4593 2891838 : q = zv_to_longptr(gel(fact, 1));
4594 2891838 : e = zv_to_longptr(gel(fact, 2));
4595 2891838 : k = lg(gel(fact, 1)) - 1;
4596 :
4597 : /* Check if the discriminant is square-free */
4598 7480639 : for (i = 0; i < k; i++)
4599 5075769 : if (e[i] > 1)
4600 486968 : break;
4601 2891838 : if (i < k)
4602 486968 : continue;
4603 :
4604 : /* L1 is initially the first factor of d if small enough, otherwise ignored. */
4605 2404870 : if (k > 1 && q[0] <= MAX_L1)
4606 1615008 : L1 = q[0];
4607 : else
4608 789862 : L1 = 0;
4609 :
4610 : /* restrict to possibly cyclic class groups */
4611 2404870 : if (k > 2)
4612 455439 : continue;
4613 :
4614 : /* Check if h(D) is too big */
4615 1949431 : DD = stoi(D);
4616 1949431 : H = classno(DD);
4617 1949431 : h = itos(H);
4618 1949431 : if (h > 2*maxh || (!L1 && h > maxh))
4619 113133 : continue;
4620 :
4621 : /* Check if ord(q) is not big enough to generate at least half the
4622 : * class group (where q is the L0-primeform). */
4623 1836298 : frm = primeform_u(DD, L0);
4624 1836298 : ordL = qfi_order(redimag(frm), H);
4625 1836298 : n = itos(ordL);
4626 1836298 : if (n < h/2 || (!L1 && n < h))
4627 262343 : continue;
4628 :
4629 : /* If q is big enough, great! Otherwise, for each potential L1,
4630 : * do a discrete log to see if it is NOT in the subgroup generated
4631 : * by L0; stop as soon as such is found. */
4632 1780273 : for (j = 0; ; j++) {
4633 1780273 : if (n == h || (L1 && ! discrete_log(primeform_u(DD, L1), frm, n))) {
4634 1492798 : dbg_printf(2)("D0=%ld good with L1=%ld\n", D, L1);
4635 1492798 : break;
4636 : }
4637 287475 : if ( ! L1) break;
4638 206318 : L1 = (j < k && k > 1 && q[j] <= MAX_L1 ? q[j] : 0);
4639 206318 : }
4640 :
4641 : /* NB: After all that, L1 is not used or saved for later. */
4642 :
4643 : /* The first bit of m indicates whether q generates a proper
4644 : * subgroup of cl(D) (hence implying that we need L1) or if q
4645 : * generates the whole class group. */
4646 1573955 : m = (n < h ? 1 : 0);
4647 : /* bits i and i+1 of m give the 2-bit number 1 + (D|p) where p is
4648 : * the ith prime. */
4649 46655696 : for (i = 1 ; i <= ((BITS_IN_LONG >> 1) - 1); i++) {
4650 45081741 : x = (ulong) (1 + kross(D, PRIMES[i]));
4651 45081741 : m |= x << (2*i);
4652 : }
4653 :
4654 : /* Insert d, h and m into the table */
4655 1573955 : tab[cnt].D = D;
4656 1573955 : tab[cnt].h = h;
4657 1573955 : tab[cnt].m = m;
4658 1573955 : cnt++;
4659 1573955 : avma = av;
4660 : }
4661 :
4662 : /* Sort the table */
4663 2568 : qsort(tab, cnt, sizeof(*tab), _qsort_cmp);
4664 2568 : *tablelen = cnt;
4665 2568 : *minD = d;
4666 2568 : return tab;
4667 : }
4668 :
4669 :
4670 : /*
4671 : * Populate Ds with discriminants (and attached data) that can be
4672 : * used to calculate the modular polynomial of level L and invariant
4673 : * inv. Return the number of discriminants found.
4674 : */
4675 : static long
4676 2566 : discriminant_with_classno_at_least(
4677 : modpoly_disc_info Ds[MODPOLY_MAX_DCNT], long L, long inv, long ignore_sparse)
4678 : {
4679 : enum { SMALL_L_BOUND = 101 };
4680 2566 : long max_max_D = 160000 * (inv ? 2 : 1);
4681 : long minD, maxD, maxh, L0, max_L1, minbits, Dcnt, flags, s, d, h, i, tmp;
4682 : D_entry *tab;
4683 : long tablen;
4684 2566 : pari_sp av = avma;
4685 2566 : double eps, best_eps = -1.0, cost, best_cost = -1.0;
4686 : modpoly_disc_info bestD[MODPOLY_MAX_DCNT];
4687 2566 : long best_cnt = 0;
4688 : pari_timer T;
4689 2566 : timer_start(&T);
4690 :
4691 2566 : s = modinv_sparse_factor(inv);
4692 2566 : d = s;
4693 : /*d = ui_ceil_ratio(L + 1, d) + 1; */
4694 2566 : tmp = (L + 1) / d;
4695 2566 : d = ((tmp * d < (L + 1)) ? tmp + 1 : tmp);
4696 2566 : d += 1;
4697 :
4698 : /* maxD of 10000 allows us to get a satisfactory discriminant in
4699 : * under 250ms in most cases. */
4700 2566 : maxD = 10000;
4701 : /* Allow the class number to overshoot L by 50%. Must be at least
4702 : * 1.1*L, and higher values don't seem to provide much benefit,
4703 : * except when L is small, in which case it's necessary to get any
4704 : * discriminant at all in some cases. */
4705 2566 : maxh = (L / s < SMALL_L_BOUND) ? 10 * L : 1.5 * L;
4706 :
4707 2566 : flags = ignore_sparse ? MODPOLY_IGNORE_SPARSE_FACTOR : 0;
4708 2566 : L0 = select_L0(L, inv, 0);
4709 2566 : max_L1 = L / 2 + 2; /* for L=11 we need L1=7 for j */
4710 2566 : minbits = modpoly_height_bound(L, inv);
4711 2566 : minD = 7;
4712 :
4713 7698 : while ( ! best_cnt) {
4714 5134 : while (maxD <= max_max_D) {
4715 : /* TODO: Find a way to re-use tab when we need multiple modpolys */
4716 2568 : tab = scanD0(&tablen, &minD, maxD, maxh, L0);
4717 2568 : dbg_printf(1)("Found %ld potential fundamental discriminants\n", tablen);
4718 :
4719 2568 : Dcnt = modpoly_pickD(Ds, L, inv, L0, max_L1, minbits, flags, tab, tablen);
4720 2568 : eps = 0.0;
4721 2568 : cost = 0.0;
4722 :
4723 2568 : if (Dcnt) {
4724 2566 : long n1 = 0;
4725 5150 : for (i = 0; i < Dcnt; ++i) {
4726 2584 : n1 = maxss(n1, Ds[i].n1);
4727 2584 : cost += Ds[i].cost;
4728 : }
4729 2566 : eps = (n1 * s - L) / (double)L;
4730 :
4731 2566 : if (best_cost < 0.0 || cost < best_cost) {
4732 2566 : (void) memcpy(bestD, Ds, Dcnt * sizeof(modpoly_disc_info));
4733 2566 : best_cost = cost;
4734 2566 : best_cnt = Dcnt;
4735 2566 : best_eps = eps;
4736 : /* We're satisfied if n1 is within 5% of L. */
4737 2566 : if (L / s <= SMALL_L_BOUND || eps < 0.05)
4738 : break;
4739 : }
4740 : } else {
4741 2 : if (log2(maxD) > BITS_IN_LONG - 2 * (log2(L) + 2))
4742 0 : pari_err(e_ARCH, "modular polynomial of given level and invariant");
4743 : }
4744 2 : maxD *= 2;
4745 2 : minD += 4;
4746 2 : if (DEBUGLEVEL > 3) {
4747 0 : err_printf(" Doubling discriminant search space (closest: %.1f%%, cost ratio: %.1f)...\n",
4748 0 : eps*100, cost/(double)(d*(L-1)));
4749 : }
4750 : }
4751 2566 : max_max_D *= 2;
4752 : }
4753 :
4754 2566 : if (DEBUGLEVEL > 3) {
4755 0 : err_printf("Found discriminant(s):\n");
4756 0 : for (i = 0; i < best_cnt; ++i) {
4757 0 : av = avma;
4758 0 : h = itos(classno(stoi(bestD[i].D1)));
4759 0 : avma = av;
4760 0 : err_printf(" D = %ld, h = %ld, u = %ld, L0 = %ld, L1 = %ld, n1 = %ld, n2 = %ld, cost = %ld\n",
4761 0 : bestD[i].D1, h, (long)sqrt((double)(bestD[i].D1 / bestD[i].D0)), bestD[i].L0, bestD[i].L1,
4762 : bestD[i].n1, bestD[i].n2, bestD[i].cost);
4763 : }
4764 0 : err_printf("(off target by %.1f%%, cost ratio: %.1f)\n",
4765 0 : best_eps*100, best_cost/(double)(d*(L-1)));
4766 : }
4767 2566 : return best_cnt;
4768 : }
|
