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