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

Generated by: LCOV version 1.13