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