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