Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - FpE.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 16962-5a32637) Lines: 977 1108 88.2 %
Date: 2014-10-29 Functions: 111 130 85.4 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 460 589 78.1 %

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

Generated by: LCOV version 1.9