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-bordeaux1.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.8.0 lcov report (development 17953-c39f2e6) Lines: 613 623 98.4 %
Date: 2015-08-31 Functions: 51 51 100.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 234 270 86.7 %

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

Generated by: LCOV version 1.9