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