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 16791-0d1274a) Lines: 302 307 98.4 %
Date: 2014-09-16 Functions: 21 21 100.0 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 82 97 84.5 %

           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                 :        625 : ellisweierstrasspoint(GEN E, GEN Q)
      20 [ +  - ][ +  + ]:        625 : { return ell_is_inf(Q) || gequal0(ec_dFdy_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                 :        510 : make_velu_curve(GEN E, GEN t, GEN w)
      28                 :            : {
      29                 :        510 :   GEN A4, A6, a1 = ell_get_a1(E), a2 = ell_get_a2(E), a3 = ell_get_a3(E);
      30                 :        510 :   A4 = gsub(ell_get_a4(E), gmulsg(5L, t));
      31                 :        510 :   A6 = gsub(ell_get_a6(E), gadd(gmul(ell_get_b2(E), t), gmulsg(7L, w)));
      32                 :        510 :   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                 :         25 : get_isog_vars(GEN phi)
      39                 :            : {
      40                 :         25 :   long vx = varn(gel(phi, 1));
      41                 :         25 :   long vy = varn(gel(phi, 2));
      42         [ +  + ]:         25 :   if (vy == vx) vy = varn(leading_term(gel(phi, 2)));
      43                 :         25 :   return mkvecsmall2(vx, vy);
      44                 :            : }
      45                 :            : 
      46                 :            : INLINE GEN
      47                 :          5 : substvec(GEN target, GEN vars, GEN in)
      48                 :            : {
      49                 :          5 :   long nvars = lg(vars) - 1, i;
      50                 :          5 :   GEN res = target;
      51         [ +  + ]:         15 :   for (i = 1; i <= nvars; ++i)
      52                 :         10 :     res = gsubst(res, vars[i], gel(in, i));
      53                 :          5 :   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                 :         10 : ellcompisog(GEN F, GEN G)
      59                 :            : {
      60                 :         10 :   pari_sp av = avma;
      61                 :            :   GEN Gh, Gh2, Gh3, f, g, h, z, numz, denz, denz2, in, comp;
      62                 :            : 
      63                 :         10 :   checkellisog(F);
      64                 :          5 :   checkellisog(G);
      65                 :          5 :   Gh = gel(G,3); Gh2 = gsqr(Gh); Gh3 = gmul(Gh, Gh2);
      66                 :          5 :   in = mkvec2(gdiv(gel(G,1), Gh2), gdiv(gel(G,2), Gh3));
      67                 :          5 :   comp = substvec(F, get_isog_vars(F), in);
      68                 :          5 :   f = gel(comp,1);
      69                 :          5 :   g = gel(comp,2);
      70                 :          5 :   z = gel(comp,3); numz = numer(z); denz = denom(z);
      71         [ -  + ]:          5 :   if (!issquareall(denom(f), &h))
      72                 :          0 :     pari_err_BUG("ellcompisog (denominator of composite abscissa not square)");
      73                 :          5 :   h  = RgX_mul(h, numz);
      74                 :          5 :   h = RgX_Rg_div(h, leading_term(h));
      75                 :            : 
      76                 :          5 :   denz2 = gsqr(denz);
      77                 :          5 :   f = RgX_mul(numer(f), denz2);
      78                 :          5 :   g = RgX_mul(numer(g), gmul(denz,denz2));
      79                 :          5 :   return gerepilecopy(av, mkvec3(f,g,h));
      80                 :            : }
      81                 :            : 
      82                 :            : /* Given an isogeny phi from ellisogeny() and a point P in the domain of phi,
      83                 :            :  * return phi(P) */
      84                 :            : GEN
      85                 :         45 : ellisogenyapply(GEN phi, GEN P)
      86                 :            : {
      87                 :         45 :   pari_sp ltop = avma;
      88                 :            :   GEN f, g, h, img_f, img_g, img_h, img_h2, img_h3, img, vars, tmp;
      89                 :            :   long vx, vy;
      90         [ +  + ]:         45 :   if (lg(P) == 4) return ellcompisog(phi,P);
      91                 :         35 :   checkellisog(phi);
      92                 :         35 :   checkellpt(P);
      93         [ +  + ]:         30 :   if (ell_is_inf(P)) return ellinf();
      94                 :         20 :   f = gel(phi, 1);
      95                 :         20 :   g = gel(phi, 2);
      96                 :         20 :   h = gel(phi, 3);
      97                 :         20 :   vars = get_isog_vars(phi);
      98                 :         20 :   vx = vars[1];
      99                 :         20 :   vy = vars[2];
     100                 :         20 :   img_h = poleval(h, gel(P, 1));
     101         [ +  + ]:         20 :   if (gequal0(img_h)) { avma = ltop; return ellinf(); }
     102                 :            : 
     103                 :         15 :   img_h2 = gsqr(img_h);
     104                 :         15 :   img_h3 = gmul(img_h, img_h2);
     105                 :         15 :   img_f = poleval(f, gel(P, 1));
     106                 :            :   /* FIXME: This calculation of g is perhaps not as efficient as it could be */
     107                 :         15 :   tmp = gsubst(g, vx, gel(P, 1));
     108                 :         15 :   img_g = gsubst(tmp, vy, gel(P, 2));
     109                 :         15 :   img = cgetg(3, t_VEC);
     110                 :         15 :   gel(img, 1) = gdiv(img_f, img_h2);
     111                 :         15 :   gel(img, 2) = gdiv(img_g, img_h3);
     112                 :         35 :   return gerepileupto(ltop, img);
     113                 :            : }
     114                 :            : 
     115                 :            : /* isog = [f, g, h] = [x, y, 1] = identity */
     116                 :            : static GEN
     117                 :        165 : isog_identity(long vx, long vy)
     118                 :        165 : { return mkvec3(pol_x(vx), pol_x(vy), pol_1(vx)); }
     119                 :            : 
     120                 :            : /* Returns an updated value for isog based (primarily) on tQ and uQ. Used to
     121                 :            :  * iteratively compute the isogeny corresponding to a subgroup generated by a
     122                 :            :  * given point. Ref: Equation 8 in Velu's paper.
     123                 :            :  * isog = NULL codes the identity */
     124                 :            : static GEN
     125                 :        360 : update_isogeny_polys(GEN isog, GEN E, GEN Q, GEN tQ, GEN uQ, long vx, long vy)
     126                 :            : {
     127                 :        360 :   pari_sp ltop = avma, av;
     128                 :        360 :   GEN xQ = gel(Q, 1), yQ = gel(Q, 2);
     129                 :        360 :   GEN rt = deg1pol_shallow(gen_1, gneg(xQ), vx);
     130                 :        360 :   GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E);
     131                 :            : 
     132                 :        360 :   GEN gQx = ec_dFdx_evalQ(E, Q);
     133                 :        360 :   GEN gQy = ec_dFdy_evalQ(E, Q);
     134                 :            :   GEN tmp1, tmp2, tmp3, tmp4, f, g, h, rt_sqr, res;
     135                 :            : 
     136                 :            :   /* g -= uQ * (2 * y + E.a1 * x + E.a3)
     137                 :            :    *   + tQ * rt * (E.a1 * rt + y - yQ)
     138                 :            :    *   + rt * (E.a1 * uQ - gQx * gQy) */
     139                 :        360 :   av = avma;
     140                 :        360 :   tmp1 = gmul(uQ, gadd(deg1pol_shallow(gen_2, gen_0, vy),
     141                 :            :                        deg1pol_shallow(a1, a3, vx)));
     142                 :        360 :   tmp1 = gerepileupto(av, tmp1);
     143                 :        360 :   av = avma;
     144                 :        360 :   tmp2 = gmul(tQ, gadd(gmul(a1, rt),
     145                 :            :                        deg1pol_shallow(gen_1, gneg(yQ), vy)));
     146                 :        360 :   tmp2 = gerepileupto(av, tmp2);
     147                 :        360 :   av = avma;
     148                 :        360 :   tmp3 = gsub(gmul(a1, uQ), gmul(gQx, gQy));
     149                 :        360 :   tmp3 = gerepileupto(av, tmp3);
     150                 :            : 
     151         [ +  + ]:        360 :   if (!isog) isog = isog_identity(vx,vy);
     152                 :        360 :   f = gel(isog, 1);
     153                 :        360 :   g = gel(isog, 2);
     154                 :        360 :   h = gel(isog, 3);
     155                 :        360 :   rt_sqr = gsqr(rt);
     156                 :        360 :   res = cgetg(4, t_VEC);
     157                 :        360 :   av = avma;
     158                 :        360 :   tmp4 = gdiv(gadd(gmul(tQ, rt), uQ), rt_sqr);
     159                 :        360 :   gel(res, 1) = gerepileupto(av, gadd(f, tmp4));
     160                 :        360 :   av = avma;
     161                 :        360 :   tmp4 = gadd(tmp1, gmul(rt, gadd(tmp2, tmp3)));
     162                 :        360 :   gel(res, 2) = gerepileupto(av, gsub(g, gdiv(tmp4, gmul(rt, rt_sqr))));
     163                 :        360 :   av = avma;
     164                 :        360 :   gel(res, 3) = gerepileupto(av, gmul(h, rt));
     165                 :        360 :   return gerepileupto(ltop, res);
     166                 :            : }
     167                 :            : 
     168                 :            : /* Given a point P on E, return the curve E/<P> and, if only_image is zero,
     169                 :            :  * the isogeny pi: E -> E/<P>. The variables vx and vy are used to describe
     170                 :            :  * the isogeny (ignored if only_image is zero) */
     171                 :            : static GEN
     172                 :        290 : isogeny_from_kernel_point(GEN E, GEN P, int only_image, long vx, long vy)
     173                 :            : {
     174                 :        290 :   pari_sp av = avma;
     175                 :            :   GEN isog, EE, f, g, h, h2, h3;
     176                 :        290 :   GEN Q = P, t = gen_0, w = gen_0;
     177                 :            :   long c;
     178         [ +  + ]:        290 :   if (!oncurve(E,P))
     179                 :          5 :     pari_err_DOMAIN("isogeny_from_kernel_point", "point", "not on", E, P);
     180         [ +  + ]:        285 :   if (ell_is_inf(P))
     181                 :            :   {
     182         [ +  + ]:         30 :     if (only_image) return E;
     183                 :         20 :     return mkvec2(E, isog_identity(vx,vy));
     184                 :            :   }
     185                 :            : 
     186                 :        255 :   isog = NULL; c = 1;
     187                 :            :   for (;;)
     188                 :            :   {
     189                 :        625 :     GEN tQ, xQ = gel(Q,1), uQ = ec_2divpol_evalx(E, xQ);
     190                 :        625 :     int stop = 0;
     191         [ +  + ]:        625 :     if (ellisweierstrasspoint(E,Q))
     192                 :            :     { /* ord(P)=2c; take Q=[c]P into consideration and stop */
     193                 :        125 :       tQ = ec_dFdx_evalQ(E, Q);
     194                 :        125 :       stop = 1;
     195                 :            :     }
     196                 :            :     else
     197                 :        500 :       tQ = ec_half_deriv_2divpol_evalx(E, xQ);
     198                 :        625 :     t = gadd(t, tQ);
     199                 :        625 :     w = gadd(w, gadd(uQ, gmul(tQ, xQ)));
     200         [ +  + ]:        625 :     if (!only_image) isog = update_isogeny_polys(isog, E, Q,tQ,uQ, vx,vy);
     201         [ +  + ]:        625 :     if (stop) break;
     202                 :            : 
     203                 :        500 :     Q = elladd(E, P, Q);
     204                 :        500 :     ++c;
     205                 :            :     /* IF x([c]P) = x([c-1]P) THEN [c]P = -[c-1]P and [2c-1]P = 0 */
     206         [ +  + ]:        500 :     if (gequal(gel(Q,1), xQ)) break;
     207         [ -  + ]:        370 :     if (gc_needed(av,1))
     208                 :            :     {
     209         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"isogeny_from_kernel_point");
     210         [ #  # ]:          0 :       gerepileall(av, isog? 4: 3, &Q, &t, &w, &isog);
     211                 :            :     }
     212                 :        370 :   }
     213                 :            : 
     214                 :        255 :   EE = make_velu_curve(E, t, w);
     215         [ +  + ]:        255 :   if (only_image) return EE;
     216                 :            : 
     217         [ -  + ]:        145 :   if (!isog) isog = isog_identity(vx,vy);
     218                 :        145 :   f = gel(isog, 1);
     219                 :        145 :   g = gel(isog, 2);
     220 [ +  - ][ -  + ]:        145 :   if ( ! (typ(f) == t_RFRAC && typ(g) == t_RFRAC))
     221                 :          0 :     pari_err_BUG("isogeny_from_kernel_point (f or g has wrong type)");
     222                 :            : 
     223                 :            :   /* Clean up the isogeny polynomials f, g and h so that the isogeny
     224                 :            :    * is given by (x,y) -> (f(x)/h(x)^2, g(x,y)/h(x)^3) */
     225                 :        145 :   h = gel(isog, 3);
     226                 :        145 :   h2 = gsqr(h);
     227                 :        145 :   h3 = gmul(h, h2);
     228                 :        145 :   f = gmul(f, h2);
     229                 :        145 :   g = gmul(g, h3);
     230 [ +  - ][ -  + ]:        145 :   if (typ(f) != t_POL || typ(g) != t_POL)
     231                 :          0 :     pari_err_BUG("isogeny_from_kernel_point (wrong denominator)");
     232                 :        285 :   return mkvec2(EE, mkvec3(f,g, gel(isog,3)));
     233                 :            : }
     234                 :            : 
     235                 :            : /* Given a t_POL x^n - s1 x^{n-1} + s2 x^{n-2} - s3 x^{n-3} + ...
     236                 :            :  * return the first three power sums (Newton's identities):
     237                 :            :  *   p1 = s1
     238                 :            :  *   p2 = s1^2 - 2 s2
     239                 :            :  *   p3 = (s1^2 - 3 s2) s1 + 3 s3 */
     240                 :            : static void
     241                 :        265 : first_three_power_sums(GEN pol, GEN *p1, GEN *p2, GEN *p3)
     242                 :            : {
     243                 :        265 :   long d = degpol(pol);
     244                 :            :   GEN s1, s2, ms3;
     245                 :            : 
     246                 :        265 :   *p1 = s1 = gneg(RgX_coeff(pol, d-1));
     247                 :            : 
     248                 :        265 :   s2 = RgX_coeff(pol, d-2);
     249                 :        265 :   *p2 = gsub(gsqr(s1), gmulsg(2L, s2));
     250                 :            : 
     251                 :        265 :   ms3 = RgX_coeff(pol, d-3);
     252                 :        265 :   *p3 = gadd(gmul(s1, gsub(*p2, s2)), gmulsg(-3L, ms3));
     253                 :        265 : }
     254                 :            : 
     255                 :            : 
     256                 :            : /* Let E and a t_POL h of degree 1 or 3 whose roots are 2-torsion points on E.
     257                 :            :  * - if only_image != 0, return [t, w] used to compute the equation of the
     258                 :            :  *   quotient by the given 2-torsion points
     259                 :            :  * - else return [t,w, f,g,h], along with the contributions f, g and
     260                 :            :  *   h to the isogeny giving the quotient by h. Variables vx and vy are used
     261                 :            :  *   to create f, g and h, or ignored if only_image is zero */
     262                 :            : 
     263                 :            : /* deg h = 1; 2-torsion contribution from Weierstrass point */
     264                 :            : static GEN
     265                 :        100 : contrib_weierstrass_pt(GEN E, GEN h, long only_image, long vx, long vy)
     266                 :            : {
     267                 :        100 :   GEN p = ellbasechar(E);
     268                 :        100 :   GEN a1 = ell_get_a1(E);
     269                 :        100 :   GEN a3 = ell_get_a3(E);
     270                 :        100 :   GEN x0 = gneg(constant_term(h)); /* h = x - x0 */
     271                 :        100 :   GEN b = gadd(gmul(a1,x0), a3);
     272                 :            :   GEN y0, Q, t, w, t1, t2, f, g;
     273                 :            : 
     274         [ +  + ]:        100 :   if (!equalis(p, 2L)) /* char(k) != 2 ==> y0 = -b/2 */
     275                 :         70 :     y0 = gmul2n(gneg(b), -1);
     276                 :            :   else
     277                 :            :   { /* char(k) = 2 ==> y0 = sqrt(f(x0)) where E is y^2 + h(x) = f(x). */
     278         [ -  + ]:         30 :     if (!gequal0(b)) pari_err_BUG("two_torsion_contrib (a1*x0+a3 != 0)");
     279                 :         30 :     y0 = gsqrt(ec_f_evalx(E, x0), 0);
     280                 :            :   }
     281                 :        100 :   Q = mkvec2(x0, y0);
     282                 :        100 :   t = ec_dFdx_evalQ(E, Q);
     283                 :        100 :   w = gmul(x0, t);
     284         [ +  + ]:        100 :   if (only_image) return mkvec2(t,w);
     285                 :            : 
     286                 :            :   /* Compute isogeny, f = (x - x0) * t */
     287                 :         50 :   f = deg1pol_shallow(t, gmul(t, gneg(x0)), vx);
     288                 :            : 
     289                 :            :   /* g = (x - x0) * t * (a1 * (x - x0) + (y - y0)) */
     290                 :         50 :   t1 = deg1pol_shallow(a1, gmul(a1, gneg(x0)), vx);
     291                 :         50 :   t2 = deg1pol_shallow(gen_1, gneg(y0), vy);
     292                 :         50 :   g = gmul(f, gadd(t1, t2));
     293                 :        100 :   return mkvec5(t, w, f, g, h);
     294                 :            : }
     295                 :            : /* deg h =3; full 2-torsion contribution. NB: assume h is monic; base field
     296                 :            :  * characteristic is odd or zero (cannot happen in char 2).*/
     297                 :            : static GEN
     298                 :         10 : contrib_full_tors(GEN E, GEN h, long only_image, long vx, long vy)
     299                 :            : {
     300                 :            :   GEN p1, p2, p3, half_b2, half_b4, t, w, f, g;
     301                 :         10 :   first_three_power_sums(h, &p1,&p2,&p3);
     302                 :         10 :   half_b2 = gmul2n(ell_get_b2(E), -1);
     303                 :         10 :   half_b4 = gmul2n(ell_get_b4(E), -1);
     304                 :            : 
     305                 :            :   /* t = 3*(p2 + b4/2) + p1 * b2/2 */
     306                 :         10 :   t = gadd(gmulsg(3L, gadd(p2, half_b4)), gmul(p1, half_b2));
     307                 :            : 
     308                 :            :   /* w = 3 * p3 + p2 * b2/2 + p1 * b4/2 */
     309                 :         10 :   w = gadd(gmulsg(3L, p3), gadd(gmul(p2, half_b2),
     310                 :            :                                 gmul(p1, half_b4)));
     311         [ +  + ]:         10 :   if (only_image) return mkvec2(t,w);
     312                 :            : 
     313                 :            :   /* Compute isogeny */
     314                 :            :   {
     315                 :          5 :     GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E), t1, t2;
     316                 :          5 :     GEN s1 = gneg(RgX_coeff(h, 2));
     317                 :          5 :     GEN dh = RgX_deriv(h);
     318                 :          5 :     GEN psi2xy = gadd(deg1pol_shallow(a1, a3, vx),
     319                 :            :                       deg1pol_shallow(gen_2, gen_0, vy));
     320                 :            : 
     321                 :            :     /* f = -3 (3 x + b2/2 + s1) h + (3 x^2 + (b2/2) x + (b4/2)) h'*/
     322                 :          5 :     t1 = RgX_mul(h, gmulsg(-3, deg1pol(stoi(3), gadd(half_b2, s1), vx)));
     323                 :          5 :     t2 = mkpoln(3, stoi(3), half_b2, half_b4);
     324                 :          5 :     setvarn(t2, vx);
     325                 :          5 :     t2 = RgX_mul(dh, t2);
     326                 :          5 :     f = RgX_add(t1, t2);
     327                 :            : 
     328                 :            :     /* 2g = psi2xy * (f'*h - f*h') - (a1*f + a3*h) * h; */
     329                 :          5 :     t1 = RgX_sub(RgX_mul(RgX_deriv(f), h), RgX_mul(f, dh));
     330                 :          5 :     t2 = RgX_mul(h, RgX_add(RgX_Rg_mul(f, a1), RgX_Rg_mul(h, a3)));
     331                 :          5 :     g = RgX_divs(gsub(gmul(psi2xy, t1), t2), 2L);
     332                 :            : 
     333                 :          5 :     f = RgX_mul(f, h);
     334                 :          5 :     g = RgX_mul(g, h);
     335                 :            :   }
     336                 :         10 :   return mkvec5(t, w, f, g, h);
     337                 :            : }
     338                 :            : 
     339                 :            : /* Given E and a t_POL T whose roots define a subgroup G of E, return the factor
     340                 :            :  * of T that corresponds to the 2-torsion points E[2] \cap G in G */
     341                 :            : INLINE GEN
     342                 :        260 : two_torsion_part(GEN E, GEN T)
     343                 :        260 : { return RgX_gcd(T, elldivpol(E, 2, varn(T))); }
     344                 :            : 
     345                 :            : /* Return the jth Hasse derivative of the polynomial f = \sum_{i=0}^n a_i x^i,
     346                 :            :  * i.e. \sum_{i=j}^n a_i \binom{i}{j} x^{i-j}. It is a derivation even when the
     347                 :            :  * coefficient ring has positive characteristic */
     348                 :            : static GEN
     349                 :         70 : derivhasse(GEN f, ulong j)
     350                 :            : {
     351                 :         70 :   ulong i, d = degpol(f);
     352                 :            :   GEN df;
     353 [ +  + ][ +  + ]:         70 :   if (gequal0(f) || d == 0) return pol_0(varn(f));
     354         [ -  + ]:         40 :   if (j == 0) return gcopy(f);
     355                 :         40 :   df = cgetg(2 + (d-j+1), t_POL);
     356                 :         40 :   df[1] = f[1];
     357         [ +  + ]:         80 :   for (i = j; i <= d; ++i) gel(df, i-j+2) = gmul(binomialuu(i,j), gel(f, i+2));
     358                 :         70 :   return normalizepol(df);
     359                 :            : }
     360                 :            : 
     361                 :            : static GEN
     362                 :        100 : non_two_torsion_abscissa(GEN E, GEN h0, GEN x)
     363                 :            : {
     364                 :            :   GEN mp1, dh0, ddh0, t, u, t1, t2, t3;
     365                 :        100 :   long m = degpol(h0);
     366                 :        100 :   mp1 = gel(h0, m + 1); /* negative of first power sum */
     367                 :        100 :   dh0 = RgX_deriv(h0);
     368                 :        100 :   ddh0 = RgX_deriv(dh0);
     369                 :        100 :   t = ec_2divpol_evalx(E, x);
     370                 :        100 :   u = ec_half_deriv_2divpol_evalx(E, x);
     371                 :        100 :   t1 = RgX_sub(RgX_sqr(dh0), RgX_mul(ddh0, h0));
     372                 :        100 :   t2 = RgX_mul(u, RgX_mul(h0, dh0));
     373                 :        100 :   t3 = RgX_mul(RgX_sqr(h0),
     374                 :        100 :                deg1pol_shallow(stoi(2*m), gmulsg(2L, mp1), varn(x)));
     375                 :            :   /* t * (dh0^2 - ddh0*h0) - u*dh0*h0 + (2*m*x - 2*s1) * h0^2); */
     376                 :        100 :   return RgX_add(RgX_sub(RgX_mul(t, t1), t2), t3);
     377                 :            : }
     378                 :            : 
     379                 :            : static GEN
     380                 :        130 : isog_abscissa(GEN E, GEN kerp, GEN h0, GEN x, GEN two_tors)
     381                 :            : {
     382                 :            :   GEN f0, f2, h2, t1, t2, t3;
     383         [ +  + ]:        130 :   f0 = (degpol(h0) > 0)? non_two_torsion_abscissa(E, h0, x): pol_0(varn(x));
     384                 :        130 :   f2 = gel(two_tors, 3);
     385                 :        130 :   h2 = gel(two_tors, 5);
     386                 :            : 
     387                 :            :   /* Combine f0 and f2 into the final abcissa of the isogeny. */
     388                 :        130 :   t1 = RgX_mul(x, RgX_sqr(kerp));
     389                 :        130 :   t2 = RgX_mul(f2, RgX_sqr(h0));
     390                 :        130 :   t3 = RgX_mul(f0, RgX_sqr(h2));
     391                 :            :   /* x * kerp^2 + f2 * h0^2 + f0 * h2^2 */
     392                 :        130 :   return RgX_add(t1, RgX_add(t2, t3));
     393                 :            : }
     394                 :            : 
     395                 :            : static GEN
     396                 :         95 : non_two_torsion_ordinate_char_not2(GEN E, GEN f, GEN h, GEN psi2)
     397                 :            : {
     398                 :         95 :   GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E);
     399                 :         95 :   GEN df = RgX_deriv(f), dh = RgX_deriv(h);
     400                 :            :   /* g = df * h * psi2/2 - f * dh * psi2
     401                 :            :    *   - (E.a1 * f + E.a3 * h^2) * h/2 */
     402                 :         95 :   GEN t1 = RgX_mul(df, RgX_mul(h, RgX_divs(psi2, 2L)));
     403                 :         95 :   GEN t2 = RgX_mul(f, RgX_mul(dh, psi2));
     404                 :         95 :   GEN t3 = RgX_mul(RgX_divs(h, 2L),
     405                 :            :                    RgX_add(RgX_Rg_mul(f, a1), RgX_Rg_mul(RgX_sqr(h), a3)));
     406                 :         95 :   return RgX_sub(RgX_sub(t1, t2), t3);
     407                 :            : }
     408                 :            : 
     409                 :            : /* h = kerq */
     410                 :            : static GEN
     411                 :         35 : non_two_torsion_ordinate_char2(GEN E, GEN h, GEN x, GEN y)
     412                 :            : {
     413                 :         35 :   GEN a1 = ell_get_a1(E), a3 = ell_get_a3(E), a4 = ell_get_a4(E);
     414                 :         35 :   GEN b2 = ell_get_b2(E), b4 = ell_get_b4(E), b6 = ell_get_b6(E);
     415                 :            :   GEN h2, dh, dh2, ddh, D2h, D2dh, H, psi2, u, t, alpha;
     416                 :            :   GEN p1, t1, t2, t3, t4;
     417                 :         35 :   long m, vx = varn(x);
     418                 :            : 
     419                 :         35 :   h2 = RgX_sqr(h);
     420                 :         35 :   dh = RgX_deriv(h);
     421                 :         35 :   dh2 = RgX_sqr(dh);
     422                 :         35 :   ddh = RgX_deriv(dh);
     423                 :         35 :   H = RgX_sub(dh2, RgX_mul(h, ddh));
     424                 :         35 :   D2h = derivhasse(h, 2);
     425                 :         35 :   D2dh = derivhasse(dh, 2);
     426                 :         35 :   psi2 = deg1pol_shallow(a1, a3, vx);
     427                 :         35 :   u = mkpoln(3, b2, gen_0, b6);
     428                 :         35 :   setvarn(u, vx);
     429                 :         35 :   t = deg1pol_shallow(b2, b4, vx);
     430                 :         35 :   alpha = mkpoln(4, a1, a3, gmul(a1, a4), gmul(a3, a4));
     431                 :         35 :   setvarn(alpha, vx);
     432                 :         35 :   m = degpol(h);
     433                 :         35 :   p1 = RgX_coeff(h, m-1); /* first power sum */
     434                 :            : 
     435                 :         35 :   t1 = gmul(gadd(gmul(a1, p1), gmulgs(a3, m)), RgX_mul(h,h2));
     436                 :            : 
     437                 :         35 :   t2 = gmul(a1, gadd(gmul(a1, gadd(y, psi2)), RgX_add(RgX_Rg_add(RgX_sqr(x), a4), t)));
     438                 :         35 :   t2 = gmul(t2, gmul(dh, h2));
     439                 :            : 
     440                 :         35 :   t3 = gadd(gmul(y, t), RgX_add(alpha, RgX_Rg_mul(u, a1)));
     441                 :         35 :   t3 = gmul(t3, RgX_mul(h, H));
     442                 :            : 
     443                 :         35 :   t4 = gmul(u, psi2);
     444                 :         35 :   t4 = gmul(t4, RgX_sub(RgX_sub(RgX_mul(h2, D2dh), RgX_mul(dh, H)),
     445                 :            :                         RgX_mul(h, RgX_mul(dh, D2h))));
     446                 :            : 
     447                 :         35 :   return gadd(t1, gadd(t2, gadd(t3, t4)));
     448                 :            : }
     449                 :            : 
     450                 :            : static GEN
     451                 :        130 : isog_ordinate(GEN E, GEN kerp, GEN kerq, GEN x, GEN y, GEN two_tors, GEN f)
     452                 :            : {
     453                 :            :   GEN g;
     454         [ +  + ]:        130 :   if (! equalis(ellbasechar(E), 2L)) {
     455                 :            :     /* FIXME: We don't use (hence don't need to calculate)
     456                 :            :      * g2 = gel(two_tors, 4) when char(k) != 2. */
     457                 :         95 :     GEN psi2 = gneg(ec_dFdy_evalQ(E, mkvec2(x, y)));
     458                 :         95 :     g = non_two_torsion_ordinate_char_not2(E, f, kerp, psi2);
     459                 :            :   } else {
     460                 :         35 :     GEN h2 = gel(two_tors, 5);
     461                 :         35 :     GEN g2 = gmul(gel(two_tors, 4), RgX_mul(kerq, RgX_sqr(kerq)));
     462                 :         35 :     GEN g0 = non_two_torsion_ordinate_char2(E, kerq, x, y);
     463                 :         35 :     g0 = gmul(g0, RgX_mul(h2, RgX_sqr(h2)));
     464                 :         35 :     g = gsub(gmul(y, RgX_mul(kerp, RgX_sqr(kerp))), gadd(g2, g0));
     465                 :            :   }
     466                 :        130 :   return g;
     467                 :            : }
     468                 :            : 
     469                 :            : /* Given an elliptic curve E and a polynomial kerp whose roots give the
     470                 :            :  * x-coordinates of a subgroup G of E, return the curve E/G and,
     471                 :            :  * if only_image is zero, the isogeny pi:E -> E/G. Variables vx and vy are
     472                 :            :  * used to describe the isogeny (and are ignored if only_image is zero). */
     473                 :            : static GEN
     474                 :        260 : isogeny_from_kernel_poly(GEN E, GEN kerp, long only_image, long vx, long vy)
     475                 :            : {
     476                 :            :   long m;
     477                 :        260 :   GEN b2 = ell_get_b2(E), b4 = ell_get_b4(E), b6 = ell_get_b6(E);
     478                 :            :   GEN p1, p2, p3, x, y, f, g, two_tors, EE, t, w;
     479                 :        260 :   GEN kerh = two_torsion_part(E, kerp);
     480                 :        260 :   GEN kerq = RgX_divrem(kerp, kerh, ONLY_DIVIDES);
     481         [ -  + ]:        260 :   if (!kerq) pari_err_BUG("isogeny_from_kernel_poly");
     482                 :            :   /* isogeny degree: 2*degpol(kerp)+1-degpol(kerh) */
     483                 :        260 :   m = degpol(kerq);
     484                 :            : 
     485                 :        260 :   kerp = RgX_Rg_div(kerp, leading_term(kerp));
     486                 :        260 :   kerq = RgX_Rg_div(kerq, leading_term(kerq));
     487                 :        260 :   kerh = RgX_Rg_div(kerh, leading_term(kerh));
     488   [ +  +  +  + ]:        260 :   switch(degpol(kerh))
     489                 :            :   {
     490                 :            :   case 0:
     491                 :        145 :     two_tors = mkvec5(gen_0, gen_0, pol_0(vx), pol_0(vx), pol_1(vx));
     492                 :        145 :     break;
     493                 :            :   case 1:
     494                 :        100 :     two_tors = contrib_weierstrass_pt(E, kerh, only_image,vx,vy);
     495                 :        100 :     break;
     496                 :            :   case 3:
     497                 :         10 :     two_tors = contrib_full_tors(E, kerh, only_image,vx,vy);
     498                 :         10 :     break;
     499                 :            :   default:
     500                 :          5 :     two_tors = NULL;
     501                 :          5 :     pari_err_DOMAIN("isogeny_from_kernel_poly", "kernel polynomial",
     502                 :            :                     "does not define a subgroup of", E, kerp);
     503                 :            :   }
     504                 :        255 :   first_three_power_sums(kerq,&p1,&p2,&p3);
     505                 :        255 :   x = pol_x(vx);
     506                 :        255 :   y = pol_x(vy);
     507                 :            : 
     508                 :            :   /* t = 6 * p2 + b2 * p1 + m * b4, */
     509                 :        255 :   t = gadd(gmulsg(6L, p2), gadd(gmul(b2, p1), gmulsg(m, b4)));
     510                 :            : 
     511                 :            :   /* w = 10 * p3 + 2 * b2 * p2 + 3 * b4 * p1 + m * b6, */
     512                 :        255 :   w = gadd(gmulsg(10L, p3),
     513                 :            :            gadd(gmul(gmulsg(2L, b2), p2),
     514                 :            :                 gadd(gmul(gmulsg(3L, b4), p1), gmulsg(m, b6))));
     515                 :            : 
     516                 :        255 :   EE = make_velu_curve(E, gadd(t, gel(two_tors, 1)),
     517                 :        255 :                           gadd(w, gel(two_tors, 2)));
     518         [ +  + ]:        255 :   if (only_image) return EE;
     519                 :            : 
     520                 :        130 :   f = isog_abscissa(E, kerp, kerq, x, two_tors);
     521                 :        130 :   g = isog_ordinate(E, kerp, kerq, x, y, two_tors, f);
     522                 :        255 :   return mkvec2(EE, mkvec3(f,g,kerp));
     523                 :            : }
     524                 :            : 
     525                 :            : /* Given an elliptic curve E and a subgroup G of E, return the curve
     526                 :            :  * E/G and, if only_image is zero, the isogeny corresponding
     527                 :            :  * to the canonical surjection pi:E -> E/G. The variables vx and
     528                 :            :  * vy are used to describe the isogeny (and are ignored if
     529                 :            :  * only_image is zero). The subgroup G may be given either as
     530                 :            :  * a generating point P on E or as a polynomial kerp whose roots are
     531                 :            :  * the x-coordinates of the points in G */
     532                 :            : GEN
     533                 :        580 : ellisogeny(GEN E, GEN G, long only_image, long vx, long vy)
     534                 :            : {
     535                 :        580 :   pari_sp av = avma;
     536                 :            :   GEN j, z;
     537                 :        580 :   checkell(E);j = ell_get_j(E);
     538         [ +  + ]:        580 :   if (vx < 0) vx = 0;
     539         [ +  + ]:        580 :   if (vy < 0) vy = fetch_user_var("y");
     540         [ +  + ]:        580 :   if (varncmp(vx, vy) >= 0) pari_err_PRIORITY("ellisogeny", pol_x(vx), "<=", vy);
     541         [ +  + ]:        575 :   if (varncmp(vy, gvar(j)) >= 0) pari_err_PRIORITY("ellisogeny", j, ">=", vy);
     542      [ +  +  + ]:        570 :   switch(typ(G))
     543                 :            :   {
     544                 :            :   case t_VEC:
     545                 :        300 :     checkellpt(G);
     546         [ +  + ]:        300 :     if (!ell_is_inf(G))
     547                 :            :     {
     548                 :        270 :       GEN x =  gel(G,1), y = gel(G,2);
     549         [ +  + ]:        270 :       if (varncmp(vy, gvar(x)) >= 0) pari_err_PRIORITY("ellisogeny", x, ">=", vy);
     550         [ +  + ]:        265 :       if (varncmp(vy, gvar(y)) >= 0) pari_err_PRIORITY("ellisogeny", y, ">=", vy);
     551                 :            :     }
     552                 :        290 :     z = isogeny_from_kernel_point(E, G, only_image, vx, vy);
     553                 :        285 :     break;
     554                 :            :   case t_POL:
     555         [ +  + ]:        265 :     if (varncmp(vy, gvar(constant_term(G))) >= 0)
     556                 :          5 :       pari_err_PRIORITY("ellisogeny", constant_term(G), ">=", vy);
     557                 :        260 :     z = isogeny_from_kernel_poly(E, G, only_image, vx, vy);
     558                 :        255 :     break;
     559                 :            :   default:
     560                 :          5 :     z = NULL;
     561                 :          5 :     pari_err_TYPE("ellisogeny", G);
     562                 :            :   }
     563                 :        540 :   return gerepilecopy(av, z);
     564                 :            : }

Generated by: LCOV version 1.9