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 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.10.0 lcov report (development 20783-cec4728) Lines: 712 736 96.7 %
Date: 2017-06-28 05:59:20 Functions: 57 57 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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : 
      16             : /* Return 1 if the point Q is a Weierstrass (2-torsion) point of the
      17             :  * curve E, return 0 otherwise */
      18             : static long
      19         903 : ellisweierstrasspoint(GEN E, GEN Q)
      20         903 : { return ell_is_inf(Q) || gequal0(ec_dmFdy_evalQ(E, Q)); }
      21             : 
      22             : 
      23             : /* Given an elliptic curve E = [a1, a2, a3, a4, a6] and t,w in the ring of
      24             :  * definition of E, return the curve
      25             :  *  E' = [a1, a2, a3, a4 - 5t, a6 - (E.b2 t + 7w)] */
      26             : static GEN
      27        1995 : make_velu_curve(GEN E, GEN t, GEN w)
      28             : {
      29        1995 :   GEN A4, A6, a1 = ell_get_a1(E), a2 = ell_get_a2(E), a3 = ell_get_a3(E);
      30        1995 :   A4 = gsub(ell_get_a4(E), gmulsg(5L, t));
      31        1995 :   A6 = gsub(ell_get_a6(E), gadd(gmul(ell_get_b2(E), t), gmulsg(7L, w)));
      32        1995 :   return mkvec5(a1,a2,a3,A4,A6);
      33             : }
      34             : 
      35             : /* If phi = (f(x)/h(x)^2, g(x,y)/h(x)^3) is an isogeny, return the
      36             :  * variables x and y in a vecsmall */
      37             : INLINE void
      38        1176 : get_isog_vars(GEN phi, long *vx, long *vy)
      39             : {
      40        1176 :   *vx = varn(gel(phi, 1));
      41        1176 :   *vy = varn(gel(phi, 2));
      42        1176 :   if (*vy == *vx) *vy = gvar2(gel(phi,2));
      43        1176 : }
      44             : 
      45             : static GEN
      46        4592 : RgX_homogenous_evalpow(GEN P, GEN A, GEN B)
      47             : {
      48        4592 :   pari_sp av = avma;
      49             :   long d, i, v;
      50             :   GEN s;
      51        4592 :   if (typ(P)!=t_POL)
      52        1435 :     return mkvec2(P, gen_1);
      53        3157 :   d = degpol(P); v = varn(A);
      54        3157 :   s = scalarpol_shallow(gel(P, d+2), v);
      55       17073 :   for (i = d-1; i >= 0; i--)
      56             :   {
      57       13916 :     s = gadd(gmul(s, A), gmul(gel(B,d+1-i), gel(P,i+2)));
      58       13916 :     if (gc_needed(av,1))
      59             :     {
      60           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"RgX_homogenous_eval(%ld)",i);
      61           0 :       s = gerepileupto(av, s);
      62             :     }
      63             :   }
      64        3157 :   s = gerepileupto(av, s);
      65        3157 :   return mkvec2(s, gel(B,d+1));
      66             : }
      67             : 
      68             : /* Given isogenies F:E' -> E and G:E'' -> E', return the composite
      69             :  * isogeny F o G:E'' -> E */
      70             : static GEN
      71        1155 : ellcompisog(GEN F, GEN G)
      72             : {
      73        1155 :   pari_sp av = avma;
      74             :   GEN Fv, Gh, Gh2, Gh3, f, g, h, h2, h3, den, num;
      75             :   GEN K, K2, K3, F0, F1, g0, g1, Gp;
      76             :   long v, vx, vy, d;
      77        1155 :   checkellisog(F);
      78        1148 :   checkellisog(G);
      79        1148 :   get_isog_vars(F, &vx, &vy);
      80        1148 :   v = fetch_var_higher();
      81        1148 :   Fv = shallowcopy(gel(F,3)); setvarn(Fv, v);
      82        1148 :   Gh = gel(G,3); Gh2 = gsqr(Gh); Gh3 = gmul(Gh, Gh2);
      83        1148 :   K = gmul(polresultant0(Fv, deg1pol(gneg(Gh2),gel(G,1), v), v, 0), Gh);
      84        1148 :   delete_var();
      85        1148 :   K = RgX_normalize(RgX_div(K, RgX_gcd(K,deriv(K,0))));
      86        1148 :   K2 = gsqr(K); K3 = gmul(K, K2);
      87        1148 :   F0 = polcoeff0(gel(F,2), 0, vy); F1 = polcoeff0(gel(F,2), 1, vy);
      88        1148 :   d = maxss(maxss(degpol(gel(F,1)),degpol(gel(F,3))),maxss(degpol(F0),degpol(F1)));
      89        1148 :   Gp = gpowers(Gh2, d);
      90        1148 :   f  = RgX_homogenous_evalpow(gel(F,1), gel(G,1), Gp);
      91        1148 :   g0 = RgX_homogenous_evalpow(F0, gel(G,1), Gp);
      92        1148 :   g1 = RgX_homogenous_evalpow(F1, gel(G,1), Gp);
      93        1148 :   h =  RgX_homogenous_evalpow(gel(F,3), gel(G,1), Gp);
      94        1148 :   h2 = mkvec2(gsqr(gel(h,1)), gsqr(gel(h,2)));
      95        1148 :   h3 = mkvec2(gmul(gel(h,1),gel(h2,1)), gmul(gel(h,2),gel(h2,2)));
      96        1148 :   f  = gdiv(gmul(gmul(K2, gel(f,1)),gel(h2,2)), gmul(gel(f,2), gel(h2,1)));
      97        1148 :   den = gmul(Gh3, gel(g1,2));
      98        1148 :   num = gadd(gmul(gel(g0,1),den), gmul(gmul(gel(G,2),gel(g1,1)),gel(g0,2)));
      99        1148 :   g = gdiv(gmul(gmul(K3,num),gel(h3,2)),gmul(gmul(gel(g0,2),den), gel(h3,1)));
     100        1148 :   return gerepilecopy(av, mkvec3(f,g,K));
     101             : }
     102             : 
     103             : /* Given an isogeny phi from ellisogeny() and a point P in the domain of phi,
     104             :  * return phi(P) */
     105             : GEN
     106        1204 : ellisogenyapply(GEN phi, GEN P)
     107             : {
     108        1204 :   pari_sp ltop = avma;
     109             :   GEN f, g, h, img_f, img_g, img_h, img_h2, img_h3, img, tmp;
     110             :   long vx, vy;
     111        1204 :   if (lg(P) == 4) return ellcompisog(phi,P);
     112          49 :   checkellisog(phi);
     113          49 :   checkellpt(P);
     114          42 :   if (ell_is_inf(P)) return ellinf();
     115          28 :   f = gel(phi, 1);
     116          28 :   g = gel(phi, 2);
     117          28 :   h = gel(phi, 3);
     118          28 :   get_isog_vars(phi, &vx, &vy);
     119          28 :   img_h = poleval(h, gel(P, 1));
     120          28 :   if (gequal0(img_h)) { avma = ltop; return ellinf(); }
     121             : 
     122          21 :   img_h2 = gsqr(img_h);
     123          21 :   img_h3 = gmul(img_h, img_h2);
     124          21 :   img_f = poleval(f, gel(P, 1));
     125             :   /* FIXME: This calculation of g is perhaps not as efficient as it could be */
     126          21 :   tmp = gsubst(g, vx, gel(P, 1));
     127          21 :   img_g = gsubst(tmp, vy, gel(P, 2));
     128          21 :   img = cgetg(3, t_VEC);
     129          21 :   gel(img, 1) = gdiv(img_f, img_h2);
     130          21 :   gel(img, 2) = gdiv(img_g, img_h3);
     131          21 :   return gerepileupto(ltop, img);
     132             : }
     133             : 
     134             : /* isog = [f, g, h] = [x, y, 1] = identity */
     135             : static GEN
     136         252 : isog_identity(long vx, long vy)
     137         252 : { return mkvec3(pol_x(vx), pol_x(vy), pol_1(vx)); }
     138             : 
     139             : /* Returns an updated value for isog based (primarily) on tQ and uQ. Used to
     140             :  * iteratively compute the isogeny corresponding to a subgroup generated by a
     141             :  * given point. Ref: Equation 8 in Velu's paper.
     142             :  * isog = NULL codes the identity */
     143             : static GEN
     144         532 : update_isogeny_polys(GEN isog, GEN E, GEN Q, GEN tQ, GEN uQ, long vx, long vy)
     145             : {
     146         532 :   pari_sp ltop = avma, av;
     147         532 :   GEN xQ = gel(Q, 1), yQ = gel(Q, 2);
     148         532 :   GEN rt = deg1pol_shallow(gen_1, gneg(xQ), vx);
     149         532 :   GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E);
     150             : 
     151         532 :   GEN gQx = ec_dFdx_evalQ(E, Q);
     152         532 :   GEN gQy = ec_dFdy_evalQ(E, Q);
     153             :   GEN tmp1, tmp2, tmp3, tmp4, f, g, h, rt_sqr, res;
     154             : 
     155             :   /* g -= uQ * (2 * y + E.a1 * x + E.a3)
     156             :    *   + tQ * rt * (E.a1 * rt + y - yQ)
     157             :    *   + rt * (E.a1 * uQ - gQx * gQy) */
     158         532 :   av = avma;
     159         532 :   tmp1 = gmul(uQ, gadd(deg1pol_shallow(gen_2, gen_0, vy),
     160             :                        deg1pol_shallow(a1, a3, vx)));
     161         532 :   tmp1 = gerepileupto(av, tmp1);
     162         532 :   av = avma;
     163         532 :   tmp2 = gmul(tQ, gadd(gmul(a1, rt),
     164             :                        deg1pol_shallow(gen_1, gneg(yQ), vy)));
     165         532 :   tmp2 = gerepileupto(av, tmp2);
     166         532 :   av = avma;
     167         532 :   tmp3 = gsub(gmul(a1, uQ), gmul(gQx, gQy));
     168         532 :   tmp3 = gerepileupto(av, tmp3);
     169             : 
     170         532 :   if (!isog) isog = isog_identity(vx,vy);
     171         532 :   f = gel(isog, 1);
     172         532 :   g = gel(isog, 2);
     173         532 :   h = gel(isog, 3);
     174         532 :   rt_sqr = gsqr(rt);
     175         532 :   res = cgetg(4, t_VEC);
     176         532 :   av = avma;
     177         532 :   tmp4 = gdiv(gadd(gmul(tQ, rt), uQ), rt_sqr);
     178         532 :   gel(res, 1) = gerepileupto(av, gadd(f, tmp4));
     179         532 :   av = avma;
     180         532 :   tmp4 = gadd(tmp1, gmul(rt, gadd(tmp2, tmp3)));
     181         532 :   gel(res, 2) = gerepileupto(av, gsub(g, gdiv(tmp4, gmul(rt, rt_sqr))));
     182         532 :   av = avma;
     183         532 :   gel(res, 3) = gerepileupto(av, gmul(h, rt));
     184         532 :   return gerepileupto(ltop, res);
     185             : }
     186             : 
     187             : /* Given a point P on E, return the curve E/<P> and, if only_image is zero,
     188             :  * the isogeny pi: E -> E/<P>. The variables vx and vy are used to describe
     189             :  * the isogeny (ignored if only_image is zero) */
     190             : static GEN
     191         427 : isogeny_from_kernel_point(GEN E, GEN P, int only_image, long vx, long vy)
     192             : {
     193         427 :   pari_sp av = avma;
     194             :   GEN isog, EE, f, g, h, h2, h3;
     195         427 :   GEN Q = P, t = gen_0, w = gen_0;
     196             :   long c;
     197         427 :   if (!oncurve(E,P))
     198           7 :     pari_err_DOMAIN("isogeny_from_kernel_point", "point", "not on", E, P);
     199         420 :   if (ell_is_inf(P))
     200             :   {
     201          42 :     if (only_image) return E;
     202          28 :     return mkvec2(E, isog_identity(vx,vy));
     203             :   }
     204             : 
     205         378 :   isog = NULL; c = 1;
     206             :   for (;;)
     207             :   {
     208         903 :     GEN tQ, xQ = gel(Q,1), uQ = ec_2divpol_evalx(E, xQ);
     209         903 :     int stop = 0;
     210         903 :     if (ellisweierstrasspoint(E,Q))
     211             :     { /* ord(P)=2c; take Q=[c]P into consideration and stop */
     212         196 :       tQ = ec_dFdx_evalQ(E, Q);
     213         196 :       stop = 1;
     214             :     }
     215             :     else
     216         707 :       tQ = ec_half_deriv_2divpol_evalx(E, xQ);
     217         903 :     t = gadd(t, tQ);
     218         903 :     w = gadd(w, gadd(uQ, gmul(tQ, xQ)));
     219         903 :     if (!only_image) isog = update_isogeny_polys(isog, E, Q,tQ,uQ, vx,vy);
     220         903 :     if (stop) break;
     221             : 
     222         707 :     Q = elladd(E, P, Q);
     223         707 :     ++c;
     224             :     /* IF x([c]P) = x([c-1]P) THEN [c]P = -[c-1]P and [2c-1]P = 0 */
     225         707 :     if (gequal(gel(Q,1), xQ)) break;
     226         525 :     if (gc_needed(av,1))
     227             :     {
     228           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"isogeny_from_kernel_point");
     229           0 :       gerepileall(av, isog? 4: 3, &Q, &t, &w, &isog);
     230             :     }
     231         525 :   }
     232             : 
     233         378 :   EE = make_velu_curve(E, t, w);
     234         378 :   if (only_image) return EE;
     235             : 
     236         224 :   if (!isog) isog = isog_identity(vx,vy);
     237         224 :   f = gel(isog, 1);
     238         224 :   g = gel(isog, 2);
     239         224 :   if ( ! (typ(f) == t_RFRAC && typ(g) == t_RFRAC))
     240           0 :     pari_err_BUG("isogeny_from_kernel_point (f or g has wrong type)");
     241             : 
     242             :   /* Clean up the isogeny polynomials f, g and h so that the isogeny
     243             :    * is given by (x,y) -> (f(x)/h(x)^2, g(x,y)/h(x)^3) */
     244         224 :   h = gel(isog, 3);
     245         224 :   h2 = gsqr(h);
     246         224 :   h3 = gmul(h, h2);
     247         224 :   f = gmul(f, h2);
     248         224 :   g = gmul(g, h3);
     249         224 :   if (typ(f) != t_POL || typ(g) != t_POL)
     250           0 :     pari_err_BUG("isogeny_from_kernel_point (wrong denominator)");
     251         224 :   return mkvec2(EE, mkvec3(f,g, gel(isog,3)));
     252             : }
     253             : 
     254             : /* Given a t_POL x^n - s1 x^{n-1} + s2 x^{n-2} - s3 x^{n-3} + ...
     255             :  * return the first three power sums (Newton's identities):
     256             :  *   p1 = s1
     257             :  *   p2 = s1^2 - 2 s2
     258             :  *   p3 = (s1^2 - 3 s2) s1 + 3 s3 */
     259             : static void
     260        1631 : first_three_power_sums(GEN pol, GEN *p1, GEN *p2, GEN *p3)
     261             : {
     262        1631 :   long d = degpol(pol);
     263             :   GEN s1, s2, ms3;
     264             : 
     265        1631 :   *p1 = s1 = gneg(RgX_coeff(pol, d-1));
     266             : 
     267        1631 :   s2 = RgX_coeff(pol, d-2);
     268        1631 :   *p2 = gsub(gsqr(s1), gmulsg(2L, s2));
     269             : 
     270        1631 :   ms3 = RgX_coeff(pol, d-3);
     271        1631 :   *p3 = gadd(gmul(s1, gsub(*p2, s2)), gmulsg(-3L, ms3));
     272        1631 : }
     273             : 
     274             : 
     275             : /* Let E and a t_POL h of degree 1 or 3 whose roots are 2-torsion points on E.
     276             :  * - if only_image != 0, return [t, w] used to compute the equation of the
     277             :  *   quotient by the given 2-torsion points
     278             :  * - else return [t,w, f,g,h], along with the contributions f, g and
     279             :  *   h to the isogeny giving the quotient by h. Variables vx and vy are used
     280             :  *   to create f, g and h, or ignored if only_image is zero */
     281             : 
     282             : /* deg h = 1; 2-torsion contribution from Weierstrass point */
     283             : static GEN
     284         896 : contrib_weierstrass_pt(GEN E, GEN h, long only_image, long vx, long vy)
     285             : {
     286         896 :   GEN p = ellbasechar(E);
     287         896 :   GEN a1 = ell_get_a1(E);
     288         896 :   GEN a3 = ell_get_a3(E);
     289         896 :   GEN x0 = gneg(constant_coeff(h)); /* h = x - x0 */
     290         896 :   GEN b = gadd(gmul(a1,x0), a3);
     291             :   GEN y0, Q, t, w, t1, t2, f, g;
     292             : 
     293         896 :   if (!equalis(p, 2L)) /* char(k) != 2 ==> y0 = -b/2 */
     294         854 :     y0 = gmul2n(gneg(b), -1);
     295             :   else
     296             :   { /* char(k) = 2 ==> y0 = sqrt(f(x0)) where E is y^2 + h(x) = f(x). */
     297          42 :     if (!gequal0(b)) pari_err_BUG("two_torsion_contrib (a1*x0+a3 != 0)");
     298          42 :     y0 = gsqrt(ec_f_evalx(E, x0), 0);
     299             :   }
     300         896 :   Q = mkvec2(x0, y0);
     301         896 :   t = ec_dFdx_evalQ(E, Q);
     302         896 :   w = gmul(x0, t);
     303         896 :   if (only_image) return mkvec2(t,w);
     304             : 
     305             :   /* Compute isogeny, f = (x - x0) * t */
     306         406 :   f = deg1pol_shallow(t, gmul(t, gneg(x0)), vx);
     307             : 
     308             :   /* g = (x - x0) * t * (a1 * (x - x0) + (y - y0)) */
     309         406 :   t1 = deg1pol_shallow(a1, gmul(a1, gneg(x0)), vx);
     310         406 :   t2 = deg1pol_shallow(gen_1, gneg(y0), vy);
     311         406 :   g = gmul(f, gadd(t1, t2));
     312         406 :   return mkvec5(t, w, f, g, h);
     313             : }
     314             : /* deg h =3; full 2-torsion contribution. NB: assume h is monic; base field
     315             :  * characteristic is odd or zero (cannot happen in char 2).*/
     316             : static GEN
     317          14 : contrib_full_tors(GEN E, GEN h, long only_image, long vx, long vy)
     318             : {
     319             :   GEN p1, p2, p3, half_b2, half_b4, t, w, f, g;
     320          14 :   first_three_power_sums(h, &p1,&p2,&p3);
     321          14 :   half_b2 = gmul2n(ell_get_b2(E), -1);
     322          14 :   half_b4 = gmul2n(ell_get_b4(E), -1);
     323             : 
     324             :   /* t = 3*(p2 + b4/2) + p1 * b2/2 */
     325          14 :   t = gadd(gmulsg(3L, gadd(p2, half_b4)), gmul(p1, half_b2));
     326             : 
     327             :   /* w = 3 * p3 + p2 * b2/2 + p1 * b4/2 */
     328          14 :   w = gadd(gmulsg(3L, p3), gadd(gmul(p2, half_b2),
     329             :                                 gmul(p1, half_b4)));
     330          14 :   if (only_image) return mkvec2(t,w);
     331             : 
     332             :   /* Compute isogeny */
     333             :   {
     334           7 :     GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E), t1, t2;
     335           7 :     GEN s1 = gneg(RgX_coeff(h, 2));
     336           7 :     GEN dh = RgX_deriv(h);
     337           7 :     GEN psi2xy = gadd(deg1pol_shallow(a1, a3, vx),
     338             :                       deg1pol_shallow(gen_2, gen_0, vy));
     339             : 
     340             :     /* f = -3 (3 x + b2/2 + s1) h + (3 x^2 + (b2/2) x + (b4/2)) h'*/
     341           7 :     t1 = RgX_mul(h, gmulsg(-3, deg1pol(stoi(3), gadd(half_b2, s1), vx)));
     342           7 :     t2 = mkpoln(3, stoi(3), half_b2, half_b4);
     343           7 :     setvarn(t2, vx);
     344           7 :     t2 = RgX_mul(dh, t2);
     345           7 :     f = RgX_add(t1, t2);
     346             : 
     347             :     /* 2g = psi2xy * (f'*h - f*h') - (a1*f + a3*h) * h; */
     348           7 :     t1 = RgX_sub(RgX_mul(RgX_deriv(f), h), RgX_mul(f, dh));
     349           7 :     t2 = RgX_mul(h, RgX_add(RgX_Rg_mul(f, a1), RgX_Rg_mul(h, a3)));
     350           7 :     g = RgX_divs(gsub(gmul(psi2xy, t1), t2), 2L);
     351             : 
     352           7 :     f = RgX_mul(f, h);
     353           7 :     g = RgX_mul(g, h);
     354             :   }
     355           7 :   return mkvec5(t, w, f, g, h);
     356             : }
     357             : 
     358             : /* Given E and a t_POL T whose roots define a subgroup G of E, return the factor
     359             :  * of T that corresponds to the 2-torsion points E[2] \cap G in G */
     360             : INLINE GEN
     361        1624 : two_torsion_part(GEN E, GEN T)
     362        1624 : { return RgX_gcd(T, elldivpol(E, 2, varn(T))); }
     363             : 
     364             : /* Return the jth Hasse derivative of the polynomial f = \sum_{i=0}^n a_i x^i,
     365             :  * i.e. \sum_{i=j}^n a_i \binom{i}{j} x^{i-j}. It is a derivation even when the
     366             :  * coefficient ring has positive characteristic */
     367             : static GEN
     368          98 : derivhasse(GEN f, ulong j)
     369             : {
     370          98 :   ulong i, d = degpol(f);
     371             :   GEN df;
     372          98 :   if (gequal0(f) || d == 0) return pol_0(varn(f));
     373          56 :   if (j == 0) return gcopy(f);
     374          56 :   df = cgetg(2 + (d-j+1), t_POL);
     375          56 :   df[1] = f[1];
     376          56 :   for (i = j; i <= d; ++i) gel(df, i-j+2) = gmul(binomialuu(i,j), gel(f, i+2));
     377          56 :   return normalizepol(df);
     378             : }
     379             : 
     380             : static GEN
     381         574 : non_two_torsion_abscissa(GEN E, GEN h0, GEN x)
     382             : {
     383             :   GEN mp1, dh0, ddh0, t, u, t1, t2, t3;
     384         574 :   long m = degpol(h0);
     385         574 :   mp1 = gel(h0, m + 1); /* negative of first power sum */
     386         574 :   dh0 = RgX_deriv(h0);
     387         574 :   ddh0 = RgX_deriv(dh0);
     388         574 :   t = ec_2divpol_evalx(E, x);
     389         574 :   u = ec_half_deriv_2divpol_evalx(E, x);
     390         574 :   t1 = RgX_sub(RgX_sqr(dh0), RgX_mul(ddh0, h0));
     391         574 :   t2 = RgX_mul(u, RgX_mul(h0, dh0));
     392         574 :   t3 = RgX_mul(RgX_sqr(h0),
     393         574 :                deg1pol_shallow(stoi(2*m), gmulsg(2L, mp1), varn(x)));
     394             :   /* t * (dh0^2 - ddh0*h0) - u*dh0*h0 + (2*m*x - 2*s1) * h0^2); */
     395         574 :   return RgX_add(RgX_sub(RgX_mul(t, t1), t2), t3);
     396             : }
     397             : 
     398             : static GEN
     399         952 : isog_abscissa(GEN E, GEN kerp, GEN h0, GEN x, GEN two_tors)
     400             : {
     401             :   GEN f0, f2, h2, t1, t2, t3;
     402         952 :   f0 = (degpol(h0) > 0)? non_two_torsion_abscissa(E, h0, x): pol_0(varn(x));
     403         952 :   f2 = gel(two_tors, 3);
     404         952 :   h2 = gel(two_tors, 5);
     405             : 
     406             :   /* Combine f0 and f2 into the final abscissa of the isogeny. */
     407         952 :   t1 = RgX_mul(x, RgX_sqr(kerp));
     408         952 :   t2 = RgX_mul(f2, RgX_sqr(h0));
     409         952 :   t3 = RgX_mul(f0, RgX_sqr(h2));
     410             :   /* x * kerp^2 + f2 * h0^2 + f0 * h2^2 */
     411         952 :   return RgX_add(t1, RgX_add(t2, t3));
     412             : }
     413             : 
     414             : static GEN
     415         903 : non_two_torsion_ordinate_char_not2(GEN E, GEN f, GEN h, GEN psi2)
     416             : {
     417         903 :   GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E);
     418         903 :   GEN df = RgX_deriv(f), dh = RgX_deriv(h);
     419             :   /* g = df * h * psi2/2 - f * dh * psi2
     420             :    *   - (E.a1 * f + E.a3 * h^2) * h/2 */
     421         903 :   GEN t1 = RgX_mul(df, RgX_mul(h, RgX_divs(psi2, 2L)));
     422         903 :   GEN t2 = RgX_mul(f, RgX_mul(dh, psi2));
     423         903 :   GEN t3 = RgX_mul(RgX_divs(h, 2L),
     424             :                    RgX_add(RgX_Rg_mul(f, a1), RgX_Rg_mul(RgX_sqr(h), a3)));
     425         903 :   return RgX_sub(RgX_sub(t1, t2), t3);
     426             : }
     427             : 
     428             : /* h = kerq */
     429             : static GEN
     430          49 : non_two_torsion_ordinate_char2(GEN E, GEN h, GEN x, GEN y)
     431             : {
     432          49 :   GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E), a4 = ell_get_a4(E);
     433          49 :   GEN b2 = ell_get_b2(E), b4 = ell_get_b4(E), b6 = ell_get_b6(E);
     434             :   GEN h2, dh, dh2, ddh, D2h, D2dh, H, psi2, u, t, alpha;
     435             :   GEN p1, t1, t2, t3, t4;
     436          49 :   long m, vx = varn(x);
     437             : 
     438          49 :   h2 = RgX_sqr(h);
     439          49 :   dh = RgX_deriv(h);
     440          49 :   dh2 = RgX_sqr(dh);
     441          49 :   ddh = RgX_deriv(dh);
     442          49 :   H = RgX_sub(dh2, RgX_mul(h, ddh));
     443          49 :   D2h = derivhasse(h, 2);
     444          49 :   D2dh = derivhasse(dh, 2);
     445          49 :   psi2 = deg1pol_shallow(a1, a3, vx);
     446          49 :   u = mkpoln(3, b2, gen_0, b6);
     447          49 :   setvarn(u, vx);
     448          49 :   t = deg1pol_shallow(b2, b4, vx);
     449          49 :   alpha = mkpoln(4, a1, a3, gmul(a1, a4), gmul(a3, a4));
     450          49 :   setvarn(alpha, vx);
     451          49 :   m = degpol(h);
     452          49 :   p1 = RgX_coeff(h, m-1); /* first power sum */
     453             : 
     454          49 :   t1 = gmul(gadd(gmul(a1, p1), gmulgs(a3, m)), RgX_mul(h,h2));
     455             : 
     456          49 :   t2 = gmul(a1, gadd(gmul(a1, gadd(y, psi2)), RgX_add(RgX_Rg_add(RgX_sqr(x), a4), t)));
     457          49 :   t2 = gmul(t2, gmul(dh, h2));
     458             : 
     459          49 :   t3 = gadd(gmul(y, t), RgX_add(alpha, RgX_Rg_mul(u, a1)));
     460          49 :   t3 = gmul(t3, RgX_mul(h, H));
     461             : 
     462          49 :   t4 = gmul(u, psi2);
     463          49 :   t4 = gmul(t4, RgX_sub(RgX_sub(RgX_mul(h2, D2dh), RgX_mul(dh, H)),
     464             :                         RgX_mul(h, RgX_mul(dh, D2h))));
     465             : 
     466          49 :   return gadd(t1, gadd(t2, gadd(t3, t4)));
     467             : }
     468             : 
     469             : static GEN
     470         952 : isog_ordinate(GEN E, GEN kerp, GEN kerq, GEN x, GEN y, GEN two_tors, GEN f)
     471             : {
     472             :   GEN g;
     473         952 :   if (! equalis(ellbasechar(E), 2L)) {
     474             :     /* FIXME: We don't use (hence don't need to calculate)
     475             :      * g2 = gel(two_tors, 4) when char(k) != 2. */
     476         903 :     GEN psi2 = ec_dmFdy_evalQ(E, mkvec2(x, y));
     477         903 :     g = non_two_torsion_ordinate_char_not2(E, f, kerp, psi2);
     478             :   } else {
     479          49 :     GEN h2 = gel(two_tors, 5);
     480          49 :     GEN g2 = gmul(gel(two_tors, 4), RgX_mul(kerq, RgX_sqr(kerq)));
     481          49 :     GEN g0 = non_two_torsion_ordinate_char2(E, kerq, x, y);
     482          49 :     g0 = gmul(g0, RgX_mul(h2, RgX_sqr(h2)));
     483          49 :     g = gsub(gmul(y, RgX_mul(kerp, RgX_sqr(kerp))), gadd(g2, g0));
     484             :   }
     485         952 :   return g;
     486             : }
     487             : 
     488             : /* Given an elliptic curve E and a polynomial kerp whose roots give the
     489             :  * x-coordinates of a subgroup G of E, return the curve E/G and,
     490             :  * if only_image is zero, the isogeny pi:E -> E/G. Variables vx and vy are
     491             :  * used to describe the isogeny (and are ignored if only_image is zero). */
     492             : static GEN
     493        1624 : isogeny_from_kernel_poly(GEN E, GEN kerp, long only_image, long vx, long vy)
     494             : {
     495             :   long m;
     496        1624 :   GEN b2 = ell_get_b2(E), b4 = ell_get_b4(E), b6 = ell_get_b6(E);
     497             :   GEN p1, p2, p3, x, y, f, g, two_tors, EE, t, w;
     498        1624 :   GEN kerh = two_torsion_part(E, kerp);
     499        1624 :   GEN kerq = RgX_divrem(kerp, kerh, ONLY_DIVIDES);
     500        1624 :   if (!kerq) pari_err_BUG("isogeny_from_kernel_poly");
     501             :   /* isogeny degree: 2*degpol(kerp)+1-degpol(kerh) */
     502        1624 :   m = degpol(kerq);
     503             : 
     504        1624 :   kerp = RgX_Rg_div(kerp, leading_coeff(kerp));
     505        1624 :   kerq = RgX_Rg_div(kerq, leading_coeff(kerq));
     506        1624 :   kerh = RgX_Rg_div(kerh, leading_coeff(kerh));
     507        1624 :   switch(degpol(kerh))
     508             :   {
     509             :   case 0:
     510        1246 :     two_tors = only_image? mkvec2(gen_0, gen_0):
     511         539 :       mkvec5(gen_0, gen_0, pol_0(vx), pol_0(vx), pol_1(vx));
     512         707 :     break;
     513             :   case 1:
     514         896 :     two_tors = contrib_weierstrass_pt(E, kerh, only_image,vx,vy);
     515         896 :     break;
     516             :   case 3:
     517          14 :     two_tors = contrib_full_tors(E, kerh, only_image,vx,vy);
     518          14 :     break;
     519             :   default:
     520           7 :     two_tors = NULL;
     521           7 :     pari_err_DOMAIN("isogeny_from_kernel_poly", "kernel polynomial",
     522             :                     "does not define a subgroup of", E, kerp);
     523             :   }
     524        1617 :   first_three_power_sums(kerq,&p1,&p2,&p3);
     525        1617 :   x = pol_x(vx);
     526        1617 :   y = pol_x(vy);
     527             : 
     528             :   /* t = 6 * p2 + b2 * p1 + m * b4, */
     529        1617 :   t = gadd(gmulsg(6L, p2), gadd(gmul(b2, p1), gmulsg(m, b4)));
     530             : 
     531             :   /* w = 10 * p3 + 2 * b2 * p2 + 3 * b4 * p1 + m * b6, */
     532        1617 :   w = gadd(gmulsg(10L, p3),
     533             :            gadd(gmul(gmulsg(2L, b2), p2),
     534             :                 gadd(gmul(gmulsg(3L, b4), p1), gmulsg(m, b6))));
     535             : 
     536        1617 :   EE = make_velu_curve(E, gadd(t, gel(two_tors, 1)),
     537        1617 :                           gadd(w, gel(two_tors, 2)));
     538        1617 :   if (only_image) return EE;
     539             : 
     540         952 :   f = isog_abscissa(E, kerp, kerq, x, two_tors);
     541         952 :   g = isog_ordinate(E, kerp, kerq, x, y, two_tors, f);
     542         952 :   return mkvec2(EE, mkvec3(f,g,kerp));
     543             : }
     544             : 
     545             : /* Given an elliptic curve E and a subgroup G of E, return the curve
     546             :  * E/G and, if only_image is zero, the isogeny corresponding
     547             :  * to the canonical surjection pi:E -> E/G. The variables vx and
     548             :  * vy are used to describe the isogeny (and are ignored if
     549             :  * only_image is zero). The subgroup G may be given either as
     550             :  * a generating point P on E or as a polynomial kerp whose roots are
     551             :  * the x-coordinates of the points in G */
     552             : GEN
     553        1092 : ellisogeny(GEN E, GEN G, long only_image, long vx, long vy)
     554             : {
     555        1092 :   pari_sp av = avma;
     556             :   GEN j, z;
     557        1092 :   checkell(E);j = ell_get_j(E);
     558        1092 :   if (vx < 0) vx = 0;
     559        1092 :   if (vy < 0) vy = 1;
     560        1092 :   if (varncmp(vx, vy) >= 0)
     561           7 :     pari_err_PRIORITY("ellisogeny", pol_x(vx), "<=", vy);
     562        1085 :   if (!only_image && varncmp(vy, gvar(j)) >= 0)
     563           7 :     pari_err_PRIORITY("ellisogeny", j, ">=", vy);
     564        1078 :   switch(typ(G))
     565             :   {
     566             :   case t_VEC:
     567         441 :     checkellpt(G);
     568         441 :     if (!ell_is_inf(G))
     569             :     {
     570         399 :       GEN x =  gel(G,1), y = gel(G,2);
     571         399 :       if (!only_image)
     572             :       {
     573         245 :         if (varncmp(vy, gvar(x)) >= 0)
     574           7 :           pari_err_PRIORITY("ellisogeny", x, ">=", vy);
     575         238 :         if (varncmp(vy, gvar(y)) >= 0)
     576           7 :           pari_err_PRIORITY("ellisogeny", y, ">=", vy);
     577             :       }
     578             :     }
     579         427 :     z = isogeny_from_kernel_point(E, G, only_image, vx, vy);
     580         420 :     break;
     581             :   case t_POL:
     582         630 :     if (!only_image && varncmp(vy, gvar(constant_coeff(G))) >= 0)
     583           7 :       pari_err_PRIORITY("ellisogeny", constant_coeff(G), ">=", vy);
     584         623 :     z = isogeny_from_kernel_poly(E, G, only_image, vx, vy);
     585         616 :     break;
     586             :   default:
     587           7 :     z = NULL;
     588           7 :     pari_err_TYPE("ellisogeny", G);
     589             :   }
     590        1036 :   return gerepilecopy(av, z);
     591             : }
     592             : 
     593             : static GEN
     594         714 : trivial_isogeny(void)
     595             : {
     596         714 :   return mkvec3(pol_x(0), scalarpol(pol_x(1), 0), pol_1(0));
     597             : }
     598             : 
     599             : static GEN
     600         224 : isogeny_a4a6(GEN E)
     601             : {
     602         224 :   GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E), b2 = ell_get_b2(E);
     603         224 :   retmkvec3(deg1pol(gen_1, gdivgs(b2, 12), 0),
     604             :             deg1pol(gdivgs(a1,2), deg1pol(gen_1, gdivgs(a3,2), 1), 0),
     605             :             pol_1(0));
     606             : }
     607             : 
     608             : static GEN
     609         224 : invisogeny_a4a6(GEN E)
     610             : {
     611         224 :   GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E), b2 = ell_get_b2(E);
     612         224 :   retmkvec3(deg1pol(gen_1, gdivgs(b2, -12), 0),
     613             :             deg1pol(gdivgs(a1,-2),
     614             :               deg1pol(gen_1, gadd(gdivgs(a3,-2), gdivgs(gmul(b2,a1), 24)), 1), 0),
     615             :             pol_1(0));
     616             : }
     617             : 
     618             : static GEN
     619         560 : RgXY_eval(GEN P, GEN x, GEN y)
     620             : {
     621         560 :   return poleval(poleval(P,x), y);
     622             : }
     623             : 
     624             : static GEN
     625         385 : twistisogeny(GEN iso, GEN d)
     626             : {
     627         385 :   GEN d2 = gsqr(d), d3 = gmul(d, d2);
     628         385 :   return mkvec3(gdiv(gel(iso,1), d2), gdiv(gel(iso,2), d3), gel(iso, 3));
     629             : }
     630             : 
     631             : static GEN
     632         616 : ellisog_by_Kohel(GEN a4, GEN a6, long n, GEN ker, GEN kert, long flag)
     633             : {
     634         616 :   GEN E = ellinit(mkvec2(a4, a6), NULL, DEFAULTPREC);
     635         616 :   GEN F = isogeny_from_kernel_poly(E, ker, flag, 0, 1);
     636         616 :   GEN Et = ellinit(flag ? F: gel(F, 1), NULL, DEFAULTPREC);
     637         616 :   GEN c4t = ell_get_c4(Et), c6t = ell_get_c6(Et), jt = ell_get_j(Et);
     638         616 :   if (!flag)
     639             :   {
     640         385 :     GEN Ft = isogeny_from_kernel_poly(Et, kert, flag, 0, 1);
     641         385 :     GEN isot = twistisogeny(gel(Ft, 2), stoi(n));
     642         385 :     return mkvec5(c4t, c6t, jt, gel(F, 2), isot);
     643             :   }
     644         231 :   else return mkvec3(c4t, c6t, jt);
     645             : }
     646             : 
     647             : static GEN
     648         483 : ellisog_by_roots(GEN a4, GEN a6, long n, GEN z, long flag)
     649             : {
     650         483 :   return ellisog_by_Kohel(a4, a6, n, deg1pol(gen_1, gneg(z), 0),
     651             :                                   deg1pol(gen_1, gmulsg(n, z), 0), flag);
     652             : }
     653             : 
     654             : static GEN
     655         707 : a4a6_divpol(GEN a4, GEN a6, long n)
     656             : {
     657         707 :   switch(n)
     658             :   {
     659             :     case 2:
     660         462 :       return mkpoln(4, gen_1, gen_0, a4, a6);
     661             :     case 3:
     662         245 :       return mkpoln(5, utoi(3), gen_0, gmulgs(a4,6) , gmulgs(a6,12),
     663             :                        gneg(gsqr(a4)));
     664             :   }
     665           0 :   return NULL;
     666             : }
     667             : 
     668             : static GEN
     669         707 : ellisograph_Kohel_iso(GEN nf, GEN e, long n, GEN z, long flag)
     670             : {
     671             :   long i, r;
     672             :   GEN R, V;
     673         707 :   GEN c4 = gel(e,1), c6 = gel(e, 2);
     674         707 :   GEN a4 = gdivgs(c4, -48), a6 = gdivgs(c6, -864);
     675         707 :   GEN P = a4a6_divpol(a4, a6, n);
     676         707 :   R = nfroots(nf, z ? RgX_div_by_X_x(P, z, NULL): P);
     677         707 :   r = lg(R);
     678         707 :   V = cgetg(r, t_VEC);
     679        1190 :   for (i=1; i < r; i++)
     680         483 :     gel(V, i) = ellisog_by_roots(a4, a6, n, gel(R, i), flag);
     681         707 :   return mkvec2(V, R);
     682             : }
     683             : 
     684             : static GEN
     685         651 : ellisograph_Kohel_r(GEN nf, GEN e, long n, GEN z, long flag)
     686             : {
     687         651 :   GEN W = ellisograph_Kohel_iso(nf, e, n, z, flag);
     688         651 :   GEN iso = gel(W, 1), R = gel(W, 2);
     689         651 :   long i, r = lg(iso);
     690         651 :   GEN V = cgetg(r, t_VEC);
     691        1078 :   for (i=1; i < r; i++)
     692         427 :     gel(V, i) = ellisograph_Kohel_r(nf, gel(iso, i), n, gmulgs(gel(R, i), -n), flag);
     693         651 :   return mkvec2(e, V);
     694             : }
     695             : 
     696             : static GEN
     697         266 : corr(GEN c4, GEN c6)
     698             : {
     699         266 :   GEN c62 = gmul2n(c6, 1);
     700         266 :   return gadd(gdiv(gsqr(c4), c62), gdiv(c62, gmulgs(c4,3)));
     701             : }
     702             : 
     703             : static GEN
     704         266 : elkies98(GEN a4, GEN a6, long l, GEN s, GEN a4t, GEN a6t)
     705             : {
     706             :   GEN C, P, S;
     707             :   long i, n, d;
     708         266 :   d = l == 2 ? 1 : l>>1;
     709         266 :   C = cgetg(d+1, t_VEC);
     710         266 :   gel(C, 1) = gdivgs(gsub(a4, a4t), 5);
     711         266 :   if (d >= 2)
     712         266 :     gel(C, 2) = gdivgs(gsub(a6, a6t), 7);
     713         266 :   if (d >= 3)
     714         210 :     gel(C, 3) = gdivgs(gsub(gsqr(gel(C, 1)), gmul(a4, gel(C, 1))), 3);
     715        2758 :   for (n = 3; n < d; ++n)
     716             :   {
     717        2492 :     GEN s = gen_0;
     718       61222 :     for (i = 1; i < n; i++)
     719       58730 :       s = gadd(s, gmul(gel(C, i), gel(C, n-i)));
     720        2492 :     gel(C, n+1) = gdivgs(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));
     721             :   }
     722         266 :   P = cgetg(d+2, t_VEC);
     723         266 :   gel(P, 1 + 0) = stoi(d);
     724         266 :   gel(P, 1 + 1) = s;
     725         266 :   if (d >= 2)
     726         266 :     gel(P, 1 + 2) = gdivgs(gsub(gel(C, 1), gmulgs(gmulsg(2, a4), d)), 6);
     727        2968 :   for (n = 2; n < d; ++n)
     728        2702 :     gel(P, 1 + n+1) = gdivgs(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);
     729         266 :   S = cgetg(d+3, t_POL);
     730         266 :   S[1] = evalsigne(1) | evalvarn(0);
     731         266 :   gel(S, 2 + d - 0) = gen_1;
     732         266 :   gel(S, 2 + d - 1) = gneg(s);
     733        3234 :   for (n = 2; n <= d; ++n)
     734             :   {
     735        2968 :     GEN s = gen_0;
     736       67844 :     for (i = 1; i <= n; ++i)
     737             :     {
     738       64876 :       GEN p = gmul(gel(P, 1+i), gel(S, 2 + d - (n-i)));
     739       64876 :       s = gadd(s, p);
     740             :     }
     741        2968 :     gel(S, 2 + d - n) = gdivgs(s, -n);
     742             :   }
     743         266 :   return S;
     744             : }
     745             : 
     746             : static GEN
     747         266 : ellisog_by_jt(GEN c4, GEN c6, GEN jt, GEN jtp, GEN s0, long n, long flag)
     748             : {
     749         266 :   GEN jtp2 = gsqr(jtp), den = gmul(jt, gsubgs(jt, 1728));
     750         266 :   GEN c4t = gdiv(jtp2, den);
     751         266 :   GEN c6t = gdiv(gmul(jtp, c4t), jt);
     752         266 :   if (flag)
     753         133 :     return mkvec3(c4t, c6t, jt);
     754             :   else
     755             :   {
     756         133 :     GEN co  = corr(c4, c6);
     757         133 :     GEN cot = corr(c4t, c6t);
     758         133 :     GEN s = gmul2n(gmulgs(gadd(gadd(s0, co), gmulgs(cot,-n)), -n), -2);
     759         133 :     GEN a4  = gdivgs(c4, -48), a6 = gdivgs(c6, -864);
     760         133 :     GEN a4t = gmul(gdivgs(c4t, -48), powuu(n,4)), a6t = gmul(gdivgs(c6t, -864), powuu(n,6));
     761         133 :     GEN ker = elkies98(a4, a6, n, s, a4t, a6t);
     762         133 :     GEN st = gmulgs(s, -n);
     763         133 :     GEN a4tt = gmul(a4,powuu(n,4)), a6tt = gmul(a6,powuu(n,6));
     764         133 :     GEN kert = elkies98(a4t, a6t, n, st, a4tt, a6tt);
     765         133 :     return ellisog_by_Kohel(a4, a6, n, ker, kert, flag);
     766             :   }
     767             : }
     768             : 
     769             : /*
     770             : Based on
     771             : RENE SCHOOF
     772             : Counting points on elliptic curves over finite fields
     773             : Journal de Theorie des Nombres de Bordeaux,
     774             : tome 7, no 1 (1995), p. 219-254.
     775             : <http://www.numdam.org/item?id=JTNB_1995__7_1_219_0>
     776             : */
     777             : 
     778             : static GEN
     779         112 : ellisog_by_j(GEN e, GEN jt, long n, GEN P, long flag)
     780             : {
     781         112 :   pari_sp av = avma;
     782         112 :   GEN c4  = gel(e,1), c6 = gel(e, 2), j = gel(e, 3);
     783         112 :   GEN Px = deriv(P, 0), Py = deriv(P, 1);
     784         112 :   GEN Pxj = RgXY_eval(Px, j, jt), Pyj = RgXY_eval(Py, j, jt);
     785         112 :   GEN Pxx  = deriv(Px, 0), Pxy = deriv(Py, 0), Pyy = deriv(Py, 1);
     786         112 :   GEN Pxxj = RgXY_eval(Pxx,j,jt);
     787         112 :   GEN Pxyj = RgXY_eval(Pxy,j,jt);
     788         112 :   GEN Pyyj = RgXY_eval(Pyy,j,jt);
     789         112 :   GEN c6c4 = gdiv(c6, c4);
     790         112 :   GEN jp = gmul(j, c6c4);
     791         112 :   GEN jtp = gdivgs(gmul(jp, gdiv(Pxj, Pyj)), -n);
     792         112 :   GEN jtpn = gmulgs(jtp, n);
     793         112 :   GEN s0 = gdiv(gadd(gadd(gmul(gsqr(jp),Pxxj),gmul(gmul(jp,jtpn),gmul2n(Pxyj,1))),
     794             :                 gmul(gsqr(jtpn),Pyyj)),gmul(jp,Pxj));
     795         112 :   GEN et = ellisog_by_jt(c4, c6, jt, jtp, s0, n, flag);
     796         112 :   return gerepilecopy(av, et);
     797             : }
     798             : 
     799             : static GEN
     800         175 : ellisograph_iso(GEN nf, GEN e, ulong p, GEN P, GEN oj, long flag)
     801             : {
     802             :   long i, r;
     803             :   GEN Pj, R, V;
     804         175 :   GEN j = gel(e, 3);
     805         175 :   Pj = poleval(P, j);
     806         175 :   R = nfroots(nf,oj ? RgX_div_by_X_x(Pj, oj, NULL):Pj);
     807         175 :   r = lg(R);
     808         175 :   V = cgetg(r, t_VEC);
     809         287 :   for (i=1; i < r; i++)
     810         112 :     gel(V, i) = ellisog_by_j(e, gel(R, i), p, P, flag);
     811         175 :   return V;
     812             : }
     813             : 
     814             : static GEN
     815         133 : ellisograph_r(GEN nf, GEN e, ulong p, GEN P, GEN oj, long flag)
     816             : {
     817         133 :   GEN iso = ellisograph_iso(nf, e, p, P, oj, flag);
     818         133 :   GEN j = gel(e, 3);
     819         133 :   long i, r = lg(iso);
     820         133 :   GEN V = cgetg(r, t_VEC);
     821         203 :   for (i=1; i < r; i++)
     822          70 :     gel(V, i) = ellisograph_r(nf, gel(iso, i), p, P, j, flag);
     823         133 :   return mkvec2(e, V);
     824             : }
     825             : 
     826             : static GEN
     827         448 : ellisograph_a4a6(GEN E, long flag)
     828             : {
     829         448 :   GEN c4 = ell_get_c4(E), c6 = ell_get_c6(E), j = ell_get_j(E);
     830         672 :   return flag ? mkvec3(c4, c6, j):
     831         224 :                 mkvec5(c4, c6, j, isogeny_a4a6(E), invisogeny_a4a6(E));
     832             : }
     833             : 
     834             : static GEN
     835         154 : ellisograph_dummy(GEN E, long n, GEN jt, GEN jtt, GEN s0, long flag)
     836             : {
     837         154 :   GEN c4 = ell_get_c4(E), c6 = ell_get_c6(E);
     838         154 :   GEN c6c4 = gdiv(c6, c4);
     839         154 :   GEN jtp = gmul(c6c4, gdivgs(gmul(jt, jtt), -n));
     840         154 :   GEN iso = ellisog_by_jt(c4, c6, jt, jtp, gmul(s0, c6c4), n, flag);
     841         154 :   GEN v = mkvec2(iso, cgetg(1, t_VEC));
     842         154 :   return mkvec2(ellisograph_a4a6(E, flag), mkvec(v));
     843             : }
     844             : 
     845             : static GEN
     846         287 : ellisograph_p(GEN nf, GEN E, ulong p, long flag)
     847             : {
     848         287 :   pari_sp av = avma;
     849         287 :   GEN iso, e = ellisograph_a4a6(E, flag);
     850         287 :   if (p > 3)
     851             :   {
     852          63 :     GEN P = polmodular_ZXX(p, 0, 0, 1);
     853          63 :     iso = ellisograph_r(nf, e, p, P, NULL, flag);
     854             :   }
     855             :   else
     856         224 :     iso = ellisograph_Kohel_r(nf, e, p, NULL, flag);
     857         287 :   return gerepilecopy(av, iso);
     858             : }
     859             : 
     860             : static long
     861        2569 : etree_nbnodes(GEN T)
     862             : {
     863        2569 :   GEN F = gel(T,2);
     864        2569 :   long n = 1, i, l = lg(F);
     865        4165 :   for (i = 1; i < l; i++)
     866        1596 :     n += etree_nbnodes(gel(F, i));
     867        2569 :   return n;
     868             : }
     869             : 
     870             : static long
     871         924 : etree_listr(GEN T, GEN V, long n, GEN u, GEN ut)
     872             : {
     873         924 :   GEN E = gel(T, 1), F = gel(T,2);
     874         924 :   long i, l = lg(F);
     875         924 :   GEN iso, isot = NULL;
     876         924 :   if (lg(E) == 6)
     877             :   {
     878         476 :     iso  = ellisogenyapply(gel(E,4), u);
     879         476 :     isot = ellisogenyapply(ut, gel(E,5));
     880         476 :     gel(V, n) = mkvec5(gel(E,1), gel(E,2), gel(E,3), iso, isot);
     881             :   } else
     882             :   {
     883         448 :     gel(V, n) = mkvec3(gel(E,1), gel(E,2), gel(E,3));
     884         448 :     iso = u;
     885             :   }
     886        1491 :   for (i = 1; i < l; i++)
     887         567 :     n = etree_listr(gel(F, i), V, n + 1, iso, isot);
     888         924 :   return n;
     889             : }
     890             : 
     891             : static GEN
     892         357 : etree_list(GEN T)
     893             : {
     894         357 :   long n = etree_nbnodes(T);
     895         357 :   GEN V = cgetg(n+1, t_VEC);
     896         357 :   (void) etree_listr(T, V, 1, trivial_isogeny(), trivial_isogeny());
     897         357 :   return V;
     898             : }
     899             : 
     900             : static long
     901         924 : etree_distmatr(GEN T, GEN M, long n)
     902             : {
     903         924 :   GEN F = gel(T,2);
     904         924 :   long i, j, lF = lg(F), m = n + 1;
     905         924 :   GEN V = cgetg(lF, t_VECSMALL);
     906         924 :   mael(M, n, n) = 0;
     907        1491 :   for(i = 1; i < lF; i++)
     908         567 :     V[i] = m = etree_distmatr(gel(F,i), M, m);
     909        1491 :   for(i = 1; i < lF; i++)
     910             :   {
     911         567 :     long mi = i==1 ? n+1: V[i-1];
     912        1288 :     for(j = mi; j < V[i]; j++)
     913             :     {
     914         721 :       mael(M,n,j) = 1 + mael(M, mi, j);
     915         721 :       mael(M,j,n) = 1 + mael(M, j, mi);
     916             :     }
     917        1414 :     for(j = 1; j < lF; j++)
     918         847 :       if (i != j)
     919             :       {
     920         280 :         long i1, j1, mj = j==1 ? n+1: V[j-1];
     921         672 :         for (i1 = mi; i1 < V[i]; i1++)
     922        1008 :           for(j1 = mj; j1 < V[j]; j1++)
     923         616 :             mael(M,i1,j1) = 2 + mael(M,mj,j1) + mael(M,i1,mi);
     924             :       }
     925             :   }
     926         924 :   return m;
     927             : }
     928             : 
     929             : static GEN
     930         357 : etree_distmat(GEN T)
     931             : {
     932         357 :   long i, n = etree_nbnodes(T);
     933         357 :   GEN M = cgetg(n+1, t_MAT);
     934        1281 :   for(i = 1; i <= n; i++)
     935         924 :     gel(M,i) = cgetg(n+1, t_VECSMALL);
     936         357 :   (void)etree_distmatr(T, M, 1);
     937         357 :   return M;
     938             : }
     939             : 
     940             : static GEN
     941         357 : list_to_crv(GEN L)
     942             : {
     943             :   long i, l;
     944         357 :   GEN V = cgetg_copy(L, &l);
     945        1463 :   for(i=1; i < l; i++)
     946             :   {
     947        1106 :     GEN Li = gel(L, i);
     948        1106 :     GEN e = mkvec2(gdivgs(gel(Li,1), -48), gdivgs(gel(Li,2), -864));
     949        1106 :     gel(V, i) = lg(Li)==6 ? mkvec3(e, gel(Li,4), gel(Li,5)): e;
     950             :   }
     951         357 :   return V;
     952             : }
     953             : 
     954             : static GEN
     955         357 : distmat_pow(GEN E, ulong p)
     956             : {
     957         357 :   long i, j, n = lg(E)-1;
     958         357 :   GEN M = cgetg(n+1, t_MAT);
     959        1281 :   for(i = 1; i <= n; i++)
     960             :   {
     961         924 :     gel(M,i) = cgetg(n+1, t_VEC);
     962        3906 :     for(j = 1; j <= n; j++)
     963        2982 :       gmael(M,i,j) = powuu(p,mael(E,i,j));
     964             :   }
     965         357 :   return M;
     966             : }
     967             : 
     968             : /* Assume there is a single p-isogeny */
     969             : 
     970             : static GEN
     971          84 : isomatdbl(GEN nf, GEN L, GEN M, ulong p, GEN T2, long flag)
     972             : {
     973          84 :   long i, j, n = lg(L) -1;
     974          84 :   GEN P = p > 3 ? polmodular_ZXX(p, 0, 0, 1): NULL;
     975          84 :   GEN V = cgetg(2*n+1, t_VEC);
     976          84 :   GEN N = cgetg(2*n+1, t_MAT);
     977         266 :   for(i=1; i <= n; i++)
     978         182 :     gel(V, i) = gel(L, i);
     979         266 :   for(i=1; i <= n; i++)
     980             :   {
     981         182 :     GEN e = gel(L, i);
     982             :     GEN F, E;
     983         182 :     if (i == 1)
     984          84 :       F = gmael(T2, 2, 1);
     985             :     else
     986             :     {
     987          98 :       if (p > 3)
     988          42 :         F = ellisograph_iso(nf, e, p, P, NULL, flag);
     989             :       else
     990          56 :         F = gel(ellisograph_Kohel_iso(nf, e, p, NULL, flag), 1);
     991          98 :       if (lg(F) != 2) pari_err_BUG("isomatdbl");
     992             :     }
     993         182 :     E = gel(F, 1);
     994         182 :     if (flag)
     995          91 :       gel(V, i+n) = mkvec3(gel(E,1), gel(E,2), gel(E,3));
     996             :     else
     997             :     {
     998          91 :       GEN iso = ellisogenyapply(gel(E,4), gel(e, 4));
     999          91 :       GEN isot = ellisogenyapply(gel(e,5), gel(E, 5));
    1000          91 :       gel(V, i+n) = mkvec5(gel(E,1), gel(E,2), gel(E,3), iso, isot);
    1001             :     }
    1002             :   }
    1003         448 :   for(i=1; i <= 2*n; i++)
    1004         364 :     gel(N, i) = cgetg(2*n+1, t_COL);
    1005         266 :   for(i=1; i <= n; i++)
    1006         588 :     for(j=1; j <= n; j++)
    1007             :     {
    1008         406 :       gcoeff(N,i,j) = gcoeff(N,i+n,j+n) = gcoeff(M,i,j);
    1009         406 :       gcoeff(N,i,j+n) = gcoeff(N,i+n,j) = muliu(gcoeff(M,i,j), p);
    1010             :     }
    1011          84 :   return mkvec2(list_to_crv(V), N);
    1012             : }
    1013             : 
    1014             : INLINE GEN
    1015         336 : mkfracss(long x, long y) { retmkfrac(stoi(x),stoi(y)); }
    1016             : 
    1017             : static ulong
    1018         329 : ellQ_exceptional_iso(GEN j, GEN *jt, GEN *jtp, GEN *s0)
    1019             : {
    1020         329 :   *jt = j; *jtp = gen_1;
    1021         329 :   if (typ(j)==t_INT)
    1022             :   {
    1023         238 :     long js = itos_or_0(j);
    1024             :     GEN j37;
    1025         238 :     if (js==-32768) { *s0 = mkfracss(-1156,539); return 11; }
    1026         224 :     if (js==-121)
    1027          14 :       { *jt = stoi(-24729001) ; *jtp = mkfracss(4973,5633);
    1028          14 :         *s0 = mkfracss(-1961682050,1204555087); return 11;}
    1029         210 :     if (js==-24729001)
    1030          14 :       { *jt = stoi(-121); *jtp = mkfracss(5633,4973);
    1031          14 :         *s0 = mkfracss(-1961682050,1063421347); return 11;}
    1032         196 :     if (js==-884736)
    1033          14 :       { *s0 = mkfracss(-1100,513); return 19; }
    1034         182 :     j37 = negi(uu32toi(37876312,1780746325));
    1035         182 :     if (js==-9317)
    1036             :     {
    1037          14 :       *jt = j37;
    1038          14 :       *jtp = mkfracss(1984136099,496260169);
    1039          14 :       *s0 = mkfrac(negi(uu32toi(457100760,4180820796UL)),
    1040             :                         uu32toi(89049913, 4077411069UL));
    1041          14 :       return 37;
    1042             :     }
    1043         168 :     if (equalii(j, j37))
    1044             :     {
    1045          14 :       *jt = stoi(-9317);
    1046          14 :       *jtp = mkfrac(utoi(496260169),utoi(1984136099UL));
    1047          14 :       *s0 = mkfrac(negi(uu32toi(41554614,2722784052UL)),
    1048             :                         uu32toi(32367030,2614994557UL));
    1049          14 :       return 37;
    1050             :     }
    1051         154 :     if (js==-884736000)
    1052          14 :     { *s0 = mkfracss(-1073708,512001); return 43; }
    1053         140 :     if (equalii(j, negi(powuu(5280,3))))
    1054          14 :     { *s0 = mkfracss(-176993228,85184001); return 67; }
    1055         126 :     if (equalii(j, negi(powuu(640320,3))))
    1056          14 :     { *s0 = mkfrac(negi(uu32toi(72512,1969695276)), uu32toi(35374,1199927297));
    1057          14 :       return 163; }
    1058             :   } else
    1059             :   {
    1060          91 :     GEN j1 = mkfracss(-297756989,2);
    1061          91 :     GEN j2 = mkfracss(-882216989,131072);
    1062          91 :     if (gequal(j, j1))
    1063             :     {
    1064          14 :       *jt = j2; *jtp = mkfracss(1503991,2878441);
    1065          14 :       *s0 = mkfrac(negi(uu32toi(121934,548114672)),uu32toi(77014,117338383));
    1066          14 :       return 17;
    1067             :     }
    1068          77 :     if (gequal(j, j2))
    1069             :     {
    1070          14 :       *jt = j1; *jtp = mkfracss(2878441,1503991);
    1071          14 :       *s0 = mkfrac(negi(uu32toi(121934,548114672)),uu32toi(40239,4202639633UL));
    1072          14 :       return 17;
    1073             :     }
    1074             :   }
    1075         175 :   return 0;
    1076             : }
    1077             : 
    1078             : static GEN
    1079         273 : mkisomat(ulong p, GEN T)
    1080             : {
    1081         273 :   pari_sp av = avma;
    1082         273 :   GEN L = list_to_crv(etree_list(T));
    1083         273 :   GEN M = distmat_pow(etree_distmat(T), p);
    1084         273 :   return gerepilecopy(av, mkvec2(L, M));
    1085             : }
    1086             : 
    1087             : static GEN
    1088          84 : mkisomatdbl(ulong p, GEN T, ulong p2, GEN T2, long flag)
    1089             : {
    1090          84 :   GEN L = etree_list(T);
    1091          84 :   GEN M = distmat_pow(etree_distmat(T), p);
    1092          84 :   return isomatdbl(NULL, L, M, p2, T2, flag);
    1093             : }
    1094             : 
    1095             : /*
    1096             : See
    1097             : M.A Kenku
    1098             : On the number of Q-isomorphism classes of elliptic curves in each Q-isogeny class
    1099             : Journal of Number Theory
    1100             : Volume 15, Issue 2, October 1982, Pages 199-202
    1101             : http://www.sciencedirect.com/science/article/pii/0022314X82900257
    1102             : */
    1103             : 
    1104             : enum { _2 = 1, _3 = 2, _5 = 4, _7 = 8, _13 = 16 };
    1105             : static ulong
    1106         175 : ellQ_goodl(GEN E)
    1107             : {
    1108             :   forprime_t T;
    1109         175 :   long i, CM = ellQ_get_CM(E);
    1110         175 :   ulong mask = 31;
    1111         175 :   GEN disc = ell_get_disc(E);
    1112         175 :   pari_sp av = avma;
    1113         175 :   u_forprime_init(&T, 17UL,ULONG_MAX);
    1114        3675 :   for(i=1; mask && i<=20; i++)
    1115             :   {
    1116        3500 :     ulong p = u_forprime_next(&T);
    1117        3500 :     if (umodiu(disc,p)==0) i--;
    1118             :     else
    1119             :     {
    1120        3500 :       long t = ellap_CM_fast(E, p, CM), D = t*t-4*p;
    1121        3500 :       if (t%2) mask &= ~_2;
    1122        3500 :       if ((mask & _3) && kross(D,3)==-1)  mask &= ~_3;
    1123        3500 :       if ((mask & _5) && kross(D,5)==-1)  mask &= ~_5;
    1124        3500 :       if ((mask & _7) && kross(D,7)==-1)  mask &= ~_7;
    1125        3500 :       if ((mask &_13) && kross(D,13)==-1) mask &= ~_13;
    1126             :     }
    1127             :   }
    1128         175 :   avma = av; return mask;
    1129             : }
    1130             : 
    1131             : static long
    1132          14 : ellQ_goodl_l(GEN E, long l)
    1133             : {
    1134             :   forprime_t T;
    1135             :   long i;
    1136          14 :   GEN disc = ell_get_disc(E);
    1137          14 :   pari_sp av = avma;
    1138          14 :   u_forprime_init(&T, 17UL,ULONG_MAX);
    1139         294 :   for(i=1; i<=20; i++)
    1140             :   {
    1141         280 :     ulong p = u_forprime_next(&T);
    1142         280 :     if (umodiu(disc,p)==0) { i--; continue; }
    1143             :     else
    1144             :     {
    1145         280 :       long t = itos(ellap(E, utoi(p)));
    1146         280 :       if (l==2)
    1147             :       {
    1148         140 :         if (t%2==1) return 0;
    1149             :       }
    1150             :       else
    1151             :       {
    1152         140 :         long D = t*t-4*p;
    1153         140 :         if (kross(D,l)==-1) return 0;
    1154             :       }
    1155         280 :       avma = av;
    1156             :     }
    1157             :   }
    1158          14 :   return 1;
    1159             : }
    1160             : 
    1161             : static long
    1162          21 : ellnf_goodl_l(GEN E, long l)
    1163             : {
    1164             :   forprime_t T;
    1165             :   long i;
    1166          21 :   GEN nf = ellnf_get_nf(E);
    1167          21 :   GEN disc = ell_get_disc(E);
    1168          21 :   pari_sp av = avma;
    1169          21 :   u_forprime_init(&T, 17UL,ULONG_MAX);
    1170         329 :   for(i=1; i<=20; i++)
    1171             :   {
    1172         315 :     ulong p = u_forprime_next(&T);
    1173         315 :     GEN pr = idealprimedec(nf, utoi(p));
    1174         315 :     long j, g = lg(pr)-1;
    1175         770 :     for (j=1; j<=g; j++)
    1176             :     {
    1177         462 :       GEN prj = gel(pr, j);
    1178         462 :       if (idealval(nf,disc,prj)>0) {i--; continue;}
    1179             :       else
    1180             :       {
    1181         455 :         long t = itos(ellap(E, prj));
    1182         455 :         if (l==2)
    1183             :         {
    1184           0 :           if (t%2==1) return 0;
    1185             :         }
    1186             :         else
    1187             :         {
    1188         455 :           GEN D = subii(sqrs(t),shifti(pr_norm(prj),2));
    1189         455 :           if (krois(D,l)==-1) return 0;
    1190             :         }
    1191             :       }
    1192             :     }
    1193         308 :     avma = av;
    1194             :   }
    1195          14 :   return 1;
    1196             : }
    1197             : 
    1198             : static GEN
    1199         329 : ellQ_isomat(GEN E, long flag)
    1200             : {
    1201         329 :   GEN K = NULL;
    1202         329 :   GEN T2 = NULL, T3 = NULL, T5, T7, T13;
    1203             :   ulong good;
    1204             :   long n2, n3, n5, n7, n13;
    1205             :   GEN jt, jtp, s0;
    1206         329 :   GEN c4 = ell_get_c4(E), c6 = ell_get_c6(E), j = ell_get_j(E);
    1207         329 :   long l = ellQ_exceptional_iso(j, &jt, &jtp, &s0);
    1208         329 :   if (l)
    1209             :   {
    1210             : #if 1
    1211         154 :     return mkisomat(l, ellisograph_dummy(E, l, jt, jtp, s0, flag));
    1212             : #else
    1213             :     return mkisomat(l, ellisograph_p(K, E, l), flag);
    1214             : #endif
    1215             :   }
    1216         175 :   good = ellQ_goodl(ellintegralmodel(E,NULL));
    1217         175 :   if (good & _2)
    1218             :   {
    1219         126 :     T2 = ellisograph_p(K, E, 2, flag);
    1220         126 :     n2 = etree_nbnodes(T2);
    1221         126 :     if (n2>4 || gequalgs(j, 1728) || gequalgs(j, 287496))
    1222          42 :       return mkisomat(2, T2);
    1223          49 :   } else n2 = 1;
    1224         133 :   if (good & _3)
    1225             :   {
    1226          84 :     T3 = ellisograph_p(K, E, 3, flag);
    1227          84 :     n3 = etree_nbnodes(T3);
    1228          84 :     if (n3>1 && n2==2) return mkisomatdbl(3,T3,2,T2, flag);
    1229          42 :     if (n3==2 && n2>1)  return mkisomatdbl(2,T2,3,T3, flag);
    1230          42 :     if (n3>2 || gequal0(j)) return mkisomat(3, T3);
    1231          49 :   } else n3 = 1;
    1232          63 :   if (good & _5)
    1233             :   {
    1234          21 :     T5 = ellisograph_p(K, E, 5, flag);
    1235          21 :     n5 = etree_nbnodes(T5);
    1236          21 :     if (n5>1 && n2>1) return mkisomatdbl(2,T2,5,T5, flag);
    1237          21 :     if (n5>1 && n3>1) return mkisomatdbl(3,T3,5,T5, flag);
    1238           7 :     if (n5>1) return mkisomat(5, T5);
    1239          42 :   } else n5 = 1;
    1240          42 :   if (good & _7)
    1241             :   {
    1242          28 :     T7 = ellisograph_p(K, E, 7, flag);
    1243          28 :     n7 = etree_nbnodes(T7);
    1244          28 :     if (n7>1 && n2>1) return mkisomatdbl(2,T2,7,T7, flag);
    1245           0 :     if (n7>1 && n3>1) return mkisomatdbl(3,T3,7,T7, flag);
    1246           0 :     if (n7>1) return mkisomat(7,T7);
    1247          14 :   } else n7 = 1;
    1248          14 :   if (n2>1) return mkisomat(2,T2);
    1249           0 :   if (n3>1) return mkisomat(3,T3);
    1250           0 :   if (good & _13)
    1251             :   {
    1252           0 :     T13 = ellisograph_p(K, E, 13, flag);
    1253           0 :     n13 = etree_nbnodes(T13);
    1254           0 :     if (n13>1) return mkisomat(13,T13);
    1255           0 :   } else n13 = 1;
    1256           0 :   if (flag)
    1257           0 :     retmkvec2(list_to_crv(mkvec(mkvec3(c4, c6, j))), matid(1));
    1258             :   else
    1259           0 :     retmkvec2(list_to_crv(mkvec(mkvec5(c4, c6, j, isogeny_a4a6(E), invisogeny_a4a6(E)))), matid(1));
    1260             : }
    1261             : 
    1262             : GEN
    1263         364 : ellisomat(GEN E, long p, long flag)
    1264             : {
    1265         364 :   pari_sp av = avma;
    1266             :   GEN r;
    1267         364 :   if (flag < 0 || flag > 1) pari_err_FLAG("ellisomat");
    1268         364 :   if (p < 0) pari_err_PRIME("ellisomat", utoi(p));
    1269         364 :   if (p == 1) { flag = 1; p = 0; } /* for backward compatibility */
    1270         364 :   checkell(E);
    1271         364 :   switch(ell_get_type(E))
    1272             :   {
    1273           0 :     default: pari_err_TYPE("ellisomat",E);
    1274           0 :       return NULL; /* NOT REACHED */
    1275             :     case t_ELL_Q:
    1276         343 :       if (!p) r = ellQ_isomat(E, flag);
    1277             :       else
    1278             :       {
    1279          14 :         if (ellQ_goodl_l(E, p))
    1280          14 :           r = mkisomat(p, ellisograph_p(NULL, E, p, flag));
    1281           0 :         else r = mkvec2(mkvec(ellisograph_a4a6(E, flag)),matid(1));
    1282             :       }
    1283         343 :       break;
    1284             :     case t_ELL_NF:
    1285          21 :       if (!p) {
    1286           0 :         pari_err_IMPL("ellisomat(E,0) for curve over number fields");
    1287           0 :         return NULL; /* NOT REACHED */ }
    1288             :       else
    1289             :       {
    1290          21 :         if (ellnf_goodl_l(E, p))
    1291          14 :           r = mkisomat(p, ellisograph_p(ellnf_get_nf(E), E, p, flag));
    1292           7 :         else r = mkvec2(mkvec(ellisograph_a4a6(E, flag)),matid(1));
    1293             :       }
    1294             :   }
    1295         364 :   return gerepilecopy(av, r);
    1296             : }

Generated by: LCOV version 1.11