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