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