Line data Source code
1 : /* Copyright (C) 2014 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : static GEN
19 217558 : pcp_get_L(GEN G) { return gmael(G,1,1)+1; }
20 : static GEN
21 217559 : pcp_get_n(GEN G) { return gmael(G,1,2)+1; }
22 : static GEN
23 217558 : pcp_get_o(GEN G) { return gmael(G,1,3)+1; }
24 : static long
25 217561 : pcp_get_L0(GEN G) { return mael(G,2,1); }
26 : static long
27 217558 : pcp_get_k(GEN G) { return mael(G,2,2); }
28 : static long
29 217558 : pcp_get_enum_cnt(GEN G) { return mael(G,2,3); }
30 :
31 : /* FIXME: Implement {ascend,descend}_volcano() in terms of the "new"
32 : * volcano traversal functions at the bottom of the file. */
33 :
34 : /* Is j = 0 or 1728 (mod p)? */
35 : INLINE int
36 359754 : is_j_exceptional(ulong j, ulong p) { return j == 0 || j == 1728 % p; }
37 :
38 : INLINE long
39 80785 : node_degree(GEN phi, long L, ulong j, ulong p, ulong pi)
40 : {
41 80785 : pari_sp av = avma;
42 80785 : long n = Flx_nbroots(Flm_Fl_polmodular_evalx(phi, L, j, p, pi), p);
43 80785 : return gc_long(av, n);
44 : }
45 :
46 : /* Given an array path = [j0, j1] of length 2, return the polynomial
47 : *
48 : * \Phi_L(X, j1) / (X - j0)
49 : *
50 : * where \Phi_L(X, Y) is the modular polynomial of level L. An error
51 : * is raised if X - j0 does not divide \Phi_L(X, j1) */
52 : INLINE GEN
53 141785 : nhbr_polynomial(ulong path[], GEN phi, ulong p, ulong pi, long L)
54 : {
55 141785 : GEN modpol = Flm_Fl_polmodular_evalx(phi, L, path[0], p, pi);
56 : ulong rem;
57 141790 : GEN nhbr_pol = Flx_div_by_X_x(modpol, path[-1], p, &rem);
58 : /* If disc End(path[0]) <= L^2, it's possible for path[0] to appear among the
59 : * roots of nhbr_pol. This should have been obviated by earlier choices */
60 141782 : if (rem) pari_err_BUG("nhbr_polynomial: invalid preceding j");
61 141783 : return nhbr_pol;
62 : }
63 :
64 : /* Path is an array with space for at least max_len+1 * elements, whose first
65 : * and second elements are the beginning of the path. I.e., the path starts
66 : * (path[0], path[1])
67 : * If the result is less than max_len, then the last element of path is on the
68 : * floor. If the result equals max_len, then it is unknown whether the last
69 : * element of path is on the floor or not */
70 : static long
71 277188 : extend_path(ulong path[], GEN phi, ulong p, ulong pi, long L, long max_len)
72 : {
73 277188 : pari_sp av = avma;
74 277188 : long d = 1;
75 357067 : for ( ; d < max_len; d++)
76 : {
77 102904 : GEN nhbr_pol = nhbr_polynomial(path + d, phi, p, pi, L);
78 102901 : ulong nhbr = Flx_oneroot_pre(nhbr_pol, p, pi);
79 102907 : set_avma(av);
80 102906 : if (nhbr == p) break; /* no root: we are on the floor. */
81 79879 : path[d+1] = nhbr;
82 : }
83 277190 : return d;
84 : }
85 :
86 : /* This is Sutherland 2009 Algorithm Ascend (p12) */
87 : ulong
88 126315 : ascend_volcano(GEN phi, ulong j, ulong p, ulong pi, long level, long L,
89 : long depth, long steps)
90 : {
91 126315 : pari_sp ltop = avma, av;
92 : /* path will never hold more than max_len < depth elements */
93 126315 : GEN path_g = cgetg(depth + 2, t_VECSMALL);
94 126315 : ulong *path = zv_to_ulongptr(path_g);
95 126315 : long max_len = depth - level;
96 126315 : int first_iter = 1;
97 :
98 126315 : if (steps <= 0 || max_len < 0) pari_err_BUG("ascend_volcano: bad params");
99 126315 : av = avma;
100 291520 : while (steps--)
101 : {
102 126315 : GEN nhbr_pol = first_iter? Flm_Fl_polmodular_evalx(phi, L, j, p, pi)
103 165201 : : nhbr_polynomial(path+1, phi, p, pi, L);
104 165204 : GEN nhbrs = Flx_roots_pre(nhbr_pol, p, pi);
105 165200 : long nhbrs_len = lg(nhbrs)-1, i;
106 165200 : pari_sp btop = avma;
107 165200 : path[0] = j;
108 165200 : first_iter = 0;
109 :
110 165200 : j = nhbrs[nhbrs_len];
111 208181 : for (i = 1; i < nhbrs_len; i++)
112 : {
113 77622 : ulong next_j = nhbrs[i], last_j;
114 : long len;
115 77622 : if (is_j_exceptional(next_j, p))
116 : {
117 : /* Fouquet & Morain, Section 4.3, if j = 0 or 1728, then it is on the
118 : * surface. So we just return it. */
119 40 : if (steps)
120 0 : pari_err_BUG("ascend_volcano: Got to the top with more steps to go!");
121 40 : j = next_j; break;
122 : }
123 77582 : path[1] = next_j;
124 77582 : len = extend_path(path, phi, p, pi, L, max_len);
125 77584 : last_j = path[len];
126 77584 : if (len == max_len
127 : /* Ended up on the surface */
128 77584 : && (is_j_exceptional(last_j, p)
129 77584 : || node_degree(phi, L, last_j, p, pi) > 1)) { j = next_j; break; }
130 42979 : set_avma(btop);
131 : }
132 165203 : path[1] = j; /* For nhbr_polynomial() at the top. */
133 :
134 165203 : max_len++; set_avma(av);
135 : }
136 126319 : return gc_long(ltop, j);
137 : }
138 :
139 : static void
140 204549 : random_distinct_neighbours_of(ulong *nhbr1, ulong *nhbr2,
141 : GEN phi, ulong j, ulong p, ulong pi, long L, long must_have_two_neighbours)
142 : {
143 204549 : pari_sp av = avma;
144 204549 : GEN modpol = Flm_Fl_polmodular_evalx(phi, L, j, p, pi);
145 : ulong rem;
146 204549 : *nhbr1 = Flx_oneroot_pre(modpol, p, pi);
147 204551 : if (*nhbr1 == p) pari_err_BUG("random_distinct_neighbours_of [no neighbour]");
148 204551 : modpol = Flx_div_by_X_x(modpol, *nhbr1, p, &rem);
149 204542 : *nhbr2 = Flx_oneroot_pre(modpol, p, pi);
150 204551 : if (must_have_two_neighbours && *nhbr2 == p)
151 0 : pari_err_BUG("random_distinct_neighbours_of [single neighbour]");
152 204551 : set_avma(av);
153 204551 : }
154 :
155 : /*
156 : * This is Sutherland 2009 Algorithm Descend (p12).
157 : */
158 : ulong
159 2937 : descend_volcano(GEN phi, ulong j, ulong p, ulong pi,
160 : long level, long L, long depth, long steps)
161 : {
162 2937 : pari_sp ltop = avma;
163 : GEN path_g;
164 : ulong *path;
165 : long max_len;
166 :
167 2937 : if (steps <= 0 || level + steps > depth) pari_err_BUG("descend_volcano");
168 2937 : max_len = depth - level;
169 2937 : path_g = cgetg(max_len + 1 + 1, t_VECSMALL);
170 2937 : path = zv_to_ulongptr(path_g);
171 2937 : path[0] = j;
172 : /* level = 0 means we're on the volcano surface... */
173 2937 : if (!level)
174 : {
175 : /* Look for any path to the floor. One of j's first three neighbours leads
176 : * to the floor, since at most two neighbours are on the surface. */
177 2652 : GEN nhbrs = Flx_roots_pre(Flm_Fl_polmodular_evalx(phi, L, j, p, pi), p, pi);
178 : long i;
179 2957 : for (i = 1; i <= 3; i++)
180 : {
181 : long len;
182 2957 : path[1] = nhbrs[i];
183 2957 : len = extend_path(path, phi, p, pi, L, max_len);
184 : /* If nhbrs[i] took us to the floor: */
185 2956 : if (len < max_len || node_degree(phi, L, path[len], p, pi) == 1) break;
186 : }
187 2652 : if (i > 3) pari_err_BUG("descend_volcano [2]");
188 : }
189 : else
190 : {
191 : ulong nhbr1, nhbr2;
192 : long len;
193 285 : random_distinct_neighbours_of(&nhbr1, &nhbr2, phi, j, p, pi, L, 1);
194 285 : path[1] = nhbr1;
195 285 : len = extend_path(path, phi, p, pi, L, max_len);
196 : /* If last j isn't on the floor */
197 285 : if (len == max_len /* Ended up on the surface. */
198 285 : && (is_j_exceptional(path[len], p)
199 245 : || node_degree(phi, L, path[len], p, pi) != 1)) {
200 : /* The other neighbour leads to the floor */
201 119 : path[1] = nhbr2;
202 119 : (void) extend_path(path, phi, p, pi, L, steps);
203 : }
204 : }
205 2937 : return gc_ulong(ltop, path[steps]);
206 : }
207 :
208 : long
209 204272 : j_level_in_volcano(
210 : GEN phi, ulong j, ulong p, ulong pi, long L, long depth)
211 : {
212 204272 : pari_sp av = avma;
213 : GEN chunk;
214 : ulong *path1, *path2;
215 : long lvl;
216 :
217 : /* Fouquet & Morain, Section 4.3, if j = 0 or 1728 then it is on the
218 : * surface. Also, if the volcano depth is zero then j has level 0 */
219 204272 : if (depth == 0 || is_j_exceptional(j, p)) return 0;
220 :
221 204266 : chunk = new_chunk(2 * (depth + 1));
222 204265 : path1 = (ulong *) &chunk[0];
223 204265 : path2 = (ulong *) &chunk[depth + 1];
224 204265 : path1[0] = path2[0] = j;
225 204265 : random_distinct_neighbours_of(&path1[1], &path2[1], phi, j, p, pi, L, 0);
226 204267 : if (path2[1] == p)
227 106139 : lvl = depth; /* Only one neighbour => j is on the floor => level = depth */
228 : else
229 : {
230 98128 : long path1_len = extend_path(path1, phi, p, pi, L, depth);
231 98127 : long path2_len = extend_path(path2, phi, p, pi, L, path1_len);
232 98129 : lvl = depth - path2_len;
233 : }
234 204268 : return gc_long(av, lvl);
235 : }
236 :
237 : INLINE GEN
238 32226777 : Flx_remove_root(GEN f, ulong a, ulong p)
239 : {
240 : ulong r;
241 32226777 : GEN g = Flx_div_by_X_x(f, a, p, &r);
242 32076020 : if (r) pari_err_BUG("Flx_remove_root");
243 32076347 : return g;
244 : }
245 :
246 : INLINE GEN
247 24299080 : get_nbrs(GEN phi, long L, ulong J, const ulong *xJ, ulong p, ulong pi)
248 : {
249 24299080 : pari_sp av = avma;
250 24299080 : GEN f = Flm_Fl_polmodular_evalx(phi, L, J, p, pi);
251 24300386 : if (xJ) f = Flx_remove_root(f, *xJ, p);
252 24208758 : return gerepileupto(av, Flx_roots_pre(f, p, pi));
253 : }
254 :
255 : /* Return a path of length n along the surface of an L-volcano of height h
256 : * starting from surface node j0. Assumes (D|L) = 1 where D = disc End(j0).
257 : *
258 : * Actually, if j0's endomorphism ring is a suborder, we return the
259 : * corresponding shorter path. W must hold space for n + h nodes.
260 : *
261 : * TODO: have two versions of this function: one that assumes J has the correct
262 : * endomorphism ring (hence avoiding several branches in the inner loop) and a
263 : * second that does not and accordingly checks for repetitions */
264 : static long
265 215219 : surface_path(
266 : ulong W[],
267 : long n, GEN phi, long L, long h, ulong J, const ulong *nJ, ulong p, ulong pi)
268 : {
269 215219 : pari_sp av = avma, bv;
270 : GEN T, v;
271 : long j, k, w, x;
272 : ulong W0;
273 :
274 215219 : W[0] = W0 = J;
275 215219 : if (n == 1) return 1;
276 :
277 215219 : T = cgetg(h+2, t_VEC); bv = avma;
278 215219 : v = get_nbrs(phi, L, J, nJ, p, pi);
279 : /* Insert known neighbour first */
280 215215 : if (nJ) v = gerepileupto(bv, vecsmall_prepend(v, *nJ));
281 215215 : gel(T,1) = v; k = lg(v)-1;
282 :
283 215215 : switch (k) {
284 0 : case 0: pari_err_BUG("surface_path"); /* We must always have neighbours */
285 8764 : case 1:
286 : /* If volcano is not flat, then we must have more than one neighbour */
287 8764 : if (h) pari_err_BUG("surface_path");
288 8764 : W[1] = uel(v, 1);
289 8764 : set_avma(av);
290 : /* Check for bad endo ring */
291 8764 : if (W[1] == W[0]) return 1;
292 8591 : return 2;
293 25982 : case 2:
294 : /* If L=2 the only way we can have 2 neighbours is if we have a double root
295 : * which can only happen for |D| <= 16 (Thm 2.2 of Fouquet-Morain)
296 : * and if it does we must have a 2-cycle. Happens for D=-15. */
297 25982 : if (L == 2)
298 : { /* The double root is the neighbour on the surface, with exactly one
299 : * neighbour other than J; the other neighbour of J has either 0 or 2
300 : * neighbours that are not J */
301 84 : GEN u = get_nbrs(phi, L, uel(v, 1), &J, p, pi);
302 84 : long n = lg(u)-1 - !!vecsmall_isin(u, J);
303 84 : W[1] = n == 1 ? uel(v,1) : uel(v,2);
304 84 : return gc_long(av, 2);
305 : }
306 : /* Volcano is not flat but found only 2 neighbours for the surface node J */
307 25898 : if (h) pari_err_BUG("surface_path");
308 :
309 25898 : W[1] = uel(v,1); /* TODO: Can we use the other root uel(v,2) somehow? */
310 4451045 : for (w = 2; w < n; w++)
311 : {
312 4425688 : v = get_nbrs(phi, L, W[w-1], &W[w-2], p, pi);
313 : /* A flat volcano must have exactly one non-previous neighbour */
314 4425651 : if (lg(v) != 2) pari_err_BUG("surface_path");
315 4425651 : W[w] = uel(v, 1);
316 : /* Detect cycle in case J doesn't have the right endo ring. */
317 4425651 : set_avma(av); if (W[w] == W0) return w;
318 : }
319 25357 : return gc_long(av, n);
320 : }
321 180469 : if (!h) pari_err_BUG("surface_path"); /* Can't have a flat volcano if k > 2 */
322 :
323 : /* At this point, each surface node has L+1 distinct neighbours, 2 of which
324 : * are on the surface */
325 180469 : w = 1;
326 6406224 : for (x = 0;; x++)
327 : {
328 : /* Get next neighbour of last known surface node to attempt to
329 : * extend the path. */
330 6406224 : W[w] = umael(T, ((w-1) % h) + 1, x + 1);
331 :
332 : /* Detect cycle in case the given J didn't have the right endo ring */
333 6406224 : if (W[w] == W0) return gc_long(av,w);
334 :
335 : /* If we have to test the last neighbour, we know it's on the
336 : * surface, and if we're done there's no need to extend. */
337 6406192 : if (x == k-1 && w == n-1) return gc_long(av,n);
338 :
339 : /* Walk forward until we hit the floor or finish. */
340 : /* NB: We don't keep the stack clean here; usage is in the order of Lh,
341 : * i.e. L roots for each level of the volcano of height h. */
342 6287614 : for (j = w;;)
343 13186400 : {
344 : long m;
345 : /* We must get 0 or L neighbours here. */
346 19474014 : v = get_nbrs(phi, L, W[j], &W[j-1], p, pi);
347 19419218 : m = lg(v)-1;
348 19419218 : if (!m) {
349 : /* We hit the floor: save the neighbours of W[w-1] and dump the rest */
350 6223528 : GEN nbrs = gel(T, ((w-1) % h) + 1);
351 6223528 : gel(T, ((w-1) % h) + 1) = gerepileupto(bv, nbrs);
352 6225755 : break;
353 : }
354 13195690 : if (m != L) pari_err_BUG("surface_path");
355 :
356 13248261 : gel(T, (j % h) + 1) = v;
357 :
358 13248261 : W[++j] = uel(v, 1);
359 : /* If we have our path by h nodes, we know W[w] is on the surface */
360 13248261 : if (j == w + h) {
361 12256507 : ++w;
362 : /* Detect cycle in case the given J didn't have the right endo ring */
363 12256507 : if (W[w] == W0) return gc_long(av,w);
364 12228756 : x = 0; k = L;
365 : }
366 13220510 : if (w == n) return gc_long(av,w);
367 : }
368 : }
369 : }
370 :
371 : long
372 146015 : next_surface_nbr(
373 : ulong *nJ,
374 : GEN phi, long L, long h, ulong J, const ulong *pJ, ulong p, ulong pi)
375 : {
376 146015 : pari_sp av = avma, bv;
377 : GEN S;
378 : ulong *P;
379 : long i, k;
380 :
381 146015 : S = get_nbrs(phi, L, J, pJ, p, pi); k = lg(S)-1;
382 : /* If there is a double root and pJ is set, then k will be zero. */
383 146002 : if (!k) return gc_long(av,0);
384 146002 : if (k == 1 || ( ! pJ && k == 2)) { *nJ = uel(S, 1); return gc_long(av,1); }
385 24090 : if (!h) pari_err_BUG("next_surface_nbr");
386 :
387 24090 : P = (ulong *) new_chunk(h + 1); bv = avma;
388 24090 : P[0] = J;
389 51390 : for (i = 0; i < k; i++)
390 : {
391 : long j;
392 51391 : P[1] = uel(S, i + 1);
393 81005 : for (j = 1; j <= h; j++)
394 : {
395 56915 : GEN T = get_nbrs(phi, L, P[j], &P[j - 1], p, pi);
396 56915 : if (lg(T) == 1) break;
397 29614 : P[j + 1] = uel(T, 1);
398 : }
399 51391 : if (j < h) pari_err_BUG("next_surface_nbr");
400 51391 : set_avma(bv); if (j > h) break;
401 : }
402 : /* TODO: We could save one get_nbrs call by iterating from i up to k-1 and
403 : * assume that the last (kth) nbr is the one we want. For now we're careful
404 : * and check that this last nbr really is on the surface */
405 24090 : if (i == k) pari_err_BUG("next_surf_nbr");
406 24090 : *nJ = uel(S, i+1); return gc_long(av,1);
407 : }
408 :
409 : /* Return the number of distinct neighbours (1 or 2) */
410 : INLINE long
411 238053 : common_nbr(ulong *nbr,
412 : ulong J1, GEN Phi1, long L1,
413 : ulong J2, GEN Phi2, long L2, ulong p, ulong pi)
414 : {
415 238053 : pari_sp av = avma;
416 : GEN d, f, g, r;
417 : long rlen;
418 :
419 238053 : g = Flm_Fl_polmodular_evalx(Phi1, L1, J1, p, pi);
420 238092 : f = Flm_Fl_polmodular_evalx(Phi2, L2, J2, p, pi);
421 238083 : d = Flx_gcd(f, g, p);
422 237935 : if (degpol(d) == 1) { *nbr = Flx_deg1_root(d, p); return gc_long(av,1); }
423 1237 : if (degpol(d) != 2) pari_err_BUG("common_neighbour");
424 1237 : r = Flx_roots_pre(d, p, pi);
425 1237 : rlen = lg(r)-1;
426 1237 : if (!rlen) pari_err_BUG("common_neighbour");
427 : /* rlen is 1 or 2 depending on whether the root is unique or not. */
428 1237 : nbr[0] = uel(r, 1);
429 1237 : nbr[1] = uel(r, rlen); return gc_long(av,rlen);
430 : }
431 :
432 : /* Return gcd(Phi1(X,J1)/(X - J0), Phi2(X,J2)). Not stack clean. */
433 : INLINE GEN
434 44097 : common_nbr_pred_poly(
435 : ulong J1, GEN Phi1, long L1,
436 : ulong J2, GEN Phi2, long L2, ulong J0, ulong p, ulong pi)
437 : {
438 : GEN f, g;
439 :
440 44097 : g = Flm_Fl_polmodular_evalx(Phi1, L1, J1, p, pi);
441 44101 : g = Flx_remove_root(g, J0, p);
442 44103 : f = Flm_Fl_polmodular_evalx(Phi2, L2, J2, p, pi);
443 44103 : return Flx_gcd(f, g, p);
444 : }
445 :
446 : /* Find common neighbour of J1 and J2, where J0 is an L1 predecessor of J1.
447 : * Return 1 if successful, 0 if not. */
448 : INLINE int
449 43024 : common_nbr_pred(ulong *nbr,
450 : ulong J1, GEN Phi1, long L1,
451 : ulong J2, GEN Phi2, long L2, ulong J0, ulong p, ulong pi)
452 : {
453 43024 : pari_sp av = avma;
454 43024 : GEN d = common_nbr_pred_poly(J1, Phi1, L1, J2, Phi2, L2, J0, p, pi);
455 43013 : int res = (degpol(d) == 1);
456 43013 : if (res) *nbr = Flx_deg1_root(d, p);
457 43019 : return gc_bool(av,res);
458 : }
459 :
460 : INLINE long
461 1073 : common_nbr_verify(ulong *nbr,
462 : ulong J1, GEN Phi1, long L1,
463 : ulong J2, GEN Phi2, long L2, ulong J0, ulong p, ulong pi)
464 : {
465 1073 : pari_sp av = avma;
466 1073 : GEN d = common_nbr_pred_poly(J1, Phi1, L1, J2, Phi2, L2, J0, p, pi);
467 :
468 1073 : if (!degpol(d)) return gc_long(av,0);
469 405 : if (degpol(d) > 1) pari_err_BUG("common_neighbour_verify");
470 405 : *nbr = Flx_deg1_root(d, p);
471 405 : return gc_long(av,1);
472 : }
473 :
474 : INLINE ulong
475 477 : Flm_Fl_polmodular_evalxy(GEN Phi, long L, ulong x, ulong y, ulong p, ulong pi)
476 : {
477 477 : pari_sp av = avma;
478 477 : GEN f = Flm_Fl_polmodular_evalx(Phi, L, x, p, pi);
479 477 : return gc_ulong(av, Flx_eval_pre(f, y, p, pi));
480 : }
481 :
482 : /* Find a common L1-neighbor of J1 and L2-neighbor of J2, given J0 an
483 : * L2-neighbor of J1 and an L1-neighbor of J2. Return 1 if successful, 0
484 : * otherwise. Will only fail if initial J-invariant had the wrong endo ring */
485 : INLINE int
486 37265 : common_nbr_corner(ulong *nbr,
487 : ulong J1, GEN Phi1, long L1, long h1,
488 : ulong J2, GEN Phi2, long L2, ulong J0, ulong p, ulong pi)
489 : {
490 : ulong nbrs[2];
491 37265 : if (common_nbr(nbrs, J1,Phi1,L1, J2,Phi2,L2, p, pi) == 2)
492 : {
493 : ulong nJ1, nJ2;
494 641 : if (!next_surface_nbr(&nJ2, Phi1, L1, h1, J2, &J0, p, pi) ||
495 355 : !next_surface_nbr(&nJ1, Phi1, L1, h1, nbrs[0], &J1, p, pi)) return 0;
496 :
497 320 : if (Flm_Fl_polmodular_evalxy(Phi2, L2, nJ1, nJ2, p, pi))
498 163 : nbrs[0] = nbrs[1];
499 314 : else if (!next_surface_nbr(&nJ1, Phi1, L1, h1, nbrs[1], &J1, p, pi) ||
500 192 : !Flm_Fl_polmodular_evalxy(Phi2, L2, nJ1, nJ2, p, pi)) return 0;
501 : }
502 37228 : *nbr = nbrs[0]; return 1;
503 : }
504 :
505 : /* Enumerate a surface L1-cycle using gcds with Phi_L2, where c_L2=c_L1^e and
506 : * |c_L1|=n, where c_a is the class of the pos def reduced primeform <a,b,c>.
507 : * Assumes n > e > 1 and roots[0],...,roots[e-1] are already present in W */
508 : static long
509 93835 : surface_gcd_cycle(
510 : ulong W[], ulong V[], long n,
511 : GEN Phi1, long L1, GEN Phi2, long L2, long e, ulong p, ulong pi)
512 : {
513 93835 : pari_sp av = avma;
514 : long i1, i2, j1, j2;
515 :
516 93835 : i1 = j2 = 0;
517 93835 : i2 = j1 = e - 1;
518 : /* If W != V we assume V actually points to an L2-isogenous parallel L1-path.
519 : * e should be 2 in this case */
520 93835 : if (W != V) { i1 = j1+1; i2 = n-1; }
521 : do {
522 : ulong t0, t1, t2, h10, h11, h20, h21;
523 : long k;
524 : GEN f, g, h1, h2;
525 :
526 4127092 : f = Flm_Fl_polmodular_evalx(Phi2, L2, V[i1], p, pi);
527 4123077 : g = Flm_Fl_polmodular_evalx(Phi1, L1, W[j1], p, pi);
528 4122676 : g = Flx_remove_root(g, W[j1 - 1], p);
529 4111347 : h1 = Flx_gcd(f, g, p);
530 4101737 : if (degpol(h1) != 1) break; /* Error */
531 4102260 : h11 = Flx_coeff(h1, 1);
532 4103356 : h10 = Flx_coeff(h1, 0); set_avma(av);
533 :
534 4104041 : f = Flm_Fl_polmodular_evalx(Phi2, L2, V[i2], p, pi);
535 4123646 : g = Flm_Fl_polmodular_evalx(Phi1, L1, W[j2], p, pi);
536 4123337 : k = j2 + 1;
537 4123337 : if (k == n) k = 0;
538 4123337 : g = Flx_remove_root(g, W[k], p);
539 4112758 : h2 = Flx_gcd(f, g, p);
540 4103499 : if (degpol(h2) != 1) break; /* Error */
541 4104099 : h21 = Flx_coeff(h2, 1);
542 4104786 : h20 = Flx_coeff(h2, 0); set_avma(av);
543 :
544 4106060 : i1++; i2--; if (i2 < 0) i2 = n-1;
545 4106060 : j1++; j2--; if (j2 < 0) j2 = n-1;
546 :
547 4106060 : t0 = Fl_mul_pre(h11, h21, p, pi);
548 4126623 : t1 = Fl_inv(t0, p);
549 4125115 : t0 = Fl_mul_pre(t1, h21, p, pi);
550 4125419 : t2 = Fl_mul_pre(t0, h10, p, pi);
551 4126120 : W[j1] = Fl_neg(t2, p);
552 4125611 : t0 = Fl_mul_pre(t1, h11, p, pi);
553 4127608 : t2 = Fl_mul_pre(t0, h20, p, pi);
554 4127665 : W[j2] = Fl_neg(t2, p);
555 4126992 : } while (j2 > j1 + 1);
556 : /* Usually the loop exits when j2 = j1 + 1, in which case we return n.
557 : * If we break early because of an error, then (j2 - (j1+1)) > 0 is the
558 : * number of elements we haven't calculated yet, and we return n minus that
559 : * quantity */
560 93735 : return gc_long(av, n - j2 + j1 + 1);
561 : }
562 :
563 : static long
564 1212 : surface_gcd_path(
565 : ulong W[], ulong V[], long n,
566 : GEN Phi1, long L1, GEN Phi2, long L2, long e, ulong p, ulong pi)
567 : {
568 1212 : pari_sp av = avma;
569 : long i, j;
570 :
571 1212 : i = 0; j = e;
572 : /* If W != V then assume V actually points to a L2-isogenous
573 : * parallel L1-path. e should be 2 in this case */
574 1212 : if (W != V) i = j;
575 4943 : while (j < n)
576 : {
577 : GEN f, g, d;
578 :
579 3731 : f = Flm_Fl_polmodular_evalx(Phi2, L2, V[i], p, pi);
580 3733 : g = Flm_Fl_polmodular_evalx(Phi1, L1, W[j - 1], p, pi);
581 3732 : g = Flx_remove_root(g, W[j - 2], p);
582 3729 : d = Flx_gcd(f, g, p);
583 3729 : if (degpol(d) != 1) break; /* Error */
584 3729 : W[j] = Flx_deg1_root(d, p);
585 3730 : i++; j++; set_avma(av);
586 : }
587 1212 : return gc_long(av, j);
588 : }
589 :
590 : /* Given a path V of length n on an L1-volcano, and W[0] L2-isogenous to V[0],
591 : * extends the path W to length n on an L1-volcano, with W[i] L2-isogenous
592 : * to V[i]. Uses gcds unless L2 is too large to make it helpful. Always uses
593 : * GCD to get W[1] to ensure consistent orientation.
594 : *
595 : * Returns the new length of W. This will almost always be n, but could be
596 : * lower if V was started with a J-invariant with bad endomorphism ring */
597 : INLINE long
598 200806 : surface_parallel_path(
599 : ulong W[], ulong V[], long n,
600 : GEN Phi1, long L1, GEN Phi2, long L2, ulong p, ulong pi, long cycle)
601 : {
602 : ulong W2, nbrs[2];
603 200806 : if (common_nbr(nbrs, W[0], Phi1, L1, V[1], Phi2, L2, p, pi) == 2)
604 : {
605 721 : if (n <= 2) return 1; /* Error: Two choices with n = 2; ambiguous */
606 721 : if (!common_nbr_verify(&W2,nbrs[0], Phi1,L1,V[2], Phi2,L2,W[0], p,pi))
607 369 : nbrs[0] = nbrs[1]; /* nbrs[1] must be the correct choice */
608 352 : else if (common_nbr_verify(&W2,nbrs[1], Phi1,L1,V[2], Phi2,L2,W[0], p,pi))
609 53 : return 1; /* Error: Both paths extend successfully */
610 : }
611 200741 : W[1] = nbrs[0];
612 200741 : if (n <= 2) return n;
613 93835 : return cycle? surface_gcd_cycle(W, V, n, Phi1, L1, Phi2, L2, 2, p, pi)
614 188892 : : surface_gcd_path (W, V, n, Phi1, L1, Phi2, L2, 2, p, pi);
615 : }
616 :
617 : GEN
618 217563 : enum_roots(ulong J0, norm_eqn_t ne, GEN fdb, GEN G, GEN vshape)
619 : { /* MAX_HEIGHT >= max_{p,n} val_p(n) where p and n are ulongs */
620 : enum { MAX_HEIGHT = BITS_IN_LONG };
621 217563 : pari_sp av, ltop = avma;
622 217563 : long s = !!pcp_get_L0(G);
623 217559 : long *n = pcp_get_n(G)+s, *L = pcp_get_L(G)+s, *o = pcp_get_o(G)+s, k = pcp_get_k(G)-s;
624 217558 : long i, t, vlen, *e, *h, *off, *poff, *M, N = pcp_get_enum_cnt(G);
625 217558 : ulong p = ne->p, pi = ne->pi, *roots;
626 : GEN Phi, vp, ve, roots_;
627 :
628 217558 : if (!k) return mkvecsmall(J0);
629 :
630 215209 : roots_ = cgetg(N + MAX_HEIGHT, t_VECSMALL);
631 215212 : roots = zv_to_ulongptr(roots_);
632 215212 : av = avma;
633 :
634 215212 : if (!vshape) vshape = factoru(ne->v);
635 215214 : vp = gel(vshape, 1); vlen = lg(vp)-1;
636 215214 : ve = gel(vshape, 2);
637 :
638 215214 : Phi = new_chunk(k);
639 215212 : e = new_chunk(k);
640 215211 : off = new_chunk(k);
641 215211 : poff = new_chunk(k);
642 : /* TODO: Surely we can work these out ahead of time? */
643 : /* h[i] is the valuation of p[i] in v */
644 215211 : h = new_chunk(k);
645 497976 : for (i = 0; i < k; ++i) {
646 282757 : h[i] = 0;
647 420146 : for (t = 1; t <= vlen; ++t)
648 326861 : if (vp[t] == L[i]) { h[i] = uel(ve, t); break; }
649 282757 : e[i] = 0;
650 282757 : off[i] = 0;
651 282757 : gel(Phi, i) = polmodular_db_getp(fdb, L[i], p);
652 : }
653 :
654 215219 : t = surface_path(roots, n[0], gel(Phi, 0), L[0], h[0], J0, NULL, p, pi);
655 215213 : if (t < n[0]) return gc_NULL(ltop); /* J0 has bad endo ring */
656 214507 : if (k == 1) { setlg(roots_, t + 1); return gc_const(av,roots_); }
657 :
658 54952 : M = new_chunk(k);
659 121319 : for (M[0] = 1, i = 1; i < k; ++i) M[i] = M[i-1] * n[i-1];
660 54952 : i = 1;
661 255701 : while (i < k) {
662 : long j, t0;
663 227577 : for (j = i + 1; j < k && ! e[j]; ++j);
664 200893 : if (j < k) {
665 80283 : if (e[i]) {
666 43020 : if (! common_nbr_pred(
667 43023 : &roots[t], roots[off[i]], gel(Phi,i), L[i],
668 43023 : roots[t - M[j]], gel(Phi, j), L[j], roots[poff[i]], p, pi)) {
669 0 : break; /* J0 has bad endo ring */
670 : }
671 37263 : } else if ( ! common_nbr_corner(
672 37260 : &roots[t], roots[off[i]], gel(Phi,i), L[i], h[i],
673 37260 : roots[t - M[j]], gel(Phi, j), L[j], roots[poff[j]], p, pi)) {
674 35 : break; /* J0 has bad endo ring */
675 : }
676 186861 : } else if ( ! next_surface_nbr(
677 120610 : &roots[t], gel(Phi,i), L[i], h[i],
678 186848 : roots[off[i]], e[i] ? &roots[poff[i]] : NULL, p, pi))
679 0 : break; /* J0 has bad endo ring */
680 200871 : if (roots[t] == roots[0]) break; /* J0 has bad endo ring */
681 :
682 200814 : poff[i] = off[i];
683 200814 : off[i] = t;
684 200814 : e[i]++;
685 238087 : for (j = i-1; j; --j) { e[j] = 0; off[j] = off[j+1]; }
686 :
687 200814 : t0 = surface_parallel_path(&roots[t], &roots[poff[i]], n[0],
688 200814 : gel(Phi, 0), L[0], gel(Phi, i), L[i], p, pi, n[0] == o[0]);
689 200802 : if (t0 < n[0]) break; /* J0 has bad endo ring */
690 :
691 : /* TODO: Do I need to check if any of the new roots is a repeat in
692 : * the case where J0 has bad endo ring? */
693 200749 : t += n[0];
694 304081 : for (i = 1; i < k && e[i] == n[i]-1; i++);
695 : }
696 54953 : if (t != N) return gc_NULL(ltop); /* J0 has wrong endo ring */
697 54808 : setlg(roots_, t + 1); return gc_const(av,roots_);
698 : }
|
