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