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 316709 : is_j_exceptional(ulong j, ulong p) { return j == 0 || j == 1728 % p; }
23 :
24 : INLINE long
25 70728 : node_degree(GEN phi, long L, ulong j, ulong p, ulong pi)
26 : {
27 70728 : pari_sp av = avma;
28 70728 : long n = Flx_nbroots(Flm_Fl_polmodular_evalx(phi, L, j, p, pi), p);
29 70727 : return gc_long(av, n);
30 : }
31 :
32 : /* Given an array path = [j0, j1] of length 2, return the polynomial
33 : *
34 : * \Phi_L(X, j1) / (X - j0)
35 : *
36 : * where \Phi_L(X, Y) is the modular polynomial of level L. An error
37 : * is raised if X - j0 does not divide \Phi_L(X, j1) */
38 : INLINE GEN
39 124586 : nhbr_polynomial(ulong path[], GEN phi, ulong p, ulong pi, long L)
40 : {
41 124586 : pari_sp ltop = avma;
42 124586 : GEN modpol = Flm_Fl_polmodular_evalx(phi, L, path[0], p, pi);
43 : ulong rem;
44 124586 : GEN nhbr_pol = Flx_div_by_X_x(modpol, path[-1], p, &rem);
45 : /* If disc End(path[0]) <= L^2, it's possible for path[0] to appear among the
46 : * roots of nhbr_pol. This should have been obviated by earlier choices */
47 124586 : if (rem) pari_err_BUG("nhbr_polynomial: invalid preceding j");
48 124586 : return gerepileupto(ltop, nhbr_pol);
49 : }
50 :
51 : /* Path is an array with space for at least max_len+1 * elements, whose first
52 : * and second elements are the beginning of the path. I.e., the path starts
53 : * (path[0], path[1])
54 : * If the result is less than max_len, then the last element of path is on the
55 : * floor. If the result equals max_len, then it is unknown whether the last
56 : * element of path is on the floor or not */
57 : static long
58 244904 : extend_path(ulong path[], GEN phi, ulong p, ulong pi, long L, long max_len)
59 : {
60 244904 : pari_sp av = avma;
61 244904 : long d = 1;
62 315430 : for ( ; d < max_len; d++)
63 : {
64 91332 : GEN nhbr_pol = nhbr_polynomial(path + d, phi, p, pi, L);
65 91332 : ulong nhbr = Flx_oneroot(nhbr_pol, p);
66 91332 : set_avma(av);
67 91332 : if (nhbr == p) break; /* no root: we are on the floor. */
68 70526 : path[d+1] = nhbr;
69 : }
70 244904 : return d;
71 : }
72 :
73 : /* This is Sutherland 2009 Algorithm Ascend (p12) */
74 : ulong
75 111173 : ascend_volcano(GEN phi, ulong j, ulong p, ulong pi, long level, long L,
76 : long depth, long steps)
77 : {
78 111173 : pari_sp ltop = avma, av;
79 : /* path will never hold more than max_len < depth elements */
80 111173 : GEN path_g = cgetg(depth + 2, t_VECSMALL);
81 111173 : ulong *path = zv_to_ulongptr(path_g);
82 111173 : long max_len = depth - level;
83 111173 : int first_iter = 1;
84 :
85 111173 : if (steps <= 0 || max_len < 0) pari_err_BUG("ascend_volcano: bad params");
86 111173 : av = avma;
87 366773 : while (steps--)
88 : {
89 144427 : GEN nhbr_pol = first_iter? Flm_Fl_polmodular_evalx(phi, L, j, p, pi)
90 144427 : : nhbr_polynomial(path+1, phi, p, pi, L);
91 144427 : GEN nhbrs = Flx_roots(nhbr_pol, p);
92 144427 : long nhbrs_len = lg(nhbrs)-1, i;
93 144427 : pari_sp btop = avma;
94 144427 : path[0] = j;
95 144427 : first_iter = 0;
96 :
97 144427 : j = nhbrs[nhbrs_len];
98 182857 : for (i = 1; i < nhbrs_len; i++)
99 : {
100 68773 : ulong next_j = nhbrs[i], last_j;
101 : long len;
102 68773 : if (is_j_exceptional(next_j, p))
103 : {
104 : /* Fouquet & Morain, Section 4.3, if j = 0 or 1728, then it is on the
105 : * surface. So we just return it. */
106 0 : if (steps)
107 0 : pari_err_BUG("ascend_volcano: Got to the top with more steps to go!");
108 0 : j = next_j; break;
109 : }
110 68773 : path[1] = next_j;
111 68773 : len = extend_path(path, phi, p, pi, L, max_len);
112 68773 : last_j = path[len];
113 68773 : if (len == max_len
114 : /* Ended up on the surface */
115 68773 : && (is_j_exceptional(last_j, p)
116 68773 : || node_degree(phi, L, last_j, p, pi) > 1)) { j = next_j; break; }
117 38430 : set_avma(btop);
118 : }
119 144427 : path[1] = j; /* For nhbr_polynomial() at the top. */
120 :
121 144427 : max_len++; set_avma(av);
122 : }
123 111173 : return gc_long(ltop, j);
124 : }
125 :
126 : static void
127 179163 : random_distinct_neighbours_of(ulong *nhbr1, ulong *nhbr2,
128 : GEN phi, ulong j, ulong p, ulong pi, long L, long must_have_two_neighbours)
129 : {
130 179163 : pari_sp av = avma;
131 179163 : GEN modpol = Flm_Fl_polmodular_evalx(phi, L, j, p, pi);
132 : ulong rem;
133 179163 : *nhbr1 = Flx_oneroot(modpol, p);
134 179163 : if (*nhbr1 == p) pari_err_BUG("random_distinct_neighbours_of [no neighbour]");
135 179163 : modpol = Flx_div_by_X_x(modpol, *nhbr1, p, &rem);
136 179163 : *nhbr2 = Flx_oneroot(modpol, p);
137 179163 : if (must_have_two_neighbours && *nhbr2 == p)
138 0 : pari_err_BUG("random_distinct_neighbours_of [single neighbour]");
139 179163 : set_avma(av);
140 179163 : }
141 :
142 :
143 : /*
144 : * This is Sutherland 2009 Algorithm Descend (p12).
145 : */
146 : ulong
147 1835 : descend_volcano(GEN phi, ulong j, ulong p, ulong pi,
148 : long level, long L, long depth, long steps)
149 : {
150 1835 : pari_sp ltop = avma;
151 : GEN path_g;
152 : ulong *path;
153 : long max_len;
154 :
155 1835 : if (steps <= 0 || level + steps > depth) pari_err_BUG("descend_volcano");
156 1835 : max_len = depth - level;
157 1835 : path_g = cgetg(max_len + 1 + 1, t_VECSMALL);
158 1836 : path = zv_to_ulongptr(path_g);
159 1836 : path[0] = j;
160 : /* level = 0 means we're on the volcano surface... */
161 1836 : if (!level)
162 : {
163 : /* Look for any path to the floor. One of j's first three neighbours leads
164 : * to the floor, since at most two neighbours are on the surface. */
165 1589 : GEN nhbrs = Flx_roots(Flm_Fl_polmodular_evalx(phi, L, j, p, pi), p);
166 : long i;
167 1734 : for (i = 1; i <= 3; i++)
168 : {
169 : long len;
170 1734 : path[1] = nhbrs[i];
171 1734 : len = extend_path(path, phi, p, pi, L, max_len);
172 : /* If nhbrs[i] took us to the floor: */
173 1734 : if (len < max_len || node_degree(phi, L, path[len], p, pi) == 1) break;
174 : }
175 1588 : if (i > 3) pari_err_BUG("descend_volcano [2]");
176 : }
177 : else
178 : {
179 : ulong nhbr1, nhbr2;
180 : long len;
181 247 : random_distinct_neighbours_of(&nhbr1, &nhbr2, phi, j, p, pi, L, 1);
182 247 : path[1] = nhbr1;
183 247 : len = extend_path(path, phi, p, pi, L, max_len);
184 : /* If last j isn't on the floor */
185 247 : if (len == max_len /* Ended up on the surface. */
186 247 : && (is_j_exceptional(path[len], p)
187 221 : || node_degree(phi, L, path[len], p, pi) != 1)) {
188 : /* The other neighbour leads to the floor */
189 96 : path[1] = nhbr2;
190 96 : (void) extend_path(path, phi, p, pi, L, steps);
191 : }
192 : }
193 1835 : return gc_ulong(ltop, path[steps]);
194 : }
195 :
196 :
197 : long
198 178916 : j_level_in_volcano(
199 : GEN phi, ulong j, ulong p, ulong pi, long L, long depth)
200 : {
201 178916 : pari_sp av = avma;
202 : GEN chunk;
203 : ulong *path1, *path2;
204 : long lvl;
205 :
206 : /* Fouquet & Morain, Section 4.3, if j = 0 or 1728 then it is on the
207 : * surface. Also, if the volcano depth is zero then j has level 0 */
208 178916 : if (depth == 0 || is_j_exceptional(j, p)) return 0;
209 :
210 178916 : chunk = new_chunk(2 * (depth + 1));
211 178916 : path1 = (ulong *) &chunk[0];
212 178916 : path2 = (ulong *) &chunk[depth + 1];
213 178916 : path1[0] = path2[0] = j;
214 178916 : random_distinct_neighbours_of(&path1[1], &path2[1], phi, j, p, pi, L, 0);
215 178916 : if (path2[1] == p)
216 91889 : lvl = depth; /* Only one neighbour => j is on the floor => level = depth */
217 : else
218 : {
219 87027 : long path1_len = extend_path(path1, phi, p, pi, L, depth);
220 87027 : long path2_len = extend_path(path2, phi, p, pi, L, path1_len);
221 87027 : lvl = depth - path2_len;
222 : }
223 178916 : return gc_long(av, lvl);
224 : }
225 :
226 : #define vecsmall_len(v) (lg(v) - 1)
227 :
228 : INLINE GEN
229 28297429 : Flx_remove_root(GEN f, ulong a, ulong p)
230 : {
231 : ulong r;
232 28297429 : GEN g = Flx_div_by_X_x(f, a, p, &r);
233 28266829 : if (r) pari_err_BUG("Flx_remove_root");
234 28268845 : return g;
235 : }
236 :
237 : INLINE GEN
238 21787182 : get_nbrs(GEN phi, long L, ulong J, const ulong *xJ, ulong p, ulong pi)
239 : {
240 21787182 : pari_sp av = avma;
241 21787182 : GEN f = Flm_Fl_polmodular_evalx(phi, L, J, p, pi);
242 21768026 : if (xJ) f = Flx_remove_root(f, *xJ, p);
243 21751041 : return gerepileupto(av, Flx_roots(f, p));
244 : }
245 :
246 : /* Return a path of length n along the surface of an L-volcano of height h
247 : * starting from surface node j0. Assumes (D|L) = 1 where D = disc End(j0).
248 : *
249 : * Actually, if j0's endomorphism ring is a suborder, we return the
250 : * corresponding shorter path. W must hold space for n + h nodes.
251 : *
252 : * TODO: have two versions of this function: one that assumes J has the correct
253 : * endomorphism ring (hence avoiding several branches in the inner loop) and a
254 : * second that does not and accordingly checks for repetitions */
255 : static long
256 182904 : surface_path(
257 : ulong W[],
258 : long n, GEN phi, long L, long h, ulong J, const ulong *nJ, ulong p, ulong pi)
259 : {
260 182904 : pari_sp av = avma, bv;
261 : GEN T, v;
262 : long j, k, w, x;
263 : ulong W0;
264 :
265 182904 : W[0] = W0 = J;
266 182904 : if (n == 1) return 1;
267 :
268 182904 : T = cgetg(h+2, t_VEC); bv = avma;
269 182903 : v = get_nbrs(phi, L, J, nJ, p, pi);
270 : /* Insert known neighbour first */
271 182898 : if (nJ) v = gerepileupto(bv, vecsmall_prepend(v, *nJ));
272 182898 : gel(T,1) = v;
273 182898 : k = vecsmall_len(v);
274 :
275 182898 : switch (k) {
276 0 : case 0: pari_err_BUG("surface_path"); /* We must always have neighbours */
277 : case 1:
278 : /* If volcano is not flat, then we must have more than one neighbour */
279 5490 : if (h) pari_err_BUG("surface_path");
280 5490 : W[1] = uel(v, 1);
281 5490 : set_avma(av);
282 : /* Check for bad endo ring */
283 5490 : if (W[1] == W[0]) return 1;
284 5411 : return 2;
285 : case 2:
286 : /* If L=2 the only way we can have 2 neighbours is if we have a double root
287 : * which can only happen for |D| <= 16 (Thm 2.2 of Fouquet-Morain)
288 : * and if it does we must have a 2-cycle. Happens for D=-15. */
289 20795 : if (L == 2)
290 : { /* The double root is the neighbour on the surface, with exactly one
291 : * neighbour other than J; the other neighbour of J has either 0 or 2
292 : * neighbours that are not J */
293 0 : GEN u = get_nbrs(phi, L, uel(v, 1), &J, p, pi);
294 0 : long n = vecsmall_len(u) - !!vecsmall_isin(u, J);
295 0 : W[1] = n == 1 ? uel(v,1) : uel(v,2);
296 0 : return gc_long(av, 2);
297 : }
298 : /* Volcano is not flat but found only 2 neighbours for the surface node J */
299 20795 : if (h) pari_err_BUG("surface_path");
300 :
301 20856 : W[1] = uel(v,1); /* TODO: Can we use the other root uel(v,2) somehow? */
302 4358018 : for (w = 2; w < n; w++)
303 : {
304 4337328 : v = get_nbrs(phi, L, W[w-1], &W[w-2], p, pi);
305 : /* A flat volcano must have exactly one non-previous neighbour */
306 4337380 : if (vecsmall_len(v) != 1) pari_err_BUG("surface_path");
307 4337380 : W[w] = uel(v, 1);
308 : /* Detect cycle in case J doesn't have the right endo ring. */
309 4337380 : set_avma(av); if (W[w] == W0) return w;
310 : }
311 20690 : return gc_long(av, n);
312 : }
313 156613 : if (!h) pari_err_BUG("surface_path"); /* Can't have a flat volcano if k > 2 */
314 :
315 : /* At this point, each surface node has L+1 distinct neighbours, 2 of which
316 : * are on the surface */
317 156781 : w = 1;
318 5619679 : for (x = 0;; x++)
319 : {
320 : /* Get next neighbour of last known surface node to attempt to
321 : * extend the path. */
322 11082577 : W[w] = umael(T, ((w-1) % h) + 1, x + 1);
323 :
324 : /* Detect cycle in case the given J didn't have the right endo ring */
325 5619679 : if (W[w] == W0) return gc_long(av,w);
326 :
327 : /* If we have to test the last neighbour, we know it's on the
328 : * surface, and if we're done there's no need to extend. */
329 5619551 : if (x == k-1 && w == n-1) return gc_long(av,n);
330 :
331 : /* Walk forward until we hit the floor or finish. */
332 : /* NB: We don't keep the stack clean here; usage is in the order of Lh,
333 : * i.e. L roots for each level of the volcano of height h. */
334 5517046 : for (j = w;;)
335 11611623 : {
336 : long m;
337 : /* We must get 0 or L neighbours here. */
338 17128669 : v = get_nbrs(phi, L, W[j], &W[j-1], p, pi);
339 17122328 : m = vecsmall_len(v);
340 17122328 : if (!m) {
341 : /* We hit the floor: save the neighbours of W[w-1] and dump the rest */
342 5460138 : GEN nbrs = gel(T, ((w-1) % h) + 1);
343 5460138 : gel(T, ((w-1) % h) + 1) = gerepileupto(bv, nbrs);
344 5462898 : break;
345 : }
346 11662190 : if (m != L) pari_err_BUG("surface_path");
347 :
348 11665610 : gel(T, (j % h) + 1) = v;
349 :
350 11665610 : W[++j] = uel(v, 1);
351 : /* If we have our path by h nodes, we know W[w] is on the surface */
352 11665610 : if (j == w + h) {
353 10709282 : ++w;
354 : /* Detect cycle in case the given J didn't have the right endo ring */
355 10709282 : if (W[w] == W0) return gc_long(av,w);
356 10684743 : x = 0; k = L;
357 : }
358 11641071 : if (w == n) return gc_long(av,w);
359 : }
360 : }
361 : }
362 :
363 : long
364 115849 : next_surface_nbr(
365 : ulong *nJ,
366 : GEN phi, long L, long h, ulong J, const ulong *pJ, ulong p, ulong pi)
367 : {
368 115849 : pari_sp av = avma, bv;
369 : GEN S;
370 : ulong *P;
371 : long i, k;
372 :
373 115849 : S = get_nbrs(phi, L, J, pJ, p, pi);
374 115851 : k = vecsmall_len(S);
375 : /* If there is a double root and pJ is set, then k will be zero. */
376 115851 : if (!k) return gc_long(av,0);
377 115851 : if (k == 1 || ( ! pJ && k == 2)) { *nJ = uel(S, 1); return gc_long(av,1); }
378 18260 : if (!h) pari_err_BUG("next_surface_nbr");
379 :
380 18260 : P = (ulong *) new_chunk(h + 1); bv = avma;
381 18260 : P[0] = J;
382 35650 : for (i = 0; i < k; i++)
383 : {
384 : long j;
385 35650 : P[1] = uel(S, i + 1);
386 58579 : for (j = 1; j <= h; j++)
387 : {
388 40319 : GEN T = get_nbrs(phi, L, P[j], &P[j - 1], p, pi);
389 40319 : if (!vecsmall_len(T)) break;
390 22929 : P[j + 1] = uel(T, 1);
391 : }
392 35650 : if (j < h) pari_err_BUG("next_surface_nbr");
393 35650 : set_avma(bv); if (j > h) break;
394 : }
395 : /* TODO: We could save one get_nbrs call by iterating from i up to k-1 and
396 : * assume that the last (kth) nbr is the one we want. For now we're careful
397 : * and check that this last nbr really is on the surface */
398 18260 : if (i == k) pari_err_BUG("next_surf_nbr");
399 18260 : *nJ = uel(S, i+1); return gc_long(av,1);
400 : }
401 :
402 : /* Return the number of distinct neighbours (1 or 2) */
403 : INLINE long
404 176824 : common_nbr(ulong *nbr,
405 : ulong J1, GEN Phi1, long L1,
406 : ulong J2, GEN Phi2, long L2, ulong p, ulong pi)
407 : {
408 176824 : pari_sp av = avma;
409 : GEN d, f, g, r;
410 : long rlen;
411 :
412 176824 : g = Flm_Fl_polmodular_evalx(Phi1, L1, J1, p, pi);
413 176823 : f = Flm_Fl_polmodular_evalx(Phi2, L2, J2, p, pi);
414 176823 : d = Flx_gcd(f, g, p);
415 176819 : if (degpol(d) == 1) { *nbr = Flx_deg1_root(d, p); return gc_long(av,1); }
416 1107 : if (degpol(d) != 2) pari_err_BUG("common_neighbour");
417 1107 : r = Flx_roots(d, p);
418 1107 : rlen = vecsmall_len(r);
419 1107 : if (!rlen) pari_err_BUG("common_neighbour");
420 : /* rlen is 1 or 2 depending on whether the root is unique or not. */
421 1107 : nbr[0] = uel(r, 1);
422 1107 : nbr[1] = uel(r, rlen); return gc_long(av,rlen);
423 : }
424 :
425 : /* Return gcd(Phi1(X,J1)/(X - J0), Phi2(X,J2)). Not stack clean. */
426 : INLINE GEN
427 42212 : common_nbr_pred_poly(
428 : ulong J1, GEN Phi1, long L1,
429 : ulong J2, GEN Phi2, long L2, ulong J0, ulong p, ulong pi)
430 : {
431 : GEN f, g;
432 :
433 42212 : g = Flm_Fl_polmodular_evalx(Phi1, L1, J1, p, pi);
434 42212 : g = Flx_remove_root(g, J0, p);
435 42212 : f = Flm_Fl_polmodular_evalx(Phi2, L2, J2, p, pi);
436 42212 : return Flx_gcd(f, g, p);
437 : }
438 :
439 : /* Find common neighbour of J1 and J2, where J0 is an L1 predecessor of J1.
440 : * Return 1 if successful, 0 if not. */
441 : INLINE int
442 41181 : common_nbr_pred(ulong *nbr,
443 : ulong J1, GEN Phi1, long L1,
444 : ulong J2, GEN Phi2, long L2, ulong J0, ulong p, ulong pi)
445 : {
446 41181 : pari_sp av = avma;
447 41181 : GEN d = common_nbr_pred_poly(J1, Phi1, L1, J2, Phi2, L2, J0, p, pi);
448 41181 : int res = (degpol(d) == 1);
449 41181 : if (res) *nbr = Flx_deg1_root(d, p);
450 41181 : return gc_bool(av,res);
451 : }
452 :
453 : INLINE long
454 1031 : common_nbr_verify(ulong *nbr,
455 : ulong J1, GEN Phi1, long L1,
456 : ulong J2, GEN Phi2, long L2, ulong J0, ulong p, ulong pi)
457 : {
458 1031 : pari_sp av = avma;
459 1031 : GEN d = common_nbr_pred_poly(J1, Phi1, L1, J2, Phi2, L2, J0, p, pi);
460 :
461 1031 : if (!degpol(d)) return gc_long(av,0);
462 377 : if (degpol(d) > 1) pari_err_BUG("common_neighbour_verify");
463 377 : *nbr = Flx_deg1_root(d, p);
464 377 : return gc_long(av,1);
465 : }
466 :
467 : INLINE ulong
468 500 : Flm_Fl_polmodular_evalxy(GEN Phi, long L, ulong x, ulong y, ulong p, ulong pi)
469 : {
470 500 : pari_sp av = avma;
471 500 : GEN f = Flm_Fl_polmodular_evalx(Phi, L, x, p, pi);
472 500 : return gc_ulong(av, Flx_eval_pre(f, y, p, pi));
473 : }
474 :
475 : /* Find a common L1-neighbor of J1 and L2-neighbor of J2, given J0 an
476 : * L2-neighbor of J1 and an L1-neighbor of J2. Return 1 if successful, 0
477 : * otherwise. Will only fail if initial J-invariant had the wrong endo ring */
478 : INLINE int
479 22080 : common_nbr_corner(ulong *nbr,
480 : ulong J1, GEN Phi1, long L1, long h1,
481 : ulong J2, GEN Phi2, long L2, ulong J0, ulong p, ulong pi)
482 : {
483 : ulong nbrs[2];
484 22080 : if (common_nbr(nbrs, J1,Phi1,L1, J2,Phi2,L2, p, pi) == 2)
485 : {
486 : ulong nJ1, nJ2;
487 606 : if (!next_surface_nbr(&nJ2, Phi1, L1, h1, J2, &J0, p, pi) ||
488 386 : !next_surface_nbr(&nJ1, Phi1, L1, h1, nbrs[0], &J1, p, pi)) return 0;
489 :
490 303 : if (Flm_Fl_polmodular_evalxy(Phi2, L2, nJ1, nJ2, p, pi))
491 106 : nbrs[0] = nbrs[1];
492 394 : else if (!next_surface_nbr(&nJ1, Phi1, L1, h1, nbrs[1], &J1, p, pi) ||
493 280 : !Flm_Fl_polmodular_evalxy(Phi2, L2, nJ1, nJ2, p, pi)) return 0;
494 : }
495 21997 : *nbr = nbrs[0]; return 1;
496 : }
497 :
498 : /* Enumerate a surface L1-cycle using gcds with Phi_L2, where c_L2=c_L1^e and
499 : * |c_L1|=n, where c_a is the class of the pos def reduced primeform <a,b,c>.
500 : * Assumes n > e > 1 and roots[0],...,roots[e-1] are already present in W */
501 : static long
502 58515 : surface_gcd_cycle(
503 : ulong W[], ulong V[], long n,
504 : GEN Phi1, long L1, GEN Phi2, long L2, long e, ulong p, ulong pi)
505 : {
506 58515 : pari_sp av = avma;
507 : long i1, i2, j1, j2;
508 :
509 58515 : i1 = j2 = 0;
510 58515 : i2 = j1 = e - 1;
511 : /* If W != V we assume V actually points to an L2-isogenous parallel L1-path.
512 : * e should be 2 in this case */
513 58515 : if (W != V) { i1 = j1+1; i2 = n-1; }
514 : do {
515 : ulong t0, t1, t2, h10, h11, h20, h21;
516 : long k;
517 : GEN f, g, h1, h2;
518 :
519 3405191 : f = Flm_Fl_polmodular_evalx(Phi2, L2, V[i1], p, pi);
520 3393287 : g = Flm_Fl_polmodular_evalx(Phi1, L1, W[j1], p, pi);
521 3397426 : g = Flx_remove_root(g, W[j1 - 1], p);
522 3392584 : h1 = Flx_gcd(f, g, p);
523 3384842 : if (degpol(h1) != 1) break; /* Error */
524 3384469 : h11 = Flx_coeff(h1, 1);
525 3386216 : h10 = Flx_coeff(h1, 0); set_avma(av);
526 :
527 3386106 : f = Flm_Fl_polmodular_evalx(Phi2, L2, V[i2], p, pi);
528 3393618 : g = Flm_Fl_polmodular_evalx(Phi1, L1, W[j2], p, pi);
529 3397820 : k = j2 + 1;
530 3397820 : if (k == n) k = 0;
531 3397820 : g = Flx_remove_root(g, W[k], p);
532 3393064 : h2 = Flx_gcd(f, g, p);
533 3385543 : if (degpol(h2) != 1) break; /* Error */
534 3385283 : h21 = Flx_coeff(h2, 1);
535 3387373 : h20 = Flx_coeff(h2, 0); set_avma(av);
536 :
537 3387077 : i1++; i2--; if (i2 < 0) i2 = n-1;
538 3387077 : j1++; j2--; if (j2 < 0) j2 = n-1;
539 :
540 3387077 : t0 = Fl_mul_pre(h11, h21, p, pi);
541 3400886 : t1 = Fl_inv(t0, p);
542 3404322 : t0 = Fl_mul_pre(t1, h21, p, pi);
543 3404484 : t2 = Fl_mul_pre(t0, h10, p, pi);
544 3405196 : W[j1] = Fl_neg(t2, p);
545 3404348 : t0 = Fl_mul_pre(t1, h11, p, pi);
546 3405593 : t2 = Fl_mul_pre(t0, h20, p, pi);
547 3405905 : W[j2] = Fl_neg(t2, p);
548 3405193 : } while (j2 > j1 + 1);
549 : /* Usually the loop exits when j2 = j1 + 1, in which case we return n.
550 : * If we break early because of an error, then (j2 - (j1+1)) > 0 is the
551 : * number of elements we haven't calculated yet, and we return n minus that
552 : * quantity */
553 58517 : return gc_long(av, n - j2 + j1 + 1);
554 : }
555 :
556 : static long
557 1176 : surface_gcd_path(
558 : ulong W[], ulong V[], long n,
559 : GEN Phi1, long L1, GEN Phi2, long L2, long e, ulong p, ulong pi)
560 : {
561 1176 : pari_sp av = avma;
562 : long i, j;
563 :
564 1176 : i = 0; j = e;
565 : /* If W != V then assume V actually points to a L2-isogenous
566 : * parallel L1-path. e should be 2 in this case */
567 1176 : if (W != V) i = j;
568 6048 : while (j < n)
569 : {
570 : GEN f, g, d;
571 :
572 3696 : f = Flm_Fl_polmodular_evalx(Phi2, L2, V[i], p, pi);
573 3696 : g = Flm_Fl_polmodular_evalx(Phi1, L1, W[j - 1], p, pi);
574 3696 : g = Flx_remove_root(g, W[j - 2], p);
575 3696 : d = Flx_gcd(f, g, p);
576 3696 : if (degpol(d) != 1) break; /* Error */
577 3696 : W[j] = Flx_deg1_root(d, p);
578 3696 : i++; j++; set_avma(av);
579 : }
580 1176 : return gc_long(av, j);
581 : }
582 :
583 : /* Given a path V of length n on an L1-volcano, and W[0] L2-isogenous to V[0],
584 : * extends the path W to length n on an L1-volcano, with W[i] L2-isogenous
585 : * to V[i]. Uses gcds unless L2 is too large to make it helpful. Always uses
586 : * GCD to get W[1] to ensure consistent orientation.
587 : *
588 : * Returns the new length of W. This will almost always be n, but could be
589 : * lower if V was started with a J-invariant with bad endomorphism ring */
590 : INLINE long
591 154743 : surface_parallel_path(
592 : ulong W[], ulong V[], long n,
593 : GEN Phi1, long L1, GEN Phi2, long L2, ulong p, ulong pi, long cycle)
594 : {
595 : ulong W2, nbrs[2];
596 154743 : if (common_nbr(nbrs, W[0], Phi1, L1, V[1], Phi2, L2, p, pi) == 2)
597 : {
598 707 : if (n <= 2) return 1; /* Error: Two choices with n = 2; ambiguous */
599 706 : if (!common_nbr_verify(&W2,nbrs[0], Phi1,L1,V[2], Phi2,L2,W[0], p,pi))
600 381 : nbrs[0] = nbrs[1]; /* nbrs[1] must be the correct choice */
601 325 : else if (common_nbr_verify(&W2,nbrs[1], Phi1,L1,V[2], Phi2,L2,W[0], p,pi))
602 52 : return 1; /* Error: Both paths extend successfully */
603 : }
604 154688 : W[1] = nbrs[0];
605 154688 : if (n <= 2) return n;
606 : return cycle? surface_gcd_cycle(W, V, n, Phi1, L1, Phi2, L2, 2, p, pi)
607 59691 : : surface_gcd_path (W, V, n, Phi1, L1, Phi2, L2, 2, p, pi);
608 : }
609 :
610 : GEN
611 183504 : enum_roots(ulong J0, norm_eqn_t ne, GEN fdb, classgp_pcp_t G)
612 : {
613 : /* MAX_HEIGHT >= max_{p,n} val_p(n) where p and n are ulongs */
614 : enum { MAX_HEIGHT = BITS_IN_LONG };
615 183504 : pari_sp av, ltop = avma;
616 183504 : long s = !!G->L0;
617 183504 : long *n = G->n + s, *L = G->L + s, *o = G->o + s, k = G->k - s;
618 183504 : long i, t, vlen, *e, *h, *off, *poff, *M, N = G->enum_cnt;
619 183504 : ulong p = ne->p, pi = ne->pi, *roots;
620 : GEN Phi, vshape, vp, ve, roots_;
621 :
622 183504 : if (!k) return mkvecsmall(J0);
623 :
624 182905 : roots_ = cgetg(N + MAX_HEIGHT, t_VECSMALL);
625 182905 : roots = zv_to_ulongptr(roots_);
626 182905 : av = avma;
627 :
628 : /* TODO: Shouldn't be factoring this every time. Store in *ne? */
629 182905 : vshape = factoru(ne->v);
630 182903 : vp = gel(vshape, 1);
631 182903 : ve = gel(vshape, 2);
632 :
633 182903 : vlen = vecsmall_len(vp);
634 182903 : Phi = new_chunk(k);
635 182904 : e = new_chunk(k);
636 182904 : off = new_chunk(k);
637 182904 : poff = new_chunk(k);
638 : /* TODO: Surely we can work these out ahead of time? */
639 : /* h[i] is the valuation of p[i] in v */
640 182904 : h = new_chunk(k);
641 415667 : for (i = 0; i < k; ++i) {
642 232764 : h[i] = 0;
643 341442 : for (t = 1; t <= vlen; ++t)
644 270968 : if (vp[t] == L[i]) { h[i] = uel(ve, t); break; }
645 232764 : e[i] = 0;
646 232764 : off[i] = 0;
647 232764 : gel(Phi, i) = polmodular_db_getp(fdb, L[i], p);
648 : }
649 :
650 182903 : M = new_chunk(k);
651 182903 : for (M[0] = 1, i = 1; i < k; ++i) M[i] = M[i-1] * n[i-1];
652 :
653 182903 : t = surface_path(roots, n[0], gel(Phi, 0), L[0], h[0], J0, NULL, p, pi);
654 : /* Error: J0 has bad endo ring */
655 182903 : if (t < n[0]) return gc_NULL(ltop);
656 182548 : if (k == 1) { set_avma(av); setlg(roots_, t + 1); return roots_; }
657 :
658 43177 : i = 1;
659 241043 : while (i < k) {
660 : long j, t0;
661 154871 : for (j = i + 1; j < k && ! e[j]; ++j);
662 154871 : if (j < k) {
663 63261 : if (e[i]) {
664 329448 : if (! common_nbr_pred(
665 164724 : &roots[t], roots[off[i]], gel(Phi,i), L[i],
666 164724 : roots[t - M[j]], gel(Phi, j), L[j], roots[poff[i]], p, pi)) {
667 0 : break; /* Error: J0 has bad endo ring */
668 : }
669 198720 : } else if ( ! common_nbr_corner(
670 110400 : &roots[t], roots[off[i]], gel(Phi,i), L[i], h[i],
671 88320 : roots[t - M[j]], gel(Phi, j), L[j], roots[poff[j]], p, pi)) {
672 83 : break; /* Error: J0 has bad endo ring */
673 : }
674 592075 : } else if ( ! next_surface_nbr(
675 366440 : &roots[t], gel(Phi,i), L[i], h[i],
676 225635 : roots[off[i]], e[i] ? &roots[poff[i]] : NULL, p, pi))
677 0 : break; /* Error: J0 has bad endo ring */
678 154790 : if (roots[t] == roots[0]) break; /* Error: J0 has bad endo ring */
679 :
680 154743 : poff[i] = off[i];
681 154743 : off[i] = t;
682 154743 : e[i]++;
683 154743 : for (j = i-1; j; --j) { e[j] = 0; off[j] = off[j+1]; }
684 :
685 464229 : t0 = surface_parallel_path(&roots[t], &roots[poff[i]], n[0],
686 464229 : gel(Phi, 0), L[0], gel(Phi, i), L[i], p, pi, n[0] == o[0]);
687 154741 : if (t0 < n[0]) break; /* Error: J0 has bad endo ring */
688 :
689 : /* TODO: Do I need to check if any of the new roots is a repeat in
690 : * the case where J0 has bad endo ring? */
691 154689 : t += n[0];
692 154689 : for (i = 1; i < k && e[i] == n[i]-1; i++);
693 : }
694 : /* Check if J0 had wrong endo ring */
695 43177 : if (t != N) return gc_NULL(ltop);
696 42995 : set_avma(av); setlg(roots_, t + 1); return roots_;
697 : }
|
