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.12.0 lcov report (development 23690-5d6e28857) Lines: 984 1067 92.2 %
Date: 2019-03-18 05:43:21 Functions: 106 114 93.0 %
Legend: Lines: hit not hit

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

Generated by: LCOV version 1.13