Line data Source code
1 : /* Copyright (C) 2014 The PARI group.
2 :
3 : This file is part of the PARI/GP package.
4 :
5 : PARI/GP is free software; you can redistribute it and/or modify it under the
6 : terms of the GNU General Public License as published by the Free Software
7 : Foundation; either version 2 of the License, or (at your option) any later
8 : version. It is distributed in the hope that it will be useful, but WITHOUT
9 : ANY WARRANTY WHATSOEVER.
10 :
11 : Check the License for details. You should have received a copy of it, along
12 : with the package; see the file 'COPYING'. If not, write to the Free Software
13 : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
14 :
15 : #include "pari.h"
16 : #include "paripriv.h"
17 :
18 : #define DEBUGLEVEL DEBUGLEVEL_ellisogeny
19 :
20 : /* Return 1 if the point Q is a Weierstrass (2-torsion) point of the
21 : * curve E, return 0 otherwise */
22 : static long
23 903 : ellisweierstrasspoint(GEN E, GEN Q)
24 903 : { return ell_is_inf(Q) || gequal0(ec_dmFdy_evalQ(E, Q)); }
25 :
26 : /* Given an elliptic curve E = [a1, a2, a3, a4, a6] and t,w in the ring of
27 : * definition of E, return the curve
28 : * E' = [a1, a2, a3, a4 - 5t, a6 - (E.b2 t + 7w)] */
29 : static GEN
30 12425 : make_velu_curve(GEN E, GEN t, GEN w)
31 : {
32 12425 : GEN A4, A6, a1 = ell_get_a1(E), a2 = ell_get_a2(E), a3 = ell_get_a3(E);
33 12425 : A4 = gsub(ell_get_a4(E), gmulsg(5L, t));
34 12425 : A6 = gsub(ell_get_a6(E), gadd(gmul(ell_get_b2(E), t), gmulsg(7L, w)));
35 12425 : return mkvec5(a1,a2,a3,A4,A6);
36 : }
37 :
38 : /* If phi = (f(x)/h(x)^2, g(x,y)/h(x)^3) is an isogeny, return the
39 : * variables x and y in a vecsmall */
40 : INLINE void
41 1764 : get_isog_vars(GEN phi, long *vx, long *vy)
42 : {
43 1764 : *vx = varn(gel(phi, 1));
44 1764 : *vy = varn(gel(phi, 2));
45 1764 : if (*vy == *vx) *vy = gvar2(gel(phi,2));
46 1764 : }
47 :
48 : /* x must be nonzero */
49 3780 : INLINE long _degree(GEN x) { return typ(x)==t_POL ? degpol(x): 0; }
50 :
51 : static GEN
52 5040 : RgH_eval(GEN P, GEN A, GEN B)
53 : {
54 5040 : if (typ(P)==t_POL)
55 : {
56 3493 : if (signe(P)==0) return mkvec2(P, gen_1);
57 3493 : return mkvec2(RgX_homogenous_evalpow(P, A, B), gel(B,degpol(P)+1));
58 : } else
59 1547 : return mkvec2(P, gen_1);
60 : }
61 :
62 : /* Given isogenies F:E' -> E and G:E'' -> E', return the composite
63 : * isogeny F o G:E'' -> E */
64 : static GEN
65 1267 : ellcompisog(GEN F, GEN G)
66 : {
67 1267 : pari_sp av = avma;
68 : GEN Gh, Gh2, Gh3, f, g, h, h2, h3, den, num;
69 : GEN K, K2, K3, F0, F1, g0, g1, Gp;
70 : long vx, vy, d;
71 1267 : checkellisog(F);
72 1260 : checkellisog(G);
73 1260 : get_isog_vars(F, &vx, &vy);
74 1260 : Gh = gel(G,3); Gh2 = gsqr(Gh); Gh3 = gmul(Gh, Gh2);
75 1260 : F0 = polcoef_i(gel(F,2), 0, vy);
76 1260 : F1 = polcoef_i(gel(F,2), 1, vy);
77 1260 : d = maxss(maxss(degpol(gel(F,1)),_degree(gel(F,3))),
78 : maxss(_degree(F0),_degree(F1)));
79 1260 : Gp = gpowers(Gh2, d);
80 1260 : f = RgH_eval(gel(F,1), gel(G,1), Gp);
81 1260 : g0 = RgH_eval(F0, gel(G,1), Gp);
82 1260 : g1 = RgH_eval(F1, gel(G,1), Gp);
83 1260 : h = RgH_eval(gel(F,3), gel(G,1), Gp);
84 1260 : K = gmul(gel(h,1), Gh);
85 1260 : K = RgX_normalize(RgX_div(K, RgX_gcd(K,RgX_deriv(K))));
86 1260 : K2 = gsqr(K); K3 = gmul(K, K2);
87 1260 : h2 = mkvec2(gsqr(gel(h,1)), gsqr(gel(h,2)));
88 1260 : h3 = mkvec2(gmul(gel(h,1),gel(h2,1)), gmul(gel(h,2),gel(h2,2)));
89 1260 : f = gdiv(gmul(gmul(K2, gel(f,1)),gel(h2,2)), gmul(gel(f,2), gel(h2,1)));
90 1260 : den = gmul(Gh3, gel(g1,2));
91 1260 : num = gadd(gmul(gel(g0,1),den), gmul(gmul(gel(G,2),gel(g1,1)),gel(g0,2)));
92 1260 : g = gdiv(gmul(gmul(K3,num),gel(h3,2)),gmul(gmul(gel(g0,2),den), gel(h3,1)));
93 1260 : return gerepilecopy(av, mkvec3(f,g,K));
94 : }
95 :
96 : static GEN
97 4760 : to_RgX(GEN P, long vx)
98 : {
99 4760 : return typ(P) == t_POL && varn(P)==vx? P: scalarpol_shallow(P, vx);
100 : }
101 :
102 : static GEN
103 476 : divy(GEN P0, GEN P1, GEN Q, GEN T, long vy)
104 : {
105 476 : GEN DP0, P0r = Q_remove_denom(P0, &DP0), P0D;
106 476 : GEN DP1, P1r = Q_remove_denom(P1, &DP1), P1D;
107 476 : GEN DQ, Qr = Q_remove_denom(Q, &DQ), P2;
108 476 : P0D = RgXQX_div(P0r, Qr, T);
109 476 : if (DP0) P0D = gdiv(P0D, DP0);
110 476 : P1D = RgXQX_div(P1r, Qr, T);
111 476 : if (DP1) P1D = gdiv(P1D, DP1);
112 476 : P2 = gadd(gmul(P1D, pol_x(vy)), P0D);
113 476 : if (DQ) P2 = gmul(P2, DQ);
114 476 : return P2;
115 : }
116 :
117 : static GEN
118 1904 : QXQH_eval(GEN P, GEN A, GEN B, GEN T)
119 : {
120 1904 : if (signe(P)==0) return mkvec2(P, pol_1(varn(P)));
121 1652 : return mkvec2(QXQX_homogenous_evalpow(P, A, B, T), gel(B,degpol(P)+1));
122 : }
123 :
124 : static GEN
125 1722 : ellnfcompisog(GEN nf, GEN F, GEN G)
126 : {
127 1722 : pari_sp av = avma;
128 : GEN Gh, Gh2, Gh3, f, g, gd, h, h21, h22, h31, h32, den;
129 : GEN K, K2, K3, F0, F1, G0, G1, g0, g1, Gp;
130 : GEN num0, num1, gn0, gn1;
131 : GEN g0d, g01, k3h32;
132 : GEN T;
133 : pari_timer ti;
134 : long vx, vy, d;
135 1722 : if (!nf) return ellcompisog(F, G);
136 476 : T = nf_get_pol(nf);
137 476 : timer_start(&ti);
138 476 : checkellisog(F); F = liftpol_shallow(F);
139 476 : checkellisog(G); G = liftpol_shallow(G);
140 476 : get_isog_vars(F, &vx, &vy);
141 476 : Gh = to_RgX(gel(G,3),vx); Gh2 = QXQX_sqr(Gh, T); Gh3 = QXQX_mul(Gh, Gh2, T);
142 476 : F0 = to_RgX(polcoef_i(gel(F,2), 0, vy), vx);
143 476 : F1 = to_RgX(polcoef_i(gel(F,2), 1, vy), vx);
144 476 : G0 = to_RgX(polcoef_i(gel(G,2), 0, vy), vx);
145 476 : G1 = to_RgX(polcoef_i(gel(G,2), 1, vy), vx);
146 476 : d = maxss(maxss(degpol(gel(F,1)),degpol(gel(F,3))),maxss(degpol(F0),degpol(F1)));
147 476 : Gp = QXQX_powers(Gh2, d, T);
148 476 : f = QXQH_eval(to_RgX(gel(F,1),vx), gel(G,1), Gp, T);
149 476 : g0 = QXQH_eval(F0, to_RgX(gel(G,1),vx), Gp, T);
150 476 : g1 = QXQH_eval(F1, to_RgX(gel(G,1),vx), Gp, T);
151 476 : h = QXQH_eval(to_RgX(gel(F,3),vx), gel(G,1), Gp, T);
152 476 : K = Q_remove_denom(QXQX_mul(to_RgX(gel(h,1),vx), Gh, T), NULL);
153 476 : K = RgXQX_div(K, nfgcd(K, RgX_deriv(K), T, NULL), T);
154 476 : K = RgX_normalize(K);
155 476 : if (DEBUGLEVEL) timer_printf(&ti,"ellnfcompisog: nfgcd");
156 476 : K2 = QXQX_sqr(K, T); K3 = QXQX_mul(K, K2, T);
157 476 : if (DEBUGLEVEL) timer_printf(&ti,"ellnfcompisog: evalpow");
158 476 : h21 = QXQX_sqr(gel(h,1),T);
159 476 : h22 = QXQX_sqr(gel(h,2),T);
160 476 : h31 = QXQX_mul(gel(h,1), h21,T);
161 476 : h32 = QXQX_mul(gel(h,2), h22,T);
162 476 : if (DEBUGLEVEL) timer_printf(&ti,"h");
163 476 : f = RgXQX_div(QXQX_mul(QXQX_mul(K2, gel(f,1), T), h22, T),
164 476 : QXQX_mul(gel(f,2), h21, T), T);
165 476 : if (DEBUGLEVEL) timer_printf(&ti,"f");
166 476 : den = QXQX_mul(Gh3, gel(g1,2), T);
167 476 : if (DEBUGLEVEL) timer_printf(&ti,"ellnfcompisog: den");
168 476 : g0d = QXQX_mul(gel(g0,1),den, T);
169 476 : g01 = QXQX_mul(gel(g1,1),gel(g0,2),T);
170 476 : num0 = RgX_add(g0d, QXQX_mul(G0,g01, T));
171 476 : num1 = QXQX_mul(G1,g01, T);
172 476 : if (DEBUGLEVEL) timer_printf(&ti,"ellnfcompisog: num");
173 476 : k3h32 = QXQX_mul(K3,h32,T);
174 476 : gn0 = QXQX_mul(num0, k3h32, T);
175 476 : gn1 = QXQX_mul(num1, k3h32, T);
176 476 : if (DEBUGLEVEL) timer_printf(&ti,"ellnfcompisog: gn");
177 476 : gd = QXQX_mul(QXQX_mul(gel(g0,2), den, T), h31, T);
178 476 : if (DEBUGLEVEL) timer_printf(&ti,"ellnfcompisog: gd");
179 476 : g = divy(gn0, gn1, gd, T, vy);
180 476 : if (DEBUGLEVEL) timer_printf(&ti,"ellnfcompisog: divy");
181 476 : return gerepilecopy(av, gmul(mkvec3(f,g,K),gmodulo(gen_1,T)));
182 : }
183 :
184 : /* Given an isogeny phi from ellisogeny() and a point P in the domain of phi,
185 : * return phi(P) */
186 : GEN
187 70 : ellisogenyapply(GEN phi, GEN P)
188 : {
189 70 : pari_sp ltop = avma;
190 : GEN f, g, h, img_f, img_g, img_h, img_h2, img_h3, img, tmp;
191 : long vx, vy;
192 70 : if (lg(P) == 4) return ellcompisog(phi,P);
193 49 : checkellisog(phi);
194 49 : if (!checkellpt_i(P)) pari_err_TYPE("ellisogenyapply",P);
195 42 : if (ell_is_inf(P)) return ellinf();
196 28 : f = gel(phi, 1);
197 28 : g = gel(phi, 2);
198 28 : h = gel(phi, 3);
199 28 : get_isog_vars(phi, &vx, &vy);
200 28 : img_h = poleval(h, gel(P, 1));
201 28 : if (gequal0(img_h)) { set_avma(ltop); return ellinf(); }
202 :
203 21 : img_h2 = gsqr(img_h);
204 21 : img_h3 = gmul(img_h, img_h2);
205 21 : img_f = poleval(f, gel(P, 1));
206 : /* FIXME: This calculation of g is perhaps not as efficient as it could be */
207 21 : tmp = gsubst(g, vx, gel(P, 1));
208 21 : img_g = gsubst(tmp, vy, gel(P, 2));
209 21 : img = cgetg(3, t_VEC);
210 21 : gel(img, 1) = gdiv(img_f, img_h2);
211 21 : gel(img, 2) = gdiv(img_g, img_h3);
212 21 : return gerepileupto(ltop, img);
213 : }
214 :
215 : /* isog = [f, g, h] = [x, y, 1] = identity */
216 : static GEN
217 252 : isog_identity(long vx, long vy)
218 252 : { return mkvec3(pol_x(vx), pol_x(vy), pol_1(vx)); }
219 :
220 : /* Returns an updated value for isog based (primarily) on tQ and uQ. Used to
221 : * iteratively compute the isogeny corresponding to a subgroup generated by a
222 : * given point. Ref: Equation 8 in Velu's paper.
223 : * isog = NULL codes the identity */
224 : static GEN
225 532 : update_isogeny_polys(GEN isog, GEN E, GEN Q, GEN tQ, GEN uQ, long vx, long vy)
226 : {
227 532 : pari_sp ltop = avma, av;
228 532 : GEN xQ = gel(Q, 1), yQ = gel(Q, 2);
229 532 : GEN rt = deg1pol_shallow(gen_1, gneg(xQ), vx);
230 532 : GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E);
231 :
232 532 : GEN gQx = ec_dFdx_evalQ(E, Q);
233 532 : GEN gQy = ec_dFdy_evalQ(E, Q);
234 : GEN tmp1, tmp2, tmp3, tmp4, f, g, h, rt_sqr, res;
235 :
236 : /* g -= uQ * (2 * y + E.a1 * x + E.a3)
237 : * + tQ * rt * (E.a1 * rt + y - yQ)
238 : * + rt * (E.a1 * uQ - gQx * gQy) */
239 532 : av = avma;
240 532 : tmp1 = gmul(uQ, gadd(deg1pol_shallow(gen_2, gen_0, vy),
241 : deg1pol_shallow(a1, a3, vx)));
242 532 : tmp1 = gerepileupto(av, tmp1);
243 532 : av = avma;
244 532 : tmp2 = gmul(tQ, gadd(gmul(a1, rt),
245 : deg1pol_shallow(gen_1, gneg(yQ), vy)));
246 532 : tmp2 = gerepileupto(av, tmp2);
247 532 : av = avma;
248 532 : tmp3 = gsub(gmul(a1, uQ), gmul(gQx, gQy));
249 532 : tmp3 = gerepileupto(av, tmp3);
250 :
251 532 : if (!isog) isog = isog_identity(vx,vy);
252 532 : f = gel(isog, 1);
253 532 : g = gel(isog, 2);
254 532 : h = gel(isog, 3);
255 532 : rt_sqr = gsqr(rt);
256 532 : res = cgetg(4, t_VEC);
257 532 : av = avma;
258 532 : tmp4 = gdiv(gadd(gmul(tQ, rt), uQ), rt_sqr);
259 532 : gel(res, 1) = gerepileupto(av, gadd(f, tmp4));
260 532 : av = avma;
261 532 : tmp4 = gadd(tmp1, gmul(rt, gadd(tmp2, tmp3)));
262 532 : gel(res, 2) = gerepileupto(av, gsub(g, gdiv(tmp4, gmul(rt, rt_sqr))));
263 532 : av = avma;
264 532 : gel(res, 3) = gerepileupto(av, gmul(h, rt));
265 532 : return gerepileupto(ltop, res);
266 : }
267 :
268 : /* Given a point P on E, return the curve E/<P> and, if only_image is zero,
269 : * the isogeny pi: E -> E/<P>. The variables vx and vy are used to describe
270 : * the isogeny (ignored if only_image is zero) */
271 : static GEN
272 427 : isogeny_from_kernel_point(GEN E, GEN P, int only_image, long vx, long vy)
273 : {
274 427 : pari_sp av = avma;
275 : GEN isog, EE, f, g, h, h2, h3;
276 427 : GEN Q = P, t = gen_0, w = gen_0;
277 427 : if (!oncurve(E,P))
278 7 : pari_err_DOMAIN("isogeny_from_kernel_point", "point", "not on", E, P);
279 420 : if (ell_is_inf(P))
280 : {
281 42 : if (only_image) return E;
282 28 : return mkvec2(E, isog_identity(vx,vy));
283 : }
284 :
285 378 : isog = NULL;
286 : for (;;)
287 525 : {
288 903 : GEN tQ, xQ = gel(Q,1), uQ = ec_2divpol_evalx(E, xQ);
289 903 : int stop = 0;
290 903 : if (ellisweierstrasspoint(E,Q))
291 : { /* ord(P)=2c; take Q=[c]P into consideration and stop */
292 196 : tQ = ec_dFdx_evalQ(E, Q);
293 196 : stop = 1;
294 : }
295 : else
296 707 : tQ = ec_half_deriv_2divpol_evalx(E, xQ);
297 903 : t = gadd(t, tQ);
298 903 : w = gadd(w, gadd(uQ, gmul(tQ, xQ)));
299 903 : if (!only_image) isog = update_isogeny_polys(isog, E, Q,tQ,uQ, vx,vy);
300 903 : if (stop) break;
301 :
302 707 : Q = elladd(E, P, Q);
303 : /* IF x([c]P) = x([c-1]P) THEN [c]P = -[c-1]P and [2c-1]P = 0 */
304 707 : if (gequal(gel(Q,1), xQ)) break;
305 525 : if (gc_needed(av,1))
306 : {
307 0 : if(DEBUGMEM>1) pari_warn(warnmem,"isogeny_from_kernel_point");
308 0 : gerepileall(av, isog? 4: 3, &Q, &t, &w, &isog);
309 : }
310 : }
311 :
312 378 : EE = make_velu_curve(E, t, w);
313 378 : if (only_image) return EE;
314 :
315 224 : if (!isog) isog = isog_identity(vx,vy);
316 224 : f = gel(isog, 1);
317 224 : g = gel(isog, 2);
318 224 : if ( ! (typ(f) == t_RFRAC && typ(g) == t_RFRAC))
319 0 : pari_err_BUG("isogeny_from_kernel_point (f or g has wrong type)");
320 :
321 : /* Clean up the isogeny polynomials f, g and h so that the isogeny
322 : * is given by (x,y) -> (f(x)/h(x)^2, g(x,y)/h(x)^3) */
323 224 : h = gel(isog, 3);
324 224 : h2 = gsqr(h);
325 224 : h3 = gmul(h, h2);
326 224 : f = gmul(f, h2);
327 224 : g = gmul(g, h3);
328 224 : if (typ(f) != t_POL || typ(g) != t_POL)
329 0 : pari_err_BUG("isogeny_from_kernel_point (wrong denominator)");
330 224 : return mkvec2(EE, mkvec3(f,g, gel(isog,3)));
331 : }
332 :
333 : /* Given a t_POL x^n - s1 x^{n-1} + s2 x^{n-2} - s3 x^{n-3} + ...
334 : * return the first three power sums (Newton's identities):
335 : * p1 = s1
336 : * p2 = s1^2 - 2 s2
337 : * p3 = (s1^2 - 3 s2) s1 + 3 s3 */
338 : static void
339 12061 : first_three_power_sums(GEN pol, GEN *p1, GEN *p2, GEN *p3)
340 : {
341 12061 : long d = degpol(pol);
342 : GEN s1, s2, ms3;
343 :
344 12061 : *p1 = s1 = gneg(RgX_coeff(pol, d-1));
345 :
346 12061 : s2 = RgX_coeff(pol, d-2);
347 12061 : *p2 = gsub(gsqr(s1), gmulsg(2L, s2));
348 :
349 12061 : ms3 = RgX_coeff(pol, d-3);
350 12061 : *p3 = gadd(gmul(s1, gsub(*p2, s2)), gmulsg(-3L, ms3));
351 12061 : }
352 :
353 : /* Let E and a t_POL h of degree 1 or 3 whose roots are 2-torsion points on E.
354 : * - if only_image != 0, return [t, w] used to compute the equation of the
355 : * quotient by the given 2-torsion points
356 : * - else return [t,w, f,g,h], along with the contributions f, g and
357 : * h to the isogeny giving the quotient by h. Variables vx and vy are used
358 : * to create f, g and h, or ignored if only_image is zero */
359 :
360 : /* deg h = 1; 2-torsion contribution from Weierstrass point */
361 : static GEN
362 7819 : contrib_weierstrass_pt(GEN E, GEN h, long only_image, long vx, long vy)
363 : {
364 7819 : GEN p = ellbasechar(E);
365 7819 : GEN a1 = ell_get_a1(E);
366 7819 : GEN a3 = ell_get_a3(E);
367 7819 : GEN x0 = gneg(constant_coeff(h)); /* h = x - x0 */
368 7819 : GEN b = gadd(gmul(a1,x0), a3);
369 : GEN y0, Q, t, w, t1, t2, f, g;
370 :
371 7819 : if (!equalis(p, 2L)) /* char(k) != 2 ==> y0 = -b/2 */
372 7777 : y0 = gmul2n(gneg(b), -1);
373 : else
374 : { /* char(k) = 2 ==> y0 = sqrt(f(x0)) where E is y^2 + h(x) = f(x). */
375 42 : if (!gequal0(b)) pari_err_BUG("two_torsion_contrib (a1*x0+a3 != 0)");
376 42 : y0 = gsqrt(ec_f_evalx(E, x0), 0);
377 : }
378 7819 : Q = mkvec2(x0, y0);
379 7819 : t = ec_dFdx_evalQ(E, Q);
380 7819 : w = gmul(x0, t);
381 7819 : if (only_image) return mkvec2(t,w);
382 :
383 : /* Compute isogeny, f = (x - x0) * t */
384 532 : f = deg1pol_shallow(t, gmul(t, gneg(x0)), vx);
385 :
386 : /* g = (x - x0) * t * (a1 * (x - x0) + (y - y0)) */
387 532 : t1 = deg1pol_shallow(a1, gmul(a1, gneg(x0)), vx);
388 532 : t2 = deg1pol_shallow(gen_1, gneg(y0), vy);
389 532 : g = gmul(f, gadd(t1, t2));
390 532 : return mkvec5(t, w, f, g, h);
391 : }
392 : /* deg h =3; full 2-torsion contribution. NB: assume h is monic; base field
393 : * characteristic is odd or zero (cannot happen in char 2).*/
394 : static GEN
395 14 : contrib_full_tors(GEN E, GEN h, long only_image, long vx, long vy)
396 : {
397 : GEN p1, p2, p3, half_b2, half_b4, t, w, f, g;
398 14 : first_three_power_sums(h, &p1,&p2,&p3);
399 14 : half_b2 = gmul2n(ell_get_b2(E), -1);
400 14 : half_b4 = gmul2n(ell_get_b4(E), -1);
401 :
402 : /* t = 3*(p2 + b4/2) + p1 * b2/2 */
403 14 : t = gadd(gmulsg(3L, gadd(p2, half_b4)), gmul(p1, half_b2));
404 :
405 : /* w = 3 * p3 + p2 * b2/2 + p1 * b4/2 */
406 14 : w = gadd(gmulsg(3L, p3), gadd(gmul(p2, half_b2),
407 : gmul(p1, half_b4)));
408 14 : if (only_image) return mkvec2(t,w);
409 :
410 : /* Compute isogeny */
411 : {
412 7 : GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E), t1, t2;
413 7 : GEN s1 = gneg(RgX_coeff(h, 2));
414 7 : GEN dh = RgX_deriv(h);
415 7 : GEN psi2xy = gadd(deg1pol_shallow(a1, a3, vx),
416 : deg1pol_shallow(gen_2, gen_0, vy));
417 :
418 : /* f = -3 (3 x + b2/2 + s1) h + (3 x^2 + (b2/2) x + (b4/2)) h'*/
419 7 : t1 = RgX_mul(h, gmulsg(-3, deg1pol(stoi(3), gadd(half_b2, s1), vx)));
420 7 : t2 = mkpoln(3, stoi(3), half_b2, half_b4);
421 7 : setvarn(t2, vx);
422 7 : t2 = RgX_mul(dh, t2);
423 7 : f = RgX_add(t1, t2);
424 :
425 : /* 2g = psi2xy * (f'*h - f*h') - (a1*f + a3*h) * h; */
426 7 : t1 = RgX_sub(RgX_mul(RgX_deriv(f), h), RgX_mul(f, dh));
427 7 : t2 = RgX_mul(h, RgX_add(RgX_Rg_mul(f, a1), RgX_Rg_mul(h, a3)));
428 7 : g = RgX_divs(gsub(gmul(psi2xy, t1), t2), 2L);
429 :
430 7 : f = RgX_mul(f, h);
431 7 : g = RgX_mul(g, h);
432 : }
433 7 : return mkvec5(t, w, f, g, h);
434 : }
435 :
436 : /* Given E and a t_POL T whose roots define a subgroup G of E, return the factor
437 : * of T that corresponds to the 2-torsion points E[2] \cap G in G */
438 : INLINE GEN
439 12054 : two_torsion_part(GEN E, GEN T)
440 12054 : { return RgX_gcd(T, elldivpol(E, 2, varn(T))); }
441 :
442 : /* Return the jth Hasse derivative of the polynomial f = \sum_{i=0}^n a_i x^i,
443 : * i.e. \sum_{i=j}^n a_i \binom{i}{j} x^{i-j}. It is a derivation even when the
444 : * coefficient ring has positive characteristic */
445 : static GEN
446 98 : derivhasse(GEN f, ulong j)
447 : {
448 98 : ulong i, d = degpol(f);
449 : GEN df;
450 98 : if (gequal0(f) || d == 0) return pol_0(varn(f));
451 56 : if (j == 0) return gcopy(f);
452 56 : df = cgetg(2 + (d-j+1), t_POL);
453 56 : df[1] = f[1];
454 112 : for (i = j; i <= d; ++i) gel(df, i-j+2) = gmul(binomialuu(i,j), gel(f, i+2));
455 56 : return normalizepol(df);
456 : }
457 :
458 : static GEN
459 812 : non_two_torsion_abscissa(GEN E, GEN h0, long vx)
460 : {
461 : GEN mp1, dh0, ddh0, t, u, t1, t2, t3;
462 812 : long m = degpol(h0);
463 812 : mp1 = gel(h0, m + 1); /* negative of first power sum */
464 812 : dh0 = RgX_deriv(h0);
465 812 : ddh0 = RgX_deriv(dh0);
466 812 : t = ec_bmodel(E, vx);
467 812 : u = ec_half_deriv_2divpol(E, vx);
468 812 : t1 = RgX_sub(RgX_sqr(dh0), RgX_mul(ddh0, h0));
469 812 : t2 = RgX_mul(u, RgX_mul(h0, dh0));
470 812 : t3 = RgX_mul(RgX_sqr(h0),
471 : deg1pol_shallow(stoi(2*m), gmulsg(2L, mp1), vx));
472 : /* t * (dh0^2 - ddh0*h0) - u*dh0*h0 + (2*m*x - 2*s1) * h0^2); */
473 812 : return RgX_add(RgX_sub(RgX_mul(t, t1), t2), t3);
474 : }
475 :
476 : static GEN
477 1316 : isog_abscissa(GEN E, GEN kerp, GEN h0, long vx, GEN two_tors)
478 : {
479 : GEN f0, f2, h2, t1, t2, t3;
480 1316 : f0 = (degpol(h0) > 0)? non_two_torsion_abscissa(E, h0, vx): pol_0(vx);
481 1316 : f2 = gel(two_tors, 3);
482 1316 : h2 = gel(two_tors, 5);
483 :
484 : /* Combine f0 and f2 into the final abscissa of the isogeny. */
485 1316 : t1 = RgX_mul(pol_x(vx), RgX_sqr(kerp));
486 1316 : t2 = RgX_mul(f2, RgX_sqr(h0));
487 1316 : t3 = RgX_mul(f0, RgX_sqr(h2));
488 : /* x * kerp^2 + f2 * h0^2 + f0 * h2^2 */
489 1316 : return RgX_add(t1, RgX_add(t2, t3));
490 : }
491 :
492 : static GEN
493 1267 : non_two_torsion_ordinate_char_not2(GEN E, GEN f, GEN h, GEN psi2)
494 : {
495 1267 : GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E);
496 1267 : GEN df = RgX_deriv(f), dh = RgX_deriv(h);
497 : /* g = df * h * psi2/2 - f * dh * psi2
498 : * - (E.a1 * f + E.a3 * h^2) * h/2 */
499 1267 : GEN t1 = RgX_mul(df, RgX_mul(h, RgX_divs(psi2, 2L)));
500 1267 : GEN t2 = RgX_mul(f, RgX_mul(dh, psi2));
501 1267 : GEN t3 = RgX_mul(RgX_divs(h, 2L),
502 : RgX_add(RgX_Rg_mul(f, a1), RgX_Rg_mul(RgX_sqr(h), a3)));
503 1267 : return RgX_sub(RgX_sub(t1, t2), t3);
504 : }
505 :
506 : /* h = kerq */
507 : static GEN
508 49 : non_two_torsion_ordinate_char2(GEN E, GEN h, GEN x, GEN y)
509 : {
510 49 : GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E), a4 = ell_get_a4(E);
511 49 : GEN b2 = ell_get_b2(E), b4 = ell_get_b4(E), b6 = ell_get_b6(E);
512 : GEN h2, dh, dh2, ddh, D2h, D2dh, H, psi2, u, t, alpha;
513 : GEN p1, t1, t2, t3, t4;
514 49 : long m, vx = varn(x);
515 :
516 49 : h2 = RgX_sqr(h);
517 49 : dh = RgX_deriv(h);
518 49 : dh2 = RgX_sqr(dh);
519 49 : ddh = RgX_deriv(dh);
520 49 : H = RgX_sub(dh2, RgX_mul(h, ddh));
521 49 : D2h = derivhasse(h, 2);
522 49 : D2dh = derivhasse(dh, 2);
523 49 : psi2 = deg1pol_shallow(a1, a3, vx);
524 49 : u = mkpoln(3, b2, gen_0, b6);
525 49 : setvarn(u, vx);
526 49 : t = deg1pol_shallow(b2, b4, vx);
527 49 : alpha = mkpoln(4, a1, a3, gmul(a1, a4), gmul(a3, a4));
528 49 : setvarn(alpha, vx);
529 49 : m = degpol(h);
530 49 : p1 = RgX_coeff(h, m-1); /* first power sum */
531 :
532 49 : t1 = gmul(gadd(gmul(a1, p1), gmulgu(a3, m)), RgX_mul(h,h2));
533 :
534 49 : t2 = gmul(a1, gadd(gmul(a1, gadd(y, psi2)), RgX_add(RgX_Rg_add(RgX_sqr(x), a4), t)));
535 49 : t2 = gmul(t2, gmul(dh, h2));
536 :
537 49 : t3 = gadd(gmul(y, t), RgX_add(alpha, RgX_Rg_mul(u, a1)));
538 49 : t3 = gmul(t3, RgX_mul(h, H));
539 :
540 49 : t4 = gmul(u, psi2);
541 49 : t4 = gmul(t4, RgX_sub(RgX_sub(RgX_mul(h2, D2dh), RgX_mul(dh, H)),
542 : RgX_mul(h, RgX_mul(dh, D2h))));
543 :
544 49 : return gadd(t1, gadd(t2, gadd(t3, t4)));
545 : }
546 :
547 : static GEN
548 1316 : isog_ordinate(GEN E, GEN kerp, GEN kerq, GEN x, GEN y, GEN two_tors, GEN f)
549 : {
550 : GEN g;
551 1316 : if (! equalis(ellbasechar(E), 2L)) {
552 : /* FIXME: We don't use (hence don't need to calculate)
553 : * g2 = gel(two_tors, 4) when char(k) != 2. */
554 1267 : GEN psi2 = ec_dmFdy_evalQ(E, mkvec2(x, y));
555 1267 : g = non_two_torsion_ordinate_char_not2(E, f, kerp, psi2);
556 : } else {
557 49 : GEN h2 = gel(two_tors, 5);
558 49 : GEN g2 = gmul(gel(two_tors, 4), RgX_mul(kerq, RgX_sqr(kerq)));
559 49 : GEN g0 = non_two_torsion_ordinate_char2(E, kerq, x, y);
560 49 : g0 = gmul(g0, RgX_mul(h2, RgX_sqr(h2)));
561 49 : g = gsub(gmul(y, RgX_mul(kerp, RgX_sqr(kerp))), gadd(g2, g0));
562 : }
563 1316 : return g;
564 : }
565 :
566 : /* Given an elliptic curve E and a polynomial kerp whose roots give the
567 : * x-coordinates of a subgroup G of E, return the curve E/G and,
568 : * if only_image is zero, the isogeny pi:E -> E/G. Variables vx and vy are
569 : * used to describe the isogeny (and are ignored if only_image is zero). */
570 : static GEN
571 12054 : isogeny_from_kernel_poly(GEN E, GEN kerp, long only_image, long vx, long vy)
572 : {
573 : long m;
574 12054 : GEN b2 = ell_get_b2(E), b4 = ell_get_b4(E), b6 = ell_get_b6(E);
575 : GEN p1, p2, p3, x, y, f, g, two_tors, EE, t, w;
576 12054 : GEN kerh = two_torsion_part(E, kerp);
577 12054 : GEN kerq = RgX_divrem(kerp, kerh, ONLY_DIVIDES);
578 12054 : if (!kerq) pari_err_BUG("isogeny_from_kernel_poly");
579 : /* isogeny degree: 2*degpol(kerp)+1-degpol(kerh) */
580 12054 : m = degpol(kerq);
581 :
582 12054 : kerp = RgX_normalize(kerp);
583 12054 : kerq = RgX_normalize(kerq);
584 12054 : kerh = RgX_normalize(kerh);
585 12054 : switch(degpol(kerh))
586 : {
587 4214 : case 0:
588 4214 : two_tors = only_image? mkvec2(gen_0, gen_0):
589 777 : mkvec5(gen_0, gen_0, pol_0(vx), pol_0(vx), pol_1(vx));
590 4214 : break;
591 7819 : case 1:
592 7819 : two_tors = contrib_weierstrass_pt(E, kerh, only_image,vx,vy);
593 7819 : break;
594 14 : case 3:
595 14 : two_tors = contrib_full_tors(E, kerh, only_image,vx,vy);
596 14 : break;
597 7 : default:
598 7 : two_tors = NULL;
599 7 : pari_err_DOMAIN("isogeny_from_kernel_poly", "kernel polynomial",
600 : "does not define a subgroup of", E, kerp);
601 : }
602 12047 : first_three_power_sums(kerq,&p1,&p2,&p3);
603 12047 : x = pol_x(vx);
604 12047 : y = pol_x(vy);
605 :
606 : /* t = 6 * p2 + b2 * p1 + m * b4, */
607 12047 : t = gadd(gmulsg(6L, p2), gadd(gmul(b2, p1), gmulsg(m, b4)));
608 :
609 : /* w = 10 * p3 + 2 * b2 * p2 + 3 * b4 * p1 + m * b6, */
610 12047 : w = gadd(gmulsg(10L, p3),
611 : gadd(gmul(gmulsg(2L, b2), p2),
612 : gadd(gmul(gmulsg(3L, b4), p1), gmulsg(m, b6))));
613 :
614 12047 : EE = make_velu_curve(E, gadd(t, gel(two_tors, 1)),
615 12047 : gadd(w, gel(two_tors, 2)));
616 12047 : if (only_image) return EE;
617 :
618 1316 : f = isog_abscissa(E, kerp, kerq, vx, two_tors);
619 1316 : g = isog_ordinate(E, kerp, kerq, x, y, two_tors, f);
620 1316 : return mkvec2(EE, mkvec3(f,g,kerp));
621 : }
622 :
623 : /* Given an elliptic curve E and a subgroup G of E, return the curve
624 : * E/G and, if only_image is zero, the isogeny corresponding
625 : * to the canonical surjection pi:E -> E/G. The variables vx and
626 : * vy are used to describe the isogeny (and are ignored if
627 : * only_image is zero). The subgroup G may be given either as
628 : * a generating point P on E or as a polynomial kerp whose roots are
629 : * the x-coordinates of the points in G */
630 : GEN
631 1120 : ellisogeny(GEN E, GEN G, long only_image, long vx, long vy)
632 : {
633 1120 : pari_sp av = avma;
634 : GEN j, z;
635 1120 : checkell(E);j = ell_get_j(E);
636 1120 : if (vx < 0) vx = 0;
637 1120 : if (vy < 0) vy = 1;
638 1120 : if (varncmp(vx, vy) >= 0)
639 7 : pari_err_PRIORITY("ellisogeny", pol_x(vx), "<=", vy);
640 1113 : if (!only_image && varncmp(vy, gvar(j)) >= 0)
641 7 : pari_err_PRIORITY("ellisogeny", j, ">=", vy);
642 1106 : switch(typ(G))
643 : {
644 441 : case t_VEC:
645 441 : if (!checkellpt_i(G)) pari_err_TYPE("ellisogeny",G);
646 441 : if (!ell_is_inf(G))
647 : {
648 399 : GEN x = gel(G,1), y = gel(G,2);
649 399 : if (!only_image)
650 : {
651 245 : if (varncmp(vy, gvar(x)) >= 0)
652 7 : pari_err_PRIORITY("ellisogeny", x, ">=", vy);
653 238 : if (varncmp(vy, gvar(y)) >= 0)
654 7 : pari_err_PRIORITY("ellisogeny", y, ">=", vy);
655 : }
656 : }
657 427 : z = isogeny_from_kernel_point(E, G, only_image, vx, vy);
658 420 : break;
659 658 : case t_POL:
660 658 : if (!only_image && varncmp(vy, gvar(constant_coeff(G))) >= 0)
661 7 : pari_err_PRIORITY("ellisogeny", constant_coeff(G), ">=", vy);
662 651 : z = isogeny_from_kernel_poly(E, G, only_image, vx, vy);
663 644 : break;
664 7 : default:
665 7 : z = NULL;
666 7 : pari_err_TYPE("ellisogeny", G);
667 : }
668 1064 : return gerepilecopy(av, z);
669 : }
670 :
671 : static GEN
672 11760 : trivial_isogeny(void)
673 : {
674 11760 : return mkvec3(pol_x(0), scalarpol(pol_x(1), 0), pol_1(0));
675 : }
676 :
677 : static GEN
678 273 : isogeny_a4a6(GEN E)
679 : {
680 273 : GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E), b2 = ell_get_b2(E);
681 273 : retmkvec3(deg1pol(gen_1, gdivgu(b2, 12), 0),
682 : deg1pol(gdivgu(a1,2), deg1pol(gen_1, gdivgu(a3,2), 1), 0),
683 : pol_1(0));
684 : }
685 :
686 : static GEN
687 273 : invisogeny_a4a6(GEN E)
688 : {
689 273 : GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E), b2 = ell_get_b2(E);
690 273 : retmkvec3(deg1pol(gen_1, gdivgs(b2, -12), 0),
691 : deg1pol(gdivgs(a1,-2),
692 : deg1pol(gen_1, gadd(gdivgs(a3,-2), gdivgu(gmul(b2,a1), 24)), 1), 0),
693 : pol_1(0));
694 : }
695 :
696 : static GEN
697 9590 : RgXY_eval(GEN P, GEN x, GEN y)
698 : {
699 9590 : return poleval(poleval(P,x), y);
700 : }
701 :
702 : static GEN
703 567 : twistisogeny(GEN iso, GEN d)
704 : {
705 567 : GEN d2 = gsqr(d), d3 = gmul(d, d2);
706 567 : return mkvec3(gdiv(gel(iso,1), d2), gdiv(gel(iso,2), d3), gel(iso, 3));
707 : }
708 :
709 : static GEN
710 10836 : ellisog_by_Kohel(GEN a4, GEN a6, long n, GEN ker, GEN kert, long flag)
711 : {
712 10836 : GEN E = ellinit(mkvec2(a4, a6), NULL, DEFAULTPREC);
713 10836 : GEN F = isogeny_from_kernel_poly(E, ker, flag, 0, 1);
714 10836 : GEN Et = ellinit(flag ? F: gel(F, 1), NULL, DEFAULTPREC);
715 10836 : GEN c4t = ell_get_c4(Et), c6t = ell_get_c6(Et), jt = ell_get_j(Et);
716 10836 : if (!flag)
717 : {
718 567 : GEN Ft = isogeny_from_kernel_poly(Et, kert, flag, 0, 1);
719 567 : GEN isot = twistisogeny(gel(Ft, 2), stoi(n));
720 567 : return mkvec5(c4t, c6t, jt, gel(F, 2), isot);
721 : }
722 10269 : else return mkvec3(c4t, c6t, jt);
723 : }
724 :
725 : static GEN
726 10668 : ellisog_by_roots(GEN a4, GEN a6, long n, GEN z, long flag)
727 : {
728 10668 : GEN k = deg1pol_shallow(gen_1, gneg(z), 0);
729 10668 : GEN kt= deg1pol_shallow(gen_1, gmulsg(n,z), 0);
730 10668 : return ellisog_by_Kohel(a4, a6, n, k, kt, flag);
731 : }
732 :
733 : /* n = 2 or 3 */
734 : static GEN
735 15421 : a4a6_divpol(GEN a4, GEN a6, long n)
736 : {
737 15421 : if (n == 2) return mkpoln(4, gen_1, gen_0, a4, a6);
738 5831 : return mkpoln(5, utoi(3), gen_0, gmulgu(a4,6) , gmulgu(a6,12),
739 : gneg(gsqr(a4)));
740 : }
741 :
742 : static GEN
743 15421 : ellisograph_Kohel_iso(GEN nf, GEN e, long n, GEN z, GEN *pR, long flag)
744 : {
745 : long i, r;
746 15421 : GEN R, V, c4 = gel(e,1), c6 = gel(e,2);
747 15421 : GEN a4 = gdivgs(c4, -48), a6 = gdivgs(c6, -864);
748 15421 : GEN P = a4a6_divpol(a4, a6, n);
749 15421 : R = nfroots(nf, z ? RgX_div_by_X_x(P, z, NULL): P);
750 15421 : if (pR) *pR = R;
751 15421 : r = lg(R); V = cgetg(r, t_VEC);
752 26089 : for (i=1; i < r; i++) gel(V,i) = ellisog_by_roots(a4, a6, n, gel(R,i), flag);
753 15421 : return V;
754 : }
755 :
756 : static GEN
757 14476 : ellisograph_Kohel_r(GEN nf, GEN e, long n, GEN z, long flag)
758 : {
759 14476 : GEN R, iso = ellisograph_Kohel_iso(nf, e, n, z, &R, flag);
760 14476 : long i, r = lg(iso);
761 14476 : GEN V = cgetg(r, t_VEC);
762 24199 : for (i=1; i < r; i++)
763 9723 : gel(V,i) = ellisograph_Kohel_r(nf, gel(iso,i), n, gmulgs(gel(R,i), -n), flag);
764 14476 : return mkvec2(e, V);
765 : }
766 :
767 : static GEN
768 336 : corr(GEN c4, GEN c6)
769 : {
770 336 : GEN c62 = gmul2n(c6, 1);
771 336 : return gadd(gdiv(gsqr(c4), c62), gdiv(c62, gmulgu(c4,3)));
772 : }
773 :
774 : static GEN
775 336 : elkies98(GEN a4, GEN a6, long l, GEN s, GEN a4t, GEN a6t)
776 : {
777 : GEN C, P, S;
778 : long i, n, d;
779 336 : d = l == 2 ? 1 : l>>1;
780 336 : C = cgetg(d+1, t_VEC);
781 336 : gel(C, 1) = gdivgu(gsub(a4, a4t), 5);
782 336 : if (d >= 2)
783 336 : gel(C, 2) = gdivgu(gsub(a6, a6t), 7);
784 336 : if (d >= 3)
785 224 : gel(C, 3) = gdivgu(gsub(gsqr(gel(C, 1)), gmul(a4, gel(C, 1))), 3);
786 2870 : for (n = 3; n < d; ++n)
787 : {
788 2534 : GEN s = gen_0;
789 61390 : for (i = 1; i < n; i++)
790 58856 : s = gadd(s, gmul(gel(C, i), gel(C, n-i)));
791 2534 : gel(C, n+1) = gdivgu(gsub(gsub(gmulsg(3, s), gmul(gmulsg((2*n-1)*(n-1), a4), gel(C, n-1))), gmul(gmulsg((2*n-2)*(n-2), a6), gel(C, n-2))), (n-1)*(2*n+5));
792 : }
793 336 : P = cgetg(d+2, t_VEC);
794 336 : gel(P, 1 + 0) = stoi(d);
795 336 : gel(P, 1 + 1) = s;
796 336 : if (d >= 2)
797 336 : gel(P, 1 + 2) = gdivgu(gsub(gel(C, 1), gmulgu(a4, 2*d)), 6);
798 3094 : for (n = 2; n < d; ++n)
799 2758 : gel(P, 1 + n+1) = gdivgu(gsub(gsub(gel(C, n), gmul(gmulsg(4*n-2, a4), gel(P, 1+n-1))), gmul(gmulsg(4*n-4, a6), gel(P, 1+n-2))), 4*n+2);
800 336 : S = cgetg(d+3, t_POL);
801 336 : S[1] = evalsigne(1) | evalvarn(0);
802 336 : gel(S, 2 + d - 0) = gen_1;
803 336 : gel(S, 2 + d - 1) = gneg(s);
804 3430 : for (n = 2; n <= d; ++n)
805 : {
806 3094 : GEN s = gen_0;
807 68362 : for (i = 1; i <= n; ++i)
808 : {
809 65268 : GEN p = gmul(gel(P, 1+i), gel(S, 2 + d - (n-i)));
810 65268 : s = gadd(s, p);
811 : }
812 3094 : gel(S, 2 + d - n) = gdivgs(s, -n);
813 : }
814 336 : return S;
815 : }
816 :
817 : static GEN
818 2072 : ellisog_by_jt(GEN c4, GEN c6, GEN jt, GEN jtp, GEN s0, long n, long flag)
819 : {
820 2072 : GEN jtp2 = gsqr(jtp), den = gmul(jt, gsubgs(jt, 1728));
821 2072 : GEN c4t = gdiv(jtp2, den);
822 2072 : GEN c6t = gdiv(gmul(jtp, c4t), jt);
823 2072 : if (flag)
824 1904 : return mkvec3(c4t, c6t, jt);
825 : else
826 : {
827 168 : GEN co = corr(c4, c6);
828 168 : GEN cot = corr(c4t, c6t);
829 168 : GEN s = gmul2n(gmulgs(gadd(gadd(s0, co), gmulgs(cot,-n)), -n), -2);
830 168 : GEN a4 = gdivgs(c4, -48), a6 = gdivgs(c6, -864);
831 168 : GEN a4t = gmul(gdivgs(c4t, -48), powuu(n,4)), a6t = gmul(gdivgs(c6t, -864), powuu(n,6));
832 168 : GEN ker = elkies98(a4, a6, n, s, a4t, a6t);
833 168 : GEN st = gmulgs(s, -n);
834 168 : GEN a4tt = gmul(a4,powuu(n,4)), a6tt = gmul(a6,powuu(n,6));
835 168 : GEN kert = elkies98(a4t, a6t, n, st, a4tt, a6tt);
836 168 : return ellisog_by_Kohel(a4, a6, n, ker, kert, flag);
837 : }
838 : }
839 :
840 : /* RENE SCHOOF, Counting points on elliptic curves over finite fields,
841 : * Journal de Theorie des Nombres de Bordeaux, tome 7, no 1 (1995), p. 219-254.
842 : * http://www.numdam.org/item?id=JTNB_1995__7_1_219_0 */
843 : static GEN
844 1918 : ellisog_by_j(GEN e, GEN jt, long n, GEN P, long flag)
845 : {
846 1918 : pari_sp av = avma;
847 1918 : GEN c4 = gel(e,1), c6 = gel(e, 2), j = gel(e, 3);
848 1918 : GEN Px = RgX_deriv(P), Py = RgXY_derivx(P);
849 1918 : GEN Pxj = RgXY_eval(Px, j, jt), Pyj = RgXY_eval(Py, j, jt);
850 1918 : GEN Pxx = RgX_deriv(Px), Pxy = RgX_deriv(Py), Pyy = RgXY_derivx(Py);
851 1918 : GEN Pxxj = RgXY_eval(Pxx,j,jt);
852 1918 : GEN Pxyj = RgXY_eval(Pxy,j,jt);
853 1918 : GEN Pyyj = RgXY_eval(Pyy,j,jt);
854 1918 : GEN c6c4 = gdiv(c6, c4);
855 1918 : GEN jp = gmul(j, c6c4);
856 1918 : GEN jtp = gdivgs(gmul(jp, gdiv(Pxj, Pyj)), -n);
857 1918 : GEN jtpn = gmulgs(jtp, n);
858 1918 : GEN s0 = gdiv(gadd(gadd(gmul(gsqr(jp),Pxxj),gmul(gmul(jp,jtpn),gmul2n(Pxyj,1))),
859 : gmul(gsqr(jtpn),Pyyj)),gmul(jp,Pxj));
860 1918 : GEN et = ellisog_by_jt(c4, c6, jt, jtp, s0, n, flag);
861 1918 : return gerepilecopy(av, et);
862 : }
863 :
864 : static GEN
865 4550 : ellisograph_iso(GEN nf, GEN e, ulong p, GEN P, GEN oj, long flag)
866 : {
867 : long i, r;
868 : GEN Pj, R, V;
869 4550 : if (!P) return ellisograph_Kohel_iso(nf, e, p, oj, NULL, flag);
870 3605 : Pj = poleval(P, gel(e,3));
871 3605 : R = nfroots(nf,oj ? RgX_div_by_X_x(Pj, oj, NULL):Pj);
872 3605 : r = lg(R);
873 3605 : V = cgetg(r, t_VEC);
874 5523 : for (i=1; i < r; i++) gel(V, i) = ellisog_by_j(e, gel(R, i), p, P, flag);
875 3605 : return V;
876 : }
877 :
878 : static GEN
879 3451 : ellisograph_r(GEN nf, GEN e, ulong p, GEN P, GEN oj, long flag)
880 : {
881 3451 : GEN j = gel(e,3), iso = ellisograph_iso(nf, e, p, P, oj, flag);
882 3451 : long i, r = lg(iso);
883 3451 : GEN V = cgetg(r, t_VEC);
884 5215 : for (i=1; i < r; i++) gel(V,i) = ellisograph_r(nf, gel(iso,i), p, P, j, flag);
885 3451 : return mkvec2(e, V);
886 : }
887 :
888 : static GEN
889 4151 : ellisograph_a4a6(GEN E, long flag)
890 : {
891 4151 : GEN c4 = ell_get_c4(E), c6 = ell_get_c6(E), j = ell_get_j(E);
892 4424 : return flag ? mkvec3(c4, c6, j):
893 273 : mkvec5(c4, c6, j, isogeny_a4a6(E), invisogeny_a4a6(E));
894 : }
895 :
896 : static GEN
897 154 : ellisograph_dummy(GEN E, long n, GEN jt, GEN jtt, GEN s0, long flag)
898 : {
899 154 : GEN c4 = ell_get_c4(E), c6 = ell_get_c6(E), c6c4 = gdiv(c6, c4);
900 154 : GEN jtp = gmul(c6c4, gdivgs(gmul(jt, jtt), -n));
901 154 : GEN iso = ellisog_by_jt(c4, c6, jt, jtp, gmul(s0, c6c4), n, flag);
902 154 : GEN v = mkvec2(iso, cgetg(1, t_VEC));
903 154 : return mkvec2(ellisograph_a4a6(E, flag), mkvec(v));
904 : }
905 :
906 : static GEN
907 6440 : isograph_p(GEN nf, GEN e, ulong p, GEN P, long flag)
908 : {
909 6440 : pari_sp av = avma;
910 : GEN iso;
911 6440 : if (P)
912 1687 : iso = ellisograph_r(nf, e, p, P, NULL, flag);
913 : else
914 4753 : iso = ellisograph_Kohel_r(nf, e, p, NULL, flag);
915 6440 : return gerepilecopy(av, iso);
916 : }
917 :
918 : static GEN
919 4088 : get_polmodular(ulong p, long vx, long vy)
920 4088 : { return p > 3 ? polmodular_ZXX(p,0,vx,vy): NULL; }
921 : static GEN
922 3164 : ellisograph_p(GEN nf, GEN E, ulong p, long flag)
923 : {
924 3164 : GEN e = ellisograph_a4a6(E, flag);
925 3164 : GEN P = get_polmodular(p, 0, 1);
926 3164 : return isograph_p(nf, e, p, P, flag);
927 : }
928 :
929 : static long
930 43652 : etree_nbnodes(GEN T)
931 : {
932 43652 : GEN F = gel(T,2);
933 43652 : long n = 1, i, l = lg(F);
934 72492 : for (i = 1; i < l; i++) n += etree_nbnodes(gel(F, i));
935 43652 : return n;
936 : }
937 :
938 : static long
939 16807 : etree_listr(GEN nf, GEN T, GEN V, long n, GEN u, GEN ut)
940 : {
941 16807 : GEN E = gel(T, 1), F = gel(T,2);
942 16807 : long i, l = lg(F);
943 16807 : GEN iso, isot = NULL;
944 16807 : if (lg(E) == 6)
945 : {
946 749 : iso = ellnfcompisog(nf,gel(E,4), u);
947 749 : isot = ellnfcompisog(nf,ut, gel(E,5));
948 749 : gel(V, n) = mkvec5(gel(E,1), gel(E,2), gel(E,3), iso, isot);
949 : } else
950 : {
951 16058 : gel(V, n) = mkvec3(gel(E,1), gel(E,2), gel(E,3));
952 16058 : iso = u;
953 : }
954 27734 : for (i = 1; i < l; i++)
955 10927 : n = etree_listr(nf, gel(F, i), V, n + 1, iso, isot);
956 16807 : return n;
957 : }
958 :
959 : static GEN
960 5880 : etree_list(GEN nf, GEN T)
961 : {
962 5880 : long n = etree_nbnodes(T);
963 5880 : GEN V = cgetg(n+1, t_VEC);
964 5880 : (void) etree_listr(nf, T, V, 1, trivial_isogeny(), trivial_isogeny());
965 5880 : return V;
966 : }
967 :
968 : static long
969 16807 : etree_distmatr(GEN T, GEN M, long n)
970 : {
971 16807 : GEN F = gel(T,2);
972 16807 : long i, j, lF = lg(F), m = n + 1;
973 16807 : GEN V = cgetg(lF, t_VECSMALL);
974 16807 : mael(M, n, n) = 0;
975 27734 : for(i = 1; i < lF; i++)
976 10927 : V[i] = m = etree_distmatr(gel(F,i), M, m);
977 27734 : for(i = 1; i < lF; i++)
978 : {
979 10927 : long mi = i==1 ? n+1: V[i-1];
980 26936 : for(j = mi; j < V[i]; j++)
981 : {
982 16009 : mael(M,n,j) = 1 + mael(M, mi, j);
983 16009 : mael(M,j,n) = 1 + mael(M, j, mi);
984 : }
985 28322 : for(j = 1; j < lF; j++)
986 17395 : if (i != j)
987 : {
988 6468 : long i1, j1, mj = j==1 ? n+1: V[j-1];
989 15211 : for (i1 = mi; i1 < V[i]; i1++)
990 20209 : for(j1 = mj; j1 < V[j]; j1++)
991 11466 : mael(M,i1,j1) = 2 + mael(M,mj,j1) + mael(M,i1,mi);
992 : }
993 : }
994 16807 : return m;
995 : }
996 :
997 : static GEN
998 5880 : etree_distmat(GEN T)
999 : {
1000 5880 : long i, n = etree_nbnodes(T);
1001 5880 : GEN M = cgetg(n+1, t_MAT);
1002 22687 : for(i = 1; i <= n; i++) gel(M,i) = cgetg(n+1, t_VECSMALL);
1003 5880 : (void)etree_distmatr(T, M, 1);
1004 5880 : return M;
1005 : }
1006 :
1007 : static GEN
1008 5880 : distmat_pow(GEN E, ulong p)
1009 : {
1010 5880 : long i, j, l = lg(E);
1011 5880 : GEN M = cgetg(l, t_MAT);
1012 22687 : for(i = 1; i < l; i++)
1013 : {
1014 16807 : gel(M,i) = cgetg(l, t_COL);
1015 77098 : for(j = 1; j < l; j++) gmael(M,i,j) = powuu(p,mael(E,i,j));
1016 : }
1017 5880 : return M;
1018 : }
1019 :
1020 : /* Assume there is a single p-isogeny */
1021 : static GEN
1022 714 : isomatdbl(GEN nf, GEN L, GEN M, ulong p, GEN T2, long flag)
1023 : {
1024 714 : long i, j, n = lg(L) -1;
1025 714 : GEN P = get_polmodular(p, 0, 1);
1026 714 : GEN V = cgetg(2*n+1, t_VEC), N = cgetg(2*n+1, t_MAT);
1027 2527 : for (i=1; i <= n; i++)
1028 : {
1029 1813 : GEN F, E, e = gel(L,i);
1030 1813 : if (i == 1)
1031 714 : F = gmael(T2, 2, 1);
1032 : else
1033 : {
1034 1099 : F = ellisograph_iso(nf, e, p, P, NULL, flag);
1035 1099 : if (lg(F) != 2) pari_err_BUG("isomatdbl");
1036 : }
1037 1813 : E = gel(F, 1);
1038 1813 : if (flag)
1039 1701 : E = mkvec3(gel(E,1), gel(E,2), gel(E,3));
1040 : else
1041 : {
1042 112 : GEN iso = ellnfcompisog(nf, gel(E,4), gel(e, 4));
1043 112 : GEN isot = ellnfcompisog(nf, gel(e,5), gel(E, 5));
1044 112 : E = mkvec5(gel(E,1), gel(E,2), gel(E,3), iso, isot);
1045 : }
1046 1813 : gel(V, i) = e;
1047 1813 : gel(V, i+n) = E;
1048 : }
1049 4340 : for (i=1; i <= 2*n; i++) gel(N, i) = cgetg(2*n+1, t_COL);
1050 2527 : for (i=1; i <= n; i++)
1051 6832 : for (j=1; j <= n; j++)
1052 : {
1053 5019 : gcoeff(N,i,j) = gcoeff(N,i+n,j+n) = gcoeff(M,i,j);
1054 5019 : gcoeff(N,i,j+n) = gcoeff(N,i+n,j) = muliu(gcoeff(M,i,j), p);
1055 : }
1056 714 : return mkvec2(V, N);
1057 : }
1058 :
1059 : static ulong
1060 2604 : ellQ_exceptional_iso(GEN j, GEN *jt, GEN *jtp, GEN *s0)
1061 : {
1062 2604 : *jt = j; *jtp = gen_1;
1063 2604 : if (typ(j)==t_INT)
1064 : {
1065 392 : long js = itos_or_0(j);
1066 : GEN j37;
1067 392 : if (js==-32768) { *s0 = mkfracss(-1156,539); return 11; }
1068 378 : if (js==-121)
1069 14 : { *jt = stoi(-24729001) ; *jtp = mkfracss(4973,5633);
1070 14 : *s0 = mkfracss(-1961682050,1204555087); return 11;}
1071 364 : if (js==-24729001)
1072 14 : { *jt = stoi(-121); *jtp = mkfracss(5633,4973);
1073 14 : *s0 = mkfracss(-1961682050,1063421347); return 11;}
1074 350 : if (js==-884736)
1075 14 : { *s0 = mkfracss(-1100,513); return 19; }
1076 336 : j37 = negi(uu32toi(37876312,1780746325));
1077 336 : if (js==-9317)
1078 : {
1079 14 : *jt = j37;
1080 14 : *jtp = mkfracss(1984136099,496260169);
1081 14 : *s0 = mkfrac(negi(uu32toi(457100760,4180820796UL)),
1082 : uu32toi(89049913, 4077411069UL));
1083 14 : return 37;
1084 : }
1085 322 : if (equalii(j, j37))
1086 : {
1087 14 : *jt = stoi(-9317);
1088 14 : *jtp = mkfrac(utoi(496260169),utoi(1984136099UL));
1089 14 : *s0 = mkfrac(negi(uu32toi(41554614,2722784052UL)),
1090 : uu32toi(32367030,2614994557UL));
1091 14 : return 37;
1092 : }
1093 308 : if (js==-884736000)
1094 14 : { *s0 = mkfracss(-1073708,512001); return 43; }
1095 294 : if (equalii(j, negi(powuu(5280,3))))
1096 14 : { *s0 = mkfracss(-176993228,85184001); return 67; }
1097 280 : if (equalii(j, negi(powuu(640320,3))))
1098 14 : { *s0 = mkfrac(negi(uu32toi(72512,1969695276)), uu32toi(35374,1199927297));
1099 14 : return 163; }
1100 : } else
1101 : {
1102 2212 : GEN j1 = mkfracss(-297756989,2);
1103 2212 : GEN j2 = mkfracss(-882216989,131072);
1104 2212 : if (gequal(j, j1))
1105 : {
1106 14 : *jt = j2; *jtp = mkfracss(1503991,2878441);
1107 14 : *s0 = mkfrac(negi(uu32toi(121934,548114672)),uu32toi(77014,117338383));
1108 14 : return 17;
1109 : }
1110 2198 : if (gequal(j, j2))
1111 : {
1112 14 : *jt = j1; *jtp = mkfracss(2878441,1503991);
1113 14 : *s0 = mkfrac(negi(uu32toi(121934,548114672)),uu32toi(40239,4202639633UL));
1114 14 : return 17;
1115 : }
1116 : }
1117 2450 : return 0;
1118 : }
1119 :
1120 : static GEN
1121 5880 : nfmkisomat(GEN nf, ulong p, GEN T)
1122 5880 : { return mkvec2(etree_list(nf,T), distmat_pow(etree_distmat(T),p)); }
1123 : static GEN
1124 2492 : mkisomat(ulong p, GEN T)
1125 2492 : { return nfmkisomat(NULL, p, T); }
1126 : static GEN
1127 714 : mkisomatdbl(ulong p, GEN T, ulong p2, GEN T2, long flag)
1128 : {
1129 714 : GEN v = mkisomat(p,T);
1130 714 : return isomatdbl(NULL, gel(v,1), gel(v,2), p2, T2, flag);
1131 : }
1132 :
1133 : /*See M.A Kenku, On the number of Q-isomorphism classes of elliptic curves in
1134 : * each Q-isogeny class, Journal of Number Theory Volume 15, Issue 2,
1135 : * October 1982, pp 199-202,
1136 : * http://www.sciencedirect.com/science/article/pii/0022314X82900257 */
1137 : enum { _2 = 1, _3 = 2, _5 = 4, _7 = 8, _13 = 16 };
1138 : static ulong
1139 2450 : ellQ_goodl(GEN E)
1140 : {
1141 : forprime_t T;
1142 2450 : long i, CM = ellQ_get_CM(E);
1143 2450 : ulong mask = 31;
1144 2450 : GEN disc = ell_get_disc(E);
1145 2450 : pari_sp av = avma;
1146 2450 : u_forprime_init(&T, 17UL,ULONG_MAX);
1147 50008 : for(i=1; mask && i<=20; i++)
1148 : {
1149 47558 : ulong p = u_forprime_next(&T);
1150 47558 : if (umodiu(disc,p)==0) i--;
1151 : else
1152 : {
1153 47180 : long t = ellap_CM_fast(E, p, CM), D = t*t - 4*p;
1154 47180 : if (t%2) mask &= ~_2;
1155 47180 : if ((mask & _3) && kross(D,3)==-1) mask &= ~_3;
1156 47180 : if ((mask & _5) && kross(D,5)==-1) mask &= ~_5;
1157 47180 : if ((mask & _7) && kross(D,7)==-1) mask &= ~_7;
1158 47180 : if ((mask &_13) && kross(D,13)==-1) mask &= ~_13;
1159 : }
1160 : }
1161 2450 : return gc_ulong(av, mask);
1162 : }
1163 :
1164 : static long
1165 182 : ellQ_goodl_l(GEN E, long l)
1166 : {
1167 : forprime_t T;
1168 : long i;
1169 182 : GEN disc = ell_get_disc(E);
1170 : pari_sp av;
1171 182 : u_forprime_init(&T, 17UL, ULONG_MAX); av = avma;
1172 2408 : for (i = 1; i <= 20; i++, set_avma(av))
1173 : {
1174 2303 : ulong p = u_forprime_next(&T);
1175 : long t;
1176 2303 : if (umodiu(disc,p)==0) { i--; continue; }
1177 2254 : t = itos(ellap(E, utoipos(p)));
1178 2254 : if (l==2)
1179 : {
1180 462 : if (odd(t)) return 0;
1181 : }
1182 : else
1183 : {
1184 1792 : long D = t*t - 4*p;
1185 1792 : if (kross(D,l) == -1) return 0;
1186 : }
1187 : }
1188 105 : return 1;
1189 : }
1190 :
1191 : static GEN
1192 112 : ellnf_goodl_l(GEN E, GEN v)
1193 : {
1194 : forprime_t T;
1195 112 : long i, lv = lg(v);
1196 112 : GEN nf = ellnf_get_nf(E), disc = ell_get_disc(E), w = const_F2v(lv-1);
1197 : pari_sp av;
1198 112 : u_forprime_init(&T, 17UL,ULONG_MAX); av = avma;
1199 2443 : for (i = 1; i <= 20; i++, set_avma(av))
1200 : {
1201 2331 : ulong p = u_forprime_next(&T);
1202 2331 : GEN pr = idealprimedec(nf, utoipos(p));
1203 2331 : long t, j, k, g = lg(pr)-1;
1204 6559 : for (j = 1; j <= g; j++)
1205 : {
1206 4228 : GEN prj = gel(pr, j);
1207 4228 : if (nfval(nf, disc, prj)) {i--; continue;}
1208 4137 : t = itos(ellap(E, prj));
1209 22288 : for(k = 1; k < lv; k++)
1210 : {
1211 18151 : GEN l = gel(v,k);
1212 18151 : if (equaliu(l,2))
1213 : {
1214 4137 : if (odd(t)) F2v_clear(w, k);
1215 : }
1216 : else
1217 : {
1218 14014 : GEN D = subii(sqrs(t), shifti(pr_norm(prj),2));
1219 14014 : if (kronecker(D,l)==-1) F2v_clear(w, k);
1220 : }
1221 : }
1222 : }
1223 : }
1224 112 : return w;
1225 : }
1226 :
1227 : static GEN
1228 6188 : ellnf_charpoly(GEN E, GEN pr)
1229 6188 : { return deg2pol_shallow(gen_1, negi(ellap(E,pr)), pr_norm(pr), 0); }
1230 :
1231 : static GEN
1232 18634 : starlaw(GEN p, GEN q)
1233 : {
1234 18634 : GEN Q = RgX_homogenize(RgX_recip(q), 1);
1235 18634 : return ZX_ZXY_resultant(p, Q);
1236 : }
1237 :
1238 : static GEN
1239 8414 : startor(GEN p, long r)
1240 : {
1241 8414 : GEN xr = pol_xn(r, 0), psir = gsub(xr, gen_1);
1242 8414 : return gsubstpol(starlaw(p, psir),xr,pol_x(0));
1243 : }
1244 :
1245 : static GEN
1246 2240 : stariter_easy(GEN E, GEN p)
1247 : {
1248 2240 : GEN nf = ellnf_get_nf(E), dec = idealprimedec(nf, p);
1249 2240 : long d = nf_get_degree(nf), l = lg(dec), i, k;
1250 2240 : GEN R, starl = deg1pol_shallow(gen_1, gen_m1, 0);
1251 6202 : for (i=1; i < l; i++)
1252 : {
1253 3962 : GEN pr = gel(dec,i), q = ellnf_charpoly(E, pr);
1254 3962 : starl = starlaw(starl, startor(q, 12*pr_get_e(pr)));
1255 : }
1256 7420 : for (k = 0, R = p; 2*k <= d; k++) R = mulii(R, poleval(starl, powiu(p,12*k)));
1257 2240 : return R;
1258 : }
1259 :
1260 : /* pr^h assumed principal */
1261 : static GEN
1262 2226 : idealgen_minpoly(GEN bnf, GEN pr, GEN h)
1263 : {
1264 2226 : GEN e = isprincipalfact(bnf, NULL, mkvec(pr), mkvec(h), nf_GEN);
1265 2226 : return minpoly(basistoalg(bnf, gel(e,2)), 0);
1266 : }
1267 :
1268 : static GEN
1269 6468 : stariter(GEN p, long r)
1270 : {
1271 6468 : GEN starl = deg1pol_shallow(gen_1, gen_m1, 0);
1272 : long i;
1273 12726 : for (i = 1; i <= r; i++) starl = starlaw(starl, p);
1274 6468 : return starl;
1275 : }
1276 :
1277 : static GEN
1278 2226 : stariter_hard(GEN E, GEN bnf, GEN pr)
1279 : {
1280 2226 : GEN nf = bnf_get_nf(bnf);
1281 2226 : long k, d = nf_get_degree(nf), d2 = d>>1;
1282 2226 : GEN h = cyc_get_expo(bnf_get_cyc(bnf)); /* could use order of pr in Cl_K */
1283 2226 : GEN pol = idealgen_minpoly(bnf, pr, h);
1284 2226 : GEN S = startor(pol,12), P = gen_1, polchar;
1285 8694 : for (k = 0; k <= d2; k++) P = gmul(P, stariter(S,k));
1286 2226 : polchar = ellnf_charpoly(E, pr);
1287 2226 : return ZX_resultant(startor(polchar,12*itou(h)), P);
1288 : }
1289 :
1290 : /* Based on a GP script by Nicolas Billerey and Theorems 2.4 and 2.8 of
1291 : * N. Billerey, Criteres d'irreductibilite pour les representations des
1292 : * courbes elliptiques, Int. J. Number Theory 7 (2011), no. 4, 1001-1032. */
1293 : static GEN
1294 49 : ellnf_prime_degree_hard(GEN E, GEN bad)
1295 : {
1296 49 : GEN nf = ellnf_get_nf(E), bnf = ellnf_get_bnf(E), R = gen_0;
1297 : forprime_t T;
1298 : long i;
1299 49 : if (!bnf) bnf = bnfinit0(nf,1,NULL,DEFAULTPREC);
1300 49 : u_forprime_init(&T, 5UL, ULONG_MAX);
1301 1078 : for (i = 1; i <= 20; i++)
1302 : {
1303 : long j, lpr;
1304 : GEN pri;
1305 1029 : ulong p = u_forprime_next(&T);
1306 1029 : if (dvdiu(bad, p)) { i--; continue; }
1307 980 : pri = idealprimedec(nf,utoi(p)); lpr = lg(pri);
1308 3206 : for(j = 1; j < lpr; j++)
1309 : {
1310 2226 : GEN R0 = stariter_hard(E,bnf,gel(pri,j)), r;
1311 2226 : R = gcdii(R, mului(p, R0));
1312 2226 : if (Z_issquareall(R, &r)) R = r;
1313 : }
1314 : }
1315 49 : return R;
1316 : }
1317 :
1318 : /* Based on a GP script by Nicolas Billerey and Theorems 2.4 and 2.8 of
1319 : * N. Billerey, Criteres d'irreductibilite pour les representations des
1320 : * courbes elliptiques, Int. J. Number Theory 7 (2011), no. 4, 1001-1032. */
1321 : static GEN
1322 112 : ellnf_prime_degree_easy(GEN E, GEN bad)
1323 : {
1324 : forprime_t T;
1325 : long i;
1326 112 : GEN B = gen_0;
1327 112 : u_forprime_init(&T, 5UL,ULONG_MAX);
1328 2576 : for (i = 1; i <= 20; i++)
1329 : {
1330 2464 : ulong p = u_forprime_next(&T);
1331 : GEN b;
1332 2464 : if (dvdiu(bad, p)) { i--; continue; }
1333 2240 : B = gcdii(B, stariter_easy(E, utoipos(p)));
1334 2240 : if (Z_issquareall(B, &b)) B = b;
1335 : }
1336 112 : return B;
1337 : }
1338 :
1339 : static GEN
1340 112 : ellnf_prime_degree(GEN E)
1341 : {
1342 112 : GEN nf = ellnf_get_nf(E), nf_bad = nf_get_ramified_primes(nf);
1343 112 : GEN g = ellglobalred(E);
1344 112 : GEN N = idealnorm(nf,gel(g,1)), Nfa = prV_primes(gmael(g,4,1));
1345 112 : GEN bad = mulii(N, nf_get_disc(nf)), P;
1346 112 : GEN B = ellnf_prime_degree_easy(E, bad);
1347 112 : if (!signe(B))
1348 : {
1349 49 : B = ellnf_prime_degree_hard(E, bad);
1350 49 : if (!signe(B))
1351 : {
1352 7 : long D = elliscm(E);
1353 7 : if (!D || isintzero(nfisincl(quadpoly(stoi(D)), ellnf_get_nf(E))))
1354 0 : pari_err_IMPL("ellisomat, very hard case");
1355 7 : return stoi(D);
1356 : }
1357 : }
1358 105 : if (!signe(B)) return NULL;
1359 105 : B = muliu(absi(B), 6);
1360 105 : P = ZV_union_shallow(ZV_union_shallow(Nfa,nf_bad), gel(Z_factor(B),1));
1361 105 : return vec_to_vecsmall(RgV_F2v_extract_shallow(P, ellnf_goodl_l(E, P)));
1362 : }
1363 :
1364 : static GEN
1365 2604 : ellQ_isomat(GEN E, long flag)
1366 : {
1367 2604 : GEN K = NULL, T2 = NULL, T3 = NULL, T5, T7, T13;
1368 : ulong good;
1369 : long n2, n3, n5, n7, n13;
1370 2604 : GEN jt, jtp, s0, j = ell_get_j(E);
1371 2604 : long l = ellQ_exceptional_iso(j, &jt, &jtp, &s0);
1372 2604 : if (l)
1373 : {
1374 : #if 1
1375 154 : return mkisomat(l, ellisograph_dummy(E, l, jt, jtp, s0, flag));
1376 : #else
1377 : return mkisomat(l, ellisograph_p(K, E, l), flag);
1378 : #endif
1379 : }
1380 2450 : good = ellQ_goodl(ellintegralmodel(E,NULL));
1381 2450 : if (good & _2)
1382 : {
1383 1834 : T2 = ellisograph_p(K, E, 2, flag);
1384 1834 : n2 = etree_nbnodes(T2);
1385 1834 : if (n2>4 || gequalgs(j, 1728) || gequalgs(j, 287496))
1386 539 : return mkisomat(2, T2);
1387 616 : } else n2 = 1;
1388 1911 : if (good & _3)
1389 : {
1390 924 : T3 = ellisograph_p(K, E, 3, flag);
1391 924 : n3 = etree_nbnodes(T3);
1392 924 : if (n3>1 && n2==2) return mkisomatdbl(3,T3,2,T2, flag);
1393 483 : if (n3==2 && n2>1) return mkisomatdbl(2,T2,3,T3, flag);
1394 364 : if (n3>2 || gequal0(j)) return mkisomat(3, T3);
1395 987 : } else n3 = 1;
1396 1099 : if (good & _5)
1397 : {
1398 224 : T5 = ellisograph_p(K, E, 5, flag);
1399 224 : n5 = etree_nbnodes(T5);
1400 224 : if (n5>1 && n2>1) return mkisomatdbl(2,T2,5,T5, flag);
1401 196 : if (n5>1 && n3>1) return mkisomatdbl(3,T3,5,T5, flag);
1402 126 : if (n5>1) return mkisomat(5, T5);
1403 875 : } else n5 = 1;
1404 875 : if (good & _7)
1405 : {
1406 70 : T7 = ellisograph_p(K, E, 7, flag);
1407 70 : n7 = etree_nbnodes(T7);
1408 70 : if (n7>1 && n2>1) return mkisomatdbl(2,T2,7,T7, flag);
1409 14 : if (n7>1 && n3>1) return mkisomatdbl(3,T3,7,T7, flag);
1410 14 : if (n7>1) return mkisomat(7,T7);
1411 805 : } else n7 = 1;
1412 805 : if (n2>1) return mkisomat(2,T2);
1413 154 : if (n3>1) return mkisomat(3,T3);
1414 112 : if (good & _13)
1415 : {
1416 0 : T13 = ellisograph_p(K, E, 13, flag);
1417 0 : n13 = etree_nbnodes(T13);
1418 0 : if (n13>1) return mkisomat(13,T13);
1419 112 : } else n13 = 1;
1420 112 : return mkvec2(mkvec(ellisograph_a4a6(E,flag)), matid(1));
1421 : }
1422 :
1423 : static long
1424 3276 : fill_LM(GEN LM, GEN L, GEN M, GEN z, long k)
1425 : {
1426 3276 : GEN Li = gel(LM,1), Mi1 = gmael(LM,2,1);
1427 3276 : long j, m = lg(Li);
1428 7616 : for (j = 2; j < m; j++)
1429 : {
1430 4340 : GEN d = gel(Mi1,j);
1431 4340 : gel(L, k) = gel(Li,j);
1432 4340 : gel(M, k) = z? mulii(d,z): d;
1433 4340 : k++;
1434 : }
1435 3276 : return k;
1436 : }
1437 : static GEN
1438 644 : ellnf_isocrv(GEN nf, GEN E, GEN v, GEN PE, long flag)
1439 : {
1440 : long i, l, lv, n, k;
1441 644 : GEN L, M, LE = cgetg_copy(v,&lv), e = ellisograph_a4a6(E, flag);
1442 2072 : for (i = n = 1; i < lv; i++)
1443 : {
1444 1428 : ulong p = uel(v,i);
1445 1428 : GEN T = isograph_p(nf, e, p, gel(PE,i), flag);
1446 1428 : GEN LM = nfmkisomat(nf, p, T);
1447 1428 : gel(LE,i) = LM;
1448 1428 : n *= lg(gel(LM,1)) - 1;
1449 : }
1450 644 : L = cgetg(n+1,t_VEC); gel(L,1) = e;
1451 644 : M = cgetg(n+1,t_COL); gel(M,1) = gen_1;
1452 2072 : for (i = 1, k = 2; i < lv; i++)
1453 : {
1454 1428 : ulong p = uel(v,i);
1455 1428 : GEN P = gel(PE,i);
1456 1428 : long kk = k;
1457 1428 : k = fill_LM(gel(LE,i), L, M, NULL, k);
1458 3276 : for (l = 2; l < kk; l++)
1459 : {
1460 1848 : GEN T = isograph_p(nf, gel(L,l), p, P, flag);
1461 1848 : GEN LMe = nfmkisomat(nf, p, T);
1462 1848 : k = fill_LM(LMe, L, M, gel(M,l), k);
1463 : }
1464 : }
1465 644 : return mkvec2(L, M);
1466 : }
1467 :
1468 : static long
1469 3654 : nfispower_quo(GEN nf, long d, GEN a, GEN b)
1470 : {
1471 3654 : if (gequal(a,b)) return 1;
1472 3248 : return nfispower(nf, nfdiv(nf, a, b), d, NULL);
1473 : }
1474 :
1475 : static long
1476 22134 : isomat_eq(GEN nf, GEN e1, GEN e2)
1477 : {
1478 22134 : if (gequal(e1,e2)) return 1;
1479 21126 : if (!gequal(gel(e1,3), gel(e2,3))) return 0;
1480 3654 : if (gequal0(gel(e1,3)))
1481 0 : return nfispower_quo(nf,6,gel(e1,2),gel(e2,2));
1482 3654 : if (gequalgs(gel(e1,3),1728))
1483 0 : return nfispower_quo(nf,4,gel(e1,1),gel(e2,1));
1484 3654 : return nfispower_quo(nf,2,gmul(gel(e1,1),gel(e2,2)),gmul(gel(e1,2),gel(e2,1)));
1485 : }
1486 :
1487 : static long
1488 4340 : isomat_find(GEN nf, GEN e, GEN L)
1489 : {
1490 4340 : long i, l = lg(L);
1491 22134 : for (i=1; i<l; i++)
1492 22134 : if (isomat_eq(nf, e, gel(L,i))) return i;
1493 : pari_err_BUG("isomat_find"); return 0; /* LCOV_EXCL_LINE */
1494 : }
1495 :
1496 : static GEN
1497 539 : isomat_perm(GEN nf, GEN E, GEN L)
1498 : {
1499 539 : long i, l = lg(E);
1500 539 : GEN v = cgetg(l, t_VECSMALL);
1501 4879 : for (i=1; i<l; i++)
1502 4340 : uel(v, i) = isomat_find(nf, gel(E,i), L);
1503 539 : return v;
1504 : }
1505 :
1506 : static GEN
1507 105 : ellnf_modpoly(GEN v, long vx, long vy)
1508 : {
1509 105 : long i, l = lg(v);
1510 105 : GEN P = cgetg(l, t_VEC);
1511 315 : for(i = 1; i < l; i++) gel(P, i) = get_polmodular(v[i],vx,vy);
1512 105 : return P;
1513 : }
1514 :
1515 : static GEN
1516 112 : ellnf_isomat(GEN E, long flag)
1517 : {
1518 112 : GEN nf = ellnf_get_nf(E);
1519 112 : GEN v = ellnf_prime_degree(E);
1520 112 : if (typ(v)!=t_INT)
1521 : {
1522 105 : long vy = fetch_var_higher();
1523 105 : long vx = fetch_var_higher();
1524 105 : GEN P = ellnf_modpoly(v,vx,vy);
1525 105 : GEN LM = ellnf_isocrv(nf, E, v, P, flag), L = gel(LM,1), M = gel(LM,2);
1526 105 : long i, l = lg(L);
1527 105 : GEN R = cgetg(l, t_MAT);
1528 105 : gel(R,1) = M;
1529 644 : for(i = 2; i < l; i++)
1530 : {
1531 539 : GEN Li = gel(L,i);
1532 539 : GEN e = mkvec2(gdivgs(gel(Li,1), -48), gdivgs(gel(Li,2), -864));
1533 539 : GEN LMi = ellnf_isocrv(nf, ellinit(e, nf, DEFAULTPREC), v, P, 1);
1534 539 : GEN LLi = gel(LMi, 1), Mi = gel(LMi, 2);
1535 539 : GEN r = isomat_perm(nf, L, LLi);
1536 539 : gel(R,i) = vecpermute(Mi, r);
1537 : }
1538 105 : delete_var(); delete_var();
1539 105 : return mkvec2(L, R);
1540 : } else
1541 7 : return v;
1542 : }
1543 :
1544 : static GEN
1545 2821 : list_to_crv(GEN L)
1546 : {
1547 : long i, l;
1548 2821 : GEN V = cgetg_copy(L, &l);
1549 14581 : for (i=1; i < l; i++)
1550 : {
1551 11760 : GEN Li = gel(L,i);
1552 11760 : GEN e = mkvec2(gdivgs(gel(Li,1), -48), gdivgs(gel(Li,2), -864));
1553 11760 : gel(V,i) = lg(Li)==6 ? mkvec3(e, gel(Li,4), gel(Li,5)): e;
1554 : }
1555 2821 : return V;
1556 : }
1557 :
1558 : GEN
1559 2919 : ellisomat(GEN E, long p, long flag)
1560 : {
1561 2919 : pari_sp av = avma;
1562 2919 : GEN r = NULL, nf = NULL;
1563 2919 : int good = 1;
1564 2919 : if (flag < 0 || flag > 1) pari_err_FLAG("ellisomat");
1565 2919 : if (p < 0) pari_err_PRIME("ellisomat", stoi(p));
1566 2919 : if (p == 1) { flag = 1; p = 0; } /* for backward compatibility */
1567 2919 : checkell(E);
1568 2919 : switch(ell_get_type(E))
1569 : {
1570 2786 : case t_ELL_Q:
1571 2786 : if (p) good = ellQ_goodl_l(E, p);
1572 2786 : break;
1573 133 : case t_ELL_NF:
1574 133 : if (p) good = F2v_coeff(ellnf_goodl_l(E, mkvecs(p)), 1);
1575 133 : nf = ellnf_get_nf(E);
1576 133 : if (nf_get_varn(nf) <= 1 - flag)
1577 14 : pari_err_PRIORITY("ellisomat", nf_get_pol(nf), "<=", 1 - flag);
1578 119 : break;
1579 0 : default: pari_err_TYPE("ellisomat",E);
1580 : }
1581 2905 : if (!good) r = mkvec2(mkvec(ellisograph_a4a6(E, flag)), matid(1));
1582 : else
1583 : {
1584 2828 : if (p)
1585 112 : r = nfmkisomat(nf, p, ellisograph_p(nf, E, p, flag));
1586 : else
1587 2716 : r = nf? ellnf_isomat(E, flag): ellQ_isomat(E, flag);
1588 2828 : if (typ(r) == t_VEC) gel(r,1) = list_to_crv(gel(r,1));
1589 : }
1590 2905 : return gerepilecopy(av, r);
1591 : }
1592 :
1593 : static GEN
1594 4361 : get_isomat(GEN v)
1595 : {
1596 : GEN M, vE, wE;
1597 : long i, l;
1598 4361 : if (typ(v) != t_VEC) return NULL;
1599 4361 : if (checkell_i(v))
1600 : {
1601 35 : if (ell_get_type(v) != t_ELL_Q) return NULL;
1602 35 : v = ellisomat(v,0,1);
1603 35 : wE = gel(v,1); l = lg(wE);
1604 35 : M = gel(v,2);
1605 : }
1606 : else
1607 : {
1608 4326 : if (lg(v) != 3) return NULL;
1609 4326 : vE = gel(v,1); l = lg(vE);
1610 4326 : M = gel(v,2);
1611 4326 : if (typ(M) != t_MAT || !RgM_is_ZM(M)) return NULL;
1612 4326 : if (typ(vE) != t_VEC || l == 1) return NULL;
1613 4326 : if (lg(gel(vE,1)) == 3) wE = shallowcopy(vE);
1614 : else
1615 : { /* [[a4,a6],f,g] */
1616 14 : wE = cgetg_copy(vE,&l);
1617 70 : for (i = 1; i < l; i++) gel(wE,i) = gel(gel(vE,i),1);
1618 : }
1619 : }
1620 : /* wE a vector of [a4,a6] */
1621 23044 : for (i = 1; i < l; i++)
1622 : {
1623 18683 : GEN e = ellinit(gel(wE,i), gen_1, DEFAULTPREC);
1624 18683 : GEN E = ellminimalmodel(e, NULL);
1625 18683 : obj_free(e); gel(wE,i) = E;
1626 : }
1627 4361 : return mkvec2(wE, M);
1628 : }
1629 :
1630 : GEN
1631 2184 : ellweilcurve(GEN E, GEN *ms)
1632 : {
1633 2184 : pari_sp av = avma;
1634 2184 : GEN vE = get_isomat(E), vL, Wx, W, XPM, Lf, Cf;
1635 : long i, l;
1636 :
1637 2184 : if (!vE) pari_err_TYPE("ellweilcurve",E);
1638 2184 : vE = gel(vE,1); l = lg(vE);
1639 2184 : Wx = msfromell(vE, 0);
1640 2184 : W = gel(Wx,1);
1641 2184 : XPM = gel(Wx,2);
1642 : /* lattice attached to the Weil curve in the isogeny class */
1643 2184 : Lf = mslattice(W, gmael(XPM,1,3));
1644 2184 : Cf = ginv(Lf); /* left-inverse */
1645 2184 : vL = cgetg(l, t_VEC);
1646 11564 : for (i=1; i < l; i++)
1647 : {
1648 9380 : GEN c, Ce, Le = gmael(XPM,i,3);
1649 9380 : Ce = Q_primitive_part(RgM_mul(Cf, Le), &c);
1650 9380 : Ce = ZM_snf(Ce);
1651 9380 : if (c) { Ce = ZC_Q_mul(Ce,c); settyp(Ce,t_VEC); }
1652 9380 : gel(vL,i) = Ce;
1653 : }
1654 11564 : for (i = 1; i < l; i++) obj_free(gel(vE,i));
1655 2184 : vE = mkvec2(vE, vL);
1656 2184 : if (!ms) return gerepilecopy(av, vE);
1657 7 : *ms = Wx; return gc_all(av, 2, &vE, ms);
1658 : }
1659 :
1660 : GEN
1661 2177 : ellisotree(GEN E)
1662 : {
1663 2177 : pari_sp av = avma;
1664 2177 : GEN L = get_isomat(E), vE, adj, M;
1665 : long i, j, n;
1666 2177 : if (!L) pari_err_TYPE("ellisotree",E);
1667 2177 : vE = gel(L,1);
1668 2177 : adj = gel(L,2);
1669 2177 : n = lg(vE)-1; L = cgetg(n+1, t_VEC);
1670 11480 : for (i = 1; i <= n; i++) gel(L,i) = ellR_area(gel(vE,i), DEFAULTPREC);
1671 2177 : M = zeromatcopy(n,n);
1672 11480 : for (i = 1; i <= n; i++)
1673 28651 : for (j = i+1; j <= n; j++)
1674 : {
1675 19348 : GEN p = gcoeff(adj,i,j);
1676 19348 : if (!isprime(p)) continue;
1677 : /* L[i] / L[j] = p or 1/p; p iff E[i].lattice \subset E[j].lattice */
1678 8022 : if (gcmp(gel(L,i), gel(L,j)) > 0)
1679 4634 : gcoeff(M,i,j) = p;
1680 : else
1681 3388 : gcoeff(M,j,i) = p;
1682 : }
1683 11480 : for (i = 1; i <= n; i++) obj_free(gel(vE,i));
1684 2177 : return gerepilecopy(av, mkvec2(vE,M));
1685 : }
1686 :
1687 : /* Shortest path between x and y on the adjencc graph of A-A~
1688 : * Dijkstra algorithm. */
1689 : static GEN
1690 2142 : shortest_path(GEN A, long x, long y)
1691 : {
1692 : GEN C, v, w;
1693 2142 : long n = lg(A)-1, i, l, z, t, m;
1694 2142 : v = const_vecsmall(n, n); v[x] = 0;
1695 2142 : w = const_vecsmall(n, 0);
1696 9170 : for (l = 1; l < n; l++)
1697 45234 : for (z = 1; z <= n; z++)
1698 : {
1699 38206 : if (v[z] != l-1) continue;
1700 55202 : for (t = 1; t <= n; t++)
1701 46501 : if (!gequal(gcoeff(A,z,t),gcoeff(A,t,z)) && v[t]>l)
1702 7028 : { v[t]=l; w[t]=z; }
1703 : }
1704 2142 : m = v[y]+1; if (m > n) return NULL;
1705 2142 : C = cgetg(m+1, t_VECSMALL);
1706 6503 : for (i = 1; i <= m; i++)
1707 : {
1708 4361 : C[m+1-i] = y;
1709 4361 : y = w[y];
1710 : }
1711 2142 : return C;
1712 : }
1713 :
1714 : static GEN
1715 2142 : path_to_manin(GEN A, long i, long k)
1716 : {
1717 2142 : GEN C = shortest_path(A, i, k);
1718 2142 : long j, lC = lg(C);
1719 2142 : GEN c = gen_1;
1720 4361 : for (j = 1; j < lC-1; j++)
1721 : {
1722 2219 : GEN t = gsub(gcoeff(A,C[j], C[j+1]), gcoeff(A,C[j+1], C[j]));
1723 2219 : if (signe(t) < 0) c = gmul(c, gneg(t));
1724 : }
1725 2142 : return c;
1726 : }
1727 :
1728 : GEN
1729 2142 : ellmaninconstant(GEN E)
1730 : {
1731 2142 : pari_sp av = avma;
1732 : GEN M, A, vS, L;
1733 2142 : long i, k, lvS, is_ell = checkell_i(E);
1734 2142 : L = is_ell ? ellisomat(E,0,1): E;
1735 2142 : A = gel(ellisotree(L), 2);
1736 2142 : vS = gel(ellweilcurve(L, NULL), 2); lvS = lg(vS);
1737 5397 : for (i = 1; i < lvS; i++)
1738 5397 : if (equali1(gmael(vS,i,1) ) && equali1(gmael(vS,i,2)))
1739 2142 : break;
1740 2142 : if (is_ell)
1741 2142 : return gerepileupto(av, path_to_manin(A, i, 1));
1742 0 : M = cgetg(lvS, t_VEC);
1743 0 : for (k = 1; k < lvS; k++)
1744 0 : gel(M,k) = path_to_manin(A, i, k);
1745 0 : return gerepileupto(av, M);
1746 : }
|
