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 4067 : make_velu_curve(GEN E, GEN t, GEN w)
31 : {
32 4067 : GEN A4, A6, a1 = ell_get_a1(E), a2 = ell_get_a2(E), a3 = ell_get_a3(E);
33 4067 : A4 = gsub(ell_get_a4(E), gmulsg(5L, t));
34 4067 : A6 = gsub(ell_get_a6(E), gadd(gmul(ell_get_b2(E), t), gmulsg(7L, w)));
35 4067 : 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 1736 : get_isog_vars(GEN phi, long *vx, long *vy)
42 : {
43 1736 : *vx = varn(gel(phi, 1));
44 1736 : *vy = varn(gel(phi, 2));
45 1736 : if (*vy == *vx) *vy = gvar2(gel(phi,2));
46 1736 : }
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 4480 : to_RgX(GEN P, long vx)
98 : {
99 4480 : return typ(P) == t_POL && varn(P)==vx? P: scalarpol_shallow(P, vx);
100 : }
101 :
102 : static GEN
103 448 : divy(GEN P0, GEN P1, GEN Q, GEN T, long vy)
104 : {
105 448 : GEN DP0, P0r = Q_remove_denom(P0, &DP0), P0D;
106 448 : GEN DP1, P1r = Q_remove_denom(P1, &DP1), P1D;
107 448 : GEN DQ, Qr = Q_remove_denom(Q, &DQ), P2;
108 448 : P0D = RgXQX_div(P0r, Qr, T);
109 448 : if (DP0) P0D = gdiv(P0D, DP0);
110 448 : P1D = RgXQX_div(P1r, Qr, T);
111 448 : if (DP1) P1D = gdiv(P1D, DP1);
112 448 : P2 = gadd(gmul(P1D, pol_x(vy)), P0D);
113 448 : if (DQ) P2 = gmul(P2, DQ);
114 448 : return P2;
115 : }
116 :
117 : static GEN
118 1792 : QXQH_eval(GEN P, GEN A, GEN B, GEN T)
119 : {
120 1792 : if (signe(P)==0) return mkvec2(P, pol_1(varn(P)));
121 1568 : return mkvec2(QXQX_homogenous_evalpow(P, A, B, T), gel(B,degpol(P)+1));
122 : }
123 :
124 : static GEN
125 1694 : ellnfcompisog(GEN nf, GEN F, GEN G)
126 : {
127 1694 : 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 1694 : if (!nf) return ellcompisog(F, G);
136 448 : T = nf_get_pol(nf);
137 448 : timer_start(&ti);
138 448 : checkellisog(F); F = liftpol_shallow(F);
139 448 : checkellisog(G); G = liftpol_shallow(G);
140 448 : get_isog_vars(F, &vx, &vy);
141 448 : Gh = to_RgX(gel(G,3),vx); Gh2 = QXQX_sqr(Gh, T); Gh3 = QXQX_mul(Gh, Gh2, T);
142 448 : F0 = to_RgX(polcoef_i(gel(F,2), 0, vy), vx);
143 448 : F1 = to_RgX(polcoef_i(gel(F,2), 1, vy), vx);
144 448 : G0 = to_RgX(polcoef_i(gel(G,2), 0, vy), vx);
145 448 : G1 = to_RgX(polcoef_i(gel(G,2), 1, vy), vx);
146 448 : d = maxss(maxss(degpol(gel(F,1)),degpol(gel(F,3))),maxss(degpol(F0),degpol(F1)));
147 448 : Gp = QXQX_powers(Gh2, d, T);
148 448 : f = QXQH_eval(to_RgX(gel(F,1),vx), gel(G,1), Gp, T);
149 448 : g0 = QXQH_eval(F0, to_RgX(gel(G,1),vx), Gp, T);
150 448 : g1 = QXQH_eval(F1, to_RgX(gel(G,1),vx), Gp, T);
151 448 : h = QXQH_eval(to_RgX(gel(F,3),vx), gel(G,1), Gp, T);
152 448 : K = Q_remove_denom(QXQX_mul(to_RgX(gel(h,1),vx), Gh, T), NULL);
153 448 : K = RgXQX_div(K, nfgcd(K, RgX_deriv(K), T, NULL), T);
154 448 : K = RgX_normalize(K);
155 448 : if (DEBUGLEVEL) timer_printf(&ti,"ellnfcompisog: nfgcd");
156 448 : K2 = QXQX_sqr(K, T); K3 = QXQX_mul(K, K2, T);
157 448 : if (DEBUGLEVEL) timer_printf(&ti,"ellnfcompisog: evalpow");
158 448 : h21 = QXQX_sqr(gel(h,1),T);
159 448 : h22 = QXQX_sqr(gel(h,2),T);
160 448 : h31 = QXQX_mul(gel(h,1), h21,T);
161 448 : h32 = QXQX_mul(gel(h,2), h22,T);
162 448 : if (DEBUGLEVEL) timer_printf(&ti,"h");
163 448 : f = RgXQX_div(QXQX_mul(QXQX_mul(K2, gel(f,1), T), h22, T),
164 448 : QXQX_mul(gel(f,2), h21, T), T);
165 448 : if (DEBUGLEVEL) timer_printf(&ti,"f");
166 448 : den = QXQX_mul(Gh3, gel(g1,2), T);
167 448 : if (DEBUGLEVEL) timer_printf(&ti,"ellnfcompisog: den");
168 448 : g0d = QXQX_mul(gel(g0,1),den, T);
169 448 : g01 = QXQX_mul(gel(g1,1),gel(g0,2),T);
170 448 : num0 = RgX_add(g0d, QXQX_mul(G0,g01, T));
171 448 : num1 = QXQX_mul(G1,g01, T);
172 448 : if (DEBUGLEVEL) timer_printf(&ti,"ellnfcompisog: num");
173 448 : k3h32 = QXQX_mul(K3,h32,T);
174 448 : gn0 = QXQX_mul(num0, k3h32, T);
175 448 : gn1 = QXQX_mul(num1, k3h32, T);
176 448 : if (DEBUGLEVEL) timer_printf(&ti,"ellnfcompisog: gn");
177 448 : gd = QXQX_mul(QXQX_mul(gel(g0,2), den, T), h31, T);
178 448 : if (DEBUGLEVEL) timer_printf(&ti,"ellnfcompisog: gd");
179 448 : g = divy(gn0, gn1, gd, T, vy);
180 448 : if (DEBUGLEVEL) timer_printf(&ti,"ellnfcompisog: divy");
181 448 : 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 : checkellpt(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 : long c;
278 427 : if (!oncurve(E,P))
279 7 : pari_err_DOMAIN("isogeny_from_kernel_point", "point", "not on", E, P);
280 420 : if (ell_is_inf(P))
281 : {
282 42 : if (only_image) return E;
283 28 : return mkvec2(E, isog_identity(vx,vy));
284 : }
285 :
286 378 : isog = NULL; c = 1;
287 : for (;;)
288 525 : {
289 903 : GEN tQ, xQ = gel(Q,1), uQ = ec_2divpol_evalx(E, xQ);
290 903 : int stop = 0;
291 903 : if (ellisweierstrasspoint(E,Q))
292 : { /* ord(P)=2c; take Q=[c]P into consideration and stop */
293 196 : tQ = ec_dFdx_evalQ(E, Q);
294 196 : stop = 1;
295 : }
296 : else
297 707 : tQ = ec_half_deriv_2divpol_evalx(E, xQ);
298 903 : t = gadd(t, tQ);
299 903 : w = gadd(w, gadd(uQ, gmul(tQ, xQ)));
300 903 : if (!only_image) isog = update_isogeny_polys(isog, E, Q,tQ,uQ, vx,vy);
301 903 : if (stop) break;
302 :
303 707 : Q = elladd(E, P, Q);
304 707 : ++c;
305 : /* IF x([c]P) = x([c-1]P) THEN [c]P = -[c-1]P and [2c-1]P = 0 */
306 707 : if (gequal(gel(Q,1), xQ)) break;
307 525 : if (gc_needed(av,1))
308 : {
309 0 : if(DEBUGMEM>1) pari_warn(warnmem,"isogeny_from_kernel_point");
310 0 : gerepileall(av, isog? 4: 3, &Q, &t, &w, &isog);
311 : }
312 : }
313 :
314 378 : EE = make_velu_curve(E, t, w);
315 378 : if (only_image) return EE;
316 :
317 224 : if (!isog) isog = isog_identity(vx,vy);
318 224 : f = gel(isog, 1);
319 224 : g = gel(isog, 2);
320 224 : if ( ! (typ(f) == t_RFRAC && typ(g) == t_RFRAC))
321 0 : pari_err_BUG("isogeny_from_kernel_point (f or g has wrong type)");
322 :
323 : /* Clean up the isogeny polynomials f, g and h so that the isogeny
324 : * is given by (x,y) -> (f(x)/h(x)^2, g(x,y)/h(x)^3) */
325 224 : h = gel(isog, 3);
326 224 : h2 = gsqr(h);
327 224 : h3 = gmul(h, h2);
328 224 : f = gmul(f, h2);
329 224 : g = gmul(g, h3);
330 224 : if (typ(f) != t_POL || typ(g) != t_POL)
331 0 : pari_err_BUG("isogeny_from_kernel_point (wrong denominator)");
332 224 : return mkvec2(EE, mkvec3(f,g, gel(isog,3)));
333 : }
334 :
335 : /* Given a t_POL x^n - s1 x^{n-1} + s2 x^{n-2} - s3 x^{n-3} + ...
336 : * return the first three power sums (Newton's identities):
337 : * p1 = s1
338 : * p2 = s1^2 - 2 s2
339 : * p3 = (s1^2 - 3 s2) s1 + 3 s3 */
340 : static void
341 3703 : first_three_power_sums(GEN pol, GEN *p1, GEN *p2, GEN *p3)
342 : {
343 3703 : long d = degpol(pol);
344 : GEN s1, s2, ms3;
345 :
346 3703 : *p1 = s1 = gneg(RgX_coeff(pol, d-1));
347 :
348 3703 : s2 = RgX_coeff(pol, d-2);
349 3703 : *p2 = gsub(gsqr(s1), gmulsg(2L, s2));
350 :
351 3703 : ms3 = RgX_coeff(pol, d-3);
352 3703 : *p3 = gadd(gmul(s1, gsub(*p2, s2)), gmulsg(-3L, ms3));
353 3703 : }
354 :
355 : /* Let E and a t_POL h of degree 1 or 3 whose roots are 2-torsion points on E.
356 : * - if only_image != 0, return [t, w] used to compute the equation of the
357 : * quotient by the given 2-torsion points
358 : * - else return [t,w, f,g,h], along with the contributions f, g and
359 : * h to the isogeny giving the quotient by h. Variables vx and vy are used
360 : * to create f, g and h, or ignored if only_image is zero */
361 :
362 : /* deg h = 1; 2-torsion contribution from Weierstrass point */
363 : static GEN
364 1792 : contrib_weierstrass_pt(GEN E, GEN h, long only_image, long vx, long vy)
365 : {
366 1792 : GEN p = ellbasechar(E);
367 1792 : GEN a1 = ell_get_a1(E);
368 1792 : GEN a3 = ell_get_a3(E);
369 1792 : GEN x0 = gneg(constant_coeff(h)); /* h = x - x0 */
370 1792 : GEN b = gadd(gmul(a1,x0), a3);
371 : GEN y0, Q, t, w, t1, t2, f, g;
372 :
373 1792 : if (!equalis(p, 2L)) /* char(k) != 2 ==> y0 = -b/2 */
374 1750 : y0 = gmul2n(gneg(b), -1);
375 : else
376 : { /* char(k) = 2 ==> y0 = sqrt(f(x0)) where E is y^2 + h(x) = f(x). */
377 42 : if (!gequal0(b)) pari_err_BUG("two_torsion_contrib (a1*x0+a3 != 0)");
378 42 : y0 = gsqrt(ec_f_evalx(E, x0), 0);
379 : }
380 1792 : Q = mkvec2(x0, y0);
381 1792 : t = ec_dFdx_evalQ(E, Q);
382 1792 : w = gmul(x0, t);
383 1792 : if (only_image) return mkvec2(t,w);
384 :
385 : /* Compute isogeny, f = (x - x0) * t */
386 518 : f = deg1pol_shallow(t, gmul(t, gneg(x0)), vx);
387 :
388 : /* g = (x - x0) * t * (a1 * (x - x0) + (y - y0)) */
389 518 : t1 = deg1pol_shallow(a1, gmul(a1, gneg(x0)), vx);
390 518 : t2 = deg1pol_shallow(gen_1, gneg(y0), vy);
391 518 : g = gmul(f, gadd(t1, t2));
392 518 : return mkvec5(t, w, f, g, h);
393 : }
394 : /* deg h =3; full 2-torsion contribution. NB: assume h is monic; base field
395 : * characteristic is odd or zero (cannot happen in char 2).*/
396 : static GEN
397 14 : contrib_full_tors(GEN E, GEN h, long only_image, long vx, long vy)
398 : {
399 : GEN p1, p2, p3, half_b2, half_b4, t, w, f, g;
400 14 : first_three_power_sums(h, &p1,&p2,&p3);
401 14 : half_b2 = gmul2n(ell_get_b2(E), -1);
402 14 : half_b4 = gmul2n(ell_get_b4(E), -1);
403 :
404 : /* t = 3*(p2 + b4/2) + p1 * b2/2 */
405 14 : t = gadd(gmulsg(3L, gadd(p2, half_b4)), gmul(p1, half_b2));
406 :
407 : /* w = 3 * p3 + p2 * b2/2 + p1 * b4/2 */
408 14 : w = gadd(gmulsg(3L, p3), gadd(gmul(p2, half_b2),
409 : gmul(p1, half_b4)));
410 14 : if (only_image) return mkvec2(t,w);
411 :
412 : /* Compute isogeny */
413 : {
414 7 : GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E), t1, t2;
415 7 : GEN s1 = gneg(RgX_coeff(h, 2));
416 7 : GEN dh = RgX_deriv(h);
417 7 : GEN psi2xy = gadd(deg1pol_shallow(a1, a3, vx),
418 : deg1pol_shallow(gen_2, gen_0, vy));
419 :
420 : /* f = -3 (3 x + b2/2 + s1) h + (3 x^2 + (b2/2) x + (b4/2)) h'*/
421 7 : t1 = RgX_mul(h, gmulsg(-3, deg1pol(stoi(3), gadd(half_b2, s1), vx)));
422 7 : t2 = mkpoln(3, stoi(3), half_b2, half_b4);
423 7 : setvarn(t2, vx);
424 7 : t2 = RgX_mul(dh, t2);
425 7 : f = RgX_add(t1, t2);
426 :
427 : /* 2g = psi2xy * (f'*h - f*h') - (a1*f + a3*h) * h; */
428 7 : t1 = RgX_sub(RgX_mul(RgX_deriv(f), h), RgX_mul(f, dh));
429 7 : t2 = RgX_mul(h, RgX_add(RgX_Rg_mul(f, a1), RgX_Rg_mul(h, a3)));
430 7 : g = RgX_divs(gsub(gmul(psi2xy, t1), t2), 2L);
431 :
432 7 : f = RgX_mul(f, h);
433 7 : g = RgX_mul(g, h);
434 : }
435 7 : return mkvec5(t, w, f, g, h);
436 : }
437 :
438 : /* Given E and a t_POL T whose roots define a subgroup G of E, return the factor
439 : * of T that corresponds to the 2-torsion points E[2] \cap G in G */
440 : INLINE GEN
441 3696 : two_torsion_part(GEN E, GEN T)
442 3696 : { return RgX_gcd(T, elldivpol(E, 2, varn(T))); }
443 :
444 : /* Return the jth Hasse derivative of the polynomial f = \sum_{i=0}^n a_i x^i,
445 : * i.e. \sum_{i=j}^n a_i \binom{i}{j} x^{i-j}. It is a derivation even when the
446 : * coefficient ring has positive characteristic */
447 : static GEN
448 98 : derivhasse(GEN f, ulong j)
449 : {
450 98 : ulong i, d = degpol(f);
451 : GEN df;
452 98 : if (gequal0(f) || d == 0) return pol_0(varn(f));
453 56 : if (j == 0) return gcopy(f);
454 56 : df = cgetg(2 + (d-j+1), t_POL);
455 56 : df[1] = f[1];
456 112 : for (i = j; i <= d; ++i) gel(df, i-j+2) = gmul(binomialuu(i,j), gel(f, i+2));
457 56 : return normalizepol(df);
458 : }
459 :
460 : static GEN
461 812 : non_two_torsion_abscissa(GEN E, GEN h0, GEN x)
462 : {
463 : GEN mp1, dh0, ddh0, t, u, t1, t2, t3;
464 812 : long m = degpol(h0);
465 812 : mp1 = gel(h0, m + 1); /* negative of first power sum */
466 812 : dh0 = RgX_deriv(h0);
467 812 : ddh0 = RgX_deriv(dh0);
468 812 : t = ec_2divpol_evalx(E, x);
469 812 : u = ec_half_deriv_2divpol_evalx(E, x);
470 812 : t1 = RgX_sub(RgX_sqr(dh0), RgX_mul(ddh0, h0));
471 812 : t2 = RgX_mul(u, RgX_mul(h0, dh0));
472 812 : t3 = RgX_mul(RgX_sqr(h0),
473 812 : deg1pol_shallow(stoi(2*m), gmulsg(2L, mp1), varn(x)));
474 : /* t * (dh0^2 - ddh0*h0) - u*dh0*h0 + (2*m*x - 2*s1) * h0^2); */
475 812 : return RgX_add(RgX_sub(RgX_mul(t, t1), t2), t3);
476 : }
477 :
478 : static GEN
479 1302 : isog_abscissa(GEN E, GEN kerp, GEN h0, GEN x, GEN two_tors)
480 : {
481 : GEN f0, f2, h2, t1, t2, t3;
482 1302 : f0 = (degpol(h0) > 0)? non_two_torsion_abscissa(E, h0, x): pol_0(varn(x));
483 1302 : f2 = gel(two_tors, 3);
484 1302 : h2 = gel(two_tors, 5);
485 :
486 : /* Combine f0 and f2 into the final abscissa of the isogeny. */
487 1302 : t1 = RgX_mul(x, RgX_sqr(kerp));
488 1302 : t2 = RgX_mul(f2, RgX_sqr(h0));
489 1302 : t3 = RgX_mul(f0, RgX_sqr(h2));
490 : /* x * kerp^2 + f2 * h0^2 + f0 * h2^2 */
491 1302 : return RgX_add(t1, RgX_add(t2, t3));
492 : }
493 :
494 : static GEN
495 1253 : non_two_torsion_ordinate_char_not2(GEN E, GEN f, GEN h, GEN psi2)
496 : {
497 1253 : GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E);
498 1253 : GEN df = RgX_deriv(f), dh = RgX_deriv(h);
499 : /* g = df * h * psi2/2 - f * dh * psi2
500 : * - (E.a1 * f + E.a3 * h^2) * h/2 */
501 1253 : GEN t1 = RgX_mul(df, RgX_mul(h, RgX_divs(psi2, 2L)));
502 1253 : GEN t2 = RgX_mul(f, RgX_mul(dh, psi2));
503 1253 : GEN t3 = RgX_mul(RgX_divs(h, 2L),
504 : RgX_add(RgX_Rg_mul(f, a1), RgX_Rg_mul(RgX_sqr(h), a3)));
505 1253 : return RgX_sub(RgX_sub(t1, t2), t3);
506 : }
507 :
508 : /* h = kerq */
509 : static GEN
510 49 : non_two_torsion_ordinate_char2(GEN E, GEN h, GEN x, GEN y)
511 : {
512 49 : GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E), a4 = ell_get_a4(E);
513 49 : GEN b2 = ell_get_b2(E), b4 = ell_get_b4(E), b6 = ell_get_b6(E);
514 : GEN h2, dh, dh2, ddh, D2h, D2dh, H, psi2, u, t, alpha;
515 : GEN p1, t1, t2, t3, t4;
516 49 : long m, vx = varn(x);
517 :
518 49 : h2 = RgX_sqr(h);
519 49 : dh = RgX_deriv(h);
520 49 : dh2 = RgX_sqr(dh);
521 49 : ddh = RgX_deriv(dh);
522 49 : H = RgX_sub(dh2, RgX_mul(h, ddh));
523 49 : D2h = derivhasse(h, 2);
524 49 : D2dh = derivhasse(dh, 2);
525 49 : psi2 = deg1pol_shallow(a1, a3, vx);
526 49 : u = mkpoln(3, b2, gen_0, b6);
527 49 : setvarn(u, vx);
528 49 : t = deg1pol_shallow(b2, b4, vx);
529 49 : alpha = mkpoln(4, a1, a3, gmul(a1, a4), gmul(a3, a4));
530 49 : setvarn(alpha, vx);
531 49 : m = degpol(h);
532 49 : p1 = RgX_coeff(h, m-1); /* first power sum */
533 :
534 49 : t1 = gmul(gadd(gmul(a1, p1), gmulgu(a3, m)), RgX_mul(h,h2));
535 :
536 49 : t2 = gmul(a1, gadd(gmul(a1, gadd(y, psi2)), RgX_add(RgX_Rg_add(RgX_sqr(x), a4), t)));
537 49 : t2 = gmul(t2, gmul(dh, h2));
538 :
539 49 : t3 = gadd(gmul(y, t), RgX_add(alpha, RgX_Rg_mul(u, a1)));
540 49 : t3 = gmul(t3, RgX_mul(h, H));
541 :
542 49 : t4 = gmul(u, psi2);
543 49 : t4 = gmul(t4, RgX_sub(RgX_sub(RgX_mul(h2, D2dh), RgX_mul(dh, H)),
544 : RgX_mul(h, RgX_mul(dh, D2h))));
545 :
546 49 : return gadd(t1, gadd(t2, gadd(t3, t4)));
547 : }
548 :
549 : static GEN
550 1302 : isog_ordinate(GEN E, GEN kerp, GEN kerq, GEN x, GEN y, GEN two_tors, GEN f)
551 : {
552 : GEN g;
553 1302 : if (! equalis(ellbasechar(E), 2L)) {
554 : /* FIXME: We don't use (hence don't need to calculate)
555 : * g2 = gel(two_tors, 4) when char(k) != 2. */
556 1253 : GEN psi2 = ec_dmFdy_evalQ(E, mkvec2(x, y));
557 1253 : g = non_two_torsion_ordinate_char_not2(E, f, kerp, psi2);
558 : } else {
559 49 : GEN h2 = gel(two_tors, 5);
560 49 : GEN g2 = gmul(gel(two_tors, 4), RgX_mul(kerq, RgX_sqr(kerq)));
561 49 : GEN g0 = non_two_torsion_ordinate_char2(E, kerq, x, y);
562 49 : g0 = gmul(g0, RgX_mul(h2, RgX_sqr(h2)));
563 49 : g = gsub(gmul(y, RgX_mul(kerp, RgX_sqr(kerp))), gadd(g2, g0));
564 : }
565 1302 : return g;
566 : }
567 :
568 : /* Given an elliptic curve E and a polynomial kerp whose roots give the
569 : * x-coordinates of a subgroup G of E, return the curve E/G and,
570 : * if only_image is zero, the isogeny pi:E -> E/G. Variables vx and vy are
571 : * used to describe the isogeny (and are ignored if only_image is zero). */
572 : static GEN
573 3696 : isogeny_from_kernel_poly(GEN E, GEN kerp, long only_image, long vx, long vy)
574 : {
575 : long m;
576 3696 : GEN b2 = ell_get_b2(E), b4 = ell_get_b4(E), b6 = ell_get_b6(E);
577 : GEN p1, p2, p3, x, y, f, g, two_tors, EE, t, w;
578 3696 : GEN kerh = two_torsion_part(E, kerp);
579 3696 : GEN kerq = RgX_divrem(kerp, kerh, ONLY_DIVIDES);
580 3696 : if (!kerq) pari_err_BUG("isogeny_from_kernel_poly");
581 : /* isogeny degree: 2*degpol(kerp)+1-degpol(kerh) */
582 3696 : m = degpol(kerq);
583 :
584 3696 : kerp = RgX_normalize(kerp);
585 3696 : kerq = RgX_normalize(kerq);
586 3696 : kerh = RgX_normalize(kerh);
587 3696 : switch(degpol(kerh))
588 : {
589 1883 : case 0:
590 1883 : two_tors = only_image? mkvec2(gen_0, gen_0):
591 777 : mkvec5(gen_0, gen_0, pol_0(vx), pol_0(vx), pol_1(vx));
592 1883 : break;
593 1792 : case 1:
594 1792 : two_tors = contrib_weierstrass_pt(E, kerh, only_image,vx,vy);
595 1792 : break;
596 14 : case 3:
597 14 : two_tors = contrib_full_tors(E, kerh, only_image,vx,vy);
598 14 : break;
599 7 : default:
600 7 : two_tors = NULL;
601 7 : pari_err_DOMAIN("isogeny_from_kernel_poly", "kernel polynomial",
602 : "does not define a subgroup of", E, kerp);
603 : }
604 3689 : first_three_power_sums(kerq,&p1,&p2,&p3);
605 3689 : x = pol_x(vx);
606 3689 : y = pol_x(vy);
607 :
608 : /* t = 6 * p2 + b2 * p1 + m * b4, */
609 3689 : t = gadd(gmulsg(6L, p2), gadd(gmul(b2, p1), gmulsg(m, b4)));
610 :
611 : /* w = 10 * p3 + 2 * b2 * p2 + 3 * b4 * p1 + m * b6, */
612 3689 : w = gadd(gmulsg(10L, p3),
613 : gadd(gmul(gmulsg(2L, b2), p2),
614 : gadd(gmul(gmulsg(3L, b4), p1), gmulsg(m, b6))));
615 :
616 3689 : EE = make_velu_curve(E, gadd(t, gel(two_tors, 1)),
617 3689 : gadd(w, gel(two_tors, 2)));
618 3689 : if (only_image) return EE;
619 :
620 1302 : f = isog_abscissa(E, kerp, kerq, x, two_tors);
621 1302 : g = isog_ordinate(E, kerp, kerq, x, y, two_tors, f);
622 1302 : return mkvec2(EE, mkvec3(f,g,kerp));
623 : }
624 :
625 : /* Given an elliptic curve E and a subgroup G of E, return the curve
626 : * E/G and, if only_image is zero, the isogeny corresponding
627 : * to the canonical surjection pi:E -> E/G. The variables vx and
628 : * vy are used to describe the isogeny (and are ignored if
629 : * only_image is zero). The subgroup G may be given either as
630 : * a generating point P on E or as a polynomial kerp whose roots are
631 : * the x-coordinates of the points in G */
632 : GEN
633 1092 : ellisogeny(GEN E, GEN G, long only_image, long vx, long vy)
634 : {
635 1092 : pari_sp av = avma;
636 : GEN j, z;
637 1092 : checkell(E);j = ell_get_j(E);
638 1092 : if (vx < 0) vx = 0;
639 1092 : if (vy < 0) vy = 1;
640 1092 : if (varncmp(vx, vy) >= 0)
641 7 : pari_err_PRIORITY("ellisogeny", pol_x(vx), "<=", vy);
642 1085 : if (!only_image && varncmp(vy, gvar(j)) >= 0)
643 7 : pari_err_PRIORITY("ellisogeny", j, ">=", vy);
644 1078 : switch(typ(G))
645 : {
646 441 : case t_VEC:
647 441 : checkellpt(G);
648 441 : if (!ell_is_inf(G))
649 : {
650 399 : GEN x = gel(G,1), y = gel(G,2);
651 399 : if (!only_image)
652 : {
653 245 : if (varncmp(vy, gvar(x)) >= 0)
654 7 : pari_err_PRIORITY("ellisogeny", x, ">=", vy);
655 238 : if (varncmp(vy, gvar(y)) >= 0)
656 7 : pari_err_PRIORITY("ellisogeny", y, ">=", vy);
657 : }
658 : }
659 427 : z = isogeny_from_kernel_point(E, G, only_image, vx, vy);
660 420 : break;
661 630 : case t_POL:
662 630 : if (!only_image && varncmp(vy, gvar(constant_coeff(G))) >= 0)
663 7 : pari_err_PRIORITY("ellisogeny", constant_coeff(G), ">=", vy);
664 623 : z = isogeny_from_kernel_poly(E, G, only_image, vx, vy);
665 616 : break;
666 7 : default:
667 7 : z = NULL;
668 7 : pari_err_TYPE("ellisogeny", G);
669 : }
670 1036 : return gerepilecopy(av, z);
671 : }
672 :
673 : static GEN
674 3388 : trivial_isogeny(void)
675 : {
676 3388 : return mkvec3(pol_x(0), scalarpol(pol_x(1), 0), pol_1(0));
677 : }
678 :
679 : static GEN
680 266 : isogeny_a4a6(GEN E)
681 : {
682 266 : GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E), b2 = ell_get_b2(E);
683 266 : retmkvec3(deg1pol(gen_1, gdivgu(b2, 12), 0),
684 : deg1pol(gdivgu(a1,2), deg1pol(gen_1, gdivgu(a3,2), 1), 0),
685 : pol_1(0));
686 : }
687 :
688 : static GEN
689 266 : invisogeny_a4a6(GEN E)
690 : {
691 266 : GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E), b2 = ell_get_b2(E);
692 266 : retmkvec3(deg1pol(gen_1, gdivgs(b2, -12), 0),
693 : deg1pol(gdivgs(a1,-2),
694 : deg1pol(gen_1, gadd(gdivgs(a3,-2), gdivgu(gmul(b2,a1), 24)), 1), 0),
695 : pol_1(0));
696 : }
697 :
698 : static GEN
699 3080 : RgXY_eval(GEN P, GEN x, GEN y)
700 : {
701 3080 : return poleval(poleval(P,x), y);
702 : }
703 :
704 : static GEN
705 560 : twistisogeny(GEN iso, GEN d)
706 : {
707 560 : GEN d2 = gsqr(d), d3 = gmul(d, d2);
708 560 : return mkvec3(gdiv(gel(iso,1), d2), gdiv(gel(iso,2), d3), gel(iso, 3));
709 : }
710 :
711 : static GEN
712 2513 : ellisog_by_Kohel(GEN a4, GEN a6, long n, GEN ker, GEN kert, long flag)
713 : {
714 2513 : GEN E = ellinit(mkvec2(a4, a6), NULL, DEFAULTPREC);
715 2513 : GEN F = isogeny_from_kernel_poly(E, ker, flag, 0, 1);
716 2513 : GEN Et = ellinit(flag ? F: gel(F, 1), NULL, DEFAULTPREC);
717 2513 : GEN c4t = ell_get_c4(Et), c6t = ell_get_c6(Et), jt = ell_get_j(Et);
718 2513 : if (!flag)
719 : {
720 560 : GEN Ft = isogeny_from_kernel_poly(Et, kert, flag, 0, 1);
721 560 : GEN isot = twistisogeny(gel(Ft, 2), stoi(n));
722 560 : return mkvec5(c4t, c6t, jt, gel(F, 2), isot);
723 : }
724 1953 : else return mkvec3(c4t, c6t, jt);
725 : }
726 :
727 : static GEN
728 2345 : ellisog_by_roots(GEN a4, GEN a6, long n, GEN z, long flag)
729 : {
730 2345 : GEN k = deg1pol_shallow(gen_1, gneg(z), 0);
731 2345 : GEN kt= deg1pol_shallow(gen_1, gmulsg(n,z), 0);
732 2345 : return ellisog_by_Kohel(a4, a6, n, k, kt, flag);
733 : }
734 :
735 : /* n = 2 or 3 */
736 : static GEN
737 3500 : a4a6_divpol(GEN a4, GEN a6, long n)
738 : {
739 3500 : if (n == 2) return mkpoln(4, gen_1, gen_0, a4, a6);
740 1848 : return mkpoln(5, utoi(3), gen_0, gmulgu(a4,6) , gmulgu(a6,12),
741 : gneg(gsqr(a4)));
742 : }
743 :
744 : static GEN
745 3500 : ellisograph_Kohel_iso(GEN nf, GEN e, long n, GEN z, GEN *pR, long flag)
746 : {
747 : long i, r;
748 3500 : GEN R, V, c4 = gel(e,1), c6 = gel(e,2);
749 3500 : GEN a4 = gdivgs(c4, -48), a6 = gdivgs(c6, -864);
750 3500 : GEN P = a4a6_divpol(a4, a6, n);
751 3500 : R = nfroots(nf, z ? RgX_div_by_X_x(P, z, NULL): P);
752 3500 : if (pR) *pR = R;
753 3500 : r = lg(R); V = cgetg(r, t_VEC);
754 5845 : for (i=1; i < r; i++) gel(V,i) = ellisog_by_roots(a4, a6, n, gel(R,i), flag);
755 3500 : return V;
756 : }
757 :
758 : static GEN
759 3311 : ellisograph_Kohel_r(GEN nf, GEN e, long n, GEN z, long flag)
760 : {
761 3311 : GEN R, iso = ellisograph_Kohel_iso(nf, e, n, z, &R, flag);
762 3311 : long i, r = lg(iso);
763 3311 : GEN V = cgetg(r, t_VEC);
764 5467 : for (i=1; i < r; i++)
765 2156 : gel(V,i) = ellisograph_Kohel_r(nf, gel(iso,i), n, gmulgs(gel(R,i), -n), flag);
766 3311 : return mkvec2(e, V);
767 : }
768 :
769 : static GEN
770 336 : corr(GEN c4, GEN c6)
771 : {
772 336 : GEN c62 = gmul2n(c6, 1);
773 336 : return gadd(gdiv(gsqr(c4), c62), gdiv(c62, gmulgu(c4,3)));
774 : }
775 :
776 : static GEN
777 336 : elkies98(GEN a4, GEN a6, long l, GEN s, GEN a4t, GEN a6t)
778 : {
779 : GEN C, P, S;
780 : long i, n, d;
781 336 : d = l == 2 ? 1 : l>>1;
782 336 : C = cgetg(d+1, t_VEC);
783 336 : gel(C, 1) = gdivgu(gsub(a4, a4t), 5);
784 336 : if (d >= 2)
785 336 : gel(C, 2) = gdivgu(gsub(a6, a6t), 7);
786 336 : if (d >= 3)
787 224 : gel(C, 3) = gdivgu(gsub(gsqr(gel(C, 1)), gmul(a4, gel(C, 1))), 3);
788 2870 : for (n = 3; n < d; ++n)
789 : {
790 2534 : GEN s = gen_0;
791 61390 : for (i = 1; i < n; i++)
792 58856 : s = gadd(s, gmul(gel(C, i), gel(C, n-i)));
793 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));
794 : }
795 336 : P = cgetg(d+2, t_VEC);
796 336 : gel(P, 1 + 0) = stoi(d);
797 336 : gel(P, 1 + 1) = s;
798 336 : if (d >= 2)
799 336 : gel(P, 1 + 2) = gdivgu(gsub(gel(C, 1), gmulgu(a4, 2*d)), 6);
800 3094 : for (n = 2; n < d; ++n)
801 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);
802 336 : S = cgetg(d+3, t_POL);
803 336 : S[1] = evalsigne(1) | evalvarn(0);
804 336 : gel(S, 2 + d - 0) = gen_1;
805 336 : gel(S, 2 + d - 1) = gneg(s);
806 3430 : for (n = 2; n <= d; ++n)
807 : {
808 3094 : GEN s = gen_0;
809 68362 : for (i = 1; i <= n; ++i)
810 : {
811 65268 : GEN p = gmul(gel(P, 1+i), gel(S, 2 + d - (n-i)));
812 65268 : s = gadd(s, p);
813 : }
814 3094 : gel(S, 2 + d - n) = gdivgs(s, -n);
815 : }
816 336 : return S;
817 : }
818 :
819 : static GEN
820 770 : ellisog_by_jt(GEN c4, GEN c6, GEN jt, GEN jtp, GEN s0, long n, long flag)
821 : {
822 770 : GEN jtp2 = gsqr(jtp), den = gmul(jt, gsubgs(jt, 1728));
823 770 : GEN c4t = gdiv(jtp2, den);
824 770 : GEN c6t = gdiv(gmul(jtp, c4t), jt);
825 770 : if (flag)
826 602 : return mkvec3(c4t, c6t, jt);
827 : else
828 : {
829 168 : GEN co = corr(c4, c6);
830 168 : GEN cot = corr(c4t, c6t);
831 168 : GEN s = gmul2n(gmulgs(gadd(gadd(s0, co), gmulgs(cot,-n)), -n), -2);
832 168 : GEN a4 = gdivgs(c4, -48), a6 = gdivgs(c6, -864);
833 168 : GEN a4t = gmul(gdivgs(c4t, -48), powuu(n,4)), a6t = gmul(gdivgs(c6t, -864), powuu(n,6));
834 168 : GEN ker = elkies98(a4, a6, n, s, a4t, a6t);
835 168 : GEN st = gmulgs(s, -n);
836 168 : GEN a4tt = gmul(a4,powuu(n,4)), a6tt = gmul(a6,powuu(n,6));
837 168 : GEN kert = elkies98(a4t, a6t, n, st, a4tt, a6tt);
838 168 : return ellisog_by_Kohel(a4, a6, n, ker, kert, flag);
839 : }
840 : }
841 :
842 : /* RENE SCHOOF, Counting points on elliptic curves over finite fields,
843 : * Journal de Theorie des Nombres de Bordeaux, tome 7, no 1 (1995), p. 219-254.
844 : * http://www.numdam.org/item?id=JTNB_1995__7_1_219_0 */
845 : static GEN
846 616 : ellisog_by_j(GEN e, GEN jt, long n, GEN P, long flag)
847 : {
848 616 : pari_sp av = avma;
849 616 : GEN c4 = gel(e,1), c6 = gel(e, 2), j = gel(e, 3);
850 616 : GEN Px = RgX_deriv(P), Py = RgXY_derivx(P);
851 616 : GEN Pxj = RgXY_eval(Px, j, jt), Pyj = RgXY_eval(Py, j, jt);
852 616 : GEN Pxx = RgX_deriv(Px), Pxy = RgX_deriv(Py), Pyy = RgXY_derivx(Py);
853 616 : GEN Pxxj = RgXY_eval(Pxx,j,jt);
854 616 : GEN Pxyj = RgXY_eval(Pxy,j,jt);
855 616 : GEN Pyyj = RgXY_eval(Pyy,j,jt);
856 616 : GEN c6c4 = gdiv(c6, c4);
857 616 : GEN jp = gmul(j, c6c4);
858 616 : GEN jtp = gdivgs(gmul(jp, gdiv(Pxj, Pyj)), -n);
859 616 : GEN jtpn = gmulgs(jtp, n);
860 616 : GEN s0 = gdiv(gadd(gadd(gmul(gsqr(jp),Pxxj),gmul(gmul(jp,jtpn),gmul2n(Pxyj,1))),
861 : gmul(gsqr(jtpn),Pyyj)),gmul(jp,Pxj));
862 616 : GEN et = ellisog_by_jt(c4, c6, jt, jtp, s0, n, flag);
863 616 : return gerepilecopy(av, et);
864 : }
865 :
866 : static GEN
867 1344 : ellisograph_iso(GEN nf, GEN e, ulong p, GEN P, GEN oj, long flag)
868 : {
869 : long i, r;
870 : GEN Pj, R, V;
871 1344 : if (!P) return ellisograph_Kohel_iso(nf, e, p, oj, NULL, flag);
872 1155 : Pj = poleval(P, gel(e,3));
873 1155 : R = nfroots(nf,oj ? RgX_div_by_X_x(Pj, oj, NULL):Pj);
874 1155 : r = lg(R);
875 1155 : V = cgetg(r, t_VEC);
876 1771 : for (i=1; i < r; i++) gel(V, i) = ellisog_by_j(e, gel(R, i), p, P, flag);
877 1155 : return V;
878 : }
879 :
880 : static GEN
881 1113 : ellisograph_r(GEN nf, GEN e, ulong p, GEN P, GEN oj, long flag)
882 : {
883 1113 : GEN j = gel(e,3), iso = ellisograph_iso(nf, e, p, P, oj, flag);
884 1113 : long i, r = lg(iso);
885 1113 : GEN V = cgetg(r, t_VEC);
886 1687 : for (i=1; i < r; i++) gel(V,i) = ellisograph_r(nf, gel(iso,i), p, P, j, flag);
887 1113 : return mkvec2(e, V);
888 : }
889 :
890 : static GEN
891 1092 : ellisograph_a4a6(GEN E, long flag)
892 : {
893 1092 : GEN c4 = ell_get_c4(E), c6 = ell_get_c6(E), j = ell_get_j(E);
894 1358 : return flag ? mkvec3(c4, c6, j):
895 266 : mkvec5(c4, c6, j, isogeny_a4a6(E), invisogeny_a4a6(E));
896 : }
897 :
898 : static GEN
899 154 : ellisograph_dummy(GEN E, long n, GEN jt, GEN jtt, GEN s0, long flag)
900 : {
901 154 : GEN c4 = ell_get_c4(E), c6 = ell_get_c6(E), c6c4 = gdiv(c6, c4);
902 154 : GEN jtp = gmul(c6c4, gdivgs(gmul(jt, jtt), -n));
903 154 : GEN iso = ellisog_by_jt(c4, c6, jt, jtp, gmul(s0, c6c4), n, flag);
904 154 : GEN v = mkvec2(iso, cgetg(1, t_VEC));
905 154 : return mkvec2(ellisograph_a4a6(E, flag), mkvec(v));
906 : }
907 :
908 : static GEN
909 1694 : isograph_p(GEN nf, GEN e, ulong p, GEN P, long flag)
910 : {
911 1694 : pari_sp av = avma;
912 : GEN iso;
913 1694 : if (P)
914 539 : iso = ellisograph_r(nf, e, p, P, NULL, flag);
915 : else
916 1155 : iso = ellisograph_Kohel_r(nf, e, p, NULL, flag);
917 1694 : return gerepilecopy(av, iso);
918 : }
919 :
920 : static GEN
921 812 : get_polmodular(ulong p, long vx, long vy)
922 812 : { return p > 3 ? polmodular_ZXX(p,0,vx,vy): NULL; }
923 : static GEN
924 560 : ellisograph_p(GEN nf, GEN E, ulong p, long flag)
925 : {
926 560 : GEN e = ellisograph_a4a6(E, flag);
927 560 : GEN P = get_polmodular(p, 0, 1);
928 560 : return isograph_p(nf, e, p, P, flag);
929 : }
930 :
931 : static long
932 10129 : etree_nbnodes(GEN T)
933 : {
934 10129 : GEN F = gel(T,2);
935 10129 : long n = 1, i, l = lg(F);
936 16415 : for (i = 1; i < l; i++) n += etree_nbnodes(gel(F, i));
937 10129 : return n;
938 : }
939 :
940 : static long
941 4424 : etree_listr(GEN nf, GEN T, GEN V, long n, GEN u, GEN ut)
942 : {
943 4424 : GEN E = gel(T, 1), F = gel(T,2);
944 4424 : long i, l = lg(F);
945 4424 : GEN iso, isot = NULL;
946 4424 : if (lg(E) == 6)
947 : {
948 735 : iso = ellnfcompisog(nf,gel(E,4), u);
949 735 : isot = ellnfcompisog(nf,ut, gel(E,5));
950 735 : gel(V, n) = mkvec5(gel(E,1), gel(E,2), gel(E,3), iso, isot);
951 : } else
952 : {
953 3689 : gel(V, n) = mkvec3(gel(E,1), gel(E,2), gel(E,3));
954 3689 : iso = u;
955 : }
956 7154 : for (i = 1; i < l; i++)
957 2730 : n = etree_listr(nf, gel(F, i), V, n + 1, iso, isot);
958 4424 : return n;
959 : }
960 :
961 : static GEN
962 1694 : etree_list(GEN nf, GEN T)
963 : {
964 1694 : long n = etree_nbnodes(T);
965 1694 : GEN V = cgetg(n+1, t_VEC);
966 1694 : (void) etree_listr(nf, T, V, 1, trivial_isogeny(), trivial_isogeny());
967 1694 : return V;
968 : }
969 :
970 : static long
971 4424 : etree_distmatr(GEN T, GEN M, long n)
972 : {
973 4424 : GEN F = gel(T,2);
974 4424 : long i, j, lF = lg(F), m = n + 1;
975 4424 : GEN V = cgetg(lF, t_VECSMALL);
976 4424 : mael(M, n, n) = 0;
977 7154 : for(i = 1; i < lF; i++)
978 2730 : V[i] = m = etree_distmatr(gel(F,i), M, m);
979 7154 : for(i = 1; i < lF; i++)
980 : {
981 2730 : long mi = i==1 ? n+1: V[i-1];
982 6188 : for(j = mi; j < V[i]; j++)
983 : {
984 3458 : mael(M,n,j) = 1 + mael(M, mi, j);
985 3458 : mael(M,j,n) = 1 + mael(M, j, mi);
986 : }
987 6706 : for(j = 1; j < lF; j++)
988 3976 : if (i != j)
989 : {
990 1246 : long i1, j1, mj = j==1 ? n+1: V[j-1];
991 2681 : for (i1 = mi; i1 < V[i]; i1++)
992 3171 : for(j1 = mj; j1 < V[j]; j1++)
993 1736 : mael(M,i1,j1) = 2 + mael(M,mj,j1) + mael(M,i1,mi);
994 : }
995 : }
996 4424 : return m;
997 : }
998 :
999 : static GEN
1000 1694 : etree_distmat(GEN T)
1001 : {
1002 1694 : long i, n = etree_nbnodes(T);
1003 1694 : GEN M = cgetg(n+1, t_MAT);
1004 6118 : for(i = 1; i <= n; i++) gel(M,i) = cgetg(n+1, t_VECSMALL);
1005 1694 : (void)etree_distmatr(T, M, 1);
1006 1694 : return M;
1007 : }
1008 :
1009 : static GEN
1010 1694 : distmat_pow(GEN E, ulong p)
1011 : {
1012 1694 : long i, j, l = lg(E);
1013 1694 : GEN M = cgetg(l, t_MAT);
1014 6118 : for(i = 1; i < l; i++)
1015 : {
1016 4424 : gel(M,i) = cgetg(l, t_COL);
1017 17500 : for(j = 1; j < l; j++) gmael(M,i,j) = powuu(p,mael(E,i,j));
1018 : }
1019 1694 : return M;
1020 : }
1021 :
1022 : /* Assume there is a single p-isogeny */
1023 : static GEN
1024 154 : isomatdbl(GEN nf, GEN L, GEN M, ulong p, GEN T2, long flag)
1025 : {
1026 154 : long i, j, n = lg(L) -1;
1027 154 : GEN P = get_polmodular(p, 0, 1);
1028 154 : GEN V = cgetg(2*n+1, t_VEC), N = cgetg(2*n+1, t_MAT);
1029 539 : for (i=1; i <= n; i++)
1030 : {
1031 385 : GEN F, E, e = gel(L,i);
1032 385 : if (i == 1)
1033 154 : F = gmael(T2, 2, 1);
1034 : else
1035 : {
1036 231 : F = ellisograph_iso(nf, e, p, P, NULL, flag);
1037 231 : if (lg(F) != 2) pari_err_BUG("isomatdbl");
1038 : }
1039 385 : E = gel(F, 1);
1040 385 : if (flag)
1041 273 : E = mkvec3(gel(E,1), gel(E,2), gel(E,3));
1042 : else
1043 : {
1044 112 : GEN iso = ellnfcompisog(nf, gel(E,4), gel(e, 4));
1045 112 : GEN isot = ellnfcompisog(nf, gel(e,5), gel(E, 5));
1046 112 : E = mkvec5(gel(E,1), gel(E,2), gel(E,3), iso, isot);
1047 : }
1048 385 : gel(V, i) = e;
1049 385 : gel(V, i+n) = E;
1050 : }
1051 924 : for (i=1; i <= 2*n; i++) gel(N, i) = cgetg(2*n+1, t_COL);
1052 539 : for (i=1; i <= n; i++)
1053 1400 : for (j=1; j <= n; j++)
1054 : {
1055 1015 : gcoeff(N,i,j) = gcoeff(N,i+n,j+n) = gcoeff(M,i,j);
1056 1015 : gcoeff(N,i,j+n) = gcoeff(N,i+n,j) = muliu(gcoeff(M,i,j), p);
1057 : }
1058 154 : return mkvec2(V, N);
1059 : }
1060 :
1061 : static ulong
1062 462 : ellQ_exceptional_iso(GEN j, GEN *jt, GEN *jtp, GEN *s0)
1063 : {
1064 462 : *jt = j; *jtp = gen_1;
1065 462 : if (typ(j)==t_INT)
1066 : {
1067 252 : long js = itos_or_0(j);
1068 : GEN j37;
1069 252 : if (js==-32768) { *s0 = mkfracss(-1156,539); return 11; }
1070 238 : if (js==-121)
1071 14 : { *jt = stoi(-24729001) ; *jtp = mkfracss(4973,5633);
1072 14 : *s0 = mkfracss(-1961682050,1204555087); return 11;}
1073 224 : if (js==-24729001)
1074 14 : { *jt = stoi(-121); *jtp = mkfracss(5633,4973);
1075 14 : *s0 = mkfracss(-1961682050,1063421347); return 11;}
1076 210 : if (js==-884736)
1077 14 : { *s0 = mkfracss(-1100,513); return 19; }
1078 196 : j37 = negi(uu32toi(37876312,1780746325));
1079 196 : if (js==-9317)
1080 : {
1081 14 : *jt = j37;
1082 14 : *jtp = mkfracss(1984136099,496260169);
1083 14 : *s0 = mkfrac(negi(uu32toi(457100760,4180820796UL)),
1084 : uu32toi(89049913, 4077411069UL));
1085 14 : return 37;
1086 : }
1087 182 : if (equalii(j, j37))
1088 : {
1089 14 : *jt = stoi(-9317);
1090 14 : *jtp = mkfrac(utoi(496260169),utoi(1984136099UL));
1091 14 : *s0 = mkfrac(negi(uu32toi(41554614,2722784052UL)),
1092 : uu32toi(32367030,2614994557UL));
1093 14 : return 37;
1094 : }
1095 168 : if (js==-884736000)
1096 14 : { *s0 = mkfracss(-1073708,512001); return 43; }
1097 154 : if (equalii(j, negi(powuu(5280,3))))
1098 14 : { *s0 = mkfracss(-176993228,85184001); return 67; }
1099 140 : if (equalii(j, negi(powuu(640320,3))))
1100 14 : { *s0 = mkfrac(negi(uu32toi(72512,1969695276)), uu32toi(35374,1199927297));
1101 14 : return 163; }
1102 : } else
1103 : {
1104 210 : GEN j1 = mkfracss(-297756989,2);
1105 210 : GEN j2 = mkfracss(-882216989,131072);
1106 210 : if (gequal(j, j1))
1107 : {
1108 14 : *jt = j2; *jtp = mkfracss(1503991,2878441);
1109 14 : *s0 = mkfrac(negi(uu32toi(121934,548114672)),uu32toi(77014,117338383));
1110 14 : return 17;
1111 : }
1112 196 : if (gequal(j, j2))
1113 : {
1114 14 : *jt = j1; *jtp = mkfracss(2878441,1503991);
1115 14 : *s0 = mkfrac(negi(uu32toi(121934,548114672)),uu32toi(40239,4202639633UL));
1116 14 : return 17;
1117 : }
1118 : }
1119 308 : return 0;
1120 : }
1121 :
1122 : static GEN
1123 1694 : nfmkisomat(GEN nf, ulong p, GEN T)
1124 1694 : { return mkvec2(etree_list(nf,T), distmat_pow(etree_distmat(T),p)); }
1125 : static GEN
1126 455 : mkisomat(ulong p, GEN T)
1127 455 : { return nfmkisomat(NULL, p, T); }
1128 : static GEN
1129 154 : mkisomatdbl(ulong p, GEN T, ulong p2, GEN T2, long flag)
1130 : {
1131 154 : GEN v = mkisomat(p,T);
1132 154 : return isomatdbl(NULL, gel(v,1), gel(v,2), p2, T2, flag);
1133 : }
1134 :
1135 : /*See M.A Kenku, On the number of Q-isomorphism classes of elliptic curves in
1136 : * each Q-isogeny class, Journal of Number Theory Volume 15, Issue 2,
1137 : * October 1982, pp 199-202,
1138 : * http://www.sciencedirect.com/science/article/pii/0022314X82900257 */
1139 : enum { _2 = 1, _3 = 2, _5 = 4, _7 = 8, _13 = 16 };
1140 : static ulong
1141 308 : ellQ_goodl(GEN E)
1142 : {
1143 : forprime_t T;
1144 308 : long i, CM = ellQ_get_CM(E);
1145 308 : ulong mask = 31;
1146 308 : GEN disc = ell_get_disc(E);
1147 308 : pari_sp av = avma;
1148 308 : u_forprime_init(&T, 17UL,ULONG_MAX);
1149 6363 : for(i=1; mask && i<=20; i++)
1150 : {
1151 6055 : ulong p = u_forprime_next(&T);
1152 6055 : if (umodiu(disc,p)==0) i--;
1153 : else
1154 : {
1155 6048 : long t = ellap_CM_fast(E, p, CM), D = t*t - 4*p;
1156 6048 : if (t%2) mask &= ~_2;
1157 6048 : if ((mask & _3) && kross(D,3)==-1) mask &= ~_3;
1158 6048 : if ((mask & _5) && kross(D,5)==-1) mask &= ~_5;
1159 6048 : if ((mask & _7) && kross(D,7)==-1) mask &= ~_7;
1160 6048 : if ((mask &_13) && kross(D,13)==-1) mask &= ~_13;
1161 : }
1162 : }
1163 308 : return gc_ulong(av, mask);
1164 : }
1165 :
1166 : static long
1167 182 : ellQ_goodl_l(GEN E, long l)
1168 : {
1169 : forprime_t T;
1170 : long i;
1171 182 : GEN disc = ell_get_disc(E);
1172 182 : pari_sp av = avma;
1173 182 : u_forprime_init(&T, 17UL, ULONG_MAX);
1174 2408 : for (i=1; i<=20; i++)
1175 : {
1176 2303 : ulong p = u_forprime_next(&T);
1177 2303 : if (umodiu(disc,p)==0) { i--; continue; }
1178 : else
1179 : {
1180 2254 : long t = itos(ellap(E, utoipos(p)));
1181 2254 : if (l==2)
1182 : {
1183 462 : if (odd(t)) return 0;
1184 : }
1185 : else
1186 : {
1187 1792 : long D = t*t - 4*p;
1188 1792 : if (kross(D,l) == -1) return 0;
1189 : }
1190 2177 : set_avma(av);
1191 : }
1192 : }
1193 105 : return 1;
1194 : }
1195 :
1196 : static ulong
1197 56 : ellnf_goodl_l(GEN E, GEN v)
1198 : {
1199 : forprime_t T;
1200 56 : GEN nf = ellnf_get_nf(E), disc = ell_get_disc(E);
1201 56 : long i, lv = lg(v);
1202 56 : ulong w = 0UL;
1203 56 : pari_sp av = avma;
1204 56 : u_forprime_init(&T, 17UL,ULONG_MAX);
1205 1260 : for (i = 1; i <= 20; i++)
1206 : {
1207 1204 : ulong p = u_forprime_next(&T);
1208 1204 : GEN pr = idealprimedec(nf, utoipos(p));
1209 1204 : long j, k, lv = lg(v), g = lg(pr)-1;
1210 2912 : for (j=1; j<=g; j++)
1211 : {
1212 1708 : GEN prj = gel(pr, j);
1213 1708 : if (idealval(nf,disc,prj) != 0) {i--; continue;}
1214 : else
1215 : {
1216 1624 : long t = itos(ellap(E, prj));
1217 7847 : for(k = 1; k < lv; k++)
1218 : {
1219 6223 : long l = v[k];
1220 6223 : if (l==2)
1221 : {
1222 1624 : if (odd(t)) w |= 1<<(k-1);
1223 : }
1224 : else
1225 : {
1226 4599 : GEN D = subii(sqrs(t), shifti(pr_norm(prj),2));
1227 4599 : if (krois(D,l)==-1) w |= 1<<(k-1);
1228 : }
1229 : }
1230 : }
1231 : }
1232 1204 : set_avma(av);
1233 : }
1234 56 : return w^((1UL<<(lv-1))-1);
1235 : }
1236 :
1237 : static GEN
1238 1743 : ellnf_charpoly(GEN E, GEN pr)
1239 1743 : { return deg2pol_shallow(gen_1, negi(ellap(E,pr)), pr_norm(pr), 0); }
1240 :
1241 : static GEN
1242 3486 : starlaw(GEN p, GEN q)
1243 : {
1244 3486 : GEN Q = RgX_homogenize(RgX_recip(q), 1);
1245 3486 : return ZX_ZXY_resultant(p, Q);
1246 : }
1247 :
1248 : static GEN
1249 1743 : startor(GEN p, long r)
1250 : {
1251 1743 : GEN xr = pol_xn(r, 0), psir = gsub(xr, gen_1);
1252 1743 : return gsubstpol(starlaw(p, psir),xr,pol_x(0));
1253 : }
1254 :
1255 : static GEN
1256 1260 : ellnf_get_degree(GEN E, GEN p)
1257 : {
1258 1260 : GEN nf = ellnf_get_nf(E), dec = idealprimedec(nf, p);
1259 1260 : long d = nf_get_degree(nf), l = lg(dec), i, k;
1260 1260 : GEN R, starl = deg1pol_shallow(gen_1, gen_m1, 0);
1261 3003 : for (i=1; i < l; i++)
1262 : {
1263 1743 : GEN pr = gel(dec,i), q = ellnf_charpoly(E, pr);
1264 1743 : starl = starlaw(starl, startor(q, 12*pr_get_e(pr)));
1265 : }
1266 3640 : for (k = 0, R = p; 2*k <= d; k++) R = mulii(R, poleval(starl, powiu(p,12*k)));
1267 1260 : return R;
1268 : }
1269 :
1270 : /* Based on a GP script by Nicolas Billerey and Theorems 2.4 and 2.8 of
1271 : * N. Billerey, Criteres d'irreductibilite pour les representations des
1272 : * courbes elliptiques, Int. J. Number Theory 7 (2011), no. 4, 1001-1032. */
1273 : static GEN
1274 63 : ellnf_prime_degree(GEN E)
1275 : {
1276 : forprime_t T;
1277 : long i;
1278 63 : GEN nf = ellnf_get_nf(E), N = nfnorm(nf, ell_get_disc(E)), B = gen_0, P;
1279 63 : GEN bad = mulii(typ(N) == t_FRAC? mulii(gel(N,1),gel(N,2)): N,
1280 : nf_get_disc(nf));
1281 63 : u_forprime_init(&T, 5UL,ULONG_MAX);
1282 1491 : for (i = 1; i <= 20; i++)
1283 : {
1284 1428 : ulong p = u_forprime_next(&T);
1285 : GEN r;
1286 1428 : if (dvdiu(bad, p)) { i--; continue; }
1287 1260 : B = gcdii(B, ellnf_get_degree(E, utoipos(p)));
1288 1260 : if (Z_issquareall(B, &r)) B = r;
1289 : }
1290 63 : if (!signe(B)) pari_err_DOMAIN("ellisomat", "E","has",strtoGENstr("CM"),E);
1291 56 : P = vec_to_vecsmall(gel(Z_factor(B),1));
1292 56 : return shallowextract(P, utoi(ellnf_goodl_l(E, P)));
1293 : }
1294 :
1295 : static GEN
1296 462 : ellQ_isomat(GEN E, long flag)
1297 : {
1298 462 : GEN K = NULL, T2 = NULL, T3 = NULL, T5, T7, T13;
1299 : ulong good;
1300 : long n2, n3, n5, n7, n13;
1301 462 : GEN jt, jtp, s0, j = ell_get_j(E);
1302 462 : long l = ellQ_exceptional_iso(j, &jt, &jtp, &s0);
1303 462 : if (l)
1304 : {
1305 : #if 1
1306 154 : return mkisomat(l, ellisograph_dummy(E, l, jt, jtp, s0, flag));
1307 : #else
1308 : return mkisomat(l, ellisograph_p(K, E, l), flag);
1309 : #endif
1310 : }
1311 308 : good = ellQ_goodl(ellintegralmodel(E,NULL));
1312 308 : if (good & _2)
1313 : {
1314 238 : T2 = ellisograph_p(K, E, 2, flag);
1315 238 : n2 = etree_nbnodes(T2);
1316 238 : if (n2>4 || gequalgs(j, 1728) || gequalgs(j, 287496))
1317 63 : return mkisomat(2, T2);
1318 70 : } else n2 = 1;
1319 245 : if (good & _3)
1320 : {
1321 161 : T3 = ellisograph_p(K, E, 3, flag);
1322 161 : n3 = etree_nbnodes(T3);
1323 161 : if (n3>1 && n2==2) return mkisomatdbl(3,T3,2,T2, flag);
1324 56 : if (n3==2 && n2>1) return mkisomatdbl(2,T2,3,T3, flag);
1325 49 : if (n3>2 || gequal0(j)) return mkisomat(3, T3);
1326 84 : } else n3 = 1;
1327 98 : if (good & _5)
1328 : {
1329 28 : T5 = ellisograph_p(K, E, 5, flag);
1330 28 : n5 = etree_nbnodes(T5);
1331 28 : if (n5>1 && n2>1) return mkisomatdbl(2,T2,5,T5, flag);
1332 28 : if (n5>1 && n3>1) return mkisomatdbl(3,T3,5,T5, flag);
1333 14 : if (n5>1) return mkisomat(5, T5);
1334 70 : } else n5 = 1;
1335 70 : if (good & _7)
1336 : {
1337 28 : T7 = ellisograph_p(K, E, 7, flag);
1338 28 : n7 = etree_nbnodes(T7);
1339 28 : if (n7>1 && n2>1) return mkisomatdbl(2,T2,7,T7, flag);
1340 0 : if (n7>1 && n3>1) return mkisomatdbl(3,T3,7,T7, flag);
1341 0 : if (n7>1) return mkisomat(7,T7);
1342 42 : } else n7 = 1;
1343 42 : if (n2>1) return mkisomat(2,T2);
1344 7 : if (n3>1) return mkisomat(3,T3);
1345 7 : if (good & _13)
1346 : {
1347 0 : T13 = ellisograph_p(K, E, 13, flag);
1348 0 : n13 = etree_nbnodes(T13);
1349 0 : if (n13>1) return mkisomat(13,T13);
1350 7 : } else n13 = 1;
1351 7 : return mkvec2(mkvec(ellisograph_a4a6(E,flag)), matid(1));
1352 : }
1353 :
1354 : static long
1355 1134 : fill_LM(GEN LM, GEN L, GEN M, GEN z, long k)
1356 : {
1357 1134 : GEN Li = gel(LM,1), Mi1 = gmael(LM,2,1);
1358 1134 : long j, m = lg(Li);
1359 2884 : for (j = 2; j < m; j++)
1360 : {
1361 1750 : GEN d = gel(Mi1,j);
1362 1750 : gel(L, k) = gel(Li,j);
1363 1750 : gel(M, k) = z? mulii(d,z): d;
1364 1750 : k++;
1365 : }
1366 1134 : return k;
1367 : }
1368 : static GEN
1369 294 : ellnf_isocrv(GEN nf, GEN E, GEN v, GEN PE, long flag)
1370 : {
1371 : long i, l, lv, n, k;
1372 294 : GEN L, M, LE = cgetg_copy(v,&lv), e = ellisograph_a4a6(E, flag);
1373 868 : for (i = n = 1; i < lv; i++)
1374 : {
1375 574 : ulong p = uel(v,i);
1376 574 : GEN T = isograph_p(nf, e, p, gel(PE,i), flag);
1377 574 : GEN LM = nfmkisomat(nf, p, T);
1378 574 : gel(LE,i) = LM;
1379 574 : n *= lg(gel(LM,1)) - 1;
1380 : }
1381 294 : L = cgetg(n+1,t_VEC); gel(L,1) = e;
1382 294 : M = cgetg(n+1,t_COL); gel(M,1) = gen_1;
1383 868 : for (i = 1, k = 2; i < lv; i++)
1384 : {
1385 574 : ulong p = uel(v,i);
1386 574 : GEN P = gel(PE,i);
1387 574 : long kk = k;
1388 574 : k = fill_LM(gel(LE,i), L, M, NULL, k);
1389 1134 : for (l = 2; l < kk; l++)
1390 : {
1391 560 : GEN T = isograph_p(nf, gel(L,l), p, P, flag);
1392 560 : GEN LMe = nfmkisomat(nf, p, T);
1393 560 : k = fill_LM(LMe, L, M, gel(M,l), k);
1394 : }
1395 : }
1396 294 : return mkvec2(L, M);
1397 : }
1398 :
1399 : static long
1400 1911 : nfispower_quo(GEN nf, long d, GEN a, GEN b)
1401 : {
1402 1911 : if (gequal(a,b)) return 1;
1403 1505 : return nfispower(nf, nfdiv(nf, a, b), d, NULL);
1404 : }
1405 :
1406 : static long
1407 8841 : isomat_eq(GEN nf, GEN e1, GEN e2)
1408 : {
1409 8841 : if (gequal(e1,e2)) return 1;
1410 8680 : if (!gequal(gel(e1,3), gel(e2,3))) return 0;
1411 1911 : if (gequal0(gel(e1,3)))
1412 0 : return nfispower_quo(nf,6,gel(e1,2),gel(e2,2));
1413 1911 : if (gequalgs(gel(e1,3),1728))
1414 0 : return nfispower_quo(nf,4,gel(e1,1),gel(e2,1));
1415 1911 : return nfispower_quo(nf,2,gmul(gel(e1,1),gel(e2,2)),gmul(gel(e1,2),gel(e2,1)));
1416 : }
1417 :
1418 : static long
1419 1750 : isomat_find(GEN nf, GEN e, GEN L)
1420 : {
1421 1750 : long i, l = lg(L);
1422 8841 : for (i=1; i<l; i++)
1423 8841 : if (isomat_eq(nf, e, gel(L,i))) return i;
1424 : pari_err_BUG("isomat_find"); return 0; /* LCOV_EXCL_LINE */
1425 : }
1426 :
1427 : static GEN
1428 238 : isomat_perm(GEN nf, GEN E, GEN L)
1429 : {
1430 238 : long i, l = lg(E);
1431 238 : GEN v = cgetg(l, t_VECSMALL);
1432 1988 : for (i=1; i<l; i++)
1433 1750 : uel(v, i) = isomat_find(nf, gel(E,i), L);
1434 238 : return v;
1435 : }
1436 :
1437 : static GEN
1438 56 : ellnf_modpoly(GEN v, long vx, long vy)
1439 : {
1440 56 : long i, l = lg(v);
1441 56 : GEN P = cgetg(l, t_VEC);
1442 154 : for(i = 1; i < l; i++) gel(P, i) = get_polmodular(v[i],vx,vy);
1443 56 : return P;
1444 : }
1445 :
1446 : static GEN
1447 63 : ellnf_isomat(GEN E, long flag)
1448 : {
1449 63 : GEN nf = ellnf_get_nf(E);
1450 63 : GEN v = ellnf_prime_degree(E);
1451 56 : long vy = fetch_var_higher();
1452 56 : long vx = fetch_var_higher();
1453 56 : GEN P = ellnf_modpoly(v,vx,vy);
1454 56 : GEN LM = ellnf_isocrv(nf, E, v, P, flag), L = gel(LM,1), M = gel(LM,2);
1455 56 : long i, l = lg(L);
1456 56 : GEN R = cgetg(l, t_MAT);
1457 56 : gel(R,1) = M;
1458 294 : for(i = 2; i < l; i++)
1459 : {
1460 238 : GEN Li = gel(L,i);
1461 238 : GEN e = mkvec2(gdivgs(gel(Li,1), -48), gdivgs(gel(Li,2), -864));
1462 238 : GEN LMi = ellnf_isocrv(nf, ellinit(e, nf, DEFAULTPREC), v, P, 1);
1463 238 : GEN LLi = gel(LMi, 1), Mi = gel(LMi, 2);
1464 238 : GEN r = isomat_perm(nf, L, LLi);
1465 238 : gel(R,i) = vecpermute(Mi, r);
1466 : }
1467 56 : delete_var(); delete_var();
1468 56 : return mkvec2(L, R);
1469 : }
1470 :
1471 : static GEN
1472 623 : list_to_crv(GEN L)
1473 : {
1474 : long i, l;
1475 623 : GEN V = cgetg_copy(L, &l);
1476 2849 : for (i=1; i < l; i++)
1477 : {
1478 2226 : GEN Li = gel(L,i);
1479 2226 : GEN e = mkvec2(gdivgs(gel(Li,1), -48), gdivgs(gel(Li,2), -864));
1480 2226 : gel(V,i) = lg(Li)==6 ? mkvec3(e, gel(Li,4), gel(Li,5)): e;
1481 : }
1482 623 : return V;
1483 : }
1484 :
1485 : GEN
1486 721 : ellisomat(GEN E, long p, long flag)
1487 : {
1488 721 : pari_sp av = avma;
1489 721 : GEN r = NULL, nf = NULL;
1490 721 : long good = 1;
1491 721 : if (flag < 0 || flag > 1) pari_err_FLAG("ellisomat");
1492 721 : if (p < 0) pari_err_PRIME("ellisomat", stoi(p));
1493 721 : if (p == 1) { flag = 1; p = 0; } /* for backward compatibility */
1494 721 : checkell(E);
1495 721 : switch(ell_get_type(E))
1496 : {
1497 644 : case t_ELL_Q:
1498 644 : if (p) good = ellQ_goodl_l(E, p);
1499 644 : break;
1500 77 : case t_ELL_NF:
1501 77 : if (p) good = ellnf_goodl_l(E, mkvecsmall(p));
1502 77 : nf = ellnf_get_nf(E);
1503 77 : if (nf_get_varn(nf) <= 0)
1504 7 : pari_err_PRIORITY("ellisomat", nf_get_pol(nf), "<=", 0);
1505 70 : if (flag==0 && nf_get_varn(nf) <= 1)
1506 7 : pari_err_PRIORITY("ellisomat", nf_get_pol(nf), "<=", 1);
1507 63 : break;
1508 0 : default: pari_err_TYPE("ellisomat",E);
1509 : }
1510 707 : if (!good) r = mkvec2(mkvec(ellisograph_a4a6(E, flag)),matid(1));
1511 : else
1512 : {
1513 630 : if (p)
1514 105 : r = nfmkisomat(nf, p, ellisograph_p(nf, E, p, flag));
1515 : else
1516 525 : r = nf? ellnf_isomat(E, flag): ellQ_isomat(E, flag);
1517 623 : gel(r,1) = list_to_crv(gel(r,1));
1518 : }
1519 700 : return gerepilecopy(av, r);
1520 : }
1521 :
1522 : static GEN
1523 77 : get_isomat(GEN v)
1524 : {
1525 : GEN M, vE, wE;
1526 : long i, l;
1527 77 : if (typ(v) != t_VEC) return NULL;
1528 77 : if (checkell_i(v))
1529 : {
1530 35 : if (ell_get_type(v) != t_ELL_Q) return NULL;
1531 35 : v = ellisomat(v,0,1);
1532 35 : wE = gel(v,1); l = lg(wE);
1533 35 : M = gel(v,2);
1534 : }
1535 : else
1536 : {
1537 42 : if (lg(v) != 3) return NULL;
1538 42 : vE = gel(v,1); l = lg(vE);
1539 42 : M = gel(v,2);
1540 42 : if (typ(M) != t_MAT || !RgM_is_ZM(M)) return NULL;
1541 42 : if (typ(vE) != t_VEC || l == 1) return NULL;
1542 42 : if (lg(gel(vE,1)) == 3) wE = shallowcopy(vE);
1543 : else
1544 : { /* [[a4,a6],f,g] */
1545 14 : wE = cgetg_copy(vE,&l);
1546 70 : for (i = 1; i < l; i++) gel(wE,i) = gel(gel(vE,i),1);
1547 : }
1548 : }
1549 : /* wE a vector of [a4,a6] */
1550 420 : for (i = 1; i < l; i++)
1551 : {
1552 343 : GEN e = ellinit(gel(wE,i), gen_1, DEFAULTPREC);
1553 343 : GEN E = ellminimalmodel(e, NULL);
1554 343 : obj_free(e); gel(wE,i) = E;
1555 : }
1556 77 : return mkvec2(wE, M);
1557 : }
1558 :
1559 : GEN
1560 42 : ellweilcurve(GEN E, GEN *ms)
1561 : {
1562 42 : pari_sp av = avma;
1563 42 : GEN vE = get_isomat(E), vL, Wx, W, XPM, Lf, Cf;
1564 : long i, l;
1565 :
1566 42 : if (!vE) pari_err_TYPE("ellweilcurve",E);
1567 42 : vE = gel(vE,1); l = lg(vE);
1568 42 : Wx = msfromell(vE, 0);
1569 42 : W = gel(Wx,1);
1570 42 : XPM = gel(Wx,2);
1571 : /* lattice attached to the Weil curve in the isogeny class */
1572 42 : Lf = mslattice(W, gmael(XPM,1,3));
1573 42 : Cf = ginv(Lf); /* left-inverse */
1574 42 : vL = cgetg(l, t_VEC);
1575 252 : for (i=1; i < l; i++)
1576 : {
1577 210 : GEN c, Ce, Le = gmael(XPM,i,3);
1578 210 : Ce = Q_primitive_part(RgM_mul(Cf, Le), &c);
1579 210 : Ce = ZM_snf(Ce);
1580 210 : if (c) { Ce = ZC_Q_mul(Ce,c); settyp(Ce,t_VEC); }
1581 210 : gel(vL,i) = Ce;
1582 : }
1583 252 : for (i = 1; i < l; i++) obj_free(gel(vE,i));
1584 42 : vE = mkvec2(vE, vL);
1585 42 : if (!ms) return gerepilecopy(av, vE);
1586 7 : *ms = Wx; return gc_all(av, 2, &vE, ms);
1587 : }
1588 :
1589 : GEN
1590 35 : ellisotree(GEN E)
1591 : {
1592 35 : pari_sp av = avma;
1593 35 : GEN L = get_isomat(E), vE, adj, M;
1594 : long i, j, n;
1595 35 : if (!L) pari_err_TYPE("ellisotree",E);
1596 35 : vE = gel(L,1);
1597 35 : adj = gel(L,2);
1598 35 : n = lg(vE)-1; L = cgetg(n+1, t_VEC);
1599 168 : for (i = 1; i <= n; i++) gel(L,i) = ellR_area(gel(vE,i), DEFAULTPREC);
1600 35 : M = zeromatcopy(n,n);
1601 168 : for (i = 1; i <= n; i++)
1602 378 : for (j = i+1; j <= n; j++)
1603 : {
1604 245 : GEN p = gcoeff(adj,i,j);
1605 245 : if (!isprime(p)) continue;
1606 : /* L[i] / L[j] = p or 1/p; p iff E[i].lattice \subset E[j].lattice */
1607 126 : if (gcmp(gel(L,i), gel(L,j)) > 0)
1608 91 : gcoeff(M,i,j) = p;
1609 : else
1610 35 : gcoeff(M,j,i) = p;
1611 : }
1612 168 : for (i = 1; i <= n; i++) obj_free(gel(vE,i));
1613 35 : return gerepilecopy(av, mkvec2(vE,M));
1614 : }
|
