Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - FpE.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.16.1 lcov report (development 28923-fbb3cf3f02) Lines: 1093 1195 91.5 %
Date: 2023-12-11 07:56:23 Functions: 121 129 93.8 %
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; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : #define DEBUGLEVEL DEBUGLEVEL_ellcard
      19             : 
      20             : /* Not so fast arithmetic with points over elliptic curves over Fp */
      21             : 
      22             : /***********************************************************************/
      23             : /**                                                                   **/
      24             : /**                              FpJ                                  **/
      25             : /**                                                                   **/
      26             : /***********************************************************************/
      27             : 
      28             : /* Arithmetic is implemented using Jacobian coordinates, representing
      29             :  * a projective point (x : y : z) on E by [z*x , z^2*y , z].  This is
      30             :  * probably not the fastest representation available for the given
      31             :  * problem, but they're easy to implement and up to 60% faster than
      32             :  * the school-book method used in FpE_mulu().
      33             :  */
      34             : 
      35             : static GEN
      36       50582 : ellinf_FpJ(void)
      37       50582 : { return mkvec3(gen_1, gen_1, gen_0); }
      38             : 
      39             : /*
      40             :  * Cost: 1M + 8S + 1*a + 10add + 1*8 + 2*2 + 1*3.
      41             :  * Source: http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#doubling-dbl-2007-bl
      42             :  */
      43             : 
      44             : GEN
      45     6641094 : FpJ_dbl(GEN P, GEN a4, GEN p)
      46             : {
      47             :   GEN X1, Y1, Z1;
      48             :   GEN XX, YY, YYYY, ZZ, S, M, T, Q;
      49             : 
      50     6641094 :   if (signe(gel(P,3)) == 0) return ellinf_FpJ();
      51             : 
      52     6630942 :   X1 = gel(P,1); Y1 = gel(P,2); Z1 = gel(P,3);
      53             : 
      54     6630942 :   XX = Fp_sqr(X1, p);
      55     6589420 :   YY = Fp_sqr(Y1, p);
      56     6584858 :   YYYY = Fp_sqr(YY, p);
      57     6581182 :   ZZ = Fp_sqr(Z1, p);
      58     6581912 :   S = Fp_double(Fp_sub(Fp_sqr(Fp_add(X1,YY,p), p), Fp_add(XX,YYYY,p), p), p);
      59     6540050 :   M = Fp_addmul(Fp_mulu(XX, 3, p), a4, Fp_sqr(ZZ, p),  p);
      60     6614088 :   T = Fp_sub(Fp_sqr(M, p), Fp_double(S, p), p);
      61     6614230 :   Q = cgetg(4, t_VEC);
      62     6598977 :   gel(Q,1) = T;
      63     6598977 :   gel(Q,2) = Fp_sub(Fp_mul(M, Fp_sub(S, T, p), p), Fp_mulu(YYYY, 8, p), p);
      64     6618026 :   gel(Q,3) = Fp_sub(Fp_sqr(Fp_add(Y1, Z1, p), p), Fp_add(YY, ZZ, p), p);
      65     6612611 :   return Q;
      66             : }
      67             : 
      68             : /*
      69             :  * Cost: 11M + 5S + 9add + 4*2.
      70             :  * Source: http://www.hyperelliptic.org/EFD/g1p/auto-shortw-jacobian.html#addition-add-2007-bl
      71             :  */
      72             : 
      73             : GEN
      74     1186489 : FpJ_add(GEN P, GEN Q, GEN a4, GEN p)
      75             : {
      76             :   GEN X1, Y1, Z1, X2, Y2, Z2;
      77             :   GEN Z1Z1, Z2Z2, U1, U2, S1, S2, H, I, J, r, V, W, R;
      78             : 
      79     1186489 :   if (signe(gel(Q,3)) == 0) return gcopy(P);
      80     1186489 :   if (signe(gel(P,3)) == 0) return gcopy(Q);
      81             : 
      82     1166636 :   X1 = gel(P,1); Y1 = gel(P,2); Z1 = gel(P,3);
      83     1166636 :   X2 = gel(Q,1); Y2 = gel(Q,2); Z2 = gel(Q,3);
      84             : 
      85     1166636 :   Z1Z1 = Fp_sqr(Z1, p);
      86     1166354 :   Z2Z2 = Fp_sqr(Z2, p);
      87     1166048 :   U1 = Fp_mul(X1, Z2Z2, p);
      88     1166177 :   U2 = Fp_mul(X2, Z1Z1, p);
      89     1166147 :   S1 = mulii(Y1, Fp_mul(Z2, Z2Z2, p));
      90     1165819 :   S2 = mulii(Y2, Fp_mul(Z1, Z1Z1, p));
      91     1166532 :   H = Fp_sub(U2, U1, p);
      92     1166793 :   r = Fp_double(Fp_sub(S2, S1, p), p);
      93             : 
      94             :   /* If points are equal we must double. */
      95     1165364 :   if (signe(H)== 0) {
      96       41620 :     if (signe(r) == 0)
      97             :       /* Points are equal so double. */
      98        1190 :       return FpJ_dbl(P, a4, p);
      99             :     else
     100       40430 :       return ellinf_FpJ();
     101             :   }
     102     1123744 :   I = Fp_sqr(Fp_double(H, p), p);
     103     1124840 :   J = Fp_mul(H, I, p);
     104     1124612 :   V = Fp_mul(U1, I, p);
     105     1124632 :   W = Fp_sub(Fp_sqr(r, p), Fp_add(J, Fp_double(V, p), p), p);
     106     1125125 :   R = cgetg(4, t_VEC);
     107     1124882 :   gel(R,1) = W;
     108     1124882 :   gel(R,2) = Fp_sub(mulii(r, subii(V, W)),
     109             :                     shifti(mulii(S1, J), 1), p);
     110     1125271 :   gel(R,3) = Fp_mul(Fp_sub(Fp_sqr(Fp_add(Z1, Z2, p), p),
     111             :                            Fp_add(Z1Z1, Z2Z2, p), p), H, p);
     112     1124632 :   return R;
     113             : }
     114             : 
     115             : GEN
     116           0 : FpJ_neg(GEN Q, GEN p)
     117             : {
     118           0 :   return mkvec3(icopy(gel(Q,1)), Fp_neg(gel(Q,2), p), icopy(gel(Q,3)));
     119             : }
     120             : 
     121             : GEN
     122      212925 : FpE_to_FpJ(GEN P)
     123             : {
     124      212925 :   return ell_is_inf(P) ? ellinf_FpJ()
     125      212925 :        : mkvec3(icopy(gel(P,1)),icopy(gel(P,2)), gen_1);
     126             : }
     127             : 
     128             : GEN
     129      212360 : FpJ_to_FpE(GEN P, GEN p)
     130             : {
     131      212360 :   if (signe(gel(P,3)) == 0) return ellinf();
     132             :   else
     133             :   {
     134      171745 :     GEN Z = Fp_inv(gel(P,3), p);
     135      171719 :     GEN Z2 = Fp_sqr(Z, p), Z3 = Fp_mul(Z, Z2, p);
     136      171719 :     retmkvec2(Fp_mul(gel(P,1), Z2, p), Fp_mul(gel(P,2), Z3, p));
     137             :   }
     138             : }
     139             : 
     140             : struct _FpE { GEN p,a4,a6; };
     141             : static GEN
     142     6638195 : _FpJ_dbl(void *E, GEN P)
     143             : {
     144     6638195 :   struct _FpE *ell = (struct _FpE *) E;
     145     6638195 :   return FpJ_dbl(P, ell->a4, ell->p);
     146             : }
     147             : static GEN
     148     1186396 : _FpJ_add(void *E, GEN P, GEN Q)
     149             : {
     150     1186396 :   struct _FpE *ell=(struct _FpE *) E;
     151     1186396 :   return FpJ_add(P, Q, ell->a4, ell->p);
     152             : }
     153             : static GEN
     154        5854 : _FpJ_mul(void *E, GEN P, GEN n)
     155             : {
     156        5854 :   pari_sp av = avma;
     157        5854 :   struct _FpE *e=(struct _FpE *) E;
     158        5854 :   long s = signe(n);
     159        5854 :   if (!s || signe(gel(P,3))==0) return ellinf_FpJ();
     160        5854 :   if (s < 0) P = FpJ_neg(P, e->p);
     161        5854 :   if (is_pm1(n)) return s > 0 ? gcopy(P): P;
     162        5854 :   return gerepilecopy(av, gen_pow_i(P, n, e, &_FpJ_dbl, &_FpJ_add));
     163             : }
     164             : 
     165             : GEN
     166        5854 : FpJ_mul(GEN P, GEN n, GEN a4, GEN p)
     167             : {
     168             :   struct _FpE E;
     169        5854 :   E.a4= a4; E.p = p;
     170        5854 :   return _FpJ_mul(&E, P, n);
     171             : }
     172             : 
     173             : /***********************************************************************/
     174             : /**                                                                   **/
     175             : /**                              FpE                                  **/
     176             : /**                                                                   **/
     177             : /***********************************************************************/
     178             : 
     179             : /* These functions deal with point over elliptic curves over Fp defined
     180             :  * by an equation of the form y^2=x^3+a4*x+a6.
     181             :  * Most of the time a6 is omitted since it can be recovered from any point
     182             :  * on the curve.
     183             :  */
     184             : 
     185             : GEN
     186        2815 : RgE_to_FpE(GEN x, GEN p)
     187             : {
     188        2815 :   if (ell_is_inf(x)) return x;
     189        2816 :   retmkvec2(Rg_to_Fp(gel(x,1),p),Rg_to_Fp(gel(x,2),p));
     190             : }
     191             : 
     192             : GEN
     193        1104 : FpE_to_mod(GEN x, GEN p)
     194             : {
     195        1104 :   if (ell_is_inf(x)) return x;
     196        1041 :   retmkvec2(Fp_to_mod(gel(x,1),p),Fp_to_mod(gel(x,2),p));
     197             : }
     198             : 
     199             : GEN
     200        1776 : FpE_changepoint(GEN P, GEN ch, GEN p)
     201             : {
     202        1776 :   pari_sp av = avma;
     203             :   GEN c, z, u, r, s, t, v, v2, v3;
     204        1776 :   if (ell_is_inf(P)) return P;
     205        1713 :   if (lgefint(p) == 3)
     206             :   {
     207         719 :     ulong pp = p[2];
     208         719 :     z = Fle_changepoint(ZV_to_Flv(P, pp), ZV_to_Flv(ch, pp), pp);
     209         719 :     return gerepileupto(av, Flv_to_ZV(z));
     210             :   }
     211         994 :   u = gel(ch,1); r = gel(ch,2); s = gel(ch,3); t = gel(ch,4);
     212         994 :   v = Fp_inv(u, p); v2 = Fp_sqr(v,p); v3 = Fp_mul(v,v2,p);
     213         994 :   c = Fp_sub(gel(P,1),r,p);
     214         994 :   z = cgetg(3,t_VEC);
     215         994 :   gel(z,1) = Fp_mul(v2, c, p);
     216         994 :   gel(z,2) = Fp_mul(v3, Fp_sub(gel(P,2), Fp_add(Fp_mul(s,c, p),t, p),p),p);
     217         994 :   return gerepileupto(av, z);
     218             : }
     219             : 
     220             : GEN
     221        2816 : FpE_changepointinv(GEN P, GEN ch, GEN p)
     222             : {
     223             :   GEN u, r, s, t, u2, u3, c, z;
     224        2816 :   if (ell_is_inf(P)) return P;
     225        2816 :   if (lgefint(p) == 3)
     226             :   {
     227        1738 :     ulong pp = p[2];
     228        1738 :     z = Fle_changepointinv(ZV_to_Flv(P, pp), ZV_to_Flv(ch, pp), pp);
     229        1738 :     return Flv_to_ZV(z);
     230             :   }
     231        1078 :   u = gel(ch,1); r = gel(ch,2); s = gel(ch,3); t = gel(ch,4);
     232        1078 :   u2 = Fp_sqr(u, p); u3 = Fp_mul(u,u2,p);
     233        1077 :   c = Fp_mul(u2, gel(P,1), p);
     234        1076 :   z = cgetg(3, t_VEC);
     235        1076 :   gel(z,1) = Fp_add(c,r,p);
     236        1076 :   gel(z,2) = Fp_add(Fp_mul(u3,gel(P,2),p), Fp_add(Fp_mul(s,c,p), t, p), p);
     237        1077 :   return z;
     238             : }
     239             : 
     240             : static GEN
     241         420 : random_nonsquare_Fp(GEN p)
     242             : {
     243         420 :   pari_sp av = avma;
     244             :   GEN a;
     245         420 :   switch(mod8(p))
     246             :   { /* easy special cases */
     247         420 :     case 3: case 5: return gen_2;
     248           0 :     case 7: return subiu(p, 1);
     249             :   }
     250             :   do
     251             :   {
     252           0 :     set_avma(av);
     253           0 :     a = randomi(p);
     254           0 :   } while (kronecker(a, p) >= 0);
     255           0 :   return a;
     256             : }
     257             : 
     258             : void
     259           0 : Fp_elltwist(GEN a4, GEN a6, GEN p, GEN *pt_a4, GEN *pt_a6)
     260             : {
     261           0 :   GEN d = random_nonsquare_Fp(p), d2 = Fp_sqr(d, p), d3 = Fp_mul(d2, d, p);
     262           0 :   *pt_a4 = Fp_mul(a4, d2, p);
     263           0 :   *pt_a6 = Fp_mul(a6, d3, p);
     264           0 : }
     265             : 
     266             : static GEN
     267      264754 : FpE_dbl_slope(GEN P, GEN a4, GEN p, GEN *slope)
     268             : {
     269             :   GEN x, y, Q;
     270      264754 :   if (ell_is_inf(P) || !signe(gel(P,2))) return ellinf();
     271      133101 :   x = gel(P,1); y = gel(P,2);
     272      133101 :   *slope = Fp_div(Fp_add(Fp_mulu(Fp_sqr(x,p), 3, p), a4, p),
     273             :                   Fp_mulu(y, 2, p), p);
     274      133101 :   Q = cgetg(3,t_VEC);
     275      133101 :   gel(Q, 1) = Fp_sub(Fp_sqr(*slope, p), Fp_mulu(x, 2, p), p);
     276      133101 :   gel(Q, 2) = Fp_sub(Fp_mul(*slope, Fp_sub(x, gel(Q, 1), p), p), y, p);
     277      133101 :   return Q;
     278             : }
     279             : 
     280             : GEN
     281      264161 : FpE_dbl(GEN P, GEN a4, GEN p)
     282             : {
     283      264161 :   pari_sp av = avma;
     284             :   GEN slope;
     285      264161 :   return gerepileupto(av, FpE_dbl_slope(P,a4,p,&slope));
     286             : }
     287             : 
     288             : static GEN
     289      916537 : FpE_add_slope(GEN P, GEN Q, GEN a4, GEN p, GEN *slope)
     290             : {
     291             :   GEN Px, Py, Qx, Qy, R;
     292      916537 :   if (ell_is_inf(P)) return Q;
     293      916047 :   if (ell_is_inf(Q)) return P;
     294      916047 :   Px = gel(P,1); Py = gel(P,2);
     295      916047 :   Qx = gel(Q,1); Qy = gel(Q,2);
     296      916047 :   if (equalii(Px, Qx))
     297             :   {
     298         574 :     if (equalii(Py, Qy))
     299         553 :       return FpE_dbl_slope(P, a4, p, slope);
     300             :     else
     301          21 :       return ellinf();
     302             :   }
     303      915473 :   *slope = Fp_div(Fp_sub(Py, Qy, p), Fp_sub(Px, Qx, p), p);
     304      915473 :   R = cgetg(3,t_VEC);
     305      915473 :   gel(R, 1) = Fp_sub(Fp_sub(Fp_sqr(*slope, p), Px, p), Qx, p);
     306      915473 :   gel(R, 2) = Fp_sub(Fp_mul(*slope, Fp_sub(Px, gel(R, 1), p), p), Py, p);
     307      915473 :   return R;
     308             : }
     309             : 
     310             : GEN
     311      916535 : FpE_add(GEN P, GEN Q, GEN a4, GEN p)
     312             : {
     313      916535 :   pari_sp av = avma;
     314             :   GEN slope;
     315      916535 :   return gerepileupto(av, FpE_add_slope(P,Q,a4,p,&slope));
     316             : }
     317             : 
     318             : static GEN
     319           0 : FpE_neg_i(GEN P, GEN p)
     320             : {
     321           0 :   if (ell_is_inf(P)) return P;
     322           0 :   return mkvec2(gel(P,1), Fp_neg(gel(P,2), p));
     323             : }
     324             : 
     325             : GEN
     326      362489 : FpE_neg(GEN P, GEN p)
     327             : {
     328      362489 :   if (ell_is_inf(P)) return ellinf();
     329      362489 :   return mkvec2(gcopy(gel(P,1)), Fp_neg(gel(P,2), p));
     330             : }
     331             : 
     332             : GEN
     333           0 : FpE_sub(GEN P, GEN Q, GEN a4, GEN p)
     334             : {
     335           0 :   pari_sp av = avma;
     336             :   GEN slope;
     337           0 :   return gerepileupto(av, FpE_add_slope(P, FpE_neg_i(Q, p), a4, p, &slope));
     338             : }
     339             : 
     340             : static GEN
     341      264161 : _FpE_dbl(void *E, GEN P)
     342             : {
     343      264161 :   struct _FpE *ell = (struct _FpE *) E;
     344      264161 :   return FpE_dbl(P, ell->a4, ell->p);
     345             : }
     346             : 
     347             : static GEN
     348      897344 : _FpE_add(void *E, GEN P, GEN Q)
     349             : {
     350      897344 :   struct _FpE *ell=(struct _FpE *) E;
     351      897344 :   return FpE_add(P, Q, ell->a4, ell->p);
     352             : }
     353             : 
     354             : static GEN
     355      916497 : _FpE_mul(void *E, GEN P, GEN n)
     356             : {
     357      916497 :   pari_sp av = avma;
     358      916497 :   struct _FpE *e=(struct _FpE *) E;
     359      916497 :   long s = signe(n);
     360             :   GEN Q;
     361      916497 :   if (!s || ell_is_inf(P)) return ellinf();
     362      916462 :   if (s<0) P = FpE_neg(P, e->p);
     363      916462 :   if (is_pm1(n)) return s>0? gcopy(P): P;
     364      476582 :   if (equalis(n,2)) return _FpE_dbl(E, P);
     365      212421 :   Q = gen_pow_i(FpE_to_FpJ(P), n, e, &_FpJ_dbl, &_FpJ_add);
     366      212360 :   return gerepileupto(av, FpJ_to_FpE(Q, e->p));
     367             : }
     368             : 
     369             : GEN
     370        1396 : FpE_mul(GEN P, GEN n, GEN a4, GEN p)
     371             : {
     372             :   struct _FpE E;
     373        1396 :   E.a4 = a4; E.p = p;
     374        1396 :   return _FpE_mul(&E, P, n);
     375             : }
     376             : 
     377             : /* Finds a random nonsingular point on E */
     378             : 
     379             : GEN
     380      194737 : random_FpE(GEN a4, GEN a6, GEN p)
     381             : {
     382      194737 :   pari_sp ltop = avma;
     383             :   GEN x, x2, y, rhs;
     384             :   do
     385             :   {
     386      343347 :     set_avma(ltop);
     387      343347 :     x   = randomi(p); /*  x^3+a4*x+a6 = x*(x^2+a4)+a6  */
     388      343347 :     x2  = Fp_sqr(x, p);
     389      343347 :     rhs = Fp_add(Fp_mul(x, Fp_add(x2, a4, p), p), a6, p);
     390       36325 :   } while ((!signe(rhs) && !signe(Fp_add(Fp_mulu(x2,3,p),a4,p)))
     391      379672 :           || kronecker(rhs, p) < 0);
     392      194737 :   y = Fp_sqrt(rhs, p);
     393      194737 :   if (!y) pari_err_PRIME("random_FpE", p);
     394      194737 :   return gerepilecopy(ltop, mkvec2(x, y));
     395             : }
     396             : 
     397             : static GEN
     398      192286 : _FpE_rand(void *E)
     399             : {
     400      192286 :   struct _FpE *e=(struct _FpE *) E;
     401      192286 :   return random_FpE(e->a4, e->a6, e->p);
     402             : }
     403             : 
     404             : static const struct bb_group FpE_group={_FpE_add,_FpE_mul,_FpE_rand,hash_GEN,ZV_equal,ell_is_inf,NULL};
     405             : 
     406             : const struct bb_group *
     407         903 : get_FpE_group(void ** pt_E, GEN a4, GEN a6, GEN p)
     408             : {
     409         903 :   struct _FpE *e = (struct _FpE *) stack_malloc(sizeof(struct _FpE));
     410         903 :   e->a4 = a4; e->a6 = a6; e->p  = p;
     411         903 :   *pt_E = (void *) e;
     412         903 :   return &FpE_group;
     413             : }
     414             : 
     415             : GEN
     416         736 : FpE_order(GEN z, GEN o, GEN a4, GEN p)
     417             : {
     418         736 :   pari_sp av = avma;
     419             :   struct _FpE e;
     420             :   GEN r;
     421         736 :   if (lgefint(p) == 3)
     422             :   {
     423         630 :     ulong pp = p[2];
     424         630 :     r = Fle_order(ZV_to_Flv(z, pp), o, umodiu(a4,pp), pp);
     425             :   }
     426             :   else
     427             :   {
     428         106 :     e.a4 = a4;
     429         106 :     e.p = p;
     430         106 :     r = gen_order(z, o, (void*)&e, &FpE_group);
     431             :   }
     432         736 :   return gerepileuptoint(av, r);
     433             : }
     434             : 
     435             : GEN
     436          49 : FpE_log(GEN a, GEN b, GEN o, GEN a4, GEN p)
     437             : {
     438          49 :   pari_sp av = avma;
     439             :   struct _FpE e;
     440             :   GEN r;
     441          49 :   if (lgefint(p) == 3)
     442             :   {
     443          49 :     ulong pp = p[2];
     444          49 :     r = Fle_log(ZV_to_Flv(a,pp), ZV_to_Flv(b,pp), o, umodiu(a4,pp), pp);
     445             :   }
     446             :   else
     447             :   {
     448           0 :     e.a4 = a4;
     449           0 :     e.p = p;
     450           0 :     r = gen_PH_log(a, b, o, (void*)&e, &FpE_group);
     451             :   }
     452          49 :   return gerepileuptoint(av, r);
     453             : }
     454             : 
     455             : /***********************************************************************/
     456             : /**                                                                   **/
     457             : /**                            Pairings                               **/
     458             : /**                                                                   **/
     459             : /***********************************************************************/
     460             : 
     461             : /* Derived from APIP from and by Jerome Milan, 2012 */
     462             : 
     463             : static GEN
     464         146 : FpE_vert(GEN P, GEN Q, GEN a4, GEN p)
     465             : {
     466         146 :   if (ell_is_inf(P))
     467          58 :     return gen_1;
     468          88 :   if (!equalii(gel(Q, 1), gel(P, 1)))
     469          80 :     return Fp_sub(gel(Q, 1), gel(P, 1), p);
     470           8 :   if (signe(gel(P,2))!=0) return gen_1;
     471           6 :   return Fp_inv(Fp_add(Fp_mulu(Fp_sqr(gel(P,1),p), 3, p), a4, p), p);
     472             : }
     473             : 
     474             : static GEN
     475          42 : FpE_Miller_line(GEN R, GEN Q, GEN slope, GEN a4, GEN p)
     476             : {
     477          42 :   GEN x = gel(Q, 1), y = gel(Q, 2);
     478          42 :   GEN tmp1 = Fp_sub(x, gel(R, 1), p);
     479          42 :   GEN tmp2 = Fp_add(Fp_mul(tmp1, slope, p), gel(R,2), p);
     480          42 :   if (!equalii(y, tmp2))
     481          37 :     return Fp_sub(y, tmp2, p);
     482           5 :   if (signe(y) == 0)
     483           3 :     return gen_1;
     484             :   else
     485             :   {
     486             :     GEN s1, s2;
     487           2 :     GEN y2i = Fp_inv(Fp_mulu(y, 2, p), p);
     488           2 :     s1 = Fp_mul(Fp_add(Fp_mulu(Fp_sqr(x, p), 3, p), a4, p), y2i, p);
     489           2 :     if (!equalii(s1, slope))
     490           2 :       return Fp_sub(s1, slope, p);
     491           0 :     s2 = Fp_mul(Fp_sub(Fp_mulu(x, 3, p), Fp_sqr(s1, p), p), y2i, p);
     492           0 :     return signe(s2)!=0 ? s2: y2i;
     493             :   }
     494             : }
     495             : 
     496             : /* Computes the equation of the line tangent to R and returns its
     497             :    evaluation at the point Q. Also doubles the point R.
     498             :  */
     499             : 
     500             : static GEN
     501          98 : FpE_tangent_update(GEN R, GEN Q, GEN a4, GEN p, GEN *pt_R)
     502             : {
     503          98 :   if (ell_is_inf(R))
     504             :   {
     505          12 :     *pt_R = ellinf();
     506          12 :     return gen_1;
     507             :   }
     508          86 :   else if (signe(gel(R,2)) == 0)
     509             :   {
     510          46 :     *pt_R = ellinf();
     511          46 :     return FpE_vert(R, Q, a4, p);
     512             :   } else {
     513             :     GEN slope;
     514          40 :     *pt_R = FpE_dbl_slope(R, a4, p, &slope);
     515          40 :     return FpE_Miller_line(R, Q, slope, a4, p);
     516             :   }
     517             : }
     518             : 
     519             : /* Computes the equation of the line through R and P, and returns its
     520             :    evaluation at the point Q. Also adds P to the point R.
     521             :  */
     522             : 
     523             : static GEN
     524           2 : FpE_chord_update(GEN R, GEN P, GEN Q, GEN a4, GEN p, GEN *pt_R)
     525             : {
     526           2 :   if (ell_is_inf(R))
     527             :   {
     528           0 :     *pt_R = gcopy(P);
     529           0 :     return FpE_vert(P, Q, a4, p);
     530             :   }
     531           2 :   else if (ell_is_inf(P))
     532             :   {
     533           0 :     *pt_R = gcopy(R);
     534           0 :     return FpE_vert(R, Q, a4, p);
     535             :   }
     536           2 :   else if (equalii(gel(P, 1), gel(R, 1)))
     537             :   {
     538           0 :     if (equalii(gel(P, 2), gel(R, 2)))
     539           0 :       return FpE_tangent_update(R, Q, a4, p, pt_R);
     540             :     else {
     541           0 :       *pt_R = ellinf();
     542           0 :       return FpE_vert(R, Q, a4, p);
     543             :     }
     544             :   } else {
     545             :     GEN slope;
     546           2 :     *pt_R = FpE_add_slope(P, R, a4, p, &slope);
     547           2 :     return FpE_Miller_line(R, Q, slope, a4, p);
     548             :   }
     549             : }
     550             : 
     551             : struct _FpE_miller { GEN p, a4, P; };
     552             : static GEN
     553          98 : FpE_Miller_dbl(void* E, GEN d)
     554             : {
     555          98 :   struct _FpE_miller *m = (struct _FpE_miller *)E;
     556          98 :   GEN p = m->p, a4 = m->a4, P = m->P;
     557             :   GEN v, line;
     558          98 :   GEN N = Fp_sqr(gel(d,1), p);
     559          98 :   GEN D = Fp_sqr(gel(d,2), p);
     560          98 :   GEN point = gel(d,3);
     561          98 :   line = FpE_tangent_update(point, P, a4, p, &point);
     562          98 :   N  = Fp_mul(N, line, p);
     563          98 :   v = FpE_vert(point, P, a4, p);
     564          98 :   D = Fp_mul(D, v, p); return mkvec3(N, D, point);
     565             : }
     566             : static GEN
     567           2 : FpE_Miller_add(void* E, GEN va, GEN vb)
     568             : {
     569           2 :   struct _FpE_miller *m = (struct _FpE_miller *)E;
     570           2 :   GEN p = m->p, a4= m->a4, P = m->P;
     571             :   GEN v, line, point;
     572           2 :   GEN na = gel(va,1), da = gel(va,2), pa = gel(va,3);
     573           2 :   GEN nb = gel(vb,1), db = gel(vb,2), pb = gel(vb,3);
     574           2 :   GEN N = Fp_mul(na, nb, p);
     575           2 :   GEN D = Fp_mul(da, db, p);
     576           2 :   line = FpE_chord_update(pa, pb, P, a4, p, &point);
     577           2 :   N = Fp_mul(N, line, p);
     578           2 :   v = FpE_vert(point, P, a4, p);
     579           2 :   D = Fp_mul(D, v, p); return mkvec3(N, D, point);
     580             : }
     581             : 
     582             : /* Returns the Miller function f_{m, Q} evaluated at the point P using
     583             :  * the standard Miller algorithm. */
     584             : static GEN
     585          46 : FpE_Miller(GEN Q, GEN P, GEN m, GEN a4, GEN p)
     586             : {
     587          46 :   pari_sp av = avma;
     588             :   struct _FpE_miller d;
     589             :   GEN v, N, D;
     590             : 
     591          46 :   d.a4 = a4; d.p = p; d.P = P;
     592          46 :   v = gen_pow_i(mkvec3(gen_1,gen_1,Q), m, (void*)&d,
     593             :                 FpE_Miller_dbl, FpE_Miller_add);
     594          46 :   N = gel(v,1); D = gel(v,2);
     595          46 :   return gerepileuptoint(av, Fp_div(N, D, p));
     596             : }
     597             : 
     598             : GEN
     599       76011 : FpE_weilpairing(GEN P, GEN Q, GEN m, GEN a4, GEN p)
     600             : {
     601       76011 :   pari_sp av = avma;
     602             :   GEN N, D, w;
     603       76011 :   if (ell_is_inf(P) || ell_is_inf(Q) || ZV_equal(P,Q)) return gen_1;
     604       51289 :   if (lgefint(p)==3 && lgefint(m)==3)
     605             :   {
     606       51266 :     ulong pp = p[2];
     607       51266 :     GEN Pp = ZV_to_Flv(P, pp), Qp = ZV_to_Flv(Q, pp);
     608       51266 :     ulong w = Fle_weilpairing(Pp, Qp, itou(m), umodiu(a4, pp), pp);
     609       51266 :     return gc_utoi(av, w);
     610             :   }
     611          23 :   N = FpE_Miller(P, Q, m, a4, p);
     612          23 :   D = FpE_Miller(Q, P, m, a4, p);
     613          23 :   w = Fp_div(N, D, p);
     614          23 :   if (mpodd(m)) w  = Fp_neg(w, p);
     615          23 :   return gerepileuptoint(av, w);
     616             : }
     617             : 
     618             : GEN
     619         203 : FpE_tatepairing(GEN P, GEN Q, GEN m, GEN a4, GEN p)
     620             : {
     621         203 :   if (ell_is_inf(P) || ell_is_inf(Q)) return gen_1;
     622         203 :   if (lgefint(p)==3 && lgefint(m)==3)
     623             :   {
     624         203 :     pari_sp av = avma;
     625         203 :     ulong pp = p[2];
     626         203 :     GEN Pp = ZV_to_Flv(P, pp), Qp = ZV_to_Flv(Q, pp);
     627         203 :     ulong w = Fle_tatepairing(Pp, Qp, itou(m), umodiu(a4, pp), pp);
     628         203 :     return gc_utoi(av,w);
     629             :   }
     630           0 :   return FpE_Miller(P, Q, m, a4, p);
     631             : }
     632             : 
     633             : /***********************************************************************/
     634             : /**                                                                   **/
     635             : /**                   CM by principal order                           **/
     636             : /**                                                                   **/
     637             : /***********************************************************************/
     638             : 
     639             : /* is jn/jd = J (mod p) */
     640             : static int
     641      652122 : is_CMj(long J, GEN jn, GEN jd, GEN p)
     642      652122 : { return dvdii(subii(mulis(jd,J), jn), p); }
     643             : #ifndef LONG_IS_64BIT
     644             : /* is jn/jd = -(2^32 a + b) (mod p) */
     645             : static int
     646       14407 : u2_is_CMj(ulong a, ulong b, GEN jn, GEN jd, GEN p)
     647             : {
     648       14407 :   GEN mJ = uu32toi(a,b);
     649       14407 :   return dvdii(addii(mulii(jd,mJ), jn), p);
     650             : }
     651             : #endif
     652             : 
     653             : static long
     654       52611 : Fp_ellj_get_CM(GEN jn, GEN jd, GEN p)
     655             : {
     656             : #define CHECK(CM,J) if (is_CMj(J,jn,jd,p)) return CM;
     657       52611 :   CHECK(-3,  0);
     658       52492 :   CHECK(-4,  1728);
     659       52383 :   CHECK(-7,  -3375);
     660       52137 :   CHECK(-8,  8000);
     661       51902 :   CHECK(-11, -32768);
     662       51653 :   CHECK(-12, 54000);
     663       51433 :   CHECK(-16, 287496);
     664       51247 :   CHECK(-19, -884736);
     665       51022 :   CHECK(-27, -12288000);
     666       50806 :   CHECK(-28, 16581375);
     667       50627 :   CHECK(-43, -884736000);
     668             : #ifdef LONG_IS_64BIT
     669       43255 :   CHECK(-67, -147197952000L);
     670       43122 :   CHECK(-163, -262537412640768000L);
     671             : #else
     672        7214 :   if (u2_is_CMj(0x00000022UL,0x45ae8000UL,jn,jd,p)) return -67;
     673        7193 :   if (u2_is_CMj(0x03a4b862UL,0xc4b40000UL,jn,jd,p)) return -163;
     674             : #endif
     675             : #undef CHECK
     676       50157 :   return 0;
     677             : }
     678             : 
     679             : /***********************************************************************/
     680             : /**                                                                   **/
     681             : /**                            issupersingular                        **/
     682             : /**                                                                   **/
     683             : /***********************************************************************/
     684             : 
     685             : /* assume x reduced mod p, monic. Return one root, or NULL if irreducible */
     686             : static GEN
     687       73563 : FqX_quad_root(GEN x, GEN T, GEN p)
     688             : {
     689       73563 :   GEN b = gel(x,3), c = gel(x,2);
     690       73563 :   GEN D = Fq_sub(Fq_sqr(b, T, p), Fq_mulu(c,4, T, p), T, p);
     691       73563 :   GEN s = Fq_sqrt(D,T, p);
     692       73563 :   if (!s) return NULL;
     693       70385 :   return Fq_halve(Fq_sub(s, b, T, p), T, p);
     694             : }
     695             : 
     696             : static GEN
     697        1229 : FpX_quad_root(GEN x, GEN p)
     698             : {
     699        1229 :   GEN s, b = gel(x,3), c = gel(x,2);
     700        1229 :   GEN D = Fp_sub(Fp_sqr(b, p), shifti(c,2), p);
     701        1229 :   if (kronecker(D,p) == -1) return NULL;
     702         781 :   s = Fp_sqrt(D,p);
     703         781 :   return Fp_halve(Fp_sub(s, b, p), p);
     704             : }
     705             : 
     706             : /*
     707             :  * pol is the modular polynomial of level 2 modulo p.
     708             :  *
     709             :  * (T, p) defines the field FF_{p^2} in which j_prev and j live.
     710             :  */
     711             : static long
     712        4878 : Fq_path_extends_to_floor(GEN j_prev, GEN j, GEN T, GEN p, GEN Phi2, long max_len)
     713             : {
     714        4878 :   pari_sp ltop = avma;
     715        4878 :   long d, i, l = lg(j);
     716             : 
     717             :   /* A path made its way to the floor if (i) its length was cut off
     718             :    * before reaching max_path_len, or (ii) it reached max_path_len but
     719             :    * only has one neighbour. */
     720       32215 :   for (d = 1; d <= max_len; ++d)
     721             :   {
     722       81746 :     for (i = 1; i < l; i++)
     723             :     {
     724       54409 :       GEN Phi2_j = FqX_div_by_X_x(FqXY_evalx(Phi2, gel(j,i), T, p), gel(j_prev,i), T, p, NULL);
     725       54409 :       GEN j_next = FqX_quad_root(Phi2_j, T, p);
     726       54409 :       if (!j_next)
     727        3178 :         return  gc_long(ltop, 1);
     728       51231 :       gel(j_prev,i) = gel(j, i); gel(j,i) = j_next;
     729             :     }
     730       27337 :     if (gc_needed(ltop, 2))
     731           0 :       gerepileall(ltop, 2, &j, &j_prev);
     732             :   }
     733        1700 :   return gc_long(ltop, 0);
     734             : }
     735             : 
     736             : static long
     737         448 : Fp_path_extends_to_floor(GEN j_prev, GEN j, GEN p, GEN Phi2, long max_len, GEN *pt_j, GEN *pt_j_prev)
     738             : {
     739         448 :   pari_sp ltop = avma;
     740         448 :   long d, i, l = lg(j);
     741             : 
     742             :   /* A path made its way to the floor if (i) its length was cut off
     743             :    * before reaching max_path_len, or (ii) it reached max_path_len but
     744             :    * only has one neighbour. */
     745         603 :   for (d = 1; d <= max_len; ++d)
     746             :   {
     747        1384 :     for (i = 1; i < l; i++)
     748             :     {
     749        1229 :       GEN Phi2_j = FpX_div_by_X_x(FpXY_evalx(Phi2, gel(j,i), p), gel(j_prev,i), p, NULL);
     750        1229 :       GEN j_next = FpX_quad_root(Phi2_j, p);
     751        1229 :       if (!j_next)
     752             :       {
     753         448 :         *pt_j = gel(j,i);
     754         448 :         *pt_j_prev = gel(j_prev,i);
     755         448 :         return 1;
     756             :       }
     757         781 :       gel(j_prev,i) = gel(j, i); gel(j,i) = j_next;
     758             :     }
     759         155 :     if (gc_needed(ltop, 2))
     760           0 :       gerepileall(ltop, 2, &j, &j_prev);
     761             :   }
     762           0 :   return gc_long(ltop, 0);
     763             : }
     764             : 
     765             : 
     766             : static int
     767        2737 : Fp_jissupersingular(GEN j, GEN p)
     768             : {
     769        2737 :   long max_path_len = expi(p)+1;
     770        2737 :   GEN Phi2 = FpXX_red(polmodular_ZXX(2,0,0,1), p);
     771        2737 :   GEN Phi2_j = FpXY_evalx(Phi2, j, p);
     772        2737 :   GEN roots = FpX_roots(Phi2_j, p);
     773        2737 :   long nbroots = lg(roots)-1;
     774        2737 :   GEN S, j_prev = NULL;
     775             : 
     776             :   /* Every node in a supersingular L-volcano has L + 1 neighbours. */
     777             :   /* Note: a multiple root only occur when j has CM by sqrt(-15). */
     778        2737 :   if (nbroots==0)
     779         665 :     return 0;
     780        2072 :   S = deg2pol_shallow(gen_1, gen_0, Fp_neg(Fp_2gener(p),p),1);
     781        2072 :   if (nbroots==1 && FpX_is_squarefree(Phi2_j, p))
     782        1624 :   { j_prev = j; j = FqX_quad_root(FpX_div_by_X_x(Phi2_j, gel(roots,1), p, NULL), S, p); }
     783             :   else
     784         448 :     if (!Fp_path_extends_to_floor(const_vec(nbroots,j), roots, p, Phi2, max_path_len, &j, &j_prev))
     785           0 :       return 1;
     786        2072 :   return !Fq_path_extends_to_floor(mkvec(j_prev), mkvec(j), S, p, Phi2, max_path_len);
     787             : }
     788             : 
     789             : static int
     790       14055 : jissupersingular(GEN j, GEN S, GEN p)
     791             : {
     792       14055 :   long max_path_len = expi(p)+1;
     793       14055 :   GEN Phi2 = FpXX_red(polmodular_ZXX(2,0,0,1), p);
     794       14055 :   GEN Phi2_j = FqXY_evalx(Phi2, j, S, p);
     795       14055 :   GEN roots = FpXQX_roots(Phi2_j, S, p);
     796       14055 :   long nbroots = lg(roots)-1;
     797             : 
     798             :   /* Every node in a supersingular L-volcano has L + 1 neighbours. */
     799             :   /* Note: a multiple root only occur when j has CM by sqrt(-15). */
     800       14055 :   if (nbroots==0 || (nbroots==1 && FqX_is_squarefree(Phi2_j, S, p)))
     801       11249 :     return 0;
     802             :   else
     803        2806 :     return !Fq_path_extends_to_floor(const_vec(nbroots,j), roots, S, p, Phi2, max_path_len);
     804             : }
     805             : 
     806             : int
     807        3711 : Fp_elljissupersingular(GEN j, GEN p)
     808             : {
     809             :   long CM;
     810        3711 :   if (abscmpiu(p, 5) <= 0) return signe(j) == 0; /* valid if p <= 5 */
     811        3571 :   CM = Fp_ellj_get_CM(j, gen_1, p);
     812        3571 :   if (CM < 0) return krosi(CM, p) < 0; /* valid if p > 3 */
     813             :   else
     814        2737 :     return Fp_jissupersingular(j, p);
     815             : }
     816             : 
     817             : /***********************************************************************/
     818             : /**                                                                   **/
     819             : /**                            Cardinal                               **/
     820             : /**                                                                   **/
     821             : /***********************************************************************/
     822             : 
     823             : /*assume a4,a6 reduced mod p odd */
     824             : static ulong
     825      466862 : Fl_elltrace_naive(ulong a4, ulong a6, ulong p)
     826             : {
     827             :   ulong i, j;
     828      466862 :   long a = 0;
     829             :   long d0, d1, d2, d3;
     830      466862 :   GEN k = const_vecsmall(p, -1);
     831      466886 :   k[1] = 0;
     832   106314335 :   for (i=1, j=1; i < p; i += 2, j = Fl_add(j, i, p))
     833   105847488 :     k[j+1] = 1;
     834      466847 :   d0 = 6%p; d1 = d0; d2 = Fl_add(a4, 1, p); d3 = a6;
     835      466845 :   for(i=0;; i++)
     836             :   {
     837   205274778 :     a -= k[1+d3];
     838   205274778 :     if (i==p-1) break;
     839   204808009 :     d3 = Fl_add(d3, d2, p);
     840   204812260 :     d2 = Fl_add(d2, d1, p);
     841   204810185 :     d1 = Fl_add(d1, d0, p);
     842             :   }
     843      466769 :   return a;
     844             : }
     845             : 
     846             : /* z1 <-- z1 + z2, with precomputed inverse */
     847             : static void
     848      305694 : FpE_add_ip(GEN z1, GEN z2, GEN a4, GEN p, GEN p2inv)
     849             : {
     850             :   GEN p1,x,x1,x2,y,y1,y2;
     851             : 
     852      305694 :   x1 = gel(z1,1); y1 = gel(z1,2);
     853      305694 :   x2 = gel(z2,1); y2 = gel(z2,2);
     854      305694 :   if (x1 == x2)
     855          67 :     p1 = Fp_add(a4, mulii(x1,mului(3,x1)), p);
     856             :   else
     857      305627 :     p1 = Fp_sub(y2,y1, p);
     858             : 
     859      305694 :   p1 = Fp_mul(p1, p2inv, p);
     860      305694 :   x = Fp_sub(sqri(p1), addii(x1,x2), p);
     861      305694 :   y = Fp_sub(mulii(p1,subii(x1,x)), y1, p);
     862      305694 :   affii(x, x1);
     863      305694 :   affii(y, y1);
     864      305694 : }
     865             : 
     866             : /* make sure *x has lgefint >= k */
     867             : static void
     868       19038 : _fix(GEN x, long k)
     869             : {
     870       19038 :   GEN y = (GEN)*x;
     871       19038 :   if (lgefint(y) < k) { GEN p1 = cgeti(k); affii(y,p1); *x = (long)p1; }
     872       19038 : }
     873             : 
     874             : /* Return the lift of a (mod b), which is closest to c */
     875             : static GEN
     876      254894 : closest_lift(GEN a, GEN b, GEN c)
     877             : {
     878      254894 :   return addii(a, mulii(b, diviiround(subii(c,a), b)));
     879             : }
     880             : 
     881             : static long
     882          78 : get_table_size(GEN pordmin, GEN B)
     883             : {
     884          78 :   pari_sp av = avma;
     885          78 :   GEN t = ceilr( sqrtr( divri(itor(pordmin, DEFAULTPREC), B) ) );
     886          78 :   if (is_bigint(t))
     887           0 :     pari_err_OVERFLOW("ellap [large prime: install the 'seadata' package]");
     888          78 :   set_avma(av);
     889          78 :   return itos(t) >> 1;
     890             : }
     891             : 
     892             : /* Find x such that kronecker(u = x^3+c4x+c6, p) is KRO.
     893             :  * Return point [x*u,u^2] on E (KRO=1) / E^twist (KRO=-1) */
     894             : static GEN
     895           0 : Fp_ellpoint(long KRO, ulong *px, GEN c4, GEN c6, GEN p)
     896             : {
     897           0 :   ulong x = *px;
     898             :   GEN u;
     899             :   for(;;)
     900             :   {
     901           0 :     x++; /* u = x^3 + c4 x + c6 */
     902           0 :     u = modii(addii(c6, mului(x, addii(c4, sqru(x)))), p);
     903           0 :     if (kronecker(u,p) == KRO) break;
     904             :   }
     905           0 :   *px = x;
     906           0 :   return mkvec2(modii(mului(x,u),p), Fp_sqr(u,p));
     907             : }
     908             : static GEN
     909        7227 : Fl_ellpoint(long KRO, ulong *px, ulong c4, ulong c6, ulong p)
     910             : {
     911        7227 :   ulong t, u, x = *px;
     912             :   for(;;)
     913             :   {
     914       14220 :     if (++x >= p) pari_err_PRIME("ellap",utoi(p));
     915       14220 :     t = Fl_add(c4, Fl_sqr(x,p), p);
     916       14220 :     u = Fl_add(c6, Fl_mul(x, t, p), p);
     917       14220 :     if (krouu(u,p) == KRO) break;
     918             :   }
     919        7227 :   *px = x;
     920        7227 :   return mkvecsmall2(Fl_mul(x,u,p), Fl_sqr(u,p));
     921             : }
     922             : 
     923             : static GEN ap_j1728(GEN a4,GEN p);
     924             : /* compute a_p using Shanks/Mestre + Montgomery's trick. Assume p > 457 */
     925             : static GEN
     926          78 : Fp_ellcard_Shanks(GEN c4, GEN c6, GEN p)
     927             : {
     928             :   pari_timer T;
     929             :   long *tx, *ty, *ti, pfinal, i, j, s, KRO, nb;
     930             :   ulong x;
     931          78 :   pari_sp av = avma, av2;
     932             :   GEN p1, P, mfh, h, F,f, fh,fg, pordmin, u, v, p1p, p2p, A, B, a4, pts;
     933          78 :   tx = NULL;
     934          78 :   ty = ti = NULL; /* gcc -Wall */
     935             : 
     936          78 :   if (!signe(c6)) {
     937           0 :     GEN ap = ap_j1728(c4, p);
     938           0 :     return gerepileuptoint(av, subii(addiu(p,1), ap));
     939             :   }
     940             : 
     941          78 :   if (DEBUGLEVEL >= 6) timer_start(&T);
     942             :   /* once #E(Fp) is know mod B >= pordmin, it is completely determined */
     943          78 :   pordmin = addiu(sqrti(gmul2n(p,4)), 1); /* ceil( 4sqrt(p) ) */
     944          78 :   p1p = addiu(p, 1);
     945          78 :   p2p = shifti(p1p, 1);
     946          78 :   x = 0; KRO = 0;
     947             :   /* how many 2-torsion points ? */
     948          78 :   switch(FpX_nbroots(mkpoln(4, gen_1, gen_0, c4, c6), p))
     949             :   {
     950           9 :     case 3:  A = gen_0; B = utoipos(4); break;
     951          31 :     case 1:  A = gen_0; B = gen_2; break;
     952          38 :     default: A = gen_1; B = gen_2; break; /* 0 */
     953             :   }
     954             :   for(;;)
     955             :   {
     956          78 :     h = closest_lift(A, B, p1p);
     957          78 :     if (!KRO) /* first time, initialize */
     958             :     {
     959          78 :       KRO = kronecker(c6,p);
     960          78 :       f = mkvec2(gen_0, Fp_sqr(c6,p));
     961             :     }
     962             :     else
     963             :     {
     964           0 :       KRO = -KRO;
     965           0 :       f = Fp_ellpoint(KRO, &x, c4,c6,p);
     966             :     }
     967             :     /* [ux, u^2] is on E_u: y^2 = x^3 + c4 u^2 x + c6 u^3
     968             :      * E_u isomorphic to E (resp. E') iff KRO = 1 (resp. -1)
     969             :      * #E(F_p) = p+1 - a_p, #E'(F_p) = p+1 + a_p
     970             :      *
     971             :      * #E_u(Fp) = A (mod B),  h is close to #E_u(Fp) */
     972          78 :     a4 = modii(mulii(c4, gel(f,2)), p); /* c4 for E_u */
     973          78 :     fh = FpE_mul(f, h, a4, p);
     974          78 :     if (ell_is_inf(fh)) goto FOUND;
     975             : 
     976          78 :     s = get_table_size(pordmin, B);
     977             :     /* look for h s.t f^h = 0 */
     978          78 :     if (!tx)
     979             :     { /* first time: initialize */
     980          78 :       tx = newblock(3*(s+1));
     981          78 :       ty = tx + (s+1);
     982          78 :       ti = ty + (s+1);
     983             :     }
     984          78 :     F = FpE_mul(f,B,a4,p);
     985          78 :     *tx = evaltyp(t_VECSMALL) | evallg(s+1);
     986             : 
     987             :     /* F = B.f */
     988          78 :     P = gcopy(fh);
     989          78 :     if (s < 3)
     990             :     { /* we're nearly done: naive search */
     991           0 :       GEN q1 = P, mF = FpE_neg(F, p); /* -F */
     992           0 :       for (i=1;; i++)
     993             :       {
     994           0 :         P = FpE_add(P,F,a4,p); /* h.f + i.F */
     995           0 :         if (ell_is_inf(P)) { h = addii(h, mului(i,B)); goto FOUND; }
     996           0 :         q1 = FpE_add(q1,mF,a4,p); /* h.f - i.F */
     997           0 :         if (ell_is_inf(q1)) { h = subii(h, mului(i,B)); goto FOUND; }
     998             :       }
     999             :     }
    1000             :     /* Baby Step/Giant Step */
    1001          78 :     nb = minss(128, s >> 1); /* > 0. Will do nb pts at a time: faster inverse */
    1002          78 :     pts = cgetg(nb+1, t_VEC);
    1003          78 :     j = lgefint(p);
    1004        9597 :     for (i=1; i<=nb; i++)
    1005             :     { /* baby steps */
    1006        9519 :       gel(pts,i) = P; /* h.f + (i-1).F */
    1007        9519 :       _fix(P+1, j); tx[i] = mod2BIL(gel(P,1));
    1008        9519 :       _fix(P+2, j); ty[i] = mod2BIL(gel(P,2));
    1009        9519 :       P = FpE_add(P,F,a4,p); /* h.f + i.F */
    1010        9519 :       if (ell_is_inf(P)) { h = addii(h, mului(i,B)); goto FOUND; }
    1011             :     }
    1012          78 :     mfh = FpE_neg(fh, p);
    1013          78 :     fg = FpE_add(P,mfh,a4,p); /* h.f + nb.F - h.f = nb.F */
    1014          78 :     if (ell_is_inf(fg)) { h = mului(nb,B); goto FOUND; }
    1015          78 :     u = cgetg(nb+1, t_VEC);
    1016          78 :     av2 = avma; /* more baby steps, nb points at a time */
    1017        1356 :     while (i <= s)
    1018             :     {
    1019             :       long maxj;
    1020      164235 :       for (j=1; j<=nb; j++) /* adding nb.F (part 1) */
    1021             :       {
    1022      162957 :         P = gel(pts,j); /* h.f + (i-nb-1+j-1).F */
    1023      162957 :         gel(u,j) = subii(gel(fg,1), gel(P,1));
    1024      162957 :         if (!signe(gel(u,j))) /* sum = 0 or doubling */
    1025             :         {
    1026           1 :           long k = i+j-2;
    1027           1 :           if (equalii(gel(P,2),gel(fg,2))) k -= 2*nb; /* fg == P */
    1028           1 :           h = addii(h, mulsi(k,B)); goto FOUND;
    1029             :         }
    1030             :       }
    1031        1278 :       v = FpV_inv(u, p);
    1032        1278 :       maxj = (i-1 + nb <= s)? nb: s % nb;
    1033      160545 :       for (j=1; j<=maxj; j++,i++) /* adding nb.F (part 2) */
    1034             :       {
    1035      159267 :         P = gel(pts,j);
    1036      159267 :         FpE_add_ip(P,fg, a4,p, gel(v,j));
    1037      159267 :         tx[i] = mod2BIL(gel(P,1));
    1038      159267 :         ty[i] = mod2BIL(gel(P,2));
    1039             :       }
    1040        1278 :       set_avma(av2);
    1041             :     }
    1042          77 :     P = FpE_add(gel(pts,j-1),mfh,a4,p); /* = (s-1).F */
    1043          77 :     if (ell_is_inf(P)) { h = mului(s-1,B); goto FOUND; }
    1044          77 :     if (DEBUGLEVEL >= 6)
    1045           0 :       timer_printf(&T, "[Fp_ellcard_Shanks] baby steps, s = %ld",s);
    1046             : 
    1047             :     /* giant steps: fg = s.F */
    1048          77 :     fg = FpE_add(P,F,a4,p);
    1049          77 :     if (ell_is_inf(fg)) { h = mului(s,B); goto FOUND; }
    1050          77 :     pfinal = mod2BIL(p); av2 = avma;
    1051             :     /* Goal of the following: sort points by increasing x-coordinate hash.
    1052             :      * Done in a complicated way to avoid allocating a large temp vector */
    1053          77 :     p1 = vecsmall_indexsort(tx); /* = permutation sorting tx */
    1054      168784 :     for (i=1; i<=s; i++) ti[i] = tx[p1[i]];
    1055             :     /* ti = tx sorted */
    1056      168784 :     for (i=1; i<=s; i++) { tx[i] = ti[i]; ti[i] = ty[p1[i]]; }
    1057             :     /* tx is sorted. ti = ty sorted */
    1058      168784 :     for (i=1; i<=s; i++) { ty[i] = ti[i]; ti[i] = p1[i]; }
    1059             :     /* ty is sorted. ti = permutation sorting tx */
    1060          77 :     if (DEBUGLEVEL >= 6) timer_printf(&T, "[Fp_ellcard_Shanks] sorting");
    1061          77 :     set_avma(av2);
    1062             : 
    1063          77 :     gaffect(fg, gel(pts,1));
    1064        9440 :     for (j=2; j<=nb; j++) /* pts[j] = j.fg = (s*j).F */
    1065             :     {
    1066        9363 :       P = FpE_add(gel(pts,j-1),fg,a4,p);
    1067        9363 :       if (ell_is_inf(P)) { h = mulii(mulss(s,j), B); goto FOUND; }
    1068        9363 :       gaffect(P, gel(pts,j));
    1069             :     }
    1070             :     /* replace fg by nb.fg since we do nb points at a time */
    1071          77 :     set_avma(av2);
    1072          77 :     fg = gcopy(gel(pts,nb)); /* copy: we modify (temporarily) pts[nb] below */
    1073          77 :     av2 = avma;
    1074             : 
    1075          77 :     for (i=1,j=1; ; i++)
    1076      152075 :     {
    1077      152152 :       GEN ftest = gel(pts,j);
    1078      152152 :       long m, l = 1, r = s+1;
    1079             :       long k, k2, j2;
    1080             : 
    1081      152152 :       set_avma(av2);
    1082      152152 :       k = mod2BIL(gel(ftest,1));
    1083     1930966 :       while (l < r)
    1084             :       {
    1085     1778814 :         m = (l+r) >> 1;
    1086     1778814 :         if (tx[m] < k) l = m+1; else r = m;
    1087             :       }
    1088      152152 :       if (r <= s && tx[r] == k)
    1089             :       {
    1090         154 :         while (r && tx[r] == k) r--;
    1091          77 :         k2 = mod2BIL(gel(ftest,2));
    1092          77 :         for (r++; r <= s && tx[r] == k; r++)
    1093          77 :           if (ty[r] == k2 || ty[r] == pfinal - k2)
    1094             :           { /* [h+j2] f == +/- ftest (= [i.s] f)? */
    1095          77 :             j2 = ti[r] - 1;
    1096          77 :             if (DEBUGLEVEL >=6)
    1097           0 :               timer_printf(&T, "[Fp_ellcard_Shanks] giant steps, i = %ld",i);
    1098          77 :             P = FpE_add(FpE_mul(F,stoi(j2),a4,p),fh,a4,p);
    1099          77 :             if (equalii(gel(P,1), gel(ftest,1)))
    1100             :             {
    1101          77 :               if (equalii(gel(P,2), gel(ftest,2))) i = -i;
    1102          77 :               h = addii(h, mulii(addis(mulss(s,i), j2), B));
    1103          77 :               goto FOUND;
    1104             :             }
    1105             :           }
    1106             :       }
    1107      152075 :       if (++j > nb)
    1108             :       { /* compute next nb points */
    1109        1149 :         long save = 0; /* gcc -Wall */;
    1110      147576 :         for (j=1; j<=nb; j++)
    1111             :         {
    1112      146427 :           P = gel(pts,j);
    1113      146427 :           gel(u,j) = subii(gel(fg,1), gel(P,1));
    1114      146427 :           if (gel(u,j) == gen_0) /* occurs once: i = j = nb, P == fg */
    1115             :           {
    1116          67 :             gel(u,j) = shifti(gel(P,2),1);
    1117          67 :             save = fg[1]; fg[1] = P[1];
    1118             :           }
    1119             :         }
    1120        1149 :         v = FpV_inv(u, p);
    1121      147576 :         for (j=1; j<=nb; j++)
    1122      146427 :           FpE_add_ip(gel(pts,j),fg,a4,p, gel(v,j));
    1123        1149 :         if (i == nb) { fg[1] = save; }
    1124        1149 :         j = 1;
    1125             :       }
    1126             :     }
    1127          78 : FOUND: /* found a point of exponent h on E_u */
    1128          78 :     h = FpE_order(f, h, a4, p);
    1129             :     /* h | #E_u(Fp) = A (mod B) */
    1130          78 :     A = Z_chinese_all(A, gen_0, B, h, &B);
    1131          78 :     if (cmpii(B, pordmin) >= 0) break;
    1132             :     /* not done: update A mod B for the _next_ curve, isomorphic to
    1133             :      * the quadratic twist of this one */
    1134           0 :     A = remii(subii(p2p,A), B); /* #E(Fp)+#E'(Fp) = 2p+2 */
    1135             :   }
    1136          78 :   if (tx) killblock(tx);
    1137          78 :   h = closest_lift(A, B, p1p);
    1138          78 :   return gerepileuptoint(av, KRO==1? h: subii(p2p,h));
    1139             : }
    1140             : 
    1141             : typedef struct
    1142             : {
    1143             :   ulong x,y,i;
    1144             : } multiple;
    1145             : 
    1146             : static int
    1147    15373710 : compare_multiples(multiple *a, multiple *b) { return a->x > b->x? 1:a->x<b->x?-1:0; }
    1148             : 
    1149             : /* find x such that h := a + b x is closest to c and return h:
    1150             :  * x = round((c-a) / b) = floor( (2(c-a) + b) / 2b )
    1151             :  * Assume 0 <= a < b < c  and b + 2c < 2^BIL */
    1152             : static ulong
    1153      261963 : uclosest_lift(ulong a, ulong b, ulong c)
    1154             : {
    1155      261963 :   ulong x = (b + ((c-a) << 1)) / (b << 1);
    1156      261963 :   return a + b * x;
    1157             : }
    1158             : 
    1159             : static long
    1160      227183 : Fle_dbl_inplace(GEN P, ulong a4, ulong p)
    1161             : {
    1162             :   ulong x, y, slope;
    1163      227183 :   if (!P[2]) return 1;
    1164      227155 :   x = P[1]; y = P[2];
    1165      227155 :   slope = Fl_div(Fl_add(Fl_triple(Fl_sqr(x,p), p), a4, p),
    1166             :                  Fl_double(y, p), p);
    1167      227155 :   P[1] = Fl_sub(Fl_sqr(slope, p), Fl_double(x, p), p);
    1168      227153 :   P[2] = Fl_sub(Fl_mul(slope, Fl_sub(x, P[1], p), p), y, p);
    1169      227153 :   return 0;
    1170             : }
    1171             : 
    1172             : static long
    1173     5797029 : Fle_add_inplace(GEN P, GEN Q, ulong a4, ulong p)
    1174             : {
    1175             :   ulong Px, Py, Qx, Qy, slope;
    1176     5797029 :   if (ell_is_inf(Q)) return 0;
    1177     5796974 :   Px = P[1]; Py = P[2];
    1178     5796974 :   Qx = Q[1]; Qy = Q[2];
    1179     5796974 :   if (Px==Qx)
    1180      238644 :     return Py==Qy ? Fle_dbl_inplace(P, a4, p): 1;
    1181     5558330 :   slope = Fl_div(Fl_sub(Py, Qy, p), Fl_sub(Px, Qx, p), p);
    1182     5558740 :   P[1] = Fl_sub(Fl_sub(Fl_sqr(slope, p), Px, p), Qx, p);
    1183     5558361 :   P[2] = Fl_sub(Fl_mul(slope, Fl_sub(Px, P[1], p), p), Py, p);
    1184     5558392 :   return 0;
    1185             : }
    1186             : 
    1187             : /* assume 99 < p < 2^(BIL-1) - 2^((BIL+1)/2) and e has good reduction at p.
    1188             :  * Should use Barett reduction + multi-inverse. See Fp_ellcard_Shanks() */
    1189             : static long
    1190      254742 : Fl_ellcard_Shanks(ulong c4, ulong c6, ulong p)
    1191             : {
    1192             :   GEN f, fh, fg, ftest, F;
    1193             :   ulong i, l, r, s, h, x, cp4, p1p, p2p, pordmin,A,B;
    1194             :   long KRO;
    1195      254742 :   pari_sp av = avma;
    1196             :   multiple *table;
    1197             : 
    1198      254742 :   if (!c6) {
    1199          14 :     GEN ap = ap_j1728(utoi(c4), utoipos(p));
    1200          14 :     return gc_long(av, p+1 - itos(ap));
    1201             :   }
    1202             : 
    1203      254728 :   pordmin = (ulong)(1 + 4*sqrt((double)p));
    1204      254728 :   p1p = p+1;
    1205      254728 :   p2p = p1p << 1;
    1206      254728 :   x = 0; KRO = 0;
    1207      254728 :   switch(Flx_nbroots(mkvecsmall5(0L, c6,c4,0L,1L), p))
    1208             :   {
    1209       51715 :     case 3:  A = 0; B = 4; break;
    1210      124407 :     case 1:  A = 0; B = 2; break;
    1211       78619 :     default: A = 1; B = 2; break; /* 0 */
    1212             :   }
    1213             :   for(;;)
    1214             :   { /* see comments in Fp_ellcard_Shanks */
    1215      261968 :     h = uclosest_lift(A, B, p1p);
    1216      261963 :     if (!KRO) /* first time, initialize */
    1217             :     {
    1218      254737 :       KRO = krouu(c6,p); /* != 0 */
    1219      254744 :       f = mkvecsmall2(0, Fl_sqr(c6,p));
    1220             :     }
    1221             :     else
    1222             :     {
    1223        7226 :       KRO = -KRO;
    1224        7226 :       f = Fl_ellpoint(KRO, &x, c4,c6,p);
    1225             :     }
    1226      261963 :     cp4 = Fl_mul(c4, f[2], p);
    1227      261963 :     fh = Fle_mulu(f, h, cp4, p);
    1228      261956 :     if (ell_is_inf(fh)) goto FOUND;
    1229             : 
    1230      255751 :     s = (ulong) (sqrt(((double)pordmin)/B) / 2);
    1231      255751 :     if (!s) s = 1;
    1232      255751 :     table = (multiple *) stack_malloc((s+1) * sizeof(multiple));
    1233      255759 :     F = Fle_mulu(f, B, cp4, p);
    1234     3346839 :     for (i=0; i < s; i++)
    1235             :     {
    1236     3102562 :       table[i].x = fh[1];
    1237     3102562 :       table[i].y = fh[2];
    1238     3102562 :       table[i].i = i;
    1239     3102562 :       if (Fle_add_inplace(fh, F, cp4, p)) { h += B*(i+1); goto FOUND; }
    1240             :     }
    1241      244277 :     qsort(table,s,sizeof(multiple),(QSCOMP)compare_multiples);
    1242      244278 :     fg = Fle_mulu(F, s, cp4, p); ftest = zv_copy(fg);
    1243      244267 :     if (ell_is_inf(ftest)) {
    1244           0 :       if (!uisprime(p)) pari_err_PRIME("ellap",utoi(p));
    1245           0 :       pari_err_BUG("ellap (f^(i*s) = 1)");
    1246             :     }
    1247     2939123 :     for (i=1; ; i++)
    1248             :     {
    1249     2939123 :       l=0; r=s;
    1250    20642547 :       while (l<r)
    1251             :       {
    1252    17703424 :         ulong m = (l+r) >> 1;
    1253    17703424 :         if (table[m].x < uel(ftest,1)) l=m+1; else r=m;
    1254             :       }
    1255     2939123 :       if (r < s && table[r].x == uel(ftest,1)) break;
    1256     2694845 :       if (Fle_add_inplace(ftest, fg, cp4, p)) pari_err_PRIME("ellap",utoi(p));
    1257             :     }
    1258      244278 :     h += table[r].i * B;
    1259      244278 :     if (table[r].y == uel(ftest,2))
    1260      126873 :       h -= s * i * B;
    1261             :     else
    1262      117405 :       h += s * i * B;
    1263      261974 : FOUND:
    1264      261974 :     h = itou(Fle_order(f, utoipos(h), cp4, p));
    1265             :     /* h | #E_u(Fp) = A (mod B) */
    1266             :     {
    1267             :       GEN C;
    1268      261953 :       A = itou( Z_chinese_all(gen_0, utoi(A), utoipos(h), utoipos(B), &C) );
    1269      261963 :       if (abscmpiu(C, pordmin) >= 0) { /* uclosest_lift could overflow */
    1270      254737 :         h = itou( closest_lift(utoi(A), C, utoipos(p1p)) );
    1271      254739 :         break;
    1272             :       }
    1273        7227 :       B = itou(C);
    1274             :     }
    1275        7227 :     A = (p2p - A) % B; set_avma(av);
    1276             :   }
    1277      254739 :   return gc_long(av, KRO==1? h: p2p-h);
    1278             : }
    1279             : 
    1280             : /** ellap from CM (original code contributed by Mark Watkins) **/
    1281             : 
    1282             : static GEN
    1283       75975 : ap_j0(GEN a6,GEN p)
    1284             : {
    1285             :   GEN a, b, e, d;
    1286       75975 :   if (umodiu(p,3) != 1) return gen_0;
    1287       37711 :   (void)cornacchia2(utoipos(27),p, &a,&b);
    1288       37797 :   if (umodiu(a, 3) == 1) a = negi(a);
    1289       37795 :   d = mulis(a6,-108);
    1290       37704 :   e = diviuexact(shifti(p,-1), 3); /* (p-1) / 6 */
    1291       37678 :   return centermod(mulii(a, Fp_pow(d, e, p)), p);
    1292             : }
    1293             : static GEN
    1294     2065805 : ap_j1728(GEN a4,GEN p)
    1295             : {
    1296             :   GEN a, b, e;
    1297     2065805 :   if (mod4(p) != 1) return gen_0;
    1298     1031891 :   (void)cornacchia2(utoipos(4),p, &a,&b);
    1299     1031891 :   if (Mod4(a)==0) a = b;
    1300     1031891 :   if (Mod2(a)==1) a = shifti(a,1);
    1301     1031891 :   if (Mod8(a)==6) a = negi(a);
    1302     1031891 :   e = shifti(p,-2); /* (p-1) / 4 */
    1303     1031891 :   return centermod(mulii(a, Fp_pow(a4, e, p)), p);
    1304             : }
    1305             : static GEN
    1306         126 : ap_j8000(GEN a6, GEN p)
    1307             : {
    1308             :   GEN a, b;
    1309         126 :   long r = mod8(p), s = 1;
    1310         126 :   if (r != 1 && r != 3) return gen_0;
    1311          49 :   (void)cornacchia2(utoipos(8),p, &a,&b);
    1312          49 :   switch(Mod16(a)) {
    1313          14 :     case 2: case 6:   if (Mod4(b)) s = -s;
    1314          14 :       break;
    1315          35 :     case 10: case 14: if (!Mod4(b)) s = -s;
    1316          35 :       break;
    1317             :   }
    1318          49 :   if (kronecker(mulis(a6, 42), p) < 0) s = -s;
    1319          49 :   return s > 0? a: negi(a);
    1320             : }
    1321             : static GEN
    1322         140 : ap_j287496(GEN a6, GEN p)
    1323             : {
    1324             :   GEN a, b;
    1325         140 :   long s = 1;
    1326         140 :   if (mod4(p) != 1) return gen_0;
    1327          70 :   (void)cornacchia2(utoipos(4),p, &a,&b);
    1328          70 :   if (Mod4(a)==0) a = b;
    1329          70 :   if (Mod2(a)==1) a = shifti(a,1);
    1330          70 :   if (Mod8(a)==6) s = -s;
    1331          70 :   if (krosi(2,p) < 0) s = -s;
    1332          70 :   if (kronecker(mulis(a6, -14), p) < 0) s = -s;
    1333          70 :   return s > 0? a: negi(a);
    1334             : }
    1335             : static GEN
    1336        1344 : ap_cm(int CM, long A6B, GEN a6, GEN p)
    1337             : {
    1338             :   GEN a, b;
    1339        1344 :   long s = 1;
    1340        1344 :   if (krosi(CM,p) < 0) return gen_0;
    1341         644 :   (void)cornacchia2(utoipos(-CM),p, &a, &b);
    1342         644 :   if ((CM&3) == 0) CM >>= 2;
    1343         644 :   if ((krois(a, -CM) > 0) ^ (CM == -7)) s = -s;
    1344         644 :   if (kronecker(mulis(a6,A6B), p) < 0) s = -s;
    1345         644 :   return s > 0? a: negi(a);
    1346             : }
    1347             : static GEN
    1348      493633 : ec_ap_cm(int CM, GEN a4, GEN a6, GEN p)
    1349             : {
    1350      493633 :   switch(CM)
    1351             :   {
    1352       29113 :     case  -3: return ap_j0(a6, p);
    1353      462910 :     case  -4: return ap_j1728(a4, p);
    1354         126 :     case  -8: return ap_j8000(a6, p);
    1355         140 :     case -16: return ap_j287496(a6, p);
    1356         154 :     case  -7: return ap_cm(CM, -2, a6, p);
    1357         147 :     case -11: return ap_cm(CM, 21, a6, p);
    1358         168 :     case -12: return ap_cm(CM, 22, a6, p);
    1359         147 :     case -19: return ap_cm(CM, 1, a6, p);
    1360         154 :     case -27: return ap_cm(CM, 253, a6, p);
    1361         140 :     case -28: return ap_cm(-7, -114, a6, p); /* yes, -7 ! */
    1362         147 :     case -43: return ap_cm(CM, 21, a6, p);
    1363         147 :     case -67: return ap_cm(CM, 217, a6, p);
    1364         140 :     case -163:return ap_cm(CM, 185801, a6, p);
    1365           0 :     default: return NULL;
    1366             :   }
    1367             : }
    1368             : 
    1369             : static GEN
    1370       49141 : Fp_ellj_nodiv(GEN a4, GEN a6, GEN p)
    1371             : {
    1372       49141 :   GEN a43 = Fp_mulu(Fp_powu(a4, 3, p), 4, p);
    1373       49148 :   GEN a62 = Fp_mulu(Fp_sqr(a6, p), 27, p);
    1374       49144 :   return mkvec2(Fp_mulu(a43, 1728, p), Fp_add(a43, a62, p));
    1375             : }
    1376             : 
    1377             : GEN
    1378          98 : Fp_ellj(GEN a4, GEN a6, GEN p)
    1379             : {
    1380          98 :   pari_sp av = avma;
    1381             :   GEN z;
    1382          98 :   if (lgefint(p) == 3)
    1383             :   {
    1384           0 :     ulong pp = p[2];
    1385           0 :     return utoi(Fl_ellj(umodiu(a4,pp), umodiu(a6,pp), pp));
    1386             :   }
    1387          98 :   z = Fp_ellj_nodiv(a4, a6, p);
    1388          98 :   return gerepileuptoint(av,Fp_div(gel(z,1),gel(z,2),p));
    1389             : }
    1390             : 
    1391             : void
    1392        1106 : Fp_ellj_to_a4a6(GEN j, GEN p, GEN *pt_a4, GEN *pt_a6)
    1393             : {
    1394        1106 :   j = modii(j, p);
    1395        1106 :   if (signe(j) == 0)    { *pt_a4 = gen_0; *pt_a6 = gen_1; }
    1396         791 :   else if (equaliu(j,umodui(1728,p))) { *pt_a4 = gen_1; *pt_a6 = gen_0; }
    1397             :   else
    1398             :   {
    1399         616 :     GEN k = Fp_sub(utoi(1728), j, p);
    1400         616 :     GEN kj = Fp_mul(k, j, p);
    1401         616 :     GEN k2j = Fp_mul(kj, k, p);
    1402         616 :     *pt_a4 = Fp_mulu(kj, 3, p);
    1403         616 :     *pt_a6 = Fp_double(k2j, p);
    1404             :   }
    1405        1106 : }
    1406             : 
    1407             : static GEN /* Only compute a mod p, so assume p>=17 */
    1408     1698753 : Fp_ellcard_CM(GEN a4, GEN a6, GEN p)
    1409             : {
    1410     1698753 :   pari_sp av = avma;
    1411             :   GEN a;
    1412     1698753 :   if (!signe(a4)) a = ap_j0(a6,p);
    1413     1651907 :   else if (!signe(a6)) a = ap_j1728(a4,p);
    1414             :   else
    1415             :   {
    1416       49026 :     GEN j = Fp_ellj_nodiv(a4, a6, p);
    1417       49045 :     long CM = Fp_ellj_get_CM(gel(j,1), gel(j,2), p);
    1418       49031 :     if (!CM) return gc_NULL(av);
    1419        1610 :     a = ec_ap_cm(CM,a4,a6,p);
    1420             :   }
    1421     1651390 :   return gerepileuptoint(av, subii(addiu(p,1),a));
    1422             : }
    1423             : 
    1424             : GEN
    1425     1959349 : Fp_ellcard(GEN a4, GEN a6, GEN p)
    1426             : {
    1427     1959349 :   long lp = expi(p);
    1428     1959475 :   ulong pp = p[2];
    1429     1959475 :   if (lp < 11)
    1430      260796 :     return utoi(pp+1 - Fl_elltrace_naive(umodiu(a4,pp), umodiu(a6,pp), pp));
    1431     1698679 :   { GEN a = Fp_ellcard_CM(a4,a6,p); if (a) return a; }
    1432       47425 :   if (lp >= 56)
    1433         868 :     return Fp_ellcard_SEA(a4, a6, p, 0);
    1434       46557 :   if (lp <= BITS_IN_LONG-2)
    1435       46479 :     return utoi(Fl_ellcard_Shanks(umodiu(a4,pp), umodiu(a6,pp), pp));
    1436          78 :   return Fp_ellcard_Shanks(a4, a6, p);
    1437             : }
    1438             : 
    1439             : long
    1440      381652 : Fl_elltrace(ulong a4, ulong a6, ulong p)
    1441             : {
    1442             :   pari_sp av;
    1443             :   long lp;
    1444             :   GEN a;
    1445      381652 :   if (p < (1<<11)) return Fl_elltrace_naive(a4, a6, p);
    1446      208254 :   lp = expu(p);
    1447      208254 :   if (lp <= minss(56, BITS_IN_LONG-2)) return p+1-Fl_ellcard_Shanks(a4, a6, p);
    1448           0 :   av = avma; a = subui(p+1, Fp_ellcard(utoi(a4), utoi(a6), utoipos(p)));
    1449           0 :   return gc_long(av, itos(a));
    1450             : }
    1451             : long
    1452      905566 : Fl_elltrace_CM(long CM, ulong a4, ulong a6, ulong p)
    1453             : {
    1454             :   pari_sp av;
    1455             :   GEN a;
    1456      905566 :   if (!CM) return Fl_elltrace(a4,a6,p);
    1457      524678 :   if (p < (1<<11)) return Fl_elltrace_naive(a4, a6, p);
    1458      492023 :   av = avma; a = ec_ap_cm(CM, utoi(a4), utoi(a6), utoipos(p));
    1459      492023 :   return gc_long(av, itos(a));
    1460             : }
    1461             : 
    1462             : static GEN
    1463       75766 : _FpE_pairorder(void *E, GEN P, GEN Q, GEN m, GEN F)
    1464             : {
    1465       75766 :   struct _FpE *e = (struct _FpE *) E;
    1466       75766 :   return  Fp_order(FpE_weilpairing(P,Q,m,e->a4,e->p), F, e->p);
    1467             : }
    1468             : 
    1469             : GEN
    1470      120715 : Fp_ellgroup(GEN a4, GEN a6, GEN N, GEN p, GEN *pt_m)
    1471             : {
    1472             :   struct _FpE e;
    1473      120715 :   e.a4=a4; e.a6=a6; e.p=p;
    1474      120715 :   return gen_ellgroup(N, subiu(p,1), pt_m, (void*)&e, &FpE_group, _FpE_pairorder);
    1475             : }
    1476             : 
    1477             : GEN
    1478         574 : Fp_ellgens(GEN a4, GEN a6, GEN ch, GEN D, GEN m, GEN p)
    1479             : {
    1480             :   GEN P;
    1481         574 :   pari_sp av = avma;
    1482             :   struct _FpE e;
    1483         574 :   e.a4=a4; e.a6=a6; e.p=p;
    1484         574 :   switch(lg(D)-1)
    1485             :   {
    1486         476 :   case 1:
    1487         476 :     P = gen_gener(gel(D,1), (void*)&e, &FpE_group);
    1488         476 :     P = mkvec(FpE_changepoint(P, ch, p));
    1489         476 :     break;
    1490          98 :   default:
    1491          98 :     P = gen_ellgens(gel(D,1), gel(D,2), m, (void*)&e, &FpE_group, _FpE_pairorder);
    1492          98 :     gel(P,1) = FpE_changepoint(gel(P,1), ch, p);
    1493          98 :     gel(P,2) = FpE_changepoint(gel(P,2), ch, p);
    1494          98 :     break;
    1495             :   }
    1496         574 :   return gerepilecopy(av, P);
    1497             : }
    1498             : 
    1499             : /* Not so fast arithmetic with points over elliptic curves over FpXQ */
    1500             : 
    1501             : /***********************************************************************/
    1502             : /**                                                                   **/
    1503             : /**                              FpXQE                                  **/
    1504             : /**                                                                   **/
    1505             : /***********************************************************************/
    1506             : 
    1507             : /* Theses functions deal with point over elliptic curves over FpXQ defined
    1508             :  * by an equation of the form y^2=x^3+a4*x+a6.
    1509             :  * Most of the time a6 is omitted since it can be recovered from any point
    1510             :  * on the curve.
    1511             :  */
    1512             : 
    1513             : GEN
    1514         976 : RgE_to_FpXQE(GEN x, GEN T, GEN p)
    1515             : {
    1516         976 :   if (ell_is_inf(x)) return x;
    1517         976 :   retmkvec2(Rg_to_FpXQ(gel(x,1),T,p),Rg_to_FpXQ(gel(x,2),T,p));
    1518             : }
    1519             : 
    1520             : GEN
    1521        1876 : FpXQE_changepoint(GEN x, GEN ch, GEN T, GEN p)
    1522             : {
    1523        1876 :   pari_sp av = avma;
    1524             :   GEN p1,z,u,r,s,t,v,v2,v3;
    1525        1876 :   if (ell_is_inf(x)) return x;
    1526         942 :   u = gel(ch,1); r = gel(ch,2);
    1527         942 :   s = gel(ch,3); t = gel(ch,4);
    1528         942 :   v = FpXQ_inv(u, T, p); v2 = FpXQ_sqr(v, T, p); v3 = FpXQ_mul(v,v2, T, p);
    1529         942 :   p1 = FpX_sub(gel(x,1),r, p);
    1530         942 :   z = cgetg(3,t_VEC);
    1531         942 :   gel(z,1) = FpXQ_mul(v2, p1, T, p);
    1532         942 :   gel(z,2) = FpXQ_mul(v3, FpX_sub(gel(x,2), FpX_add(FpXQ_mul(s,p1, T, p),t, p), p), T, p);
    1533         942 :   return gerepileupto(av, z);
    1534             : }
    1535             : 
    1536             : GEN
    1537         976 : FpXQE_changepointinv(GEN x, GEN ch, GEN T, GEN p)
    1538             : {
    1539             :   GEN u, r, s, t, X, Y, u2, u3, u2X, z;
    1540         976 :   if (ell_is_inf(x)) return x;
    1541         976 :   X = gel(x,1); Y = gel(x,2);
    1542         976 :   u = gel(ch,1); r = gel(ch,2);
    1543         976 :   s = gel(ch,3); t = gel(ch,4);
    1544         976 :   u2 = FpXQ_sqr(u, T, p); u3 = FpXQ_mul(u,u2, T, p);
    1545         976 :   u2X = FpXQ_mul(u2,X, T, p);
    1546         976 :   z = cgetg(3, t_VEC);
    1547         976 :   gel(z,1) = FpX_add(u2X,r, p);
    1548         976 :   gel(z,2) = FpX_add(FpXQ_mul(u3,Y, T, p), FpX_add(FpXQ_mul(s,u2X, T, p), t, p), p);
    1549         976 :   return z;
    1550             : }
    1551             : 
    1552             : static GEN
    1553         840 : random_nonsquare_FpXQ(GEN T, GEN p)
    1554             : {
    1555         840 :   pari_sp av = avma;
    1556         840 :   long n = degpol(T), v = varn(T);
    1557             :   GEN a;
    1558         840 :   if (odd(n))
    1559             :   {
    1560         420 :     GEN z = cgetg(3, t_POL);
    1561         420 :     z[1] = evalsigne(1) | evalvarn(v);
    1562         420 :     gel(z,2) = random_nonsquare_Fp(p); return z;
    1563             :   }
    1564             :   do
    1565             :   {
    1566         777 :     set_avma(av);
    1567         777 :     a = random_FpX(n, v, p);
    1568         777 :   } while (FpXQ_issquare(a, T, p));
    1569         420 :   return a;
    1570             : }
    1571             : 
    1572             : void
    1573         840 : FpXQ_elltwist(GEN a4, GEN a6, GEN T, GEN p, GEN *pt_a4, GEN *pt_a6)
    1574             : {
    1575         840 :   GEN d = random_nonsquare_FpXQ(T, p);
    1576         840 :   GEN d2 = FpXQ_sqr(d, T, p), d3 = FpXQ_mul(d2, d, T, p);
    1577         840 :   *pt_a4 = FpXQ_mul(a4, d2, T, p);
    1578         840 :   *pt_a6 = FpXQ_mul(a6, d3, T, p);
    1579         840 : }
    1580             : 
    1581             : static GEN
    1582      340967 : FpXQE_dbl_slope(GEN P, GEN a4, GEN T, GEN p, GEN *slope)
    1583             : {
    1584             :   GEN x, y, Q;
    1585      340967 :   if (ell_is_inf(P) || !signe(gel(P,2))) return ellinf();
    1586      339320 :   x = gel(P,1); y = gel(P,2);
    1587      339320 :   *slope = FpXQ_div(FpX_add(FpX_mulu(FpXQ_sqr(x, T, p), 3, p), a4, p),
    1588             :                             FpX_mulu(y, 2, p), T, p);
    1589      339320 :   Q = cgetg(3,t_VEC);
    1590      339320 :   gel(Q, 1) = FpX_sub(FpXQ_sqr(*slope, T, p), FpX_mulu(x, 2, p), p);
    1591      339320 :   gel(Q, 2) = FpX_sub(FpXQ_mul(*slope, FpX_sub(x, gel(Q, 1), p), T, p), y, p);
    1592      339320 :   return Q;
    1593             : }
    1594             : 
    1595             : GEN
    1596      327821 : FpXQE_dbl(GEN P, GEN a4, GEN T, GEN p)
    1597             : {
    1598      327821 :   pari_sp av = avma;
    1599             :   GEN slope;
    1600      327821 :   return gerepileupto(av, FpXQE_dbl_slope(P,a4,T,p,&slope));
    1601             : }
    1602             : 
    1603             : static GEN
    1604      268891 : FpXQE_add_slope(GEN P, GEN Q, GEN a4, GEN T, GEN p, GEN *slope)
    1605             : {
    1606             :   GEN Px, Py, Qx, Qy, R;
    1607      268891 :   if (ell_is_inf(P)) return Q;
    1608      268877 :   if (ell_is_inf(Q)) return P;
    1609      268877 :   Px = gel(P,1); Py = gel(P,2);
    1610      268877 :   Qx = gel(Q,1); Qy = gel(Q,2);
    1611      268877 :   if (ZX_equal(Px, Qx))
    1612             :   {
    1613        1228 :     if (ZX_equal(Py, Qy))
    1614           7 :       return FpXQE_dbl_slope(P, a4, T, p, slope);
    1615             :     else
    1616        1221 :       return ellinf();
    1617             :   }
    1618      267649 :   *slope = FpXQ_div(FpX_sub(Py, Qy, p), FpX_sub(Px, Qx, p), T, p);
    1619      267649 :   R = cgetg(3,t_VEC);
    1620      267649 :   gel(R, 1) = FpX_sub(FpX_sub(FpXQ_sqr(*slope, T, p), Px, p), Qx, p);
    1621      267649 :   gel(R, 2) = FpX_sub(FpXQ_mul(*slope, FpX_sub(Px, gel(R, 1), p), T, p), Py, p);
    1622      267649 :   return R;
    1623             : }
    1624             : 
    1625             : GEN
    1626      267106 : FpXQE_add(GEN P, GEN Q, GEN a4, GEN T, GEN p)
    1627             : {
    1628      267106 :   pari_sp av = avma;
    1629             :   GEN slope;
    1630      267106 :   return gerepileupto(av, FpXQE_add_slope(P,Q,a4,T,p,&slope));
    1631             : }
    1632             : 
    1633             : static GEN
    1634           0 : FpXQE_neg_i(GEN P, GEN p)
    1635             : {
    1636           0 :   if (ell_is_inf(P)) return P;
    1637           0 :   return mkvec2(gel(P,1), FpX_neg(gel(P,2), p));
    1638             : }
    1639             : 
    1640             : GEN
    1641       73329 : FpXQE_neg(GEN P, GEN T, GEN p)
    1642             : {
    1643             :   (void) T;
    1644       73329 :   if (ell_is_inf(P)) return ellinf();
    1645       73329 :   return mkvec2(gcopy(gel(P,1)), FpX_neg(gel(P,2), p));
    1646             : }
    1647             : 
    1648             : GEN
    1649           0 : FpXQE_sub(GEN P, GEN Q, GEN a4, GEN T, GEN p)
    1650             : {
    1651           0 :   pari_sp av = avma;
    1652             :   GEN slope;
    1653           0 :   return gerepileupto(av, FpXQE_add_slope(P, FpXQE_neg_i(Q, p), a4, T, p, &slope));
    1654             : }
    1655             : 
    1656             : struct _FpXQE { GEN a4,a6,T,p; };
    1657             : static GEN
    1658      327821 : _FpXQE_dbl(void *E, GEN P)
    1659             : {
    1660      327821 :   struct _FpXQE *ell = (struct _FpXQE *) E;
    1661      327821 :   return FpXQE_dbl(P, ell->a4, ell->T, ell->p);
    1662             : }
    1663             : static GEN
    1664      267106 : _FpXQE_add(void *E, GEN P, GEN Q)
    1665             : {
    1666      267106 :   struct _FpXQE *ell=(struct _FpXQE *) E;
    1667      267106 :   return FpXQE_add(P, Q, ell->a4, ell->T, ell->p);
    1668             : }
    1669             : static GEN
    1670       83136 : _FpXQE_mul(void *E, GEN P, GEN n)
    1671             : {
    1672       83136 :   pari_sp av = avma;
    1673       83136 :   struct _FpXQE *e=(struct _FpXQE *) E;
    1674       83136 :   long s = signe(n);
    1675       83136 :   if (!s || ell_is_inf(P)) return ellinf();
    1676       83136 :   if (s<0) P = FpXQE_neg(P, e->T, e->p);
    1677       83136 :   if (is_pm1(n)) return s>0? gcopy(P): P;
    1678        9575 :   return gerepilecopy(av, gen_pow_i(P, n, e, &_FpXQE_dbl, &_FpXQE_add));
    1679             : }
    1680             : 
    1681             : GEN
    1682         934 : FpXQE_mul(GEN P, GEN n, GEN a4, GEN T, GEN p)
    1683             : {
    1684             :   struct _FpXQE E;
    1685         934 :   E.a4= a4; E.T = T; E.p = p;
    1686         934 :   return _FpXQE_mul(&E, P, n);
    1687             : }
    1688             : 
    1689             : /* Finds a random nonsingular point on E */
    1690             : 
    1691             : GEN
    1692        1203 : random_FpXQE(GEN a4, GEN a6, GEN T, GEN p)
    1693             : {
    1694        1203 :   pari_sp ltop = avma;
    1695             :   GEN x, x2, y, rhs;
    1696        1203 :   long v = get_FpX_var(T), d = get_FpX_degree(T);
    1697             :   do
    1698             :   {
    1699        2376 :     set_avma(ltop);
    1700        2376 :     x   = random_FpX(d,v,p); /*  x^3+a4*x+a6 = x*(x^2+a4)+a6  */
    1701        2376 :     x2  = FpXQ_sqr(x, T, p);
    1702        2376 :     rhs = FpX_add(FpXQ_mul(x, FpX_add(x2, a4, p), T, p), a6, p);
    1703           0 :   } while ((!signe(rhs) && !signe(FpX_add(FpX_mulu(x2,3,p), a4, p)))
    1704        2376 :           || !FpXQ_issquare(rhs, T, p));
    1705        1203 :   y = FpXQ_sqrt(rhs, T, p);
    1706        1203 :   if (!y) pari_err_PRIME("random_FpE", p);
    1707        1203 :   return gerepilecopy(ltop, mkvec2(x, y));
    1708             : }
    1709             : 
    1710             : static GEN
    1711         269 : _FpXQE_rand(void *E)
    1712             : {
    1713         269 :   struct _FpXQE *e=(struct _FpXQE *) E;
    1714         269 :   return random_FpXQE(e->a4, e->a6, e->T, e->p);
    1715             : }
    1716             : 
    1717             : static const struct bb_group FpXQE_group={_FpXQE_add,_FpXQE_mul,_FpXQE_rand,hash_GEN,ZXV_equal,ell_is_inf};
    1718             : 
    1719             : const struct bb_group *
    1720          16 : get_FpXQE_group(void ** pt_E, GEN a4, GEN a6, GEN T, GEN p)
    1721             : {
    1722          16 :   struct _FpXQE *e = (struct _FpXQE *) stack_malloc(sizeof(struct _FpXQE));
    1723          16 :   e->a4 = a4; e->a6 = a6; e->T = T; e->p = p;
    1724          16 :   *pt_E = (void *) e;
    1725          16 :   return &FpXQE_group;
    1726             : }
    1727             : 
    1728             : GEN
    1729          14 : FpXQE_order(GEN z, GEN o, GEN a4, GEN T, GEN p)
    1730             : {
    1731          14 :   pari_sp av = avma;
    1732             :   struct _FpXQE e;
    1733          14 :   e.a4=a4; e.T=T; e.p=p;
    1734          14 :   return gerepileuptoint(av, gen_order(z, o, (void*)&e, &FpXQE_group));
    1735             : }
    1736             : 
    1737             : GEN
    1738           0 : FpXQE_log(GEN a, GEN b, GEN o, GEN a4, GEN T, GEN p)
    1739             : {
    1740           0 :   pari_sp av = avma;
    1741             :   struct _FpXQE e;
    1742           0 :   e.a4=a4; e.T=T; e.p=p;
    1743           0 :   return gerepileuptoint(av, gen_PH_log(a, b, o, (void*)&e, &FpXQE_group));
    1744             : }
    1745             : 
    1746             : /***********************************************************************/
    1747             : /**                                                                   **/
    1748             : /**                            Pairings                               **/
    1749             : /**                                                                   **/
    1750             : /***********************************************************************/
    1751             : 
    1752             : /* Derived from APIP from and by Jerome Milan, 2012 */
    1753             : 
    1754             : static GEN
    1755       15372 : FpXQE_vert(GEN P, GEN Q, GEN a4, GEN T, GEN p)
    1756             : {
    1757       15372 :   long vT = get_FpX_var(T);
    1758       15372 :   if (ell_is_inf(P))
    1759         245 :     return pol_1(get_FpX_var(T));
    1760       15127 :   if (!ZX_equal(gel(Q, 1), gel(P, 1)))
    1761       15127 :     return FpX_sub(gel(Q, 1), gel(P, 1), p);
    1762           0 :   if (signe(gel(P,2))!=0) return pol_1(vT);
    1763           0 :   return FpXQ_inv(FpX_add(FpX_mulu(FpXQ_sqr(gel(P,1), T, p), 3, p),
    1764             :                   a4, p), T, p);
    1765             : }
    1766             : 
    1767             : static GEN
    1768       14924 : FpXQE_Miller_line(GEN R, GEN Q, GEN slope, GEN a4, GEN T, GEN p)
    1769             : {
    1770       14924 :   long vT = get_FpX_var(T);
    1771       14924 :   GEN x = gel(Q, 1), y = gel(Q, 2);
    1772       14924 :   GEN tmp1  = FpX_sub(x, gel(R, 1), p);
    1773       14924 :   GEN tmp2  = FpX_add(FpXQ_mul(tmp1, slope, T, p), gel(R, 2), p);
    1774       14924 :   if (!ZX_equal(y, tmp2))
    1775       14924 :     return FpX_sub(y, tmp2, p);
    1776           0 :   if (signe(y) == 0)
    1777           0 :     return pol_1(vT);
    1778             :   else
    1779             :   {
    1780             :     GEN s1, s2;
    1781           0 :     GEN y2i = FpXQ_inv(FpX_mulu(y, 2, p), T, p);
    1782           0 :     s1 = FpXQ_mul(FpX_add(FpX_mulu(FpXQ_sqr(x, T, p), 3, p), a4, p), y2i, T, p);
    1783           0 :     if (!ZX_equal(s1, slope))
    1784           0 :       return FpX_sub(s1, slope, p);
    1785           0 :     s2 = FpXQ_mul(FpX_sub(FpX_mulu(x, 3, p), FpXQ_sqr(s1, T, p), p), y2i, T, p);
    1786           0 :     return signe(s2)!=0 ? s2: y2i;
    1787             :   }
    1788             : }
    1789             : 
    1790             : /* Computes the equation of the line tangent to R and returns its
    1791             :    evaluation at the point Q. Also doubles the point R.
    1792             :  */
    1793             : 
    1794             : static GEN
    1795       13314 : FpXQE_tangent_update(GEN R, GEN Q, GEN a4, GEN T, GEN p, GEN *pt_R)
    1796             : {
    1797       13314 :   if (ell_is_inf(R))
    1798             :   {
    1799          42 :     *pt_R = ellinf();
    1800          42 :     return pol_1(get_FpX_var(T));
    1801             :   }
    1802       13272 :   else if (!signe(gel(R,2)))
    1803             :   {
    1804         133 :     *pt_R = ellinf();
    1805         133 :     return FpXQE_vert(R, Q, a4, T, p);
    1806             :   } else {
    1807             :     GEN slope;
    1808       13139 :     *pt_R = FpXQE_dbl_slope(R, a4, T, p, &slope);
    1809       13139 :     return FpXQE_Miller_line(R, Q, slope, a4, T, p);
    1810             :   }
    1811             : }
    1812             : 
    1813             : /* Computes the equation of the line through R and P, and returns its
    1814             :    evaluation at the point Q. Also adds P to the point R.
    1815             :  */
    1816             : 
    1817             : static GEN
    1818        1855 : FpXQE_chord_update(GEN R, GEN P, GEN Q, GEN a4, GEN T, GEN p, GEN *pt_R)
    1819             : {
    1820        1855 :   if (ell_is_inf(R))
    1821             :   {
    1822           0 :     *pt_R = gcopy(P);
    1823           0 :     return FpXQE_vert(P, Q, a4, T, p);
    1824             :   }
    1825        1855 :   else if (ell_is_inf(P))
    1826             :   {
    1827           0 :     *pt_R = gcopy(R);
    1828           0 :     return FpXQE_vert(R, Q, a4, T, p);
    1829             :   }
    1830        1855 :   else if (ZX_equal(gel(P, 1), gel(R, 1)))
    1831             :   {
    1832          70 :     if (ZX_equal(gel(P, 2), gel(R, 2)))
    1833           0 :       return FpXQE_tangent_update(R, Q, a4, T, p, pt_R);
    1834             :     else
    1835             :     {
    1836          70 :       *pt_R = ellinf();
    1837          70 :       return FpXQE_vert(R, Q, a4, T, p);
    1838             :     }
    1839             :   } else {
    1840             :     GEN slope;
    1841        1785 :     *pt_R = FpXQE_add_slope(P, R, a4, T, p, &slope);
    1842        1785 :     return FpXQE_Miller_line(R, Q, slope, a4, T, p);
    1843             :   }
    1844             : }
    1845             : 
    1846             : struct _FpXQE_miller { GEN p, T, a4, P; };
    1847             : static GEN
    1848       13314 : FpXQE_Miller_dbl(void* E, GEN d)
    1849             : {
    1850       13314 :   struct _FpXQE_miller *m = (struct _FpXQE_miller *)E;
    1851       13314 :   GEN p  = m->p;
    1852       13314 :   GEN T = m->T, a4 = m->a4, P = m->P;
    1853             :   GEN v, line;
    1854       13314 :   GEN N = FpXQ_sqr(gel(d,1), T, p);
    1855       13314 :   GEN D = FpXQ_sqr(gel(d,2), T, p);
    1856       13314 :   GEN point = gel(d,3);
    1857       13314 :   line = FpXQE_tangent_update(point, P, a4, T, p, &point);
    1858       13314 :   N = FpXQ_mul(N, line, T, p);
    1859       13314 :   v = FpXQE_vert(point, P, a4, T, p);
    1860       13314 :   D = FpXQ_mul(D, v, T, p); return mkvec3(N, D, point);
    1861             : }
    1862             : 
    1863             : static GEN
    1864        1855 : FpXQE_Miller_add(void* E, GEN va, GEN vb)
    1865             : {
    1866        1855 :   struct _FpXQE_miller *m = (struct _FpXQE_miller *)E;
    1867        1855 :   GEN p = m->p;
    1868        1855 :   GEN T = m->T, a4 = m->a4, P = m->P;
    1869             :   GEN v, line, point;
    1870        1855 :   GEN na = gel(va,1), da = gel(va,2), pa = gel(va,3);
    1871        1855 :   GEN nb = gel(vb,1), db = gel(vb,2), pb = gel(vb,3);
    1872        1855 :   GEN N = FpXQ_mul(na, nb, T, p);
    1873        1855 :   GEN D = FpXQ_mul(da, db, T, p);
    1874        1855 :   line = FpXQE_chord_update(pa, pb, P, a4, T, p, &point);
    1875        1855 :   N = FpXQ_mul(N, line, T, p);
    1876        1855 :   v = FpXQE_vert(point, P, a4, T, p);
    1877        1855 :   D = FpXQ_mul(D, v, T, p); return mkvec3(N, D, point);
    1878             : }
    1879             : 
    1880             : /* Returns the Miller function f_{m, Q} evaluated at the point P using
    1881             :  * the standard Miller algorithm. */
    1882             : static GEN
    1883         203 : FpXQE_Miller(GEN Q, GEN P, GEN m, GEN a4, GEN T, GEN p)
    1884             : {
    1885         203 :   pari_sp av = avma;
    1886             :   struct _FpXQE_miller d;
    1887             :   GEN v, N, D, g1;
    1888             : 
    1889         203 :   d.a4 = a4; d.T = T; d.p = p; d.P = P;
    1890         203 :   g1 = pol_1(get_FpX_var(T));
    1891         203 :   v = gen_pow_i(mkvec3(g1,g1,Q), m, (void*)&d,
    1892             :                 FpXQE_Miller_dbl, FpXQE_Miller_add);
    1893         203 :   N = gel(v,1); D = gel(v,2);
    1894         203 :   return gerepileupto(av, FpXQ_div(N, D, T, p));
    1895             : }
    1896             : 
    1897             : GEN
    1898          98 : FpXQE_weilpairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, GEN p)
    1899             : {
    1900          98 :   pari_sp av = avma;
    1901             :   GEN N, D, w;
    1902          98 :   if (ell_is_inf(P) || ell_is_inf(Q) || ZXV_equal(P,Q))
    1903           0 :     return pol_1(get_FpX_var(T));
    1904          98 :   N = FpXQE_Miller(P, Q, m, a4, T, p);
    1905          98 :   D = FpXQE_Miller(Q, P, m, a4, T, p);
    1906          98 :   w = FpXQ_div(N, D, T, p);
    1907          98 :   if (mpodd(m)) w = FpX_neg(w, p);
    1908          98 :   return gerepileupto(av, w);
    1909             : }
    1910             : 
    1911             : GEN
    1912           7 : FpXQE_tatepairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, GEN p)
    1913             : {
    1914           7 :   if (ell_is_inf(P) || ell_is_inf(Q)) return pol_1(get_FpX_var(T));
    1915           7 :   return FpXQE_Miller(P, Q, m, a4, T, p);
    1916             : }
    1917             : 
    1918             : /***********************************************************************/
    1919             : /**                                                                   **/
    1920             : /**                           issupersingular                         **/
    1921             : /**                                                                   **/
    1922             : /***********************************************************************/
    1923             : 
    1924             : GEN
    1925        1718 : FpXQ_ellj(GEN a4, GEN a6, GEN T, GEN p)
    1926             : {
    1927        1718 :   if (absequaliu(p,3)) return pol_0(get_FpX_var(T));
    1928             :   else
    1929             :   {
    1930        1718 :     pari_sp av=avma;
    1931        1718 :     GEN a43 = FpXQ_mul(a4,FpXQ_sqr(a4,T,p),T,p);
    1932        1718 :     GEN a62 = FpXQ_sqr(a6,T,p);
    1933        1718 :     GEN num = FpX_mulu(a43,6912,p);
    1934        1718 :     GEN den = FpX_add(FpX_mulu(a43,4,p),FpX_mulu(a62,27,p),p);
    1935        1718 :     return gerepileuptoleaf(av, FpXQ_div(num, den, T, p));
    1936             :   }
    1937             : }
    1938             : 
    1939             : static GEN
    1940       33530 : FpXQ_is_quad(GEN x, GEN T, GEN p)
    1941             : {
    1942       33530 :   pari_sp av = avma;
    1943             :   GEN K;
    1944       33530 :   long d = degpol(T);
    1945       33530 :   x = FpXQ_red(x,T,p);
    1946       33530 :   if (lgpol(x)<=1) return NULL;
    1947       33530 :   if (d==2) return FpXQ_minpoly(x, T, p);
    1948       33530 :   if (odd(degpol(T))) return NULL;
    1949       33530 :   K = FpM_ker(FpXQ_matrix_pow(x, d, 3, T, p), p);
    1950       33530 :   if (lg(K)!=2) return gc_NULL(av);
    1951         588 :   return RgV_to_RgX(gel(K,1), get_FpX_var(T));
    1952             : }
    1953             : 
    1954             : int
    1955      165515 : FpXQ_elljissupersingular(GEN j, GEN T, GEN p)
    1956             : {
    1957      165515 :   pari_sp ltop = avma;
    1958             : 
    1959             :   /* All supersingular j-invariants are in FF_{p^2}, so we first check
    1960             :    * whether j is in FF_{p^2}.  If d is odd, then FF_{p^2} is not a
    1961             :    * subfield of FF_{p^d} so the j-invariants are all in FF_p.  Hence
    1962             :    * the j-invariants are in FF_{p^{2 - e}}. */
    1963      165515 :   ulong d = get_FpX_degree(T);
    1964             :   GEN S;
    1965      165515 :   if (degpol(j) <= 0) return Fp_elljissupersingular(constant_coeff(j), p);
    1966      164660 :   j = FpXQ_red(j, T, p);
    1967      164660 :   if (degpol(j) <= 0) return gc_bool(ltop, Fp_elljissupersingular(constant_coeff(j), p));
    1968             :   /* Now j is not in F_p */
    1969      164660 :   if (abscmpiu(p, 5) <= 0) return gc_bool(ltop,0); /* j != 0*/
    1970      164653 :   if (odd(d)) return 0;
    1971             :   /* Set S so that FF_p[T]/(S) is isomorphic to FF_{p^2}: */
    1972       46997 :   if (d == 2)
    1973       13467 :     S = T;
    1974             :   else /* d > 2 */
    1975             :   {
    1976       33530 :     S = FpXQ_is_quad(j, T, p);
    1977       33530 :     if (!S) return gc_bool(ltop,0);
    1978         588 :     j = pol_x(varn(S));
    1979             :   }
    1980       14055 :   return gc_bool(ltop, jissupersingular(j,S,p));
    1981             : }
    1982             : 
    1983             : int
    1984        1050 : Fq_elljissupersingular(GEN j, GEN T, GEN p)
    1985         959 : { return typ(j)==t_INT? Fp_elljissupersingular(j, p)
    1986        2009 :                       : FpXQ_elljissupersingular(j, T, p); }
    1987             : 
    1988             : /* p > 5 prime; return d such that (-d/p) = -1 */
    1989             : static ulong
    1990        1183 : find_inert_disc(GEN p)
    1991             : {
    1992        1183 :   long s = mod4(p) == 1? -1: 1; /* - (-1/p) */
    1993        1183 :   ulong d = 3;
    1994             :   while(1)
    1995             :   {
    1996        1190 :     if (kroui(d,p) == s) return d; /* = 3 mod (16) */
    1997         595 :     d++;
    1998         595 :     if (kroui(d>>2,p) == s) return d; /* = 4 mod (16) */
    1999         266 :     d += 3;
    2000         266 :     if (kroui(d,p) == s) return d; /* = 7 mod (16) */
    2001         105 :     d++;
    2002         105 :     if (kroui(d>>2,p) == s) return d; /* = 8 mod (16) */
    2003          35 :     d += 3;
    2004          35 :     if (kroui(d,p) == s) return d; /* = 11 mod (16) */
    2005           7 :     d += 4;
    2006           7 :     if (kroui(d,p) == s) return d; /* = 15 mod (16) */
    2007           7 :     d += 4;
    2008             :   }
    2009             : }
    2010             : 
    2011             : /* p > 5 */
    2012             : static GEN
    2013        1183 : ellsupersingularj_easy_FpXQ(GEN T, GEN p)
    2014             : {
    2015        1183 :   long d = find_inert_disc(p);
    2016        1183 :   GEN R = FpXQX_roots(polclass(stoi(-d), 0, 0), T, p);
    2017        1183 :   return gel(R,1);
    2018             : }
    2019             : 
    2020             : GEN
    2021        1204 : ellsupersingularj_FpXQ(GEN T, GEN p)
    2022             : {
    2023             :   GEN j, j2, R, Phi2;
    2024             :   long i, ep, lp;
    2025        1204 :   if (cmpiu(p, 5) <= 0) return pol_0(get_FpX_var(T));
    2026        1183 :   j2 = ellsupersingularj_easy_FpXQ(T, p);
    2027        1183 :   Phi2 = polmodular_ZXX(2,0,0,1);
    2028        1183 :   R = FpXQX_roots(FqXY_evalx(Phi2, j2, T, p), T, p);
    2029        1183 :   j = gel(R,1+random_Fl(lg(R)-1));
    2030        1183 :   ep = expi(p); lp = ep + random_Fl(ep);
    2031       18713 :   for (i = 1; i <= lp; i++)
    2032             :   {
    2033       17530 :     GEN Phi2_j = FqX_div_by_X_x(FqXY_evalx(Phi2, j, T, p), j2, T, p, NULL);
    2034       17530 :     R = FqX_quad_root(Phi2_j, T, p);
    2035       17530 :     if (!R) pari_err_PRIME("ellsupersingularj",p);
    2036       17530 :     j2 = j; j = random_bits(1) ? R: Fq_neg(Fq_add(gel(Phi2_j,3), R, T, p), T, p);
    2037             :   }
    2038        1183 :   return j;
    2039             : }
    2040             : 
    2041             : /***********************************************************************/
    2042             : /**                                                                   **/
    2043             : /**                           Point counting                          **/
    2044             : /**                                                                   **/
    2045             : /***********************************************************************/
    2046             : 
    2047             : GEN
    2048       15484 : elltrace_extension(GEN t, long n, GEN q)
    2049             : {
    2050       15484 :   pari_sp av = avma;
    2051       15484 :   GEN v = RgX_to_RgC(RgXQ_powu(pol_x(0), n, mkpoln(3,gen_1,negi(t),q)),2);
    2052       15484 :   GEN te = addii(shifti(gel(v,1),1), mulii(t,gel(v,2)));
    2053       15484 :   return gerepileuptoint(av, te);
    2054             : }
    2055             : 
    2056             : GEN
    2057       14721 : Fp_ffellcard(GEN a4, GEN a6, GEN q, long n, GEN p)
    2058             : {
    2059       14721 :   pari_sp av = avma;
    2060       14721 :   GEN ap = subii(addiu(p, 1), Fp_ellcard(a4, a6, p));
    2061       14721 :   GEN te = elltrace_extension(ap, n, p);
    2062       14721 :   return gerepileuptoint(av, subii(addiu(q, 1), te));
    2063             : }
    2064             : 
    2065             : static GEN
    2066        1687 : FpXQ_ellcardj(GEN a4, GEN a6, GEN j, GEN T, GEN q, GEN p, long n)
    2067             : {
    2068        1687 :   GEN q1 = addiu(q,1);
    2069        1687 :   if (signe(j)==0)
    2070             :   {
    2071             :     GEN W, w, t, N;
    2072         560 :     if (umodiu(q,6)!=1) return q1;
    2073         420 :     N = Fp_ffellcard(gen_0,gen_1,q,n,p);
    2074         420 :     t = subii(q1, N);
    2075         420 :     W = FpXQ_pow(a6,diviuexact(shifti(q,-1), 3),T,p);
    2076         420 :     if (degpol(W)>0) /*p=5 mod 6*/
    2077         126 :       return ZX_equal1(FpXQ_powu(W,3,T,p)) ? addii(q1,shifti(t,-1)):
    2078          42 :                                              subii(q1,shifti(t,-1));
    2079         336 :     w = modii(gel(W,2),p);
    2080         336 :     if (equali1(w))  return N;
    2081         266 :     if (equalii(w,subiu(p,1))) return addii(q1,t);
    2082             :     else /*p=1 mod 6*/
    2083             :     {
    2084         196 :       GEN u = shifti(t,-1), v = sqrtint(diviuexact(subii(q,sqri(u)),3));
    2085         196 :       GEN a = addii(u,v), b = shifti(v,1);
    2086         196 :       if (equali1(Fp_powu(w,3,p)))
    2087             :       {
    2088          98 :         if (dvdii(addmulii(a, w, b), p))
    2089          21 :           return subii(q1,subii(shifti(b,1),a));
    2090             :         else
    2091          77 :           return addii(q1,addii(a,b));
    2092             :       }
    2093             :       else
    2094             :       {
    2095          98 :         if (dvdii(submulii(a, w, b), p))
    2096          21 :           return subii(q1,subii(a,shifti(b,1)));
    2097             :         else
    2098          77 :           return subii(q1,addii(a,b));
    2099             :       }
    2100             :     }
    2101        1127 :   } else if (equalii(j,modsi(1728,p)))
    2102             :   {
    2103             :     GEN w, W, N, t;
    2104         567 :     if (mod4(q)==3) return q1;
    2105         427 :     W = FpXQ_pow(a4,shifti(q,-2),T,p);
    2106         427 :     if (degpol(W)>0) return q1; /*p=3 mod 4*/
    2107         315 :     w = modii(gel(W,2),p);
    2108         315 :     N = Fp_ffellcard(gen_1,gen_0,q,n,p);
    2109         315 :     if (equali1(w)) return N;
    2110         238 :     t = subii(q1, N);
    2111         238 :     if (equalii(w,subiu(p,1))) return addii(q1,t);
    2112             :     else /*p=1 mod 4*/
    2113             :     {
    2114         154 :       GEN u = shifti(t,-1), v = sqrtint(subii(q,sqri(u)));
    2115         154 :       if (dvdii(addmulii(u, w, v), p))
    2116          77 :         return subii(q1,shifti(v,1));
    2117             :       else
    2118          77 :         return addii(q1,shifti(v,1));
    2119             :     }
    2120             :   } else
    2121             :   {
    2122         560 :     GEN g = Fp_div(j, Fp_sub(utoi(1728), j, p), p);
    2123         560 :     GEN l = FpXQ_div(FpX_mulu(a6,3,p),FpX_mulu(a4,2,p),T,p);
    2124         560 :     GEN N = Fp_ffellcard(Fp_mulu(g,3,p),Fp_double(g,p),q,n,p);
    2125         560 :     if (FpXQ_issquare(l,T,p)) return N;
    2126         280 :     return subii(shifti(q1,1),N);
    2127             :   }
    2128             : }
    2129             : 
    2130             : static GEN
    2131           8 : FpXQ_ffellcard(GEN a4, GEN a6, GEN M, GEN q, GEN T, GEN p, long n)
    2132             : {
    2133           8 :   long m = degpol(M);
    2134           8 :   GEN j = pol_x(get_FpX_var(T));
    2135           8 :   GEN g = FpXQ_div(j, Fp_FpX_sub(utoi(1728), j, p), M, p);
    2136           8 :   GEN N = FpXQ_ellcard(FpX_mulu(g,3,p),FpX_mulu(g,2,p),M,p);
    2137           8 :   GEN qm = powiu(p, m), q1 = addiu(q, 1), qm1 = addiu(qm, 1);
    2138           8 :   GEN l = FpXQ_mul(FpX_mulu(a6,3,p),FpX_mulu(a4,2,p),T,p);
    2139           8 :   GEN te = elltrace_extension(subii(qm1, N), n/m, qm);
    2140           8 :   return FpXQ_issquare(l,T,p) ? subii(q1, te): addii(q1, te);
    2141             : }
    2142             : 
    2143             : static int
    2144           7 : FpXQ_is4power(GEN x, GEN T, GEN p)
    2145             : {
    2146           7 :   long d = get_FpX_degree(T);
    2147           7 :   if (lg(x) == 2 || absequalui(2, p)) return 1;
    2148           7 :   if (Mod4(p)==1)
    2149           7 :     return equali1(Fp_pow(FpXQ_norm(x,T,p),shifti(p,-2), p));
    2150           0 :   if (odd(d))
    2151           0 :     return FpXQ_issquare(x, T, p);
    2152           0 :   return ZX_equal1(FpXQ_pow(x, shifti(powiu(p, d),-2), T, p));
    2153             : }
    2154             : 
    2155             : /* http://www.numdam.org/article/ASENS_1969_4_2_4_521_0.pdf */
    2156             : 
    2157             : GEN
    2158           7 : FpXQ_ellcard_supersingular(GEN a4, GEN a6, GEN T, GEN p)
    2159             : {
    2160           7 :   pari_sp av = avma;
    2161           7 :   long d = get_FpX_degree(T);
    2162             :   GEN r;
    2163           7 :   if (equaliu(p,3))
    2164           0 :     r = Flxq_ellcard(ZX_to_Flx(a4,3), ZX_to_Flx(a6,3), ZXT_to_FlxT(T,3), 3);
    2165           7 :   else if (signe(a4)==0)
    2166           0 :     r = FpXQ_ellcardj(a4, a6, gen_0, T, powiu(p, d), p, d);
    2167           7 :   else if (signe(a6)==0)
    2168           0 :     r = FpXQ_ellcardj(a4, a6, modsi(1728,p), T, powiu(p, d), p, d);
    2169             :   else
    2170             :   {
    2171             :     GEN q, q2, t, D;
    2172           7 :     long qm4 = (odd(d>>1) && Mod4(p)==3);
    2173           7 :     if (odd(d)) return gen_0;
    2174           7 :     q2 = powiu(p, d>>1); q = sqri(q2);
    2175           7 :     t = shifti(q2, 1);
    2176           7 :     D = FpX_sub(FpX_Fp_mul(FpXQ_powu(a4,3,T,p), stoi(-4), p),
    2177             :                 FpX_mulu(FpXQ_sqr(a6,T,p), 27, p), p);
    2178          14 :     r = qm4 ^ FpXQ_is4power(D, T, p) ? subii(addiu(q, 1), t)
    2179           7 :                                      : addii(addiu(q, 1), t);
    2180             :   }
    2181           7 :   return gerepileuptoint(av, r);
    2182             : }
    2183             : 
    2184             : GEN
    2185          21 : Fq_ellcard_supersingular(GEN a4, GEN a6, GEN T, GEN p)
    2186          21 : { return T ? FpXQ_ellcard_supersingular(a4, a6, T, p) : addiu(p, 1); }
    2187             : 
    2188             : static GEN
    2189        8571 : FpXQ_ellcard_i(GEN a4, GEN a6, GEN T, GEN p)
    2190             : {
    2191        8571 :   long n = get_FpX_degree(T);
    2192        8571 :   GEN q = powiu(p, n);
    2193        8571 :   if (degpol(a4)<=0 && degpol(a6)<=0)
    2194         833 :     return Fp_ffellcard(constant_coeff(a4),constant_coeff(a6),q,n,p);
    2195        7738 :   if (lgefint(p)==3)
    2196             :   {
    2197        6020 :     ulong pp = p[2];
    2198        6020 :     return Flxq_ellcard(ZX_to_Flx(a4,pp),ZX_to_Flx(a6,pp),ZX_to_Flx(T,pp),pp);
    2199             :   }
    2200             :   else
    2201             :   {
    2202        1718 :     GEN J = FpXQ_ellj(a4,a6,T,p), M;
    2203        1718 :     if (degpol(J) <= 0)
    2204        1687 :       return FpXQ_ellcardj(a4,a6,constant_coeff(J),T,q,p,n);
    2205          31 :     M = FpXQ_minpoly(J,T,p);
    2206          31 :     if (degpol(M) < degpol(T))
    2207           8 :       return FpXQ_ffellcard(a4, a6, M, q, T, p, n);
    2208          23 :     return Fq_ellcard_SEA(a4, a6, q, T, p, 0);
    2209             :   }
    2210             : }
    2211             : 
    2212             : GEN
    2213        8571 : FpXQ_ellcard(GEN a4, GEN a6, GEN T, GEN p)
    2214             : {
    2215        8571 :   pari_sp av = avma;
    2216        8571 :   return gerepileuptoint(av, FpXQ_ellcard_i(a4, a6, T, p));
    2217             : }
    2218             : 
    2219             : static GEN
    2220          91 : _FpXQE_pairorder(void *E, GEN P, GEN Q, GEN m, GEN F)
    2221             : {
    2222          91 :   struct _FpXQE *e = (struct _FpXQE *) E;
    2223          91 :   return  FpXQ_order(FpXQE_weilpairing(P,Q,m,e->a4,e->T,e->p), F, e->T, e->p);
    2224             : }
    2225             : 
    2226             : GEN
    2227          15 : FpXQ_ellgroup(GEN a4, GEN a6, GEN N, GEN T, GEN p, GEN *pt_m)
    2228             : {
    2229             :   struct _FpXQE e;
    2230          15 :   GEN q = powiu(p, get_FpX_degree(T));
    2231          15 :   e.a4=a4; e.a6=a6; e.T=T; e.p=p;
    2232          15 :   return gen_ellgroup(N, subiu(q,1), pt_m, (void*)&e, &FpXQE_group, _FpXQE_pairorder);
    2233             : }
    2234             : 
    2235             : GEN
    2236           8 : FpXQ_ellgens(GEN a4, GEN a6, GEN ch, GEN D, GEN m, GEN T, GEN p)
    2237             : {
    2238             :   GEN P;
    2239           8 :   pari_sp av = avma;
    2240             :   struct _FpXQE e;
    2241           8 :   e.a4=a4; e.a6=a6; e.T=T; e.p=p;
    2242           8 :   switch(lg(D)-1)
    2243             :   {
    2244           8 :   case 1:
    2245           8 :     P = gen_gener(gel(D,1), (void*)&e, &FpXQE_group);
    2246           8 :     P = mkvec(FpXQE_changepoint(P, ch, T, p));
    2247           8 :     break;
    2248           0 :   default:
    2249           0 :     P = gen_ellgens(gel(D,1), gel(D,2), m, (void*)&e, &FpXQE_group, _FpXQE_pairorder);
    2250           0 :     gel(P,1) = FpXQE_changepoint(gel(P,1), ch, T, p);
    2251           0 :     gel(P,2) = FpXQE_changepoint(gel(P,2), ch, T, p);
    2252           0 :     break;
    2253             :   }
    2254           8 :   return gerepilecopy(av, P);
    2255             : }

Generated by: LCOV version 1.14