Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - FpE.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.10.0 lcov report (development 19825-b77c7f8) Lines: 911 997 91.4 %
Date: 2016-12-04 05:49:01 Functions: 97 105 92.4 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2009  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             : #include "paripriv.h"
      16             : 
      17             : /* Not so fast arithmetic with points over elliptic curves over Fp */
      18             : 
      19             : /***********************************************************************/
      20             : /**                                                                   **/
      21             : /**                              FpE                                  **/
      22             : /**                                                                   **/
      23             : /***********************************************************************/
      24             : 
      25             : /* These functions deal with point over elliptic curves over Fp defined
      26             :  * by an equation of the form y^2=x^3+a4*x+a6.
      27             :  * Most of the time a6 is omitted since it can be recovered from any point
      28             :  * on the curve.
      29             :  */
      30             : 
      31             : GEN
      32        1283 : RgE_to_FpE(GEN x, GEN p)
      33             : {
      34        1283 :   if (ell_is_inf(x)) return x;
      35        1284 :   retmkvec2(Rg_to_Fp(gel(x,1),p),Rg_to_Fp(gel(x,2),p));
      36             : }
      37             : 
      38             : GEN
      39         440 : FpE_to_mod(GEN x, GEN p)
      40             : {
      41         440 :   if (ell_is_inf(x)) return x;
      42         384 :   retmkvec2(Fp_to_mod(gel(x,1),p),Fp_to_mod(gel(x,2),p));
      43             : }
      44             : 
      45             : GEN
      46        1112 : FpE_changepoint(GEN x, GEN ch, GEN p)
      47             : {
      48        1112 :   pari_sp av = avma;
      49             :   GEN p1,z,u,r,s,t,v,v2,v3;
      50        1112 :   if (ell_is_inf(x)) return x;
      51        1056 :   u = gel(ch,1); r = gel(ch,2);
      52        1056 :   s = gel(ch,3); t = gel(ch,4);
      53        1056 :   v = Fp_inv(u, p); v2 = Fp_sqr(v,p); v3 = Fp_mul(v,v2,p);
      54        1056 :   p1 = Fp_sub(gel(x,1),r,p);
      55        1056 :   z = cgetg(3,t_VEC);
      56        1056 :   gel(z,1) = Fp_mul(v2, p1, p);
      57        1056 :   gel(z,2) = Fp_mul(v3, Fp_sub(gel(x,2), Fp_add(Fp_mul(s,p1, p),t, p),p),p);
      58        1056 :   return gerepileupto(av, z);
      59             : }
      60             : 
      61             : GEN
      62        2180 : FpE_changepointinv(GEN x, GEN ch, GEN p)
      63             : {
      64             :   GEN u, r, s, t, X, Y, u2, u3, u2X, z;
      65        2180 :   if (ell_is_inf(x)) return x;
      66        2180 :   X = gel(x,1); Y = gel(x,2);
      67        2180 :   u = gel(ch,1); r = gel(ch,2);
      68        2180 :   s = gel(ch,3); t = gel(ch,4);
      69        2180 :   u2 = Fp_sqr(u, p); u3 = Fp_mul(u,u2,p);
      70        2179 :   u2X = Fp_mul(u2,X, p);
      71        2179 :   z = cgetg(3, t_VEC);
      72        2180 :   gel(z,1) = Fp_add(u2X,r,p);
      73        2182 :   gel(z,2) = Fp_add(Fp_mul(u3,Y,p), Fp_add(Fp_mul(s,u2X,p), t, p), p);
      74        2179 :   return z;
      75             : }
      76             : 
      77             : static GEN
      78         420 : nonsquare_Fp(GEN p)
      79             : {
      80         420 :   pari_sp av = avma;
      81             :   GEN a;
      82             :   do
      83             :   {
      84         833 :     avma = av;
      85         833 :     a = randomi(p);
      86         833 :   } while (kronecker(a, p) >= 0);
      87         420 :   return a;
      88             : }
      89             : 
      90             : void
      91           0 : Fp_elltwist(GEN a4, GEN a6, GEN p, GEN *pt_a4, GEN *pt_a6)
      92             : {
      93           0 :   GEN d = nonsquare_Fp(p), d2 = Fp_sqr(d, p), d3 = Fp_mul(d2, d, p);
      94           0 :   *pt_a4 = Fp_mul(a4, d2, p);
      95           0 :   *pt_a6 = Fp_mul(a6, d3, p);
      96           0 : }
      97             : 
      98             : static GEN
      99     2021548 : FpE_dbl_slope(GEN P, GEN a4, GEN p, GEN *slope)
     100             : {
     101             :   GEN x, y, Q;
     102     2021548 :   if (ell_is_inf(P) || !signe(gel(P,2))) return ellinf();
     103     2004865 :   x = gel(P,1); y = gel(P,2);
     104     2004865 :   *slope = Fp_div(Fp_add(Fp_mulu(Fp_sqr(x,p), 3, p), a4, p),
     105             :                   Fp_mulu(y, 2, p), p);
     106     1980474 :   Q = cgetg(3,t_VEC);
     107     1983838 :   gel(Q, 1) = Fp_sub(Fp_sqr(*slope, p), Fp_mulu(x, 2, p), p);
     108     1994591 :   gel(Q, 2) = Fp_sub(Fp_mul(*slope, Fp_sub(x, gel(Q, 1), p), p), y, p);
     109     1991232 :   return Q;
     110             : }
     111             : 
     112             : GEN
     113     2005526 : FpE_dbl(GEN P, GEN a4, GEN p)
     114             : {
     115     2005526 :   pari_sp av = avma;
     116             :   GEN slope;
     117     2005526 :   return gerepileupto(av, FpE_dbl_slope(P,a4,p,&slope));
     118             : }
     119             : 
     120             : static GEN
     121     1281239 : FpE_add_slope(GEN P, GEN Q, GEN a4, GEN p, GEN *slope)
     122             : {
     123             :   GEN Px, Py, Qx, Qy, R;
     124     1281239 :   if (ell_is_inf(P)) return Q;
     125     1279717 :   if (ell_is_inf(Q)) return P;
     126     1279795 :   Px = gel(P,1); Py = gel(P,2);
     127     1279795 :   Qx = gel(Q,1); Qy = gel(Q,2);
     128     1279795 :   if (equalii(Px, Qx))
     129             :   {
     130        6080 :     if (equalii(Py, Qy))
     131         622 :       return FpE_dbl_slope(P, a4, p, slope);
     132             :     else
     133        5458 :       return ellinf();
     134             :   }
     135     1273559 :   *slope = Fp_div(Fp_sub(Py, Qy, p), Fp_sub(Px, Qx, p), p);
     136     1273340 :   R = cgetg(3,t_VEC);
     137     1273380 :   gel(R, 1) = Fp_sub(Fp_sub(Fp_sqr(*slope, p), Px, p), Qx, p);
     138     1273497 :   gel(R, 2) = Fp_sub(Fp_mul(*slope, Fp_sub(Px, gel(R, 1), p), p), Py, p);
     139     1273560 :   return R;
     140             : }
     141             : 
     142             : GEN
     143     1278527 : FpE_add(GEN P, GEN Q, GEN a4, GEN p)
     144             : {
     145     1278527 :   pari_sp av = avma;
     146             :   GEN slope;
     147     1278527 :   return gerepileupto(av, FpE_add_slope(P,Q,a4,p,&slope));
     148             : }
     149             : 
     150             : static GEN
     151           0 : FpE_neg_i(GEN P, GEN p)
     152             : {
     153           0 :   if (ell_is_inf(P)) return P;
     154           0 :   return mkvec2(gel(P,1), Fp_neg(gel(P,2), p));
     155             : }
     156             : 
     157             : GEN
     158      372582 : FpE_neg(GEN P, GEN p)
     159             : {
     160      372582 :   if (ell_is_inf(P)) return ellinf();
     161      372582 :   return mkvec2(gcopy(gel(P,1)), Fp_neg(gel(P,2), p));
     162             : }
     163             : 
     164             : GEN
     165           0 : FpE_sub(GEN P, GEN Q, GEN a4, GEN p)
     166             : {
     167           0 :   pari_sp av = avma;
     168             :   GEN slope;
     169           0 :   return gerepileupto(av, FpE_add_slope(P, FpE_neg_i(Q, p), a4, p, &slope));
     170             : }
     171             : 
     172             : struct _FpE
     173             : {
     174             :   GEN a4,a6;
     175             :   GEN p;
     176             : };
     177             : 
     178             : static GEN
     179     2006275 : _FpE_dbl(void *E, GEN P)
     180             : {
     181     2006275 :   struct _FpE *ell = (struct _FpE *) E;
     182     2006275 :   return FpE_dbl(P, ell->a4, ell->p);
     183             : }
     184             : 
     185             : static GEN
     186     1259437 : _FpE_add(void *E, GEN P, GEN Q)
     187             : {
     188     1259437 :   struct _FpE *ell=(struct _FpE *) E;
     189     1259437 :   return FpE_add(P, Q, ell->a4, ell->p);
     190             : }
     191             : 
     192             : static GEN
     193      488016 : _FpE_mul(void *E, GEN P, GEN n)
     194             : {
     195      488016 :   pari_sp av = avma;
     196      488016 :   struct _FpE *e=(struct _FpE *) E;
     197      488016 :   long s = signe(n);
     198      488016 :   if (!s || ell_is_inf(P)) return ellinf();
     199      487988 :   if (s<0) P = FpE_neg(P, e->p);
     200      487988 :   if (is_pm1(n)) return s>0? gcopy(P): P;
     201       90926 :   return gerepileupto(av, gen_pow(P, n, e, &_FpE_dbl, &_FpE_add));
     202             : }
     203             : 
     204             : GEN
     205         688 : FpE_mul(GEN P, GEN n, GEN a4, GEN p)
     206             : {
     207             :   struct _FpE E;
     208         688 :   E.a4= a4; E.p = p;
     209         688 :   return _FpE_mul(&E, P, n);
     210             : }
     211             : 
     212             : /* Finds a random non-singular point on E */
     213             : 
     214             : GEN
     215       29029 : random_FpE(GEN a4, GEN a6, GEN p)
     216             : {
     217       29029 :   pari_sp ltop = avma;
     218             :   GEN x, x2, y, rhs;
     219             :   do
     220             :   {
     221       51812 :     avma= ltop;
     222       51812 :     x   = randomi(p); /*  x^3+a4*x+a6 = x*(x^2+a4)+a6  */
     223       51812 :     x2  = Fp_sqr(x, p);
     224       51812 :     rhs = Fp_add(Fp_mul(x, Fp_add(x2, a4, p), p), a6, p);
     225       60425 :   } while ((!signe(rhs) && !signe(Fp_add(Fp_mulu(x2,3,p),a4,p)))
     226      103624 :           || kronecker(rhs, p) < 0);
     227       29029 :   y = Fp_sqrt(rhs, p);
     228       29029 :   if (!y) pari_err_PRIME("random_FpE", p);
     229       29029 :   return gerepilecopy(ltop, mkvec2(x, y));
     230             : }
     231             : 
     232             : static GEN
     233       29001 : _FpE_rand(void *E)
     234             : {
     235       29001 :   struct _FpE *e=(struct _FpE *) E;
     236       29001 :   return random_FpE(e->a4, e->a6, e->p);
     237             : }
     238             : 
     239             : static const struct bb_group FpE_group={_FpE_add,_FpE_mul,_FpE_rand,hash_GEN,ZV_equal,ell_is_inf,NULL};
     240             : 
     241             : const struct bb_group *
     242         840 : get_FpE_group(void ** pt_E, GEN a4, GEN a6, GEN p)
     243             : {
     244         840 :   struct _FpE *e = (struct _FpE *) stack_malloc(sizeof(struct _FpE));
     245         840 :   e->a4 = a4; e->a6 = a6; e->p  = p;
     246         840 :   *pt_E = (void *) e;
     247         840 :   return &FpE_group;
     248             : }
     249             : 
     250             : GEN
     251         819 : FpE_order(GEN z, GEN o, GEN a4, GEN p)
     252             : {
     253         819 :   pari_sp av = avma;
     254             :   struct _FpE e;
     255             :   GEN r;
     256         819 :   if (lgefint(p) == 3)
     257             :   {
     258         713 :     ulong pp = p[2];
     259         713 :     r = Fle_order(ZV_to_Flv(z, pp), o, umodiu(a4,pp), pp);
     260             :   }
     261             :   else
     262             :   {
     263         106 :     e.a4 = a4;
     264         106 :     e.p = p;
     265         106 :     r = gen_order(z, o, (void*)&e, &FpE_group);
     266             :   }
     267         819 :   return gerepileuptoint(av, r);
     268             : }
     269             : 
     270             : GEN
     271          42 : FpE_log(GEN a, GEN b, GEN o, GEN a4, GEN p)
     272             : {
     273          42 :   pari_sp av = avma;
     274             :   struct _FpE e;
     275             :   GEN r;
     276          42 :   if (lgefint(p) == 3)
     277             :   {
     278          42 :     ulong pp = p[2];
     279          42 :     r = Fle_log(ZV_to_Flv(a,pp), ZV_to_Flv(b,pp), o, umodiu(a4,pp), pp);
     280             :   }
     281             :   else
     282             :   {
     283           0 :     e.a4 = a4;
     284           0 :     e.p = p;
     285           0 :     r = gen_PH_log(a, b, o, (void*)&e, &FpE_group);
     286             :   }
     287          42 :   return gerepileuptoint(av, r);
     288             : }
     289             : 
     290             : /***********************************************************************/
     291             : /**                                                                   **/
     292             : /**                            Pairings                               **/
     293             : /**                                                                   **/
     294             : /***********************************************************************/
     295             : 
     296             : /* Derived from APIP from and by Jerome Milan, 2012 */
     297             : 
     298             : static GEN
     299       56323 : FpE_vert(GEN P, GEN Q, GEN a4, GEN p)
     300             : {
     301       56323 :   if (ell_is_inf(P))
     302       19906 :     return gen_1;
     303       36417 :   if (!equalii(gel(Q, 1), gel(P, 1)))
     304       33834 :     return Fp_sub(gel(Q, 1), gel(P, 1), p);
     305        2583 :   if (signe(gel(P,2))!=0) return gen_1;
     306        2212 :   return Fp_inv(Fp_add(Fp_mulu(Fp_sqr(gel(P,1),p), 3, p), a4, p), p);
     307             : }
     308             : 
     309             : static GEN
     310       20159 : FpE_Miller_line(GEN R, GEN Q, GEN slope, GEN a4, GEN p)
     311             : {
     312       20159 :   GEN x = gel(Q, 1), y = gel(Q, 2);
     313       20159 :   GEN tmp1 = Fp_sub(x, gel(R, 1), p);
     314       20159 :   GEN tmp2 = Fp_add(Fp_mul(tmp1, slope, p), gel(R,2), p);
     315       20159 :   if (!equalii(y, tmp2))
     316       18731 :     return Fp_sub(y, tmp2, p);
     317        1428 :   if (signe(y) == 0)
     318        1092 :     return gen_1;
     319             :   else
     320             :   {
     321             :     GEN s1, s2;
     322         336 :     GEN y2i = Fp_inv(Fp_mulu(y, 2, p), p);
     323         336 :     s1 = Fp_mul(Fp_add(Fp_mulu(Fp_sqr(x, p), 3, p), a4, p), y2i, p);
     324         336 :     if (!equalii(s1, slope))
     325         196 :       return Fp_sub(s1, slope, p);
     326         140 :     s2 = Fp_mul(Fp_sub(Fp_mulu(x, 3, p), Fp_sqr(s1, p), p), y2i, p);
     327         140 :     return signe(s2)!=0 ? s2: y2i;
     328             :   }
     329             : }
     330             : 
     331             : /* Computes the equation of the line tangent to R and returns its
     332             :    evaluation at the point Q. Also doubles the point R.
     333             :  */
     334             : 
     335             : static GEN
     336       35571 : FpE_tangent_update(GEN R, GEN Q, GEN a4, GEN p, GEN *pt_R)
     337             : {
     338       35571 :   if (ell_is_inf(R))
     339             :   {
     340        4166 :     *pt_R = ellinf();
     341        4166 :     return gen_1;
     342             :   }
     343       31405 :   else if (signe(gel(R,2)) == 0)
     344             :   {
     345       13941 :     *pt_R = ellinf();
     346       13941 :     return FpE_vert(R, Q, a4, p);
     347             :   } else {
     348             :     GEN slope;
     349       17464 :     *pt_R = FpE_dbl_slope(R, a4, p, &slope);
     350       17464 :     return FpE_Miller_line(R, Q, slope, a4, p);
     351             :   }
     352             : }
     353             : 
     354             : /* Computes the equation of the line through R and P, and returns its
     355             :    evaluation at the point Q. Also adds P to the point R.
     356             :  */
     357             : 
     358             : static GEN
     359        4753 : FpE_chord_update(GEN R, GEN P, GEN Q, GEN a4, GEN p, GEN *pt_R)
     360             : {
     361        4753 :   if (ell_is_inf(R))
     362             :   {
     363         259 :     *pt_R = gcopy(P);
     364         259 :     return FpE_vert(P, Q, a4, p);
     365             :   }
     366        4494 :   else if (ell_is_inf(P))
     367             :   {
     368           0 :     *pt_R = gcopy(R);
     369           0 :     return FpE_vert(R, Q, a4, p);
     370             :   }
     371        4494 :   else if (equalii(gel(P, 1), gel(R, 1)))
     372             :   {
     373        1799 :     if (equalii(gel(P, 2), gel(R, 2)))
     374           0 :       return FpE_tangent_update(R, Q, a4, p, pt_R);
     375             :     else {
     376        1799 :       *pt_R = ellinf();
     377        1799 :       return FpE_vert(R, Q, a4, p);
     378             :     }
     379             :   } else {
     380             :     GEN slope;
     381        2695 :     *pt_R = FpE_add_slope(P, R, a4, p, &slope);
     382        2695 :     return FpE_Miller_line(R, Q, slope, a4, p);
     383             :   }
     384             : }
     385             : 
     386             : /* Returns the Miller function f_{m, Q} evaluated at the point P using
     387             :    the standard Miller algorithm.
     388             :  */
     389             : 
     390             : struct _FpE_miller
     391             : {
     392             :   GEN p, a4, P;
     393             : };
     394             : 
     395             : static GEN
     396       35571 : FpE_Miller_dbl(void* E, GEN d)
     397             : {
     398       35571 :   struct _FpE_miller *m = (struct _FpE_miller *)E;
     399       35571 :   GEN p = m->p, a4 = m->a4, P = m->P;
     400             :   GEN v, line;
     401       35571 :   GEN num = Fp_sqr(gel(d,1), p);
     402       35571 :   GEN denom = Fp_sqr(gel(d,2), p);
     403       35571 :   GEN point = gel(d,3);
     404       35571 :   line = FpE_tangent_update(point, P, a4, p, &point);
     405       35571 :   num  = Fp_mul(num, line, p);
     406       35571 :   v = FpE_vert(point, P, a4, p);
     407       35571 :   denom = Fp_mul(denom, v, p);
     408       35571 :   return mkvec3(num, denom, point);
     409             : }
     410             : 
     411             : static GEN
     412        4753 : FpE_Miller_add(void* E, GEN va, GEN vb)
     413             : {
     414        4753 :   struct _FpE_miller *m = (struct _FpE_miller *)E;
     415        4753 :   GEN p = m->p, a4= m->a4, P = m->P;
     416             :   GEN v, line, point;
     417        4753 :   GEN na = gel(va,1), da = gel(va,2), pa = gel(va,3);
     418        4753 :   GEN nb = gel(vb,1), db = gel(vb,2), pb = gel(vb,3);
     419        4753 :   GEN num   = Fp_mul(na, nb, p);
     420        4753 :   GEN denom = Fp_mul(da, db, p);
     421        4753 :   line = FpE_chord_update(pa, pb, P, a4, p, &point);
     422        4753 :   num  = Fp_mul(num, line, p);
     423        4753 :   v = FpE_vert(point, P, a4, p);
     424        4753 :   denom = Fp_mul(denom, v, p);
     425        4753 :   return mkvec3(num, denom, point);
     426             : }
     427             : 
     428             : static GEN
     429       15481 : FpE_Miller(GEN Q, GEN P, GEN m, GEN a4, GEN p)
     430             : {
     431       15481 :   pari_sp ltop = avma;
     432             :   struct _FpE_miller d;
     433             :   GEN v, num, denom;
     434             : 
     435       15481 :   d.a4 = a4; d.p = p; d.P = P;
     436       15481 :   v = gen_pow(mkvec3(gen_1,gen_1,Q), m, (void*)&d, FpE_Miller_dbl, FpE_Miller_add);
     437       15481 :   num = gel(v,1); denom = gel(v,2);
     438       15481 :   return gerepileupto(ltop, Fp_div(num, denom, p));
     439             : }
     440             : 
     441             : GEN
     442       10660 : FpE_weilpairing(GEN P, GEN Q, GEN m, GEN a4, GEN p)
     443             : {
     444       10660 :   pari_sp ltop = avma;
     445             :   GEN num, denom, result;
     446       10660 :   if (ell_is_inf(P) || ell_is_inf(Q) || ZV_equal(P,Q))
     447        3021 :     return gen_1;
     448        7639 :   num    = FpE_Miller(P, Q, m, a4, p);
     449        7639 :   denom  = FpE_Miller(Q, P, m, a4, p);
     450        7639 :   result = Fp_div(num, denom, p);
     451        7639 :   if (mpodd(m))
     452         658 :     result  = Fp_neg(result, p);
     453        7639 :   return gerepileupto(ltop, result);
     454             : }
     455             : 
     456             : GEN
     457         203 : FpE_tatepairing(GEN P, GEN Q, GEN m, GEN a4, GEN p)
     458             : {
     459         203 :   if (ell_is_inf(P) || ell_is_inf(Q))
     460           0 :     return gen_1;
     461         203 :   return FpE_Miller(P, Q, m, a4, p);
     462             : }
     463             : 
     464             : /***********************************************************************/
     465             : /**                                                                   **/
     466             : /**                   CM by principal order                           **/
     467             : /**                                                                   **/
     468             : /***********************************************************************/
     469             : 
     470             : /* is jn/jd = J (mod p) */
     471             : static int
     472      927675 : is_CMj(long J, GEN jn, GEN jd, GEN p)
     473      927675 : { return remii(subii(mulis(jd,J), jn), p) == gen_0; }
     474             : #ifndef LONG_IS_64BIT
     475             : /* is jn/jd = -(2^32 a + b) (mod p) */
     476             : static int
     477       20538 : u2_is_CMj(ulong a, ulong b, GEN jn, GEN jd, GEN p)
     478             : {
     479       20538 :   GEN mJ = uu32toi(a,b);
     480       20538 :   return remii(addii(mulii(jd,mJ), jn), p) == gen_0;
     481             : }
     482             : #endif
     483             : 
     484             : static long
     485       73836 : Fp_ellj_get_CM(GEN jn, GEN jd, GEN p)
     486             : {
     487             : #define CHECK(CM,J) if (is_CMj(J,jn,jd,p)) return CM;
     488       73836 :   CHECK(-3,  0);
     489       73794 :   CHECK(-4,  1728);
     490       73745 :   CHECK(-7,  -3375);
     491       73584 :   CHECK(-8,  8000);
     492       73416 :   CHECK(-11, -32768);
     493       73255 :   CHECK(-12, 54000);
     494       72961 :   CHECK(-16, 287496);
     495       72800 :   CHECK(-19, -884736);
     496       72625 :   CHECK(-27, -12288000);
     497       72317 :   CHECK(-28, 16581375);
     498       72114 :   CHECK(-43, -884736000);
     499             : #ifdef LONG_IS_64BIT
     500       61668 :   CHECK(-67, -147197952000L);
     501       61560 :   CHECK(-163, -262537412640768000L);
     502             : #else
     503       10278 :   if (u2_is_CMj(0x00000022UL,0x45ae8000UL,jn,jd,p)) return -67;
     504       10260 :   if (u2_is_CMj(0x03a4b862UL,0xc4b40000UL,jn,jd,p)) return -163;
     505             : #endif
     506             : #undef CHECK
     507       71659 :   return 0;
     508             : }
     509             : 
     510             : /***********************************************************************/
     511             : /**                                                                   **/
     512             : /**                            issupersingular                        **/
     513             : /**                                                                   **/
     514             : /***********************************************************************/
     515             : 
     516             : /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
     517             : static GEN
     518        5691 : FqX_quad_root(GEN x, GEN T, GEN p)
     519             : {
     520        5691 :   GEN b = gel(x,3), c = gel(x,2);
     521        5691 :   GEN D = Fq_sub(Fq_sqr(b, T, p), Fq_mulu(c,4, T, p), T, p);
     522        5691 :   GEN s = Fq_sqrt(D,T, p);
     523        5691 :   if (!s) return NULL;
     524        3374 :   return Fq_Fp_mul(Fq_sub(s, b, T, p), shifti(addis(p, 1),-1),T, p);
     525             : }
     526             : 
     527             : /*
     528             :  * pol is the modular polynomial of level 2 modulo p.
     529             :  *
     530             :  * (T, p) defines the field FF_{p^2} in which j_prev and j live.
     531             :  */
     532             : static long
     533        2583 : path_extends_to_floor(GEN j_prev, GEN j, GEN T, GEN p, GEN Phi2, ulong max_len)
     534             : {
     535        2583 :   pari_sp ltop = avma;
     536             :   GEN Phi2_j;
     537             :   ulong mult, d;
     538             : 
     539             :   /* A path made its way to the floor if (i) its length was cut off
     540             :    * before reaching max_path_len, or (ii) it reached max_path_len but
     541             :    * only has one neighbour. */
     542        5957 :   for (d = 1; d < max_len; ++d) {
     543             :     GEN j_next;
     544             : 
     545        5691 :     Phi2_j = FqX_div_by_X_x(FqXY_evalx(Phi2, j, T, p), j_prev, T, p, NULL);
     546        5691 :     j_next = FqX_quad_root(Phi2_j, T, p);
     547        5691 :     if (!j_next)
     548             :     { /* j is on the floor */
     549        2317 :       avma = ltop;
     550        2317 :       return 1;
     551             :     }
     552             : 
     553        3374 :     j_prev = j; j = j_next;
     554        3374 :     if (gc_needed(ltop, 2))
     555           0 :       gerepileall(ltop, 2, &j, &j_prev);
     556             :   }
     557             : 
     558             :   /* Check that we didn't end up at the floor on the last step (j will
     559             :    * point to the last element in the path. */
     560         266 :   Phi2_j = FqX_div_by_X_x(FqXY_evalx(Phi2, j, T, p), j_prev, T, p, NULL);
     561         266 :   mult = FqX_nbroots(Phi2_j, T, p);
     562         266 :   avma = ltop;
     563         266 :   return mult == 0;
     564             : }
     565             : 
     566             : static int
     567       13860 : jissupersingular(GEN j, GEN S, GEN p)
     568             : {
     569       13860 :   long max_path_len = expi(p)+1;
     570       13860 :   GEN Phi2 = FpXX_red(polmodular_ZXX(2,0,0,1), p);
     571       13860 :   GEN Phi2_j = FqXY_evalx(Phi2, j, S, p);
     572       13860 :   GEN roots = FpXQX_roots(Phi2_j, S, p);
     573       13860 :   long nbroots = lg(roots)-1;
     574       13860 :   int res = 1;
     575             : 
     576             :   /* Every node in a supersingular L-volcano has L + 1 neighbours. */
     577             :   /* Note: a multiple root only occur when j has CM by sqrt(-15). */
     578       13860 :   if (nbroots==0 || (nbroots==1 && FqX_is_squarefree(Phi2_j, S, p)))
     579       11431 :     res = 0;
     580             :   else {
     581        2429 :     long i, l = lg(roots);
     582        2597 :     for (i = 1; i < l; ++i) {
     583        2583 :       if (path_extends_to_floor(j, gel(roots, i), S, p, Phi2, max_path_len)) {
     584        2415 :         res = 0;
     585        2415 :         break;
     586             :       }
     587             :     }
     588             :   }
     589             :   /* If none of the paths reached the floor, then the j-invariant is
     590             :    * supersingular. */
     591       13860 :   return res;
     592             : }
     593             : 
     594             : int
     595        1050 : Fp_elljissupersingular(GEN j, GEN p)
     596             : {
     597        1050 :   pari_sp ltop = avma;
     598             :   long CM;
     599        1050 :   if (abscmpiu(p, 5) <= 0) return signe(j) == 0; /* valid if p <= 5 */
     600         938 :   CM = Fp_ellj_get_CM(j, gen_1, p);
     601         938 :   if (CM < 0) return krosi(CM, p) < 0; /* valid if p > 3 */
     602             :   else
     603             :   {
     604         609 :     GEN S = init_Fq(p, 2, fetch_var());
     605         609 :     int res = jissupersingular(j, S, p);
     606         609 :     (void)delete_var(); avma = ltop; return res;
     607             :   }
     608             : }
     609             : 
     610             : /***********************************************************************/
     611             : /**                                                                   **/
     612             : /**                            Cardinal                               **/
     613             : /**                                                                   **/
     614             : /***********************************************************************/
     615             : 
     616             : /*assume a4,a6 reduced mod p odd */
     617             : static ulong
     618      208666 : Fl_elltrace_naive(ulong a4, ulong a6, ulong p)
     619             : {
     620             :   ulong i, j;
     621      208666 :   long a = 0;
     622             :   long d0, d1, d2, d3;
     623      208666 :   GEN k = const_vecsmall(p, -1);
     624      208666 :   k[1] = 0;
     625    61920368 :   for (i=1, j=1; i < p; i += 2, j = Fl_add(j, i, p))
     626    61711702 :     k[j+1] = 1;
     627      208666 :   d0 = 6%p; d1 = d0; d2 = Fl_add(a4, 1, p); d3 = a6;
     628   123632070 :   for(i=0;; i++)
     629             :   {
     630   123632070 :     a -= k[1+d3];
     631   123632070 :     if (i==p-1) break;
     632   123423404 :     d3 = Fl_add(d3, d2, p);
     633   123423404 :     d2 = Fl_add(d2, d1, p);
     634   123423404 :     d1 = Fl_add(d1, d0, p);
     635   123423404 :   }
     636      208666 :   return a;
     637             : }
     638             : 
     639             : /* z1 <-- z1 + z2, with precomputed inverse */
     640             : static void
     641      305362 : FpE_add_ip(GEN z1, GEN z2, GEN a4, GEN p, GEN p2inv)
     642             : {
     643             :   GEN p1,x,x1,x2,y,y1,y2;
     644             : 
     645      305362 :   x1 = gel(z1,1); y1 = gel(z1,2);
     646      305362 :   x2 = gel(z2,1); y2 = gel(z2,2);
     647      305362 :   if (x1 == x2)
     648          66 :     p1 = Fp_add(a4, mulii(x1,mului(3,x1)), p);
     649             :   else
     650      305296 :     p1 = Fp_sub(y2,y1, p);
     651             : 
     652      305362 :   p1 = Fp_mul(p1, p2inv, p);
     653      305362 :   x = Fp_sub(sqri(p1), addii(x1,x2), p);
     654      305362 :   y = Fp_sub(mulii(p1,subii(x1,x)), y1, p);
     655      305362 :   affii(x, x1);
     656      305362 :   affii(y, y1);
     657      305362 : }
     658             : 
     659             : /* make sure *x has lgefint >= k */
     660             : static void
     661       18872 : _fix(GEN x, long k)
     662             : {
     663       18872 :   GEN y = (GEN)*x;
     664       18872 :   if (lgefint(y) < k) { GEN p1 = cgeti(k); affii(y,p1); *x = (long)p1; }
     665       18872 : }
     666             : 
     667             : /* Return the lift of a (mod b), which is closest to c */
     668             : static GEN
     669      410200 : closest_lift(GEN a, GEN b, GEN c)
     670             : {
     671      410200 :   return addii(a, mulii(b, diviiround(subii(c,a), b)));
     672             : }
     673             : 
     674             : static long
     675          77 : get_table_size(GEN pordmin, GEN B)
     676             : {
     677          77 :   pari_sp av = avma;
     678          77 :   GEN t = ceilr( sqrtr( divri(itor(pordmin, DEFAULTPREC), B) ) );
     679          77 :   if (is_bigint(t))
     680           0 :     pari_err_OVERFLOW("ellap [large prime: install the 'seadata' package]");
     681          77 :   avma = av;
     682          77 :   return itos(t) >> 1;
     683             : }
     684             : 
     685             : /* Find x such that kronecker(u = x^3+c4x+c6, p) is KRO.
     686             :  * Return point [x*u,u^2] on E (KRO=1) / E^twist (KRO=-1) */
     687             : static GEN
     688           0 : Fp_ellpoint(long KRO, ulong *px, GEN c4, GEN c6, GEN p)
     689             : {
     690           0 :   ulong x = *px;
     691             :   GEN u;
     692             :   for(;;)
     693             :   {
     694           0 :     x++; /* u = x^3 + c4 x + c6 */
     695           0 :     u = modii(addii(c6, mului(x, addii(c4, sqru(x)))), p);
     696           0 :     if (kronecker(u,p) == KRO) break;
     697           0 :   }
     698           0 :   *px = x;
     699           0 :   return mkvec2(modii(mului(x,u),p), Fp_sqr(u,p));
     700             : }
     701             : static GEN
     702        7021 : Fl_ellpoint(long KRO, ulong *px, ulong c4, ulong c6, ulong p)
     703             : {
     704        7021 :   ulong t, u, x = *px;
     705             :   for(;;)
     706             :   {
     707       13335 :     if (++x >= p) pari_err_PRIME("ellap",utoi(p));
     708       13335 :     t = Fl_add(c4, Fl_sqr(x,p), p);
     709       13335 :     u = Fl_add(c6, Fl_mul(x, t, p), p);
     710       13335 :     if (krouu(u,p) == KRO) break;
     711        6314 :   }
     712        7021 :   *px = x;
     713        7021 :   return mkvecsmall2(Fl_mul(x,u,p), Fl_sqr(u,p));
     714             : }
     715             : 
     716             : static GEN ap_j1728(GEN a4,GEN p);
     717             : /* compute a_p using Shanks/Mestre + Montgomery's trick. Assume p > 457 */
     718             : static GEN
     719          77 : Fp_ellcard_Shanks(GEN c4, GEN c6, GEN p)
     720             : {
     721             :   pari_timer T;
     722             :   long *tx, *ty, *ti, pfinal, i, j, s, KRO, nb;
     723             :   ulong x;
     724          77 :   pari_sp av = avma, av2;
     725             :   GEN p1, P, mfh, h, F,f, fh,fg, pordmin, u, v, p1p, p2p, A, B, a4, pts;
     726          77 :   tx = NULL;
     727          77 :   ty = ti = NULL; /* gcc -Wall */
     728             : 
     729          77 :   if (!signe(c6)) {
     730           0 :     GEN ap = ap_j1728(c4, p);
     731           0 :     return gerepileuptoint(av, subii(addiu(p,1), ap));
     732             :   }
     733             : 
     734          77 :   if (DEBUGLEVEL >= 6) timer_start(&T);
     735             :   /* once #E(Fp) is know mod B >= pordmin, it is completely determined */
     736          77 :   pordmin = addis(sqrti(gmul2n(p,4)), 1); /* ceil( 4sqrt(p) ) */
     737          77 :   p1p = addsi(1, p);
     738          77 :   p2p = shifti(p1p, 1);
     739          77 :   x = 0; KRO = 0;
     740             :   /* how many 2-torsion points ? */
     741          77 :   switch(FpX_nbroots(mkpoln(4, gen_1, gen_0, c4, c6), p))
     742             :   {
     743           9 :     case 3:  A = gen_0; B = utoipos(4); break;
     744          31 :     case 1:  A = gen_0; B = gen_2; break;
     745          37 :     default: A = gen_1; B = gen_2; break; /* 0 */
     746             :   }
     747             :   for(;;)
     748             :   {
     749          77 :     h = closest_lift(A, B, p1p);
     750          77 :     if (!KRO) /* first time, initialize */
     751             :     {
     752          77 :       KRO = kronecker(c6,p);
     753          77 :       f = mkvec2(gen_0, Fp_sqr(c6,p));
     754             :     }
     755             :     else
     756             :     {
     757           0 :       KRO = -KRO;
     758           0 :       f = Fp_ellpoint(KRO, &x, c4,c6,p);
     759             :     }
     760             :     /* [ux, u^2] is on E_u: y^2 = x^3 + c4 u^2 x + c6 u^3
     761             :      * E_u isomorphic to E (resp. E') iff KRO = 1 (resp. -1)
     762             :      * #E(F_p) = p+1 - a_p, #E'(F_p) = p+1 + a_p
     763             :      *
     764             :      * #E_u(Fp) = A (mod B),  h is close to #E_u(Fp) */
     765          77 :     a4 = modii(mulii(c4, gel(f,2)), p); /* c4 for E_u */
     766          77 :     fh = FpE_mul(f, h, a4, p);
     767          77 :     if (ell_is_inf(fh)) goto FOUND;
     768             : 
     769          77 :     s = get_table_size(pordmin, B);
     770             :     /* look for h s.t f^h = 0 */
     771          77 :     if (!tx)
     772             :     { /* first time: initialize */
     773          77 :       tx = newblock(3*(s+1));
     774          77 :       ty = tx + (s+1);
     775          77 :       ti = ty + (s+1);
     776             :     }
     777          77 :     F = FpE_mul(f,B,a4,p);
     778          77 :     *tx = evaltyp(t_VECSMALL) | evallg(s+1);
     779             : 
     780             :     /* F = B.f */
     781          77 :     P = gcopy(fh);
     782          77 :     if (s < 3)
     783             :     { /* we're nearly done: naive search */
     784           0 :       GEN q1 = P, mF = FpE_neg(F, p); /* -F */
     785           0 :       for (i=1;; i++)
     786             :       {
     787           0 :         P = FpE_add(P,F,a4,p); /* h.f + i.F */
     788           0 :         if (ell_is_inf(P)) { h = addii(h, mului(i,B)); goto FOUND; }
     789           0 :         q1 = FpE_add(q1,mF,a4,p); /* h.f - i.F */
     790           0 :         if (ell_is_inf(q1)) { h = subii(h, mului(i,B)); goto FOUND; }
     791           0 :       }
     792             :     }
     793             :     /* Baby Step/Giant Step */
     794          77 :     nb = minss(128, s >> 1); /* > 0. Will do nb pts at a time: faster inverse */
     795          77 :     pts = cgetg(nb+1, t_VEC);
     796          77 :     j = lgefint(p);
     797        9513 :     for (i=1; i<=nb; i++)
     798             :     { /* baby steps */
     799        9436 :       gel(pts,i) = P; /* h.f + (i-1).F */
     800        9436 :       _fix(P+1, j); tx[i] = mod2BIL(gel(P,1));
     801        9436 :       _fix(P+2, j); ty[i] = mod2BIL(gel(P,2));
     802        9436 :       P = FpE_add(P,F,a4,p); /* h.f + i.F */
     803        9436 :       if (ell_is_inf(P)) { h = addii(h, mului(i,B)); goto FOUND; }
     804             :     }
     805          77 :     mfh = FpE_neg(fh, p);
     806          77 :     fg = FpE_add(P,mfh,a4,p); /* h.f + nb.F - h.f = nb.F */
     807          77 :     if (ell_is_inf(fg)) { h = mului(nb,B); goto FOUND; }
     808          77 :     u = cgetg(nb+1, t_VEC);
     809          77 :     av2 = avma; /* more baby steps, nb points at a time */
     810        1431 :     while (i <= s)
     811             :     {
     812             :       long maxj;
     813      164151 :       for (j=1; j<=nb; j++) /* adding nb.F (part 1) */
     814             :       {
     815      162874 :         P = gel(pts,j); /* h.f + (i-nb-1+j-1).F */
     816      162874 :         gel(u,j) = subii(gel(fg,1), gel(P,1));
     817      162874 :         if (!signe(gel(u,j))) /* sum = 0 or doubling */
     818             :         {
     819           1 :           long k = i+j-2;
     820           1 :           if (equalii(gel(P,2),gel(fg,2))) k -= 2*nb; /* fg == P */
     821           1 :           h = addii(h, mulsi(k,B)); goto FOUND;
     822             :         }
     823             :       }
     824        1277 :       v = FpV_inv(u, p);
     825        1277 :       maxj = (i-1 + nb <= s)? nb: s % nb;
     826      160461 :       for (j=1; j<=maxj; j++,i++) /* adding nb.F (part 2) */
     827             :       {
     828      159184 :         P = gel(pts,j);
     829      159184 :         FpE_add_ip(P,fg, a4,p, gel(v,j));
     830      159184 :         tx[i] = mod2BIL(gel(P,1));
     831      159184 :         ty[i] = mod2BIL(gel(P,2));
     832             :       }
     833        1277 :       avma = av2;
     834             :     }
     835          76 :     P = FpE_add(gel(pts,j-1),mfh,a4,p); /* = (s-1).F */
     836          76 :     if (ell_is_inf(P)) { h = mului(s-1,B); goto FOUND; }
     837          76 :     if (DEBUGLEVEL >= 6)
     838           0 :       timer_printf(&T, "[Fp_ellcard_Shanks] baby steps, s = %ld",s);
     839             : 
     840             :     /* giant steps: fg = s.F */
     841          76 :     fg = FpE_add(P,F,a4,p);
     842          76 :     if (ell_is_inf(fg)) { h = mului(s,B); goto FOUND; }
     843          76 :     pfinal = mod2BIL(p); av2 = avma;
     844             :     /* Goal of the following: sort points by increasing x-coordinate hash.
     845             :      * Done in a complicated way to avoid allocating a large temp vector */
     846          76 :     p1 = vecsmall_indexsort(tx); /* = permutation sorting tx */
     847          76 :     for (i=1; i<=s; i++) ti[i] = tx[p1[i]];
     848             :     /* ti = tx sorted */
     849          76 :     for (i=1; i<=s; i++) { tx[i] = ti[i]; ti[i] = ty[p1[i]]; }
     850             :     /* tx is sorted. ti = ty sorted */
     851          76 :     for (i=1; i<=s; i++) { ty[i] = ti[i]; ti[i] = p1[i]; }
     852             :     /* ty is sorted. ti = permutation sorting tx */
     853          76 :     if (DEBUGLEVEL >= 6) timer_printf(&T, "[Fp_ellcard_Shanks] sorting");
     854          76 :     avma = av2;
     855             : 
     856          76 :     gaffect(fg, gel(pts,1));
     857        9357 :     for (j=2; j<=nb; j++) /* pts[j] = j.fg = (s*j).F */
     858             :     {
     859        9281 :       P = FpE_add(gel(pts,j-1),fg,a4,p);
     860        9281 :       if (ell_is_inf(P)) { h = mulii(mulss(s,j), B); goto FOUND; }
     861        9281 :       gaffect(P, gel(pts,j));
     862             :     }
     863             :     /* replace fg by nb.fg since we do nb points at a time */
     864          76 :     avma = av2;
     865          76 :     fg = gcopy(gel(pts,nb)); /* copy: we modify (temporarily) pts[nb] below */
     866          76 :     av2 = avma;
     867             : 
     868      151888 :     for (i=1,j=1; ; i++)
     869             :     {
     870      151888 :       GEN ftest = gel(pts,j);
     871      151888 :       long m, l = 1, r = s+1;
     872             :       long k, k2, j2;
     873             : 
     874      151888 :       avma = av2;
     875      151888 :       k = mod2BIL(gel(ftest,1));
     876     2080615 :       while (l < r)
     877             :       {
     878     1776839 :         m = (l+r) >> 1;
     879     1776839 :         if (tx[m] < k) l = m+1; else r = m;
     880             :       }
     881      151888 :       if (r <= s && tx[r] == k)
     882             :       {
     883          76 :         while (r && tx[r] == k) r--;
     884          76 :         k2 = mod2BIL(gel(ftest,2));
     885          76 :         for (r++; r <= s && tx[r] == k; r++)
     886          76 :           if (ty[r] == k2 || ty[r] == pfinal - k2)
     887             :           { /* [h+j2] f == +/- ftest (= [i.s] f)? */
     888          76 :             j2 = ti[r] - 1;
     889          76 :             if (DEBUGLEVEL >=6)
     890           0 :               timer_printf(&T, "[Fp_ellcard_Shanks] giant steps, i = %ld",i);
     891          76 :             P = FpE_add(FpE_mul(F,stoi(j2),a4,p),fh,a4,p);
     892          76 :             if (equalii(gel(P,1), gel(ftest,1)))
     893             :             {
     894          76 :               if (equalii(gel(P,2), gel(ftest,2))) i = -i;
     895          76 :               h = addii(h, mulii(addis(mulss(s,i), j2), B));
     896          76 :               goto FOUND;
     897             :             }
     898             :           }
     899             :       }
     900      151812 :       if (++j > nb)
     901             :       { /* compute next nb points */
     902        1146 :         long save = 0; /* gcc -Wall */;
     903      147324 :         for (j=1; j<=nb; j++)
     904             :         {
     905      146178 :           P = gel(pts,j);
     906      146178 :           gel(u,j) = subii(gel(fg,1), gel(P,1));
     907      146178 :           if (gel(u,j) == gen_0) /* occurs once: i = j = nb, P == fg */
     908             :           {
     909          66 :             gel(u,j) = shifti(gel(P,2),1);
     910          66 :             save = fg[1]; fg[1] = P[1];
     911             :           }
     912             :         }
     913        1146 :         v = FpV_inv(u, p);
     914      147324 :         for (j=1; j<=nb; j++)
     915      146178 :           FpE_add_ip(gel(pts,j),fg,a4,p, gel(v,j));
     916        1146 :         if (i == nb) { fg[1] = save; }
     917        1146 :         j = 1;
     918             :       }
     919      151812 :     }
     920             : FOUND: /* found a point of exponent h on E_u */
     921          77 :     h = FpE_order(f, h, a4, p);
     922             :     /* h | #E_u(Fp) = A (mod B) */
     923          77 :     A = Z_chinese_all(A, gen_0, B, h, &B);
     924          77 :     if (cmpii(B, pordmin) >= 0) break;
     925             :     /* not done: update A mod B for the _next_ curve, isomorphic to
     926             :      * the quadratic twist of this one */
     927           0 :     A = remii(subii(p2p,A), B); /* #E(Fp)+#E'(Fp) = 2p+2 */
     928           0 :   }
     929          77 :   if (tx) killblock(tx);
     930          77 :   h = closest_lift(A, B, p1p);
     931          77 :   return gerepileuptoint(av, KRO==1? h: subii(p2p,h));
     932             : }
     933             : 
     934             : typedef struct
     935             : {
     936             :   ulong x,y,i;
     937             : } multiple;
     938             : 
     939             : static int
     940    21609093 : compare_multiples(multiple *a, multiple *b) { return a->x > b->x? 1:a->x<b->x?-1:0; }
     941             : 
     942             : /* find x such that h := a + b x is closest to c and return h:
     943             :  * x = round((c-a) / b) = floor( (2(c-a) + b) / 2b )
     944             :  * Assume 0 <= a < b < c  and b + 2c < 2^BIL */
     945             : static ulong
     946      417067 : uclosest_lift(ulong a, ulong b, ulong c)
     947             : {
     948      417067 :   ulong x = (b + ((c-a) << 1)) / (b << 1);
     949      417067 :   return a + b * x;
     950             : }
     951             : 
     952             : static long
     953      374942 : Fle_dbl_inplace(GEN P, ulong a4, ulong p)
     954             : {
     955             :   ulong x, y, slope;
     956      374942 :   if (!P[2]) return 1;
     957      374907 :   x = P[1]; y = P[2];
     958      374907 :   slope = Fl_div(Fl_add(Fl_triple(Fl_sqr(x,p), p), a4, p),
     959             :                  Fl_double(y, p), p);
     960      374907 :   P[1] = Fl_sub(Fl_sqr(slope, p), Fl_double(x, p), p);
     961      374907 :   P[2] = Fl_sub(Fl_mul(slope, Fl_sub(x, P[1], p), p), y, p);
     962      374907 :   return 0;
     963             : }
     964             : 
     965             : static long
     966    10097820 : Fle_add_inplace(GEN P, GEN Q, ulong a4, ulong p)
     967             : {
     968             :   ulong Px, Py, Qx, Qy, slope;
     969    10097820 :   if (ell_is_inf(Q)) return 0;
     970    10097820 :   Px = P[1]; Py = P[2];
     971    10097820 :   Qx = Q[1]; Qy = Q[2];
     972    10097820 :   if (Px==Qx)
     973      389068 :     return Py==Qy ? Fle_dbl_inplace(P, a4, p): 1;
     974     9708752 :   slope = Fl_div(Fl_sub(Py, Qy, p), Fl_sub(Px, Qx, p), p);
     975     9708752 :   P[1] = Fl_sub(Fl_sub(Fl_sqr(slope, p), Px, p), Qx, p);
     976     9708752 :   P[2] = Fl_sub(Fl_mul(slope, Fl_sub(Px, P[1], p), p), Py, p);
     977     9708752 :   return 0;
     978             : }
     979             : 
     980             : /* assume 99 < p < 2^(BIL-1) - 2^((BIL+1)/2) and e has good reduction at p.
     981             :  * Should use Barett reduction + multi-inverse. See Fp_ellcard_Shanks() */
     982             : static long
     983      410060 : Fl_ellcard_Shanks(ulong c4, ulong c6, ulong p)
     984             : {
     985             :   GEN f, fh, fg, ftest, F;
     986             :   ulong i, l, r, s, h, x, cp4, p1p, p2p, pordmin,A,B;
     987             :   long KRO;
     988      410060 :   pari_sp av = avma;
     989             :   multiple *table;
     990             : 
     991      410060 :   if (!c6) {
     992          14 :     GEN ap = ap_j1728(utoi(c4), utoipos(p));
     993          14 :     avma = av; return p+1 - itos(ap);
     994             :   }
     995             : 
     996      410046 :   pordmin = (ulong)(1 + 4*sqrt((double)p));
     997      410046 :   p1p = p+1;
     998      410046 :   p2p = p1p << 1;
     999      410046 :   x = 0; KRO = 0;
    1000      410046 :   switch(Flx_nbroots(mkvecsmall5(0L, c6,c4,0L,1L), p))
    1001             :   {
    1002       70481 :     case 3:  A = 0; B = 4; break;
    1003      199245 :     case 1:  A = 0; B = 2; break;
    1004      140320 :     default: A = 1; B = 2; break; /* 0 */
    1005             :   }
    1006             :   for(;;)
    1007             :   { /* see comments in Fp_ellcard_Shanks */
    1008      417067 :     h = uclosest_lift(A, B, p1p);
    1009      417067 :     if (!KRO) /* first time, initialize */
    1010             :     {
    1011      410046 :       KRO = krouu(c6,p); /* != 0 */
    1012      410046 :       f = mkvecsmall2(0, Fl_sqr(c6,p));
    1013             :     }
    1014             :     else
    1015             :     {
    1016        7021 :       KRO = -KRO;
    1017        7021 :       f = Fl_ellpoint(KRO, &x, c4,c6,p);
    1018             :     }
    1019      417067 :     cp4 = Fl_mul(c4, f[2], p);
    1020      417067 :     fh = Fle_mulu(f, h, cp4, p);
    1021      417067 :     if (ell_is_inf(fh)) goto FOUND;
    1022             : 
    1023      411166 :     s = (ulong) (sqrt(((double)pordmin)/B) / 2);
    1024      411166 :     if (!s) s = 1;
    1025      411166 :     table = (multiple *) stack_malloc((s+1) * sizeof(multiple));
    1026      411166 :     F = Fle_mulu(f, B, cp4, p);
    1027     5823606 :     for (i=0; i < s; i++)
    1028             :     {
    1029     5426601 :       table[i].x = fh[1];
    1030     5426601 :       table[i].y = fh[2];
    1031     5426601 :       table[i].i = i;
    1032     5426601 :       if (Fle_add_inplace(fh, F, cp4, p)) { h += B*(i+1); goto FOUND; }
    1033             :     }
    1034      397005 :     qsort(table,s,sizeof(multiple),(QSCOMP)compare_multiples);
    1035      397005 :     fg = Fle_mulu(F, s, cp4, p); ftest = zv_copy(fg);
    1036      397005 :     if (ell_is_inf(ftest)) {
    1037           0 :       if (!uisprime(p)) pari_err_PRIME("ellap",utoi(p));
    1038           0 :       pari_err_BUG("ellap (f^(i*s) = 1)");
    1039             :     }
    1040     5068224 :     for (i=1; ; i++)
    1041             :     {
    1042     5068224 :       l=0; r=s;
    1043    36294009 :       while (l<r)
    1044             :       {
    1045    26157561 :         ulong m = (l+r) >> 1;
    1046    26157561 :         if (table[m].x < uel(ftest,1)) l=m+1; else r=m;
    1047             :       }
    1048     5068224 :       if (r < s && table[r].x == uel(ftest,1)) break;
    1049     4671219 :       if (Fle_add_inplace(ftest, fg, cp4, p))
    1050           0 :         pari_err_PRIME("ellap",utoi(p));
    1051     4671219 :     }
    1052      397005 :     h += table[r].i * B;
    1053      397005 :     if (table[r].y == uel(ftest,2))
    1054      204489 :       h -= s * i * B;
    1055             :     else
    1056      192516 :       h += s * i * B;
    1057             : FOUND:
    1058      417067 :     h = itou(Fle_order(f, utoipos(h), cp4, p));
    1059             :     /* h | #E_u(Fp) = A (mod B) */
    1060             :     {
    1061             :       GEN C;
    1062      417067 :       A = itou( Z_chinese_all(gen_0, utoi(A), utoipos(h), utoipos(B), &C) );
    1063      417067 :       if (abscmpiu(C, pordmin) >= 0) { /* uclosest_lift could overflow */
    1064      410046 :         h = itou( closest_lift(utoi(A), C, utoipos(p1p)) );
    1065      410046 :         break;
    1066             :       }
    1067        7021 :       B = itou(C);
    1068             :     }
    1069        7021 :     A = (p2p - A) % B; avma = av;
    1070        7021 :   }
    1071      610682 :   avma = av; return KRO==1? h: p2p-h;
    1072             : }
    1073             : 
    1074             : /** ellap from CM (original code contributed by Mark Watkins) **/
    1075             : 
    1076             : static GEN
    1077       43386 : ap_j0(GEN a6,GEN p)
    1078             : {
    1079             :   GEN a, b, e, d;
    1080       43386 :   if (umodiu(p,3) != 1) return gen_0;
    1081       21609 :   (void)cornacchia2(utoipos(27),p, &a,&b);
    1082       21609 :   if (umodiu(a, 3) == 1) a = negi(a);
    1083       21609 :   d = mulis(a6,-108);
    1084       21609 :   e = diviuexact(shifti(p,-1), 3); /* (p-1) / 6 */
    1085       21609 :   return centermod(mulii(a, Fp_pow(d, e, p)), p);
    1086             : }
    1087             : static GEN
    1088     2617825 : ap_j1728(GEN a4,GEN p)
    1089             : {
    1090             :   GEN a, b, e;
    1091     2617825 :   if (mod4(p) != 1) return gen_0;
    1092     1307922 :   (void)cornacchia2(utoipos(4),p, &a,&b);
    1093     1307922 :   if (Mod4(a)==0) a = b;
    1094     1307922 :   if (Mod2(a)==1) a = shifti(a,1);
    1095     1307922 :   if (Mod8(a)==6) a = negi(a);
    1096     1307922 :   e = shifti(p,-2); /* (p-1) / 4 */
    1097     1307922 :   return centermod(mulii(a, Fp_pow(a4, e, p)), p);
    1098             : }
    1099             : static GEN
    1100         140 : ap_j8000(GEN a6, GEN p)
    1101             : {
    1102             :   GEN a, b;
    1103         140 :   long r = mod8(p), s = 1;
    1104         140 :   if (r != 1 && r != 3) return gen_0;
    1105          63 :   (void)cornacchia2(utoipos(8),p, &a,&b);
    1106          63 :   switch(Mod16(a)) {
    1107          28 :     case 2: case 6:   if (Mod4(b)) s = -s;
    1108          28 :       break;
    1109          35 :     case 10: case 14: if (!Mod4(b)) s = -s;
    1110          35 :       break;
    1111             :   }
    1112          63 :   if (kronecker(mulis(a6, 42), p) < 0) s = -s;
    1113          63 :   return s > 0? a: negi(a);
    1114             : }
    1115             : static GEN
    1116         161 : ap_j287496(GEN a6, GEN p)
    1117             : {
    1118             :   GEN a, b;
    1119         161 :   long s = 1;
    1120         161 :   if (mod4(p) != 1) return gen_0;
    1121          70 :   (void)cornacchia2(utoipos(4),p, &a,&b);
    1122          70 :   if (Mod4(a)==0) a = b;
    1123          70 :   if (Mod2(a)==1) a = shifti(a,1);
    1124          70 :   if (Mod8(a)==6) s = -s;
    1125          70 :   if (krosi(2,p) < 0) s = -s;
    1126          70 :   if (kronecker(mulis(a6, -14), p) < 0) s = -s;
    1127          70 :   return s > 0? a: negi(a);
    1128             : }
    1129             : static GEN
    1130        1547 : ap_cm(int CM, long A6B, GEN a6, GEN p)
    1131             : {
    1132             :   GEN a, b;
    1133        1547 :   long s = 1;
    1134        1547 :   if (krosi(CM,p) < 0) return gen_0;
    1135         777 :   (void)cornacchia2(utoipos(-CM),p, &a, &b);
    1136         777 :   if ((CM&3) == 0) CM >>= 2;
    1137         777 :   if ((krois(a, -CM) > 0) ^ (CM == -7)) s = -s;
    1138         777 :   if (kronecker(mulis(a6,A6B), p) < 0) s = -s;
    1139         777 :   return s > 0? a: negi(a);
    1140             : }
    1141             : static GEN
    1142       11522 : ec_ap_cm(int CM, GEN a4, GEN a6, GEN p)
    1143             : {
    1144       11522 :   switch(CM)
    1145             :   {
    1146           0 :     case  -3: return ap_j0(a6, p);
    1147        9674 :     case  -4: return ap_j1728(a4, p);
    1148         140 :     case  -8: return ap_j8000(a6, p);
    1149         161 :     case -16: return ap_j287496(a6, p);
    1150         140 :     case  -7: return ap_cm(CM, -2, a6, p);
    1151         147 :     case -11: return ap_cm(CM, 21, a6, p);
    1152         252 :     case -12: return ap_cm(CM, 22, a6, p);
    1153         126 :     case -19: return ap_cm(CM, 1, a6, p);
    1154         266 :     case -27: return ap_cm(CM, 253, a6, p);
    1155         182 :     case -28: return ap_cm(-7, -114, a6, p); /* yes, -7 ! */
    1156         161 :     case -43: return ap_cm(CM, 21, a6, p);
    1157         126 :     case -67: return ap_cm(CM, 217, a6, p);
    1158         147 :     case -163:return ap_cm(CM, 185801, a6, p);
    1159           0 :     default: return NULL;
    1160             :   }
    1161             : }
    1162             : 
    1163             : static GEN
    1164       72898 : Fp_ellj_nodiv(GEN a4, GEN a6, GEN p)
    1165             : {
    1166       72898 :   GEN a43 = Fp_mulu(Fp_powu(a4, 3, p), 4, p);
    1167       72898 :   GEN a62 = Fp_mulu(Fp_sqr(a6, p), 27, p);
    1168       72898 :   return mkvec2(Fp_mulu(a43, 1728, p), Fp_add(a43, a62, p));
    1169             : }
    1170             : 
    1171             : GEN
    1172           0 : Fp_ellj(GEN a4, GEN a6, GEN p)
    1173             : {
    1174           0 :   pari_sp av=avma;
    1175           0 :   GEN z = Fp_ellj_nodiv(a4, a6, p);
    1176           0 :   return gerepileuptoint(av,Fp_div(gel(z,1),gel(z,2),p));
    1177             : }
    1178             : 
    1179             : static GEN /* Only compute a mod p, so assume p>=17 */
    1180     2724421 : Fp_ellcard_CM(GEN a4, GEN a6, GEN p)
    1181             : {
    1182     2724421 :   pari_sp av = avma;
    1183             :   GEN a;
    1184     2724421 :   if (!signe(a4)) a = ap_j0(a6,p);
    1185     2681035 :   else if (!signe(a6)) a = ap_j1728(a4,p);
    1186             :   else
    1187             :   {
    1188       72898 :     GEN j = Fp_ellj_nodiv(a4, a6, p);
    1189       72898 :     long CM = Fp_ellj_get_CM(gel(j,1), gel(j,2), p);
    1190       72898 :     if (!CM) { avma = av; return NULL; }
    1191        1848 :     a = ec_ap_cm(CM,a4,a6,p);
    1192             :   }
    1193     2653371 :   return gerepileuptoint(av, subii(addis(p,1),a));
    1194             : }
    1195             : 
    1196             : GEN
    1197     2840299 : Fp_ellcard(GEN a4, GEN a6, GEN p)
    1198             : {
    1199     2840299 :   long lp = expi(p);
    1200     2840299 :   ulong pp = p[2];
    1201     2840299 :   if (lp < 11)
    1202      115878 :     return utoi(pp+1 - Fl_elltrace_naive(umodiu(a4,pp), umodiu(a6,pp), pp));
    1203     2724421 :   { GEN a = Fp_ellcard_CM(a4,a6,p); if (a) return a; }
    1204       71050 :   if (lp >= 56)
    1205         854 :     return Fp_ellcard_SEA(a4, a6, p, 0);
    1206       70196 :   if (lp <= BITS_IN_LONG-2)
    1207       70119 :     return utoi(Fl_ellcard_Shanks(umodiu(a4,pp), umodiu(a6,pp), pp));
    1208          77 :   return Fp_ellcard_Shanks(a4, a6, p);
    1209             : }
    1210             : 
    1211             : long
    1212      410775 : Fl_elltrace(ulong a4, ulong a6, ulong p)
    1213             : {
    1214             :   pari_sp av;
    1215             :   long lp;
    1216             :   GEN a;
    1217      410775 :   if (p < (1<<11)) return Fl_elltrace_naive(a4, a6, p);
    1218      339941 :   lp = expu(p);
    1219      339941 :   if (lp <= minss(56, BITS_IN_LONG-2)) return p+1-Fl_ellcard_Shanks(a4, a6, p);
    1220           0 :   av = avma; a = subui(p+1, Fp_ellcard(utoi(a4), utoi(a6), utoipos(p)));
    1221           0 :   avma = av; return itos(a);
    1222             : }
    1223             : long
    1224      442361 : Fl_elltrace_CM(long CM, ulong a4, ulong a6, ulong p)
    1225             : {
    1226             :   pari_sp av;
    1227             :   GEN a;
    1228      442361 :   if (!CM) return Fl_elltrace(a4,a6,p);
    1229       31628 :   if (p < (1<<11)) return Fl_elltrace_naive(a4, a6, p);
    1230        9674 :   av = avma; a = ec_ap_cm(CM, utoi(a4), utoi(a6), utoipos(p));
    1231        9674 :   avma = av; return itos(a);
    1232             : }
    1233             : 
    1234             : static GEN
    1235       10415 : _FpE_pairorder(void *E, GEN P, GEN Q, GEN m, GEN F)
    1236             : {
    1237       10415 :   struct _FpE *e = (struct _FpE *) E;
    1238       10415 :   return  Fp_order(FpE_weilpairing(P,Q,m,e->a4,e->p), F, e->p);
    1239             : }
    1240             : 
    1241             : GEN
    1242       21525 : Fp_ellgroup(GEN a4, GEN a6, GEN N, GEN p, GEN *pt_m)
    1243             : {
    1244             :   struct _FpE e;
    1245       21525 :   e.a4=a4; e.a6=a6; e.p=p;
    1246       21525 :   return gen_ellgroup(N, subis(p, 1), pt_m, (void*)&e, &FpE_group, _FpE_pairorder);
    1247             : }
    1248             : 
    1249             : GEN
    1250         574 : Fp_ellgens(GEN a4, GEN a6, GEN ch, GEN D, GEN m, GEN p)
    1251             : {
    1252             :   GEN P;
    1253         574 :   pari_sp av = avma;
    1254             :   struct _FpE e;
    1255         574 :   e.a4=a4; e.a6=a6; e.p=p;
    1256         574 :   switch(lg(D)-1)
    1257             :   {
    1258             :   case 1:
    1259         476 :     P = gen_gener(gel(D,1), (void*)&e, &FpE_group);
    1260         476 :     P = mkvec(FpE_changepoint(P, ch, p));
    1261         476 :     break;
    1262             :   default:
    1263          98 :     P = gen_ellgens(gel(D,1), gel(D,2), m, (void*)&e, &FpE_group, _FpE_pairorder);
    1264          98 :     gel(P,1) = FpE_changepoint(gel(P,1), ch, p);
    1265          98 :     gel(P,2) = FpE_changepoint(gel(P,2), ch, p);
    1266          98 :     break;
    1267             :   }
    1268         574 :   return gerepilecopy(av, P);
    1269             : }
    1270             : 
    1271             : /* Not so fast arithmetic with points over elliptic curves over FpXQ */
    1272             : 
    1273             : /***********************************************************************/
    1274             : /**                                                                   **/
    1275             : /**                              FpXQE                                  **/
    1276             : /**                                                                   **/
    1277             : /***********************************************************************/
    1278             : 
    1279             : /* Theses functions deal with point over elliptic curves over FpXQ defined
    1280             :  * by an equation of the form y^2=x^3+a4*x+a6.
    1281             :  * Most of the time a6 is omitted since it can be recovered from any point
    1282             :  * on the curve.
    1283             :  */
    1284             : 
    1285             : GEN
    1286         896 : RgE_to_FpXQE(GEN x, GEN T, GEN p)
    1287             : {
    1288         896 :   if (ell_is_inf(x)) return x;
    1289         896 :   retmkvec2(Rg_to_FpXQ(gel(x,1),T,p),Rg_to_FpXQ(gel(x,2),T,p));
    1290             : }
    1291             : 
    1292             : GEN
    1293        1715 : FpXQE_changepoint(GEN x, GEN ch, GEN T, GEN p)
    1294             : {
    1295        1715 :   pari_sp av = avma;
    1296             :   GEN p1,z,u,r,s,t,v,v2,v3;
    1297        1715 :   if (ell_is_inf(x)) return x;
    1298         861 :   u = gel(ch,1); r = gel(ch,2);
    1299         861 :   s = gel(ch,3); t = gel(ch,4);
    1300         861 :   v = FpXQ_inv(u, T, p); v2 = FpXQ_sqr(v, T, p); v3 = FpXQ_mul(v,v2, T, p);
    1301         861 :   p1 = FpX_sub(gel(x,1),r, p);
    1302         861 :   z = cgetg(3,t_VEC);
    1303         861 :   gel(z,1) = FpXQ_mul(v2, p1, T, p);
    1304         861 :   gel(z,2) = FpXQ_mul(v3, FpX_sub(gel(x,2), FpX_add(FpXQ_mul(s,p1, T, p),t, p), p), T, p);
    1305         861 :   return gerepileupto(av, z);
    1306             : }
    1307             : 
    1308             : GEN
    1309         896 : FpXQE_changepointinv(GEN x, GEN ch, GEN T, GEN p)
    1310             : {
    1311             :   GEN u, r, s, t, X, Y, u2, u3, u2X, z;
    1312         896 :   if (ell_is_inf(x)) return x;
    1313         896 :   X = gel(x,1); Y = gel(x,2);
    1314         896 :   u = gel(ch,1); r = gel(ch,2);
    1315         896 :   s = gel(ch,3); t = gel(ch,4);
    1316         896 :   u2 = FpXQ_sqr(u, T, p); u3 = FpXQ_mul(u,u2, T, p);
    1317         896 :   u2X = FpXQ_mul(u2,X, T, p);
    1318         896 :   z = cgetg(3, t_VEC);
    1319         896 :   gel(z,1) = FpX_add(u2X,r, p);
    1320         896 :   gel(z,2) = FpX_add(FpXQ_mul(u3,Y, T, p), FpX_add(FpXQ_mul(s,u2X, T, p), t, p), p);
    1321         896 :   return z;
    1322             : }
    1323             : 
    1324             : static GEN
    1325         840 : nonsquare_FpXQ(GEN T, GEN p)
    1326             : {
    1327         840 :   pari_sp av = avma;
    1328         840 :   long n = degpol(T), v = varn(T);
    1329             :   GEN a;
    1330         840 :   if (odd(n))
    1331             :   {
    1332         420 :     GEN z = cgetg(3, t_POL);
    1333         420 :     z[1] = evalsigne(1) | evalvarn(v);
    1334         420 :     gel(z,2) = nonsquare_Fp(p); return z;
    1335             :   }
    1336             :   do
    1337             :   {
    1338         784 :     avma = av;
    1339         784 :     a = random_FpX(n, v, p);
    1340         784 :   } while (FpXQ_issquare(a, T, p));
    1341         420 :   return a;
    1342             : }
    1343             : 
    1344             : void
    1345         840 : FpXQ_elltwist(GEN a4, GEN a6, GEN T, GEN p, GEN *pt_a4, GEN *pt_a6)
    1346             : {
    1347         840 :   GEN d = nonsquare_FpXQ(T, p);
    1348         840 :   GEN d2 = FpXQ_sqr(d, T, p), d3 = FpXQ_mul(d2, d, T, p);
    1349         840 :   *pt_a4 = FpXQ_mul(a4, d2, T, p);
    1350         840 :   *pt_a6 = FpXQ_mul(a6, d3, T, p);
    1351         840 : }
    1352             : 
    1353             : static GEN
    1354      186144 : FpXQE_dbl_slope(GEN P, GEN a4, GEN T, GEN p, GEN *slope)
    1355             : {
    1356             :   GEN x, y, Q;
    1357      186144 :   if (ell_is_inf(P) || !signe(gel(P,2))) return ellinf();
    1358      184883 :   x = gel(P,1); y = gel(P,2);
    1359      184883 :   *slope = FpXQ_div(FpX_add(FpX_mulu(FpXQ_sqr(x, T, p), 3, p), a4, p),
    1360             :                             FpX_mulu(y, 2, p), T, p);
    1361      184883 :   Q = cgetg(3,t_VEC);
    1362      184883 :   gel(Q, 1) = FpX_sub(FpXQ_sqr(*slope, T, p), FpX_mulu(x, 2, p), p);
    1363      184883 :   gel(Q, 2) = FpX_sub(FpXQ_mul(*slope, FpX_sub(x, gel(Q, 1), p), T, p), y, p);
    1364      184883 :   return Q;
    1365             : }
    1366             : 
    1367             : GEN
    1368      181188 : FpXQE_dbl(GEN P, GEN a4, GEN T, GEN p)
    1369             : {
    1370      181188 :   pari_sp av = avma;
    1371             :   GEN slope;
    1372      181188 :   return gerepileupto(av, FpXQE_dbl_slope(P,a4,T,p,&slope));
    1373             : }
    1374             : 
    1375             : static GEN
    1376       35527 : FpXQE_add_slope(GEN P, GEN Q, GEN a4, GEN T, GEN p, GEN *slope)
    1377             : {
    1378             :   GEN Px, Py, Qx, Qy, R;
    1379       35527 :   if (ell_is_inf(P)) return Q;
    1380       35527 :   if (ell_is_inf(Q)) return P;
    1381       35527 :   Px = gel(P,1); Py = gel(P,2);
    1382       35527 :   Qx = gel(Q,1); Qy = gel(Q,2);
    1383       35527 :   if (ZX_equal(Px, Qx))
    1384             :   {
    1385         636 :     if (ZX_equal(Py, Qy))
    1386           7 :       return FpXQE_dbl_slope(P, a4, T, p, slope);
    1387             :     else
    1388         629 :       return ellinf();
    1389             :   }
    1390       34891 :   *slope = FpXQ_div(FpX_sub(Py, Qy, p), FpX_sub(Px, Qx, p), T, p);
    1391       34891 :   R = cgetg(3,t_VEC);
    1392       34891 :   gel(R, 1) = FpX_sub(FpX_sub(FpXQ_sqr(*slope, T, p), Px, p), Qx, p);
    1393       34891 :   gel(R, 2) = FpX_sub(FpXQ_mul(*slope, FpX_sub(Px, gel(R, 1), p), T, p), Py, p);
    1394       34891 :   return R;
    1395             : }
    1396             : 
    1397             : GEN
    1398       34715 : FpXQE_add(GEN P, GEN Q, GEN a4, GEN T, GEN p)
    1399             : {
    1400       34715 :   pari_sp av = avma;
    1401             :   GEN slope;
    1402       34715 :   return gerepileupto(av, FpXQE_add_slope(P,Q,a4,T,p,&slope));
    1403             : }
    1404             : 
    1405             : static GEN
    1406           0 : FpXQE_neg_i(GEN P, GEN p)
    1407             : {
    1408           0 :   if (ell_is_inf(P)) return P;
    1409           0 :   return mkvec2(gel(P,1), FpX_neg(gel(P,2), p));
    1410             : }
    1411             : 
    1412             : GEN
    1413         749 : FpXQE_neg(GEN P, GEN T, GEN p)
    1414             : {
    1415             :   (void) T;
    1416         749 :   if (ell_is_inf(P)) return ellinf();
    1417         749 :   return mkvec2(gcopy(gel(P,1)), FpX_neg(gel(P,2), p));
    1418             : }
    1419             : 
    1420             : GEN
    1421           0 : FpXQE_sub(GEN P, GEN Q, GEN a4, GEN T, GEN p)
    1422             : {
    1423           0 :   pari_sp av = avma;
    1424             :   GEN slope;
    1425           0 :   return gerepileupto(av, FpXQE_add_slope(P, FpXQE_neg_i(Q, p), a4, T, p, &slope));
    1426             : }
    1427             : 
    1428             : struct _FpXQE
    1429             : {
    1430             :   GEN a4,a6;
    1431             :   GEN T,p;
    1432             : };
    1433             : 
    1434             : static GEN
    1435      181188 : _FpXQE_dbl(void *E, GEN P)
    1436             : {
    1437      181188 :   struct _FpXQE *ell = (struct _FpXQE *) E;
    1438      181188 :   return FpXQE_dbl(P, ell->a4, ell->T, ell->p);
    1439             : }
    1440             : 
    1441             : static GEN
    1442       34715 : _FpXQE_add(void *E, GEN P, GEN Q)
    1443             : {
    1444       34715 :   struct _FpXQE *ell=(struct _FpXQE *) E;
    1445       34715 :   return FpXQE_add(P, Q, ell->a4, ell->T, ell->p);
    1446             : }
    1447             : 
    1448             : static GEN
    1449        2825 : _FpXQE_mul(void *E, GEN P, GEN n)
    1450             : {
    1451        2825 :   pari_sp av = avma;
    1452        2825 :   struct _FpXQE *e=(struct _FpXQE *) E;
    1453        2825 :   long s = signe(n);
    1454        2825 :   if (!s || ell_is_inf(P)) return ellinf();
    1455        2825 :   if (s<0) P = FpXQE_neg(P, e->T, e->p);
    1456        2825 :   if (is_pm1(n)) return s>0? gcopy(P): P;
    1457        1971 :   return gerepileupto(av, gen_pow(P, n, e, &_FpXQE_dbl, &_FpXQE_add));
    1458             : }
    1459             : 
    1460             : GEN
    1461         854 : FpXQE_mul(GEN P, GEN n, GEN a4, GEN T, GEN p)
    1462             : {
    1463             :   struct _FpXQE E;
    1464         854 :   E.a4= a4; E.T = T; E.p = p;
    1465         854 :   return _FpXQE_mul(&E, P, n);
    1466             : }
    1467             : 
    1468             : /* Finds a random non-singular point on E */
    1469             : 
    1470             : GEN
    1471         985 : random_FpXQE(GEN a4, GEN a6, GEN T, GEN p)
    1472             : {
    1473         985 :   pari_sp ltop = avma;
    1474             :   GEN x, x2, y, rhs;
    1475         985 :   long v = get_FpX_var(T), d = get_FpX_degree(T);
    1476             :   do
    1477             :   {
    1478        1858 :     avma= ltop;
    1479        1858 :     x   = random_FpX(d,v,p); /*  x^3+a4*x+a6 = x*(x^2+a4)+a6  */
    1480        1858 :     x2  = FpXQ_sqr(x, T, p);
    1481        1858 :     rhs = FpX_add(FpXQ_mul(x, FpX_add(x2, a4, p), T, p), a6, p);
    1482        1858 :   } while ((!signe(rhs) && !signe(FpX_add(FpX_mulu(x2,3,p), a4, p)))
    1483        3716 :           || !FpXQ_issquare(rhs, T, p));
    1484         985 :   y = FpXQ_sqrt(rhs, T, p);
    1485         985 :   if (!y) pari_err_PRIME("random_FpE", p);
    1486         985 :   return gerepilecopy(ltop, mkvec2(x, y));
    1487             : }
    1488             : 
    1489             : static GEN
    1490         131 : _FpXQE_rand(void *E)
    1491             : {
    1492         131 :   struct _FpXQE *e=(struct _FpXQE *) E;
    1493         131 :   return random_FpXQE(e->a4, e->a6, e->T, e->p);
    1494             : }
    1495             : 
    1496             : static const struct bb_group FpXQE_group={_FpXQE_add,_FpXQE_mul,_FpXQE_rand,hash_GEN,ZXV_equal,ell_is_inf};
    1497             : 
    1498             : const struct bb_group *
    1499           8 : get_FpXQE_group(void ** pt_E, GEN a4, GEN a6, GEN T, GEN p)
    1500             : {
    1501           8 :   struct _FpXQE *e = (struct _FpXQE *) stack_malloc(sizeof(struct _FpXQE));
    1502           8 :   e->a4 = a4; e->a6 = a6; e->T = T; e->p = p;
    1503           8 :   *pt_E = (void *) e;
    1504           8 :   return &FpXQE_group;
    1505             : }
    1506             : 
    1507             : GEN
    1508          14 : FpXQE_order(GEN z, GEN o, GEN a4, GEN T, GEN p)
    1509             : {
    1510          14 :   pari_sp av = avma;
    1511             :   struct _FpXQE e;
    1512          14 :   e.a4=a4; e.T=T; e.p=p;
    1513          14 :   return gerepileuptoint(av, gen_order(z, o, (void*)&e, &FpXQE_group));
    1514             : }
    1515             : 
    1516             : GEN
    1517           0 : FpXQE_log(GEN a, GEN b, GEN o, GEN a4, GEN T, GEN p)
    1518             : {
    1519           0 :   pari_sp av = avma;
    1520             :   struct _FpXQE e;
    1521           0 :   e.a4=a4; e.T=T; e.p=p;
    1522           0 :   return gerepileuptoint(av, gen_PH_log(a, b, o, (void*)&e, &FpXQE_group));
    1523             : }
    1524             : 
    1525             : 
    1526             : /***********************************************************************/
    1527             : /**                                                                   **/
    1528             : /**                            Pairings                               **/
    1529             : /**                                                                   **/
    1530             : /***********************************************************************/
    1531             : 
    1532             : /* Derived from APIP from and by Jerome Milan, 2012 */
    1533             : 
    1534             : static GEN
    1535        5936 : FpXQE_vert(GEN P, GEN Q, GEN a4, GEN T, GEN p)
    1536             : {
    1537        5936 :   long vT = get_FpX_var(T);
    1538        5936 :   if (ell_is_inf(P))
    1539          98 :     return pol_1(get_FpX_var(T));
    1540        5838 :   if (!ZX_equal(gel(Q, 1), gel(P, 1)))
    1541        5838 :     return FpX_sub(gel(Q, 1), gel(P, 1), p);
    1542           0 :   if (signe(gel(P,2))!=0) return pol_1(vT);
    1543           0 :   return FpXQ_inv(FpX_add(FpX_mulu(FpXQ_sqr(gel(P,1), T, p), 3, p),
    1544             :                   a4, p), T, p);
    1545             : }
    1546             : 
    1547             : static GEN
    1548        5761 : FpXQE_Miller_line(GEN R, GEN Q, GEN slope, GEN a4, GEN T, GEN p)
    1549             : {
    1550        5761 :   long vT = get_FpX_var(T);
    1551        5761 :   GEN x = gel(Q, 1), y = gel(Q, 2);
    1552        5761 :   GEN tmp1  = FpX_sub(x, gel(R, 1), p);
    1553        5761 :   GEN tmp2  = FpX_add(FpXQ_mul(tmp1, slope, T, p), gel(R, 2), p);
    1554        5761 :   if (!ZX_equal(y, tmp2))
    1555        5761 :     return FpX_sub(y, tmp2, p);
    1556           0 :   if (signe(y) == 0)
    1557           0 :     return pol_1(vT);
    1558             :   else
    1559             :   {
    1560             :     GEN s1, s2;
    1561           0 :     GEN y2i = FpXQ_inv(FpX_mulu(y, 2, p), T, p);
    1562           0 :     s1 = FpXQ_mul(FpX_add(FpX_mulu(FpXQ_sqr(x, T, p), 3, p), a4, p), y2i, T, p);
    1563           0 :     if (!ZX_equal(s1, slope))
    1564           0 :       return FpX_sub(s1, slope, p);
    1565           0 :     s2 = FpXQ_mul(FpX_sub(FpX_mulu(x, 3, p), FpXQ_sqr(s1, T, p), p), y2i, T, p);
    1566           0 :     return signe(s2)!=0 ? s2: y2i;
    1567             :   }
    1568             : }
    1569             : 
    1570             : /* Computes the equation of the line tangent to R and returns its
    1571             :    evaluation at the point Q. Also doubles the point R.
    1572             :  */
    1573             : 
    1574             : static GEN
    1575        5026 : FpXQE_tangent_update(GEN R, GEN Q, GEN a4, GEN T, GEN p, GEN *pt_R)
    1576             : {
    1577        5026 :   if (ell_is_inf(R))
    1578             :   {
    1579          21 :     *pt_R = ellinf();
    1580          21 :     return pol_1(get_FpX_var(T));
    1581             :   }
    1582        5005 :   else if (!signe(gel(R,2)))
    1583             :   {
    1584          56 :     *pt_R = ellinf();
    1585          56 :     return FpXQE_vert(R, Q, a4, T, p);
    1586             :   } else {
    1587             :     GEN slope;
    1588        4949 :     *pt_R = FpXQE_dbl_slope(R, a4, T, p, &slope);
    1589        4949 :     return FpXQE_Miller_line(R, Q, slope, a4, T, p);
    1590             :   }
    1591             : }
    1592             : 
    1593             : /* Computes the equation of the line through R and P, and returns its
    1594             :    evaluation at the point Q. Also adds P to the point R.
    1595             :  */
    1596             : 
    1597             : static GEN
    1598         833 : FpXQE_chord_update(GEN R, GEN P, GEN Q, GEN a4, GEN T, GEN p, GEN *pt_R)
    1599             : {
    1600         833 :   if (ell_is_inf(R))
    1601             :   {
    1602           0 :     *pt_R = gcopy(P);
    1603           0 :     return FpXQE_vert(P, Q, a4, T, p);
    1604             :   }
    1605         833 :   else if (ell_is_inf(P))
    1606             :   {
    1607           0 :     *pt_R = gcopy(R);
    1608           0 :     return FpXQE_vert(R, Q, a4, T, p);
    1609             :   }
    1610         833 :   else if (ZX_equal(gel(P, 1), gel(R, 1)))
    1611             :   {
    1612          21 :     if (ZX_equal(gel(P, 2), gel(R, 2)))
    1613           0 :       return FpXQE_tangent_update(R, Q, a4, T, p, pt_R);
    1614             :     else
    1615             :     {
    1616          21 :       *pt_R = ellinf();
    1617          21 :       return FpXQE_vert(R, Q, a4, T, p);
    1618             :     }
    1619             :   } else {
    1620             :     GEN slope;
    1621         812 :     *pt_R = FpXQE_add_slope(P, R, a4, T, p, &slope);
    1622         812 :     return FpXQE_Miller_line(R, Q, slope, a4, T, p);
    1623             :   }
    1624             : }
    1625             : 
    1626             : /* Returns the Miller function f_{m, Q} evaluated at the point P using
    1627             :    the standard Miller algorithm.
    1628             :  */
    1629             : 
    1630             : struct _FpXQE_miller
    1631             : {
    1632             :   GEN p;
    1633             :   GEN T, a4, P;
    1634             : };
    1635             : 
    1636             : static GEN
    1637        5026 : FpXQE_Miller_dbl(void* E, GEN d)
    1638             : {
    1639        5026 :   struct _FpXQE_miller *m = (struct _FpXQE_miller *)E;
    1640        5026 :   GEN p  = m->p;
    1641        5026 :   GEN T = m->T, a4 = m->a4, P = m->P;
    1642             :   GEN v, line;
    1643        5026 :   GEN num = FpXQ_sqr(gel(d,1), T, p);
    1644        5026 :   GEN denom = FpXQ_sqr(gel(d,2), T, p);
    1645        5026 :   GEN point = gel(d,3);
    1646        5026 :   line = FpXQE_tangent_update(point, P, a4, T, p, &point);
    1647        5026 :   num  = FpXQ_mul(num, line, T, p);
    1648        5026 :   v = FpXQE_vert(point, P, a4, T, p);
    1649        5026 :   denom = FpXQ_mul(denom, v, T, p);
    1650        5026 :   return mkvec3(num, denom, point);
    1651             : }
    1652             : 
    1653             : static GEN
    1654         833 : FpXQE_Miller_add(void* E, GEN va, GEN vb)
    1655             : {
    1656         833 :   struct _FpXQE_miller *m = (struct _FpXQE_miller *)E;
    1657         833 :   GEN p = m->p;
    1658         833 :   GEN T = m->T, a4 = m->a4, P = m->P;
    1659             :   GEN v, line, point;
    1660         833 :   GEN na = gel(va,1), da = gel(va,2), pa = gel(va,3);
    1661         833 :   GEN nb = gel(vb,1), db = gel(vb,2), pb = gel(vb,3);
    1662         833 :   GEN num   = FpXQ_mul(na, nb, T, p);
    1663         833 :   GEN denom = FpXQ_mul(da, db, T, p);
    1664         833 :   line = FpXQE_chord_update(pa, pb, P, a4, T, p, &point);
    1665         833 :   num  = FpXQ_mul(num, line, T, p);
    1666         833 :   v = FpXQE_vert(point, P, a4, T, p);
    1667         833 :   denom = FpXQ_mul(denom, v, T, p);
    1668         833 :   return mkvec3(num, denom, point);
    1669             : }
    1670             : 
    1671             : static GEN
    1672          77 : FpXQE_Miller(GEN Q, GEN P, GEN m, GEN a4, GEN T, GEN p)
    1673             : {
    1674          77 :   pari_sp ltop = avma;
    1675             :   struct _FpXQE_miller d;
    1676             :   GEN v, num, denom, g1;
    1677             : 
    1678          77 :   d.a4 = a4; d.T = T; d.p = p; d.P = P;
    1679          77 :   g1 = pol_1(get_FpX_var(T));
    1680          77 :   v = gen_pow(mkvec3(g1,g1,Q), m, (void*)&d, FpXQE_Miller_dbl, FpXQE_Miller_add);
    1681          77 :   num = gel(v,1); denom = gel(v,2);
    1682          77 :   return gerepileupto(ltop, FpXQ_div(num, denom, T, p));
    1683             : }
    1684             : 
    1685             : GEN
    1686          35 : FpXQE_weilpairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, GEN p)
    1687             : {
    1688          35 :   pari_sp ltop = avma;
    1689             :   GEN num, denom, result;
    1690          35 :   if (ell_is_inf(P) || ell_is_inf(Q) || ZXV_equal(P,Q))
    1691           0 :     return pol_1(get_FpX_var(T));
    1692          35 :   num    = FpXQE_Miller(P, Q, m, a4, T, p);
    1693          35 :   denom  = FpXQE_Miller(Q, P, m, a4, T, p);
    1694          35 :   result = FpXQ_div(num, denom, T, p);
    1695          35 :   if (mpodd(m))
    1696           0 :     result  = FpX_neg(result, p);
    1697          35 :   return gerepileupto(ltop, result);
    1698             : }
    1699             : 
    1700             : GEN
    1701           7 : FpXQE_tatepairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, GEN p)
    1702             : {
    1703           7 :   if (ell_is_inf(P) || ell_is_inf(Q))
    1704           0 :     return pol_1(get_FpX_var(T));
    1705           7 :   return FpXQE_Miller(P, Q, m, a4, T, p);
    1706             : }
    1707             : 
    1708             : /***********************************************************************/
    1709             : /**                                                                   **/
    1710             : /**                           issupersingular                         **/
    1711             : /**                                                                   **/
    1712             : /***********************************************************************/
    1713             : 
    1714             : GEN
    1715        1695 : FpXQ_ellj(GEN a4, GEN a6, GEN T, GEN p)
    1716             : {
    1717        1695 :   if (absequaliu(p,3)) return pol_0(get_FpX_var(T));
    1718             :   else
    1719             :   {
    1720        1695 :     pari_sp av=avma;
    1721        1695 :     GEN a43 = FpXQ_mul(a4,FpXQ_sqr(a4,T,p),T,p);
    1722        1695 :     GEN a62 = FpXQ_sqr(a6,T,p);
    1723        1695 :     GEN num = FpX_mulu(a43,6912,p);
    1724        1695 :     GEN den = FpX_add(FpX_mulu(a43,4,p),FpX_mulu(a62,27,p),p);
    1725        1695 :     return gerepileuptoleaf(av, FpXQ_div(num, den, T, p));
    1726             :   }
    1727             : }
    1728             : 
    1729             : int
    1730      164227 : FpXQ_elljissupersingular(GEN j, GEN T, GEN p)
    1731             : {
    1732      164227 :   pari_sp ltop = avma;
    1733             : 
    1734             :   /* All supersingular j-invariants are in FF_{p^2}, so we first check
    1735             :    * whether j is in FF_{p^2}.  If d is odd, then FF_{p^2} is not a
    1736             :    * subfield of FF_{p^d} so the j-invariants are all in FF_p.  Hence
    1737             :    * the j-invariants are in FF_{p^{2 - e}}. */
    1738      164227 :   ulong d = get_FpX_degree(T);
    1739             :   GEN S;
    1740             :   int res;
    1741             : 
    1742      164227 :   if (degpol(j) <= 0) return Fp_elljissupersingular(constant_coeff(j), p);
    1743      163786 :   if (abscmpiu(p, 5) <= 0) return 0; /* j != 0*/
    1744             : 
    1745             :   /* Set S so that FF_p[T]/(S) is isomorphic to FF_{p^2}: */
    1746      163779 :   if (d == 2)
    1747       12663 :     S = T;
    1748             :   else { /* d > 2 */
    1749             :     /* We construct FF_{p^2} = FF_p[t]/((T - j)(T - j^p)) which
    1750             :      * injects into FF_{p^d} via the map T |--> j. */
    1751      151116 :     GEN j_pow_p = FpXQ_pow(j, p, T, p);
    1752      151116 :     GEN j_sum = FpX_add(j, j_pow_p, p), j_prod;
    1753      151116 :     long var = varn(T);
    1754      151116 :     if (degpol(j_sum) > 0) { avma = ltop; return 0; /* j not in Fp^2 */ }
    1755         588 :     j_prod = FpXQ_mul(j, j_pow_p, T, p);
    1756         588 :     if (degpol(j_prod) > 0 ) { avma = ltop; return 0; /* j not in Fp^2 */ }
    1757         588 :     j_sum = constant_coeff(j_sum); j_prod = constant_coeff(j_prod);
    1758         588 :     S = mkpoln(3, gen_1, Fp_neg(j_sum, p), j_prod);
    1759         588 :     setvarn(S, var);
    1760         588 :     j = pol_x(var);
    1761             :   }
    1762       13251 :   res = jissupersingular(j, S, p);
    1763       13251 :   avma = ltop;
    1764       13251 :   return res;
    1765             : }
    1766             : 
    1767             : /***********************************************************************/
    1768             : /**                                                                   **/
    1769             : /**                           Point counting                          **/
    1770             : /**                                                                   **/
    1771             : /***********************************************************************/
    1772             : 
    1773             : GEN
    1774       13678 : elltrace_extension(GEN t, long n, GEN q)
    1775             : {
    1776       13678 :   pari_sp av = avma;
    1777       13678 :   GEN v = RgX_to_RgC(RgXQ_powu(pol_x(0), n, mkpoln(3,gen_1,negi(t),q)),2);
    1778       13678 :   GEN te = addii(shifti(gel(v,1),1), mulii(t,gel(v,2)));
    1779       13678 :   return gerepileuptoint(av, te);
    1780             : }
    1781             : 
    1782             : GEN
    1783       13083 : Fp_ffellcard(GEN a4, GEN a6, GEN q, long n, GEN p)
    1784             : {
    1785       13083 :   pari_sp av = avma;
    1786       13083 :   GEN ap = subii(addis(p, 1), Fp_ellcard(a4, a6, p));
    1787       13083 :   GEN te = elltrace_extension(ap, n, p);
    1788       13083 :   return gerepileuptoint(av, subii(addis(q, 1), te));
    1789             : }
    1790             : 
    1791             : static GEN
    1792        1687 : FpXQ_ellcardj(GEN a4, GEN a6, GEN j, GEN T, GEN q, GEN p, long n)
    1793             : {
    1794        1687 :   GEN q1 = addis(q,1);
    1795        1687 :   if (signe(j)==0)
    1796             :   {
    1797             :     GEN W, w, t, N;
    1798         560 :     if (umodiu(q,6)!=1) return q1;
    1799         420 :     N = Fp_ffellcard(gen_0,gen_1,q,n,p);
    1800         420 :     t = subii(q1, N);
    1801         420 :     W = FpXQ_pow(a6,diviuexact(shifti(q,-1), 3),T,p);
    1802         420 :     if (degpol(W)>0) /*p=5 mod 6*/
    1803         168 :       return ZX_equal1(FpXQ_powu(W,3,T,p)) ? addii(q1,shifti(t,-1)):
    1804          56 :                                              subii(q1,shifti(t,-1));
    1805         308 :     w = modii(gel(W,2),p);
    1806         308 :     if (equali1(w))  return N;
    1807         238 :     if (equalii(w,subiu(p,1))) return addii(q1,t);
    1808             :     else /*p=1 mod 6*/
    1809             :     {
    1810         168 :       GEN u = shifti(t,-1), v = sqrtint(diviuexact(subii(q,sqri(u)),3));
    1811         168 :       GEN a = addii(u,v), b = shifti(v,1);
    1812         168 :       if (equali1(Fp_powu(w,3,p)))
    1813             :       {
    1814          84 :         if (signe(Fp_add(modii(a,p),Fp_mul(w,modii(b,p),p),p))==0)
    1815          49 :           return subii(q1,subii(shifti(b,1),a));
    1816             :         else
    1817          35 :           return addii(q1,addii(a,b));
    1818             :       }
    1819             :       else
    1820             :       {
    1821          84 :         if (signe(Fp_sub(modii(a,p),Fp_mul(w,modii(b,p),p),p))==0)
    1822          49 :           return subii(q1,subii(a,shifti(b,1)));
    1823             :         else
    1824          35 :           return subii(q1,addii(a,b));
    1825             :       }
    1826             :     }
    1827        1127 :   } else if (equalii(j,modsi(1728,p)))
    1828             :   {
    1829             :     GEN w, W, N, t;
    1830         567 :     if (mod4(q)==3) return q1;
    1831         427 :     W = FpXQ_pow(a4,shifti(q,-2),T,p);
    1832         427 :     if (degpol(W)>0) return q1; /*p=3 mod 4*/
    1833         385 :     w = modii(gel(W,2),p);
    1834         385 :     N = Fp_ffellcard(gen_1,gen_0,q,n,p);
    1835         385 :     if (equali1(w)) return N;
    1836         245 :     t = subii(q1, N);
    1837         245 :     if (equalii(w,subiu(p,1))) return addii(q1,t);
    1838             :     else /*p=1 mod 4*/
    1839             :     {
    1840          98 :       GEN u = shifti(t,-1), v = sqrtint(subii(q,sqri(u)));
    1841          98 :       if (signe(Fp_add(modii(u,p),Fp_mul(w,modii(v,p),p),p))==0)
    1842          49 :         return subii(q1,shifti(v,1));
    1843             :       else
    1844          49 :         return addii(q1,shifti(v,1));
    1845             :     }
    1846             :   } else
    1847             :   {
    1848         560 :     GEN g = Fp_div(j, Fp_sub(utoi(1728), j, p), p);
    1849         560 :     GEN l = FpXQ_div(FpX_mulu(a6,3,p),FpX_mulu(a4,2,p),T,p);
    1850         560 :     GEN N = Fp_ffellcard(Fp_mulu(g,3,p),Fp_mulu(g,2,p),q,n,p);
    1851         560 :     if (FpXQ_issquare(l,T,p)) return N;
    1852         280 :     return subii(shifti(q1,1),N);
    1853             :   }
    1854             : }
    1855             : 
    1856             : GEN
    1857        3662 : FpXQ_ellcard(GEN a4, GEN a6, GEN T, GEN p)
    1858             : {
    1859        3662 :   pari_sp av = avma;
    1860        3662 :   long n = get_FpX_degree(T);
    1861        3662 :   GEN q = powiu(p, n), r, J;
    1862        3662 :   if (degpol(a4)<=0 && degpol(a6)<=0)
    1863         133 :     r = Fp_ffellcard(constant_coeff(a4),constant_coeff(a6),q,n,p);
    1864        3529 :   else if (lgefint(p)==3)
    1865             :   {
    1866        1834 :     ulong pp = p[2];
    1867        1834 :     r =  Flxq_ellcard(ZX_to_Flx(a4,pp),ZX_to_Flx(a6,pp),ZX_to_Flx(T,pp),pp);
    1868             :   }
    1869        1695 :   else if (degpol(J=FpXQ_ellj(a4,a6,T,p))<=0)
    1870        1687 :     r = FpXQ_ellcardj(a4,a6,constant_coeff(J),T,q,p,n);
    1871             :   else
    1872           8 :     r = Fq_ellcard_SEA(a4, a6, q, T, p, 0);
    1873        3662 :   return gerepileuptoint(av, r);
    1874             : }
    1875             : 
    1876             : static GEN
    1877          28 : _FpXQE_pairorder(void *E, GEN P, GEN Q, GEN m, GEN F)
    1878             : {
    1879          28 :   struct _FpXQE *e = (struct _FpXQE *) E;
    1880          28 :   return  FpXQ_order(FpXQE_weilpairing(P,Q,m,e->a4,e->T,e->p), F, e->T, e->p);
    1881             : }
    1882             : 
    1883             : GEN
    1884          14 : FpXQ_ellgroup(GEN a4, GEN a6, GEN N, GEN T, GEN p, GEN *pt_m)
    1885             : {
    1886             :   struct _FpXQE e;
    1887          14 :   GEN q = powiu(p, get_FpX_degree(T));
    1888          14 :   e.a4=a4; e.a6=a6; e.T=T; e.p=p;
    1889          14 :   return gen_ellgroup(N, subis(q,1), pt_m, (void*)&e, &FpXQE_group, _FpXQE_pairorder);
    1890             : }
    1891             : 
    1892             : GEN
    1893           7 : FpXQ_ellgens(GEN a4, GEN a6, GEN ch, GEN D, GEN m, GEN T, GEN p)
    1894             : {
    1895             :   GEN P;
    1896           7 :   pari_sp av = avma;
    1897             :   struct _FpXQE e;
    1898           7 :   e.a4=a4; e.a6=a6; e.T=T; e.p=p;
    1899           7 :   switch(lg(D)-1)
    1900             :   {
    1901             :   case 1:
    1902           7 :     P = gen_gener(gel(D,1), (void*)&e, &FpXQE_group);
    1903           7 :     P = mkvec(FpXQE_changepoint(P, ch, T, p));
    1904           7 :     break;
    1905             :   default:
    1906           0 :     P = gen_ellgens(gel(D,1), gel(D,2), m, (void*)&e, &FpXQE_group, _FpXQE_pairorder);
    1907           0 :     gel(P,1) = FpXQE_changepoint(gel(P,1), ch, T, p);
    1908           0 :     gel(P,2) = FpXQE_changepoint(gel(P,2), ch, T, p);
    1909           0 :     break;
    1910             :   }
    1911           7 :   return gerepilecopy(av, P);
    1912             : }
    1913             : 
    1914             : 

Generated by: LCOV version 1.11