Code coverage tests

This page documents the degree to which the PARI/GP source code is tested by our public test suite, distributed with the source distribution in directory src/test/. This is measured by the gcov utility; we then process gcov output using the lcov frond-end.

We test a few variants depending on Configure flags on the pari.math.u-bordeaux.fr machine (x86_64 architecture), and agregate them in the final report:

The target is to exceed 90% coverage for all mathematical modules (given that branches depending on DEBUGLEVEL or DEBUGMEM are not covered). This script is run to produce the results below.

LCOV - code coverage report
Current view: top level - basemath - ellisog.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.2 lcov report (development 29115-f22e516b23) Lines: 978 991 98.7 %
Date: 2024-02-22 08:05:55 Functions: 85 85 100.0 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.14