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