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 : /* FIXME: Implement {ascend,descend}_volcano() in terms of the "new"
18 : * volcano traversal functions at the bottom of the file. */
19 :
20 : /* Is j = 0 or 1728 (mod p)? */
21 : INLINE int
22 315229 : is_j_exceptional(ulong j, ulong p)
23 : {
24 315229 : return j == 0 || j == 1728 % p;
25 : }
26 :
27 :
28 : INLINE long
29 70483 : node_degree(GEN phi, long L, ulong j, ulong p, ulong pi)
30 : {
31 70483 : pari_sp av = avma;
32 70483 : long n = Flx_nbroots(Flm_Fl_polmodular_evalx(phi, L, j, p, pi), p);
33 70483 : avma = av;
34 70483 : return n;
35 : }
36 :
37 :
38 : /*
39 : * Given an array path = [j0, j1] of length 2, return the polynomial
40 : *
41 : * \Phi_L(X, j1) / (X - j0)
42 : *
43 : * where \Phi_L(X, Y) is the modular polynomial of level L. An error
44 : * is raised if X - j0 does not divide \Phi_L(X, j1).
45 : */
46 : INLINE GEN
47 124325 : nhbr_polynomial(ulong path[], GEN phi, ulong p, ulong pi, long L)
48 : {
49 124325 : pari_sp ltop = avma;
50 124325 : GEN modpol = Flm_Fl_polmodular_evalx(phi, L, path[0], p, pi);
51 :
52 : /* Note that, if the discriminant of End(path[0]) is less than L^2,
53 : * then it's possible for path[0] to appear among the roots of
54 : * nhbr_pol. This case should have been obviated by earlier
55 : * choices. */
56 124325 : ulong rem = 0;
57 124325 : GEN nhbr_pol = Flx_div_by_X_x(modpol, path[-1], p, &rem);
58 124325 : if (rem)
59 0 : pari_err_BUG("nhbr_polynomial: invalid preceding j");
60 124325 : return gerepileupto(ltop, nhbr_pol);
61 : }
62 :
63 :
64 : /*
65 : * Assumes path is an array with space for at least max_len + 1
66 : * elements, whose first and second elements are the beginning of the
67 : * path. I.e., the path starts
68 : *
69 : * (path[0], path[1])
70 : *
71 : * If the result is less than max_len, then the last element of path
72 : * is definitely on the floor. If the result equals max_len, then in
73 : * general it is unknown whether the last element of path is on the
74 : * floor or not.
75 : */
76 : static long
77 244545 : extend_path(
78 : ulong path[], GEN phi, ulong p, ulong pi, long L, long max_len)
79 : {
80 244545 : pari_sp av = avma;
81 244545 : long d = 1;
82 314980 : for ( ; d < max_len; ++d) {
83 : ulong nhbr;
84 : GEN nhbr_pol;
85 :
86 91081 : nhbr_pol = nhbr_polynomial(path + d, phi, p, pi, L);
87 91081 : nhbr = Flx_oneroot(nhbr_pol, p);
88 91081 : avma = av;
89 : /* Flx_oneroot didn't find a root; so we must be on the floor. */
90 91081 : if (nhbr == p)
91 20646 : break;
92 70435 : path[d + 1] = nhbr;
93 : }
94 244545 : return d;
95 : }
96 :
97 :
98 : /*
99 : * This is Sutherland 2009 Algorithm Ascend (p12).
100 : */
101 : ulong
102 110072 : ascend_volcano(
103 : GEN phi, ulong j, ulong p, ulong pi, long level, long L,
104 : long depth, long steps)
105 : {
106 110072 : pari_sp ltop = avma, av;
107 : /* path will never hold more than max_len elements, and max_len <
108 : * depth always. */
109 110072 : GEN path_g = cgetg(depth + 2, t_VECSMALL);
110 110072 : ulong *path = zv_to_ulongptr(path_g);
111 110072 : long max_len = depth - level;
112 110072 : int first_iter = 1;
113 :
114 110072 : if (steps <= 0 || max_len < 0)
115 0 : pari_err_BUG("ascend_volcano: bad params");
116 :
117 110072 : av = avma;
118 363460 : while (steps--) {
119 143316 : GEN nhbr_pol = first_iter
120 : ? Flm_Fl_polmodular_evalx(phi, L, j, p, pi)
121 143316 : : nhbr_polynomial(path + 1, phi, p, pi, L);
122 143316 : GEN nhbrs = Flx_roots(nhbr_pol, p);
123 143316 : long nhbrs_len = lg(nhbrs) - 1, i;
124 143316 : pari_sp btop = avma;
125 143316 : path[0] = j;
126 143316 : first_iter = 0;
127 :
128 143316 : j = nhbrs[nhbrs_len];
129 181637 : for (i = 1; i < nhbrs_len; ++i) {
130 68554 : ulong next_j = nhbrs[i], last_j;
131 : long len;
132 68554 : if (is_j_exceptional(next_j, p)) {
133 : /* According to Fouquet & Morain, Section 4.3, if j = 0 or
134 : * 1728, then it is necessarily on the surface. So we just
135 : * return it. */
136 0 : if (steps) {
137 0 : pari_err_BUG("ascend_volcano: "
138 : "Got to the top with more steps to go!");
139 : }
140 0 : j = next_j;
141 0 : break;
142 : }
143 68554 : path[1] = next_j;
144 68554 : len = extend_path(path, phi, p, pi, L, max_len);
145 68554 : last_j = path[len];
146 68554 : if (len == max_len
147 : /* Ended up on the surface */
148 68554 : && (is_j_exceptional(last_j, p)
149 68554 : || node_degree(phi, L, last_j, p, pi) > 1)) {
150 30233 : j = next_j;
151 30233 : break;
152 : }
153 38321 : avma = btop;
154 : }
155 143316 : path[1] = j; /* For nhbr_polynomial() at the top. */
156 :
157 143316 : ++max_len;
158 143316 : avma = av;
159 : }
160 110072 : avma = ltop;
161 110072 : return j;
162 : }
163 :
164 :
165 : static void
166 178121 : random_distinct_neighbours_of(
167 : ulong *nhbr1, ulong *nhbr2,
168 : GEN phi, ulong j, ulong p, ulong pi, long L,
169 : long must_have_two_neighbours)
170 : {
171 178121 : pari_sp ltop = avma;
172 178121 : GEN modpol = Flm_Fl_polmodular_evalx(phi, L, j, p, pi);
173 : ulong rem;
174 178121 : *nhbr1 = Flx_oneroot(modpol, p);
175 178121 : if (*nhbr1 == p) {
176 : /* Didn't even find one root! */
177 0 : char *err = stack_sprintf("random_distinct_neighbours_of: "
178 : "No neighbours for j = %lu (mod %lu) "
179 : "in %lu-volcano.", j, p, L);
180 0 : pari_err_BUG(err);
181 : }
182 178121 : modpol = Flx_div_by_X_x(modpol, *nhbr1, p, &rem);
183 178121 : *nhbr2 = Flx_oneroot(modpol, p);
184 178121 : if (must_have_two_neighbours && *nhbr2 == p) {
185 : /* Didn't find distinct root! */
186 0 : char *err = stack_sprintf("random_distinct_neighbours_of: "
187 : "Only one neighbour for j = %lu (mod %lu) "
188 : "in %lu-volcano.", j, p, L);
189 0 : pari_err_BUG(err);
190 : }
191 178121 : avma = ltop;
192 178121 : }
193 :
194 :
195 : /*
196 : * This is Sutherland 2009 Algorithm Descend (p12).
197 : */
198 : ulong
199 1810 : descend_volcano(
200 : GEN phi, ulong j, ulong p, ulong pi,
201 : long level, long L, long depth, long steps)
202 : {
203 1810 : pari_sp ltop = avma;
204 : GEN path_g;
205 : ulong *path, res;
206 : long max_len;
207 :
208 1810 : if (steps <= 0 || level + steps > depth)
209 0 : pari_err_BUG("descend_volcano: bad params");
210 :
211 1810 : max_len = depth - level;
212 1810 : path_g = cgetg(max_len + 1 + 1, t_VECSMALL);
213 1810 : path = zv_to_ulongptr(path_g);
214 1810 : path[0] = j;
215 : /* level = 0 means we're on the volcano surface... */
216 1810 : if ( ! level) {
217 : /* Look for any path to the floor. One of j's first three
218 : * neighbours must lead to the floor, since at most two neighbours
219 : * are on the surface. */
220 1616 : GEN nhbrs = Flx_roots(Flm_Fl_polmodular_evalx(phi, L, j, p, pi), p);
221 : long i;
222 1761 : for (i = 1; i <= 3; ++i) {
223 : long len;
224 1761 : path[1] = nhbrs[i];
225 1761 : len = extend_path(path, phi, p, pi, L, max_len);
226 : /* If nhbrs[i] took us to the floor: */
227 1761 : if (len < max_len || node_degree(phi, L, path[len], p, pi) == 1)
228 : break;
229 : }
230 :
231 1616 : if (i > 3) {
232 0 : pari_err_BUG("descend_volcano: "
233 : "None of three neighbours lead to the floor");
234 : }
235 : } else {
236 : ulong nhbr1, nhbr2;
237 : long len;
238 194 : random_distinct_neighbours_of(&nhbr1, &nhbr2, phi, j, p, pi, L, 1);
239 194 : path[1] = nhbr1;
240 194 : len = extend_path(path, phi, p, pi, L, max_len);
241 : /* If last j isn't on the floor */
242 194 : if (len == max_len /* Ended up on the surface. */
243 194 : && (is_j_exceptional(path[len], p)
244 168 : || node_degree(phi, L, path[len], p, pi) != 1)) {
245 : /* The other neighbour leads to the floor */
246 86 : path[1] = nhbr2;
247 86 : (void) extend_path(path, phi, p, pi, L, steps);
248 : }
249 : }
250 1810 : res = path[steps];
251 1810 : avma = ltop;
252 1810 : return res;
253 : }
254 :
255 :
256 : long
257 177927 : j_level_in_volcano(
258 : GEN phi, ulong j, ulong p, ulong pi, long L, long depth)
259 : {
260 177927 : pari_sp av = avma;
261 : GEN chunk;
262 : ulong *path1, *path2;
263 : long lvl;
264 :
265 177927 : if (depth == 0 || is_j_exceptional(j, p)) {
266 : /* According to Fouquet & Morain, Section 4.3, if j = 0 or 1728,
267 : * then it is necessarily on the surface. Also, if the volcano
268 : * depth is zero then j necessarily has level 0. */
269 0 : return 0;
270 : }
271 :
272 177927 : chunk = new_chunk(2 * (depth + 1));
273 177927 : path1 = (ulong *) &chunk[0];
274 177927 : path2 = (ulong *) &chunk[depth + 1];
275 :
276 177927 : path1[0] = path2[0] = j;
277 :
278 177927 : random_distinct_neighbours_of(&path1[1], &path2[1],
279 : phi, j, p, pi, L, 0);
280 177927 : if (path2[1] == p) {
281 : /* Only found one neighbour, hence j is on the floor, hence
282 : * level == depth. */
283 90952 : lvl = depth;
284 : } else {
285 86975 : long path1_len = extend_path(path1, phi, p, pi, L, depth);
286 86975 : long path2_len = extend_path(path2, phi, p, pi, L, path1_len);
287 86975 : lvl = depth - path2_len;
288 : }
289 177927 : avma = av;
290 177927 : return lvl;
291 : }
292 :
293 :
294 : #define vecsmall_len(v) (lg(v) - 1)
295 :
296 : INLINE GEN
297 28407507 : Flx_remove_root(GEN f, ulong a, ulong p)
298 : {
299 : ulong r;
300 28407507 : GEN g = Flx_div_by_X_x(f, a, p, &r);
301 28424587 : if (r) pari_err_BUG("Flx_remove_root");
302 28432562 : return g;
303 : }
304 :
305 : INLINE GEN
306 21842508 : get_nbrs(GEN phi, long L, ulong J, const ulong *xJ, ulong p, ulong pi)
307 : {
308 21842508 : pari_sp av = avma;
309 21842508 : GEN f = Flm_Fl_polmodular_evalx(phi, L, J, p, pi);
310 21827574 : if (xJ)
311 21580817 : f = Flx_remove_root(f, *xJ, p);
312 21830695 : return gerepileupto(av, Flx_roots(f, p));
313 : }
314 :
315 : /* Return a path of length n along the surface of an L-volcano of
316 : * height h starting from surface node j0. Assumes (D|L) = 1 where D
317 : * is the discriminant of End(j0).
318 : *
319 : * Actually, if j0's endomorphism ring is a suborder, we return
320 : * the corresponding shorter path.
321 : *
322 : * W must hold space for n + h nodes.
323 : *
324 : * TODO: It could be worth having two versions of this function: one
325 : * that assumes J has the correct endomorphism ring (hence avoiding
326 : * several branches in the inner loop) and a second that doesn't
327 : * assume J has the correct endo ring and accordingly checks for
328 : * repetitions. */
329 : static long
330 181864 : surface_path(
331 : ulong W[],
332 : long n, GEN phi, long L, long h, ulong J, const ulong *nJ, ulong p, ulong pi)
333 : {
334 181864 : pari_sp av = avma, bv;
335 : GEN T, v;
336 : long j, k, w, x;
337 : ulong W0;
338 :
339 181864 : W[0] = W0 = J;
340 181864 : if (n == 1)
341 0 : return 1;
342 :
343 181864 : T = cgetg(h + 2, t_VEC);
344 181864 : bv = avma;
345 :
346 181864 : v = get_nbrs(phi, L, J, nJ, p, pi);
347 : /* Insert known neighbour first */
348 181863 : if (nJ)
349 0 : v = gerepileupto(bv, vecsmall_prepend(v, *nJ));
350 181865 : gel(T, 1) = v;
351 181865 : k = vecsmall_len(v);
352 :
353 181865 : switch (k) {
354 : case 0:
355 : /* We must always have neighbours */
356 0 : pari_err_BUG("surface_path");
357 : case 1:
358 : /* If volcano is not flat, then we must have more than one
359 : * neighbour */
360 5247 : if (h) pari_err_BUG("surface_path");
361 5247 : W[1] = uel(v, 1);
362 5247 : avma = av;
363 : /* Check for bad endo ring */
364 5247 : if (W[1] == W[0])
365 79 : return 1;
366 5168 : return 2;
367 : case 2:
368 : /* If L=2 the only way we can have 2 neighbours is if we have a
369 : * double root. This can only happen for |D| <= 16 (by Thm 2.2 of
370 : * Fouquet-Morain), and if it does we must have a 2-cycle. This
371 : * case happens for D=-15. */
372 20683 : if (L == 2) {
373 : /* The double root is the neighbour on the surface, who will
374 : * have exactly one neighbour other than J; the other neighbour
375 : * of J has either 0 or 2 neighbours that are not J. */
376 0 : GEN u = get_nbrs(phi, L, uel(v, 1), &J, p, pi);
377 0 : long n = vecsmall_len(u) - !!vecsmall_isin(u, J);
378 0 : W[1] = n == 1 ? uel(v, 1) : uel(v, 2);
379 0 : avma = av;
380 0 : return 2;
381 : }
382 :
383 : /* Volcano is not flat, but we only found 2 neighbours for the
384 : * surface node J */
385 20683 : if (h) pari_err_BUG("surface_path");
386 :
387 20861 : W[1] = uel(v, 1);
388 : /* TODO: Can we use the other root -- uel(v,2) -- somehow? */
389 4359714 : for (w = 2; w < n; ++w) {
390 4339135 : v = get_nbrs(phi, L, W[w - 1], &W[w - 2], p, pi);
391 : /* A flat volcano must have exactly one non-previous
392 : * neighbour */
393 4338485 : if (vecsmall_len(v) != 1) pari_err_BUG("surface_path");
394 4338957 : W[w] = uel(v, 1);
395 4338957 : avma = av;
396 : /* Detect cycle in case J doesn't have the right endo ring. */
397 4338957 : if (W[w] == W0)
398 104 : return w;
399 : }
400 20579 : avma = av;
401 20579 : return n;
402 : }
403 :
404 : /* Can't have a flat volcano if k > 2. */
405 155935 : if (!h) pari_err_BUG("surface_path");
406 :
407 : /* At this point, each surface node must have L + 1 distinct
408 : * neighbours, 2 of which are on the surface. */
409 156295 : w = 1;
410 5631411 : for (x = 0; ; ++x) {
411 : /* Get next neighbour of last known surface node to attempt to
412 : * extend the path. */
413 5631411 : W[w] = umael(T, ((w - 1) % h) + 1, x + 1);
414 :
415 : /* Detect cycle in case the given J didn't have the right endo
416 : * ring. */
417 5631411 : if (W[w] == W0) {
418 128 : avma = av;
419 128 : return w;
420 : }
421 :
422 : /* If we have to test the last neighbour, we know it's on the
423 : * surface, and if we're done there's no need to extend. */
424 5631283 : if (x == k - 1 && w == n - 1) {
425 102181 : avma = av;
426 102181 : return n;
427 : }
428 :
429 : /* Walk forward until we hit the floor or finish. */
430 : /* NB: We don't keep the stack completely clean here; usage is
431 : * determined by the size of T, which will be in the order of Lh,
432 : * i.e. L roots for each level of the volcano of height h. */
433 5529102 : for (j = w; ; ) {
434 : long m;
435 : /* We must get 0 or L neighbours here. */
436 17178239 : v = get_nbrs(phi, L, W[j], &W[j - 1], p, pi);
437 17172019 : m = vecsmall_len(v);
438 17172019 : if ( ! m) {
439 : /* We hit the floor; save the neighbours of W[w-1] and throw
440 : * away the rest. */
441 5473606 : GEN nbrs = gel(T, ((w - 1) % h) + 1);
442 5473606 : gel(T, ((w - 1) % h) + 1) = gerepileupto(bv, nbrs);
443 5475116 : break;
444 : }
445 11698413 : if (m != L) pari_err_BUG("surface_path");
446 :
447 11702767 : gel(T, (j % h) + 1) = v;
448 :
449 11702767 : W[++j] = uel(v, 1);
450 : /* If we have our path by h nodes, we know W[w] is on the
451 : * surface */
452 11702767 : if (j == w + h) {
453 10746701 : ++w;
454 : /* Detect cycle in case the given J didn't have the right endo
455 : * ring. */
456 10746701 : if (W[w] == W0) {
457 24477 : avma = av;
458 24477 : return w;
459 : }
460 10722224 : x = 0;
461 10722224 : k = L;
462 : }
463 11678290 : if (w == n) {
464 29153 : avma = av;
465 29153 : return w;
466 : }
467 11649137 : }
468 5475116 : }
469 : return -1; /*LCOV_EXCL_LINE*/
470 : }
471 :
472 : long
473 107795 : next_surface_nbr(
474 : ulong *nJ,
475 : GEN phi, long L, long h, ulong J, const ulong *pJ, ulong p, ulong pi)
476 : {
477 107795 : pari_sp av = avma, bv;
478 : GEN S;
479 : ulong *P;
480 : long i, k;
481 :
482 107795 : S = get_nbrs(phi, L, J, pJ, p, pi);
483 107794 : k = vecsmall_len(S);
484 : /* If there is a double root and pJ is set, then k will be zero. */
485 107794 : if ( ! k) {
486 0 : avma = av;
487 0 : return 0;
488 : }
489 107794 : if (k == 1 || ( ! pJ && k == 2)) {
490 89807 : *nJ = uel(S, 1);
491 89807 : avma = av;
492 89807 : return 1;
493 : }
494 17987 : if ( ! h) pari_err_BUG("next_surface_nbr");
495 :
496 17987 : P = (ulong *) new_chunk(h + 1);
497 17987 : P[0] = J;
498 17987 : bv = avma;
499 35226 : for (i = 0; i < k; ++i) {
500 : long j;
501 35226 : P[1] = uel(S, i + 1);
502 57835 : for (j = 1; j <= h; ++j) {
503 39848 : GEN T = get_nbrs(phi, L, P[j], &P[j - 1], p, pi);
504 39848 : if ( ! vecsmall_len(T))
505 17239 : break;
506 22609 : P[j + 1] = uel(T, 1);
507 : }
508 35226 : if (j < h) pari_err_BUG("next_surface_nbr");
509 35226 : avma = bv;
510 35226 : if (j > h)
511 17987 : break;
512 : }
513 : /* TODO: We could save one get_nbrs call by iterating from i up to
514 : * k-1 and assume that the last (kth) nbr is the one we want. For
515 : * now we're careful and check that this last nbr really is on the
516 : * surface. */
517 17987 : if (i == k) pari_err_BUG("next_surf_nbr");
518 17987 : *nJ = uel(S, i + 1);
519 17987 : avma = av;
520 17987 : return 1;
521 : }
522 :
523 : /* If flag is set, requires a unique common neighor. Otherwise it
524 : * will allow 2 common neighbors and return both. */
525 : INLINE long
526 169451 : common_nbr(
527 : ulong *nbr,
528 : ulong J1, GEN Phi1, long L1,
529 : ulong J2, GEN Phi2, long L2,
530 : ulong p, ulong pi, long flag)
531 : {
532 169451 : pari_sp av = avma;
533 : GEN d, f, g, r;
534 : long rlen;
535 :
536 169451 : g = Flm_Fl_polmodular_evalx(Phi1, L1, J1, p, pi);
537 169452 : f = Flm_Fl_polmodular_evalx(Phi2, L2, J2, p, pi);
538 169451 : d = Flx_gcd(f, g, p);
539 169452 : if (degpol(d) == 1) {
540 168437 : *nbr = Flx_deg1_root(d, p);
541 168438 : avma = av;
542 168438 : return 1;
543 : }
544 1015 : if (flag || degpol(d) != 2) pari_err_BUG("common_neighbour");
545 1015 : r = Flx_roots(d, p);
546 1015 : rlen = vecsmall_len(r);
547 1015 : if ( ! rlen) pari_err_BUG("common_neighbour");
548 1015 : nbr[0] = uel(r, 1);
549 : /* rlen is 1 or 2 depending on whether the root is unique or not. */
550 1015 : nbr[1] = uel(r, rlen);
551 1015 : avma = av;
552 1015 : return 2;
553 : }
554 :
555 : /* Return gcd(Phi1(X,J1)/(X - J0), Phi2(X,J2)). Not stack clean. */
556 : INLINE GEN
557 42204 : common_nbr_pred_poly(
558 : ulong J1, GEN Phi1, long L1,
559 : ulong J2, GEN Phi2, long L2,
560 : ulong J0, ulong p, ulong pi)
561 : {
562 : GEN f, g;
563 :
564 42204 : g = Flm_Fl_polmodular_evalx(Phi1, L1, J1, p, pi);
565 42204 : g = Flx_remove_root(g, J0, p);
566 42204 : f = Flm_Fl_polmodular_evalx(Phi2, L2, J2, p, pi);
567 42204 : return Flx_gcd(f, g, p);
568 : }
569 :
570 : /* Find common neighbour of J1 and J2, where J0 is an L1 predecessor
571 : * of J1. Return 1 if successful, 0 if not. */
572 : INLINE int
573 41181 : common_nbr_pred(
574 : ulong *nbr,
575 : ulong J1, GEN Phi1, long L1,
576 : ulong J2, GEN Phi2, long L2,
577 : ulong J0, ulong p, ulong pi)
578 : {
579 41181 : pari_sp av = avma;
580 : int res;
581 : GEN d;
582 :
583 41181 : d = common_nbr_pred_poly(J1, Phi1, L1, J2, Phi2, L2, J0, p, pi);
584 41181 : res = (degpol(d) == 1);
585 41181 : if (res)
586 41181 : *nbr = Flx_deg1_root(d, p);
587 41181 : avma = av;
588 41181 : return res;
589 : }
590 :
591 : INLINE long
592 1023 : common_nbr_verify(
593 : ulong *nbr,
594 : ulong J1, GEN Phi1, long L1,
595 : ulong J2, GEN Phi2, long L2,
596 : ulong J0, ulong p, ulong pi)
597 : {
598 1023 : pari_sp av = avma;
599 : GEN d;
600 :
601 1023 : d = common_nbr_pred_poly(J1, Phi1, L1, J2, Phi2, L2, J0, p, pi);
602 1023 : if ( ! degpol(d)) {
603 632 : avma = av;
604 632 : return 0;
605 : }
606 391 : if (degpol(d) > 1) pari_err_BUG("common_neighbour_verify");
607 391 : *nbr = Flx_deg1_root(d, p);
608 391 : avma = av;
609 391 : return 1;
610 : }
611 :
612 : INLINE ulong
613 360 : Flm_Fl_polmodular_evalxy(
614 : GEN Phi, long L, ulong x, ulong y, ulong p, ulong pi)
615 : {
616 360 : pari_sp av = avma;
617 360 : GEN f = Flm_Fl_polmodular_evalx(Phi, L, x, p, pi);
618 360 : ulong r = Flx_eval_pre(f, y, p, pi);
619 360 : avma = av;
620 360 : return r;
621 : }
622 :
623 : /* Find a common L1-neighbor of J1 and L2-neighbor of J2, given J0 an
624 : * L2-neighbor of J1 and an L1-neighbor of J2. Return 1 if successful,
625 : * 0 otherwise. Will only fail if the initial J-invariant had the
626 : * wrong endo ring. */
627 : INLINE int
628 21870 : common_nbr_corner(
629 : ulong *nbr,
630 : ulong J1, GEN Phi1, long L1, long h1,
631 : ulong J2, GEN Phi2, long L2,
632 : ulong J0, ulong p, ulong pi)
633 : {
634 : ulong nbrs[2];
635 :
636 21870 : if (common_nbr(nbrs, J1, Phi1, L1, J2, Phi2, L2, p, pi, 0) == 2
637 233 : && nbrs[0] != nbrs[1]) {
638 : ulong nJ1, nJ2;
639 :
640 233 : if ( ! next_surface_nbr(&nJ2, Phi1, L1, h1, J2, &J0, p, pi))
641 13 : return 0; /* Error */
642 233 : if ( ! next_surface_nbr(&nJ1, Phi1, L1, h1, nbrs[0], &J1, p, pi))
643 0 : return 0; /* Error */
644 :
645 233 : if (Flm_Fl_polmodular_evalxy(Phi2, L2, nJ1, nJ2, p, pi))
646 106 : nbrs[0] = nbrs[1];
647 : else {
648 127 : if ( ! next_surface_nbr(&nJ1, Phi1, L1, h1, nbrs[1], &J1, p, pi))
649 0 : return 0; /* Error */
650 127 : if ( ! Flm_Fl_polmodular_evalxy(Phi2, L2, nJ1, nJ2, p, pi))
651 13 : return 0; /* Error */
652 : }
653 : }
654 21857 : *nbr = nbrs[0];
655 21857 : return 1;
656 : }
657 :
658 : /*
659 : * Enumerates a surface L1-cycle using gcds with Phi_L2, where
660 : * alpha_L2=alpha_L1^e and |alpha_L1|=n (Here alpha_a denotes the
661 : * class represented by the pos def reduced primeform <a,b,c>).
662 : * Assumes n > e > 1 and that roots[0],...,roots[e-1] are already
663 : * present in W.
664 : */
665 : static long
666 58404 : surface_gcd_cycle(
667 : ulong W[], ulong V[], long n,
668 : GEN Phi1, long L1, GEN Phi2, long L2, long e, ulong p, ulong pi)
669 : {
670 58404 : pari_sp av = avma;
671 : long i1, i2, j1, j2;
672 :
673 58404 : i1 = j2 = 0;
674 58404 : i2 = j1 = e - 1;
675 : /* If W != V we assume V actually points to an L2-isogenous
676 : * parallel L1-path. e should be 2 in this case */
677 58404 : if (W != V) {
678 58406 : i1 = j1 + 1;
679 58406 : i2 = n - 1;
680 : }
681 : do {
682 : ulong t0, t1, t2, h10, h11, h20, h21;
683 : long k;
684 : GEN f, g, h1, h2;
685 :
686 3409110 : f = Flm_Fl_polmodular_evalx(Phi2, L2, V[i1], p, pi);
687 3400835 : g = Flm_Fl_polmodular_evalx(Phi1, L1, W[j1], p, pi);
688 3404400 : g = Flx_remove_root(g, W[j1 - 1], p);
689 3405320 : h1 = Flx_gcd(f, g, p);
690 3402629 : if (degpol(h1) != 1) break; /* Error */
691 3402593 : h11 = Flx_coeff(h1, 1);
692 3403014 : h10 = Flx_coeff(h1, 0);
693 3404086 : avma = av;
694 :
695 3404086 : f = Flm_Fl_polmodular_evalx(Phi2, L2, V[i2], p, pi);
696 3401595 : g = Flm_Fl_polmodular_evalx(Phi1, L1, W[j2], p, pi);
697 3404451 : k = j2 + 1;
698 3404451 : if (k == n)
699 52824 : k = 0;
700 3404451 : g = Flx_remove_root(g, W[k], p);
701 3405246 : h2 = Flx_gcd(f, g, p);
702 3402674 : if (degpol(h2) != 1) break; /* Error */
703 3402530 : h21 = Flx_coeff(h2, 1);
704 3403349 : h20 = Flx_coeff(h2, 0);
705 3403866 : avma = av;
706 :
707 3403866 : ++i1;
708 3403866 : --i2;
709 3403866 : if (i2 < 0)
710 0 : i2 = n - 1;
711 3403866 : ++j1;
712 3403866 : --j2;
713 3403866 : if (j2 < 0)
714 58406 : j2 = n - 1;
715 :
716 3403866 : t0 = Fl_mul_pre(h11, h21, p, pi);
717 3406678 : t1 = Fl_inv(t0, p);
718 3407959 : t0 = Fl_mul_pre(t1, h21, p, pi);
719 3409223 : t2 = Fl_mul_pre(t0, h10, p, pi);
720 3407821 : W[j1] = Fl_neg(t2, p);
721 3407375 : t0 = Fl_mul_pre(t1, h11, p, pi);
722 3409541 : t2 = Fl_mul_pre(t0, h20, p, pi);
723 3409612 : W[j2] = Fl_neg(t2, p);
724 3409113 : } while (j2 > j1 + 1);
725 58407 : avma = av;
726 : /* Under normal circumstances, the loop exits when j2 == j1 + 1, in
727 : * which case we return n. If we break early because of an error,
728 : * then j2 > j1 + 1, so (j2 - (j1 + 1)) > 0 (which is the number of
729 : * elements we haven't calculated yet), and we return n minus that
730 : * quantity. */
731 58407 : return n - j2 + (j1 + 1);
732 : }
733 :
734 : static long
735 1176 : surface_gcd_path(
736 : ulong W[], ulong V[], long n,
737 : GEN Phi1, long L1, GEN Phi2, long L2, long e, ulong p, ulong pi)
738 : {
739 1176 : pari_sp av = avma;
740 : long i, j;
741 :
742 1176 : i = 0;
743 1176 : j = e;
744 : /* If W != V then assume V actually points to a L2-isogenous
745 : * parallel L1-path. e should be 2 in this case */
746 1176 : if (W != V)
747 1176 : i = j;
748 6048 : while (j < n) {
749 : GEN f, g, d;
750 :
751 3696 : f = Flm_Fl_polmodular_evalx(Phi2, L2, V[i], p, pi);
752 3696 : g = Flm_Fl_polmodular_evalx(Phi1, L1, W[j - 1], p, pi);
753 3696 : g = Flx_remove_root(g, W[j - 2], p);
754 3696 : d = Flx_gcd(f, g, p);
755 3696 : if (degpol(d) != 1) break; /* Error */
756 3696 : W[j] = Flx_deg1_root(d, p);
757 3696 : i++;
758 3696 : j++;
759 3696 : avma = av;
760 : }
761 1176 : avma = av;
762 1176 : return j;
763 : }
764 :
765 : /*
766 : * Given a path V of length n on an L1-volcano, and W[0] L2-isogenous
767 : * to V[0], extends the path W to length n on an L1-volcano, with W[i]
768 : * L2-isogenous to V[i].
769 : *
770 : * Uses gcds unless L2 is too large to make it helpful. Always uses
771 : * GCD to get W[1] to ensure consistent orientation.
772 : *
773 : * Returns the new length of W. This will of course almost always be
774 : * n, but could be lower if V was started with a J-invariant with bad
775 : * endomorphism ring.
776 : */
777 : INLINE long
778 147580 : surface_parallel_path(
779 : ulong W[], ulong V[], long n,
780 : GEN Phi1, long L1, GEN Phi2, long L2, ulong p, ulong pi, long cycle)
781 : {
782 : ulong W2, nbrs[2];
783 :
784 147580 : if (common_nbr(nbrs, W[0], Phi1, L1, V[1], Phi2, L2, p, pi, 0) == 2
785 782 : && nbrs[0] != nbrs[1]) {
786 684 : if (n > 2) {
787 1368 : if ( ! common_nbr_verify(&W2, nbrs[0], Phi1, L1,
788 684 : V[2], Phi2, L2, W[0], p, pi)) {
789 : /* nbrs[1] must be the correct choice */
790 345 : nbrs[0] = nbrs[1];
791 678 : } else if (common_nbr_verify(&W2, nbrs[1], Phi1, L1,
792 339 : V[2], Phi2, L2, W[0], p, pi)) {
793 : /* Error: Both paths extend successfully */
794 52 : return 1;
795 : }
796 : /* Error: Two distinct choices with n = 2; ambiguous. */
797 0 : } else return 1;
798 : }
799 147531 : W[1] = nbrs[0];
800 147531 : if (n > 2) {
801 59582 : return cycle
802 : ? surface_gcd_cycle(W, V, n, Phi1, L1, Phi2, L2, 2, p, pi)
803 59583 : : surface_gcd_path(W, V, n, Phi1, L1, Phi2, L2, 2, p, pi);
804 : }
805 87948 : return n;
806 : }
807 :
808 : GEN
809 182467 : enum_roots(ulong J0, norm_eqn_t ne, GEN fdb, classgp_pcp_t G)
810 : {
811 : /* MAX_HEIGHT just needs to be bigger than the biggest value of
812 : * val_p(n) where p and n are ulongs. */
813 : enum { MAX_HEIGHT = BITS_IN_LONG };
814 182467 : pari_sp av, ltop = avma;
815 182467 : long s = !!G->L0;
816 182467 : long *n = G->n + s, *L = G->L + s, *o = G->o + s, k = G->k - s;
817 182467 : long i, t, vlen, *e, *h, *off, *poff, *M, N = G->enum_cnt;
818 182467 : ulong p = ne->p, pi = ne->pi, *roots;
819 : GEN vshape, vp, ve, roots_;
820 : GEN Phi;
821 :
822 182467 : if ( ! k) {
823 599 : roots_ = cgetg(2, t_VECSMALL);
824 599 : uel(roots_, 1) = J0;
825 599 : return roots_;
826 : }
827 :
828 181868 : roots_ = cgetg(N + MAX_HEIGHT, t_VECSMALL);
829 181866 : roots = zv_to_ulongptr(roots_);
830 181866 : av = avma;
831 :
832 : /* TODO: Shouldn't be factoring this every time. Store in *ne? */
833 181866 : vshape = factoru(ne->v);
834 181859 : vp = gel(vshape, 1);
835 181859 : ve = gel(vshape, 2);
836 :
837 181859 : vlen = vecsmall_len(vp);
838 181859 : Phi = new_chunk(k);
839 181859 : e = new_chunk(k);
840 181863 : off = new_chunk(k);
841 181865 : poff = new_chunk(k);
842 : /* TODO: Surely we can work these out ahead of time? */
843 : /* h[i] is the valuation of p[i] in v */
844 181864 : h = new_chunk(k);
845 412994 : for (i = 0; i < k; ++i) {
846 231129 : h[i] = 0;
847 338975 : for (t = 1; t <= vlen; ++t) {
848 269453 : if (vp[t] == L[i]) {
849 161607 : h[i] = uel(ve, t);
850 161607 : break;
851 : }
852 : }
853 231129 : e[i] = 0;
854 231129 : off[i] = 0;
855 231129 : gel(Phi, i) = polmodular_db_getp(fdb, L[i], p);
856 : }
857 :
858 181865 : M = new_chunk(k);
859 231131 : for (M[0] = 1, i = 1; i < k; ++i)
860 49266 : M[i] = M[i - 1] * n[i - 1];
861 :
862 181865 : t = surface_path(roots, n[0], gel(Phi, 0), L[0], h[0], J0, NULL, p, pi);
863 181870 : if (t < n[0]) {
864 : /* Error: J0 has bad endo ring */
865 355 : avma = ltop;
866 355 : return NULL;
867 : }
868 181515 : if (k == 1) {
869 138720 : avma = av;
870 138720 : setlg(roots_, t + 1);
871 138719 : return roots_;
872 : }
873 :
874 42795 : i = 1;
875 233120 : while (i < k) {
876 : long j, t0;
877 147643 : for (j = i + 1; j < k && ! e[j]; ++j);
878 147643 : if (j < k) {
879 63051 : if (e[i]) {
880 329448 : if (! common_nbr_pred(
881 164724 : &roots[t], roots[off[i]], gel(Phi, i), L[i],
882 164724 : roots[t - M[j]], gel(Phi, j), L[j], roots[poff[i]], p, pi)) {
883 0 : break; /* Error: J0 has bad endo ring */
884 : }
885 196830 : } else if ( ! common_nbr_corner(
886 109350 : &roots[t], roots[off[i]], gel(Phi, i), L[i], h[i],
887 87480 : roots[t - M[j]], gel(Phi, j), L[j], roots[poff[j]], p, pi)) {
888 13 : break; /* Error: J0 has bad endo ring */
889 : }
890 543471 : } else if ( ! next_surface_nbr(
891 338368 : &roots[t], gel(Phi, i), L[i], h[i],
892 205103 : roots[off[i]], e[i] ? &roots[poff[i]] : NULL, p, pi)) {
893 0 : break; /* Error: J0 has bad endo ring */
894 : }
895 :
896 147628 : if (roots[t] == roots[0]) break; /* Error: J0 has bad endo ring */
897 :
898 147581 : poff[i] = off[i];
899 147581 : off[i] = t;
900 147581 : ++e[i];
901 169451 : for (j = i - 1; j; --j) {
902 21870 : e[j] = 0;
903 21870 : off[j] = off[j + 1];
904 : }
905 :
906 442743 : t0 = surface_parallel_path(&roots[t], &roots[poff[i]], n[0],
907 442743 : gel(Phi, 0), L[0], gel(Phi, i), L[i], p, pi, n[0] == o[0]);
908 147582 : if (t0 < n[0]) break; /* Error: J0 has bad endo ring */
909 :
910 : /* TODO: Do I need to check if any of the new roots is a repeat in
911 : * the case where J0 has bad endo ring? */
912 147530 : t += n[0];
913 147530 : for (i = 1; i < k && e[i] == n[i] - 1; ++i);
914 : }
915 : /* Check if J0 had wrong endo ring */
916 42794 : if (t != N) { avma = ltop; return NULL; }
917 :
918 42682 : avma = av;
919 42682 : setlg(roots_, t + 1);
920 42682 : return roots_;
921 : }
|
