Code coverage tests

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

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

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

LCOV - code coverage report
Current view: top level - basemath - FlxqE.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.14.0 lcov report (development 27775-aca467eab2) Lines: 797 819 97.3 %
Date: 2022-07-03 07:33:15 Functions: 79 80 98.8 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2012  The PARI group.
       2             : 
       3             : This file is part of the PARI/GP package.
       4             : 
       5             : PARI/GP is free software; you can redistribute it and/or modify it under the
       6             : terms of the GNU General Public License as published by the Free Software
       7             : Foundation; either version 2 of the License, or (at your option) any later
       8             : version. It is distributed in the hope that it will be useful, but WITHOUT
       9             : ANY WARRANTY WHATSOEVER.
      10             : 
      11             : Check the License for details. You should have received a copy of it, along
      12             : with the package; see the file 'COPYING'. If not, write to the Free Software
      13             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      14             : 
      15             : #include "pari.h"
      16             : #include "paripriv.h"
      17             : 
      18             : #define DEBUGLEVEL DEBUGLEVEL_ellcard
      19             : 
      20             : /* Not so fast arithmetic with points over elliptic curves over Fq,
      21             : small characteristic. */
      22             : 
      23             : /***********************************************************************/
      24             : /**                                                                   **/
      25             : /**                              FlxqE                                **/
      26             : /**                                                                   **/
      27             : /***********************************************************************/
      28             : 
      29             : /* Theses functions deal with point over elliptic curves over Fq defined
      30             :  * by an equation of the form y^2=x^3+a4*x+a6.
      31             :  * Most of the time a6 is omitted since it can be recovered from any point
      32             :  * on the curve.
      33             :  */
      34             : 
      35             : GEN
      36       63966 : RgE_to_FlxqE(GEN x, GEN T, ulong p)
      37             : {
      38       63966 :   if (ell_is_inf(x)) return x;
      39       63966 :   retmkvec2(Rg_to_Flxq(gel(x,1),T,p),Rg_to_Flxq(gel(x,2),T,p));
      40             : }
      41             : 
      42             : GEN
      43      154174 : FlxqE_changepoint(GEN x, GEN ch, GEN T, ulong p)
      44             : {
      45      154174 :   pari_sp av = avma;
      46             :   GEN p1,z,u,r,s,t,v,v2,v3;
      47      154174 :   if (ell_is_inf(x)) return x;
      48       91986 :   u = gel(ch,1); r = gel(ch,2);
      49       91986 :   s = gel(ch,3); t = gel(ch,4);
      50       91986 :   v = Flxq_inv(u, T, p); v2 = Flxq_sqr(v, T, p); v3 = Flxq_mul(v,v2, T, p);
      51       91986 :   p1 = Flx_sub(gel(x,1),r, p);
      52       91986 :   z = cgetg(3,t_VEC);
      53       91986 :   gel(z,1) = Flxq_mul(v2, p1, T, p);
      54       91986 :   gel(z,2) = Flxq_mul(v3, Flx_sub(gel(x,2), Flx_add(Flxq_mul(s, p1, T, p),t, p), p), T, p);
      55       91986 :   return gerepileupto(av, z);
      56             : }
      57             : 
      58             : GEN
      59       63966 : FlxqE_changepointinv(GEN x, GEN ch, GEN T, ulong p)
      60             : {
      61             :   GEN u, r, s, t, X, Y, u2, u3, u2X, z;
      62       63966 :   if (ell_is_inf(x)) return x;
      63       63966 :   X = gel(x,1); Y = gel(x,2);
      64       63966 :   u = gel(ch,1); r = gel(ch,2);
      65       63966 :   s = gel(ch,3); t = gel(ch,4);
      66       63966 :   u2 = Flxq_sqr(u, T, p); u3 = Flxq_mul(u,u2, T, p);
      67       63966 :   u2X = Flxq_mul(u2,X, T, p);
      68       63966 :   z = cgetg(3, t_VEC);
      69       63966 :   gel(z,1) = Flx_add(u2X,r, p);
      70       63966 :   gel(z,2) = Flx_add(Flxq_mul(u3,Y, T, p), Flx_add(Flxq_mul(s,u2X, T, p), t, p), p);
      71       63966 :   return z;
      72             : }
      73             : 
      74             : static GEN
      75       22834 : nonsquare_Flxq(GEN T, ulong p)
      76             : {
      77       22834 :   pari_sp av = avma;
      78       22834 :   long n = degpol(T), vs = T[1];
      79             :   GEN a;
      80       22834 :   if (odd(n))
      81        7686 :     return mkvecsmall2(vs, nonsquare_Fl(p));
      82             :   do
      83             :   {
      84       30821 :     set_avma(av);
      85       30821 :     a = random_Flx(n, vs, p);
      86       30821 :   } while (Flxq_issquare(a, T, p));
      87       15148 :   return a;
      88             : }
      89             : 
      90             : void
      91       22834 : Flxq_elltwist(GEN a, GEN a6, GEN T, ulong p, GEN *pt_a, GEN *pt_a6)
      92             : {
      93       22834 :   GEN d = nonsquare_Flxq(T, p);
      94       22834 :   GEN d2 = Flxq_sqr(d, T, p), d3 = Flxq_mul(d2, d, T, p);
      95       22834 :   if (typ(a)==t_VECSMALL)
      96             :   {
      97       15232 :     *pt_a  = Flxq_mul(a,  d2, T, p);
      98       15232 :     *pt_a6 = Flxq_mul(a6, d3, T, p);
      99             :   } else
     100             :   {
     101        7602 :     *pt_a  = mkvec(Flxq_mul(gel(a,1), d, T, p));
     102        7602 :     *pt_a6 = Flxq_mul(a6, d3, T, p);
     103             :   }
     104       22834 : }
     105             : 
     106             : static GEN
     107     1288706 : FlxqE_dbl_slope(GEN P, GEN a4, GEN T, ulong p, GEN *slope)
     108             : {
     109             :   GEN x, y, Q;
     110     1288706 :   if (ell_is_inf(P) || !lgpol(gel(P,2))) return ellinf();
     111     1186009 :   x = gel(P,1); y = gel(P,2);
     112     1186009 :   if (p==3UL)
     113      534415 :     *slope = typ(a4)==t_VEC ? Flxq_div(Flxq_mul(x, gel(a4, 1), T, p), y, T, p)
     114      534415 :                             : Flxq_div(a4, Flx_neg(y, p), T, p);
     115             :   else
     116             :   {
     117      651594 :     GEN sx = Flx_add(Flx_triple(Flxq_sqr(x, T, p), p), a4, p);
     118      651594 :     *slope = Flxq_div(sx, Flx_double(y, p), T, p);
     119             :   }
     120     1186009 :   Q = cgetg(3,t_VEC);
     121     1186009 :   gel(Q, 1) = Flx_sub(Flxq_sqr(*slope, T, p), Flx_double(x, p), p);
     122     1186009 :   if (typ(a4)==t_VEC) gel(Q, 1) = Flx_sub(gel(Q, 1), gel(a4, 1), p);
     123     1186009 :   gel(Q, 2) = Flx_sub(Flxq_mul(*slope, Flx_sub(x, gel(Q, 1), p), T, p), y, p);
     124     1186009 :   return Q;
     125             : }
     126             : 
     127             : GEN
     128     1259733 : FlxqE_dbl(GEN P, GEN a4, GEN T, ulong p)
     129             : {
     130     1259733 :   pari_sp av = avma;
     131             :   GEN slope;
     132     1259733 :   return gerepileupto(av, FlxqE_dbl_slope(P,a4, T, p,&slope));
     133             : }
     134             : 
     135             : static GEN
     136      538237 : FlxqE_add_slope(GEN P, GEN Q, GEN a4, GEN T, ulong p, GEN *slope)
     137             : {
     138             :   GEN Px, Py, Qx, Qy, R;
     139      538237 :   if (ell_is_inf(P)) return Q;
     140      534575 :   if (ell_is_inf(Q)) return P;
     141      534281 :   Px = gel(P,1); Py = gel(P,2);
     142      534281 :   Qx = gel(Q,1); Qy = gel(Q,2);
     143      534281 :   if (Flx_equal(Px, Qx))
     144             :   {
     145       47986 :     if (Flx_equal(Py, Qy))
     146        1413 :       return FlxqE_dbl_slope(P, a4, T, p, slope);
     147             :     else
     148       46573 :       return ellinf();
     149             :   }
     150      486295 :   *slope = Flxq_div(Flx_sub(Py, Qy, p), Flx_sub(Px, Qx, p), T, p);
     151      486295 :   R = cgetg(3,t_VEC);
     152      486295 :   gel(R, 1) = Flx_sub(Flx_sub(Flxq_sqr(*slope, T, p), Px, p), Qx, p);
     153      486295 :   if (typ(a4)==t_VEC) gel(R, 1) = Flx_sub(gel(R, 1),gel(a4, 1), p);
     154      486295 :   gel(R, 2) = Flx_sub(Flxq_mul(*slope, Flx_sub(Px, gel(R, 1), p), T, p), Py, p);
     155      486295 :   return R;
     156             : }
     157             : 
     158             : GEN
     159      533819 : FlxqE_add(GEN P, GEN Q, GEN a4, GEN T, ulong p)
     160             : {
     161      533819 :   pari_sp av = avma;
     162             :   GEN slope;
     163      533819 :   return gerepileupto(av, FlxqE_add_slope(P,Q,a4, T, p,&slope));
     164             : }
     165             : 
     166             : static GEN
     167        1391 : FlxqE_neg_i(GEN P, ulong p)
     168             : {
     169        1391 :   if (ell_is_inf(P)) return P;
     170        1391 :   return mkvec2(gel(P,1), Flx_neg(gel(P,2), p));
     171             : }
     172             : 
     173             : GEN
     174         490 : FlxqE_neg(GEN P, GEN T, ulong p)
     175             : {
     176             :   (void) T;
     177         490 :   if (ell_is_inf(P)) return ellinf();
     178         490 :   return mkvec2(gcopy(gel(P,1)), Flx_neg(gel(P,2), p));
     179             : }
     180             : 
     181             : GEN
     182        1391 : FlxqE_sub(GEN P, GEN Q, GEN a4, GEN T, ulong p)
     183             : {
     184        1391 :   pari_sp av = avma;
     185             :   GEN slope;
     186        1391 :   return gerepileupto(av, FlxqE_add_slope(P, FlxqE_neg_i(Q, p), a4, T, p, &slope));
     187             : }
     188             : 
     189             : struct _FlxqE
     190             : {
     191             :   GEN a4, a6;
     192             :   GEN T;
     193             :   ulong p;
     194             : };
     195             : 
     196             : static GEN
     197     1259733 : _FlxqE_dbl(void *E, GEN P)
     198             : {
     199     1259733 :   struct _FlxqE *ell = (struct _FlxqE *) E;
     200     1259733 :   return FlxqE_dbl(P, ell->a4, ell->T, ell->p);
     201             : }
     202             : 
     203             : static GEN
     204      524988 : _FlxqE_add(void *E, GEN P, GEN Q)
     205             : {
     206      524988 :   struct _FlxqE *ell=(struct _FlxqE *) E;
     207      524988 :   return FlxqE_add(P, Q, ell->a4, ell->T, ell->p);
     208             : }
     209             : 
     210             : static GEN
     211      228296 : _FlxqE_mul(void *E, GEN P, GEN n)
     212             : {
     213      228296 :   pari_sp av = avma;
     214      228296 :   struct _FlxqE *e=(struct _FlxqE *) E;
     215      228296 :   long s = signe(n);
     216      228296 :   if (!s || ell_is_inf(P)) return ellinf();
     217      228078 :   if (s<0) P = FlxqE_neg(P, e->T, e->p);
     218      228078 :   if (is_pm1(n)) return s>0? gcopy(P): P;
     219      221464 :   return gerepilecopy(av, gen_pow_i(P, n, e, &_FlxqE_dbl, &_FlxqE_add));
     220             : }
     221             : 
     222             : GEN
     223       65437 : FlxqE_mul(GEN P, GEN n, GEN a4, GEN T, ulong p)
     224             : {
     225             :   struct _FlxqE E;
     226       65437 :   E.a4= a4; E.T = T; E.p = p;
     227       65437 :   return _FlxqE_mul(&E, P, n);
     228             : }
     229             : 
     230             : /* 3*x^2+2*a2*x = -a2*x, and a2!=0 */
     231             : 
     232             : /* Finds a random nonsingular point on E */
     233             : static GEN
     234       77546 : random_F3xqE(GEN a2, GEN a6, GEN T)
     235             : {
     236       77546 :   pari_sp ltop = avma;
     237             :   GEN x, y, rhs;
     238       77546 :   const ulong p = 3;
     239             :   do
     240             :   {
     241      154728 :     set_avma(ltop);
     242      154728 :     x   = random_Flx(get_Flx_degree(T),get_Flx_var(T),p);
     243      154728 :     rhs = Flx_add(Flxq_mul(Flxq_sqr(x, T, p), Flx_add(x, a2, p), T, p), a6, p);
     244      154728 :   } while ((!lgpol(rhs) && !lgpol(x)) || !Flxq_issquare(rhs, T, p));
     245       77546 :   y = Flxq_sqrt(rhs, T, p);
     246       77546 :   if (!y) pari_err_PRIME("random_F3xqE", T);
     247       77546 :   return gerepilecopy(ltop, mkvec2(x, y));
     248             : }
     249             : 
     250             : /* Finds a random nonsingular point on E */
     251             : GEN
     252      145336 : random_FlxqE(GEN a4, GEN a6, GEN T, ulong p)
     253             : {
     254      145336 :   pari_sp ltop = avma;
     255             :   GEN x, x2, y, rhs;
     256      145336 :   if (typ(a4)==t_VEC)
     257       77546 :     return random_F3xqE(gel(a4,1), a6, T);
     258             :   do
     259             :   {
     260      132842 :     set_avma(ltop);
     261      132842 :     x   = random_Flx(get_Flx_degree(T),get_Flx_var(T),p);
     262      132842 :     x2  = Flxq_sqr(x, T, p); /*  x^3+a4*x+a6 = x*(x^2+a4)+a6  */
     263      132842 :     rhs = Flx_add(Flxq_mul(x, Flx_add(x2, a4, p), T, p), a6, p);
     264      134809 :   } while ((!lgpol(rhs) && !lgpol(Flx_add(Flx_triple(x2, p), a4, p)))
     265      134767 :           || !Flxq_issquare(rhs, T, p));
     266       67790 :   y = Flxq_sqrt(rhs, T, p);
     267       67790 :   if (!y) pari_err_PRIME("random_FlxqE", T);
     268       67790 :   return gerepilecopy(ltop, mkvec2(x, y));
     269             : }
     270             : 
     271             : static GEN
     272       70212 : _FlxqE_rand(void *E)
     273             : {
     274       70212 :   struct _FlxqE *ell=(struct _FlxqE *) E;
     275       70212 :   return random_FlxqE(ell->a4, ell->a6, ell->T, ell->p);
     276             : }
     277             : 
     278             : static const struct bb_group FlxqE_group={_FlxqE_add,_FlxqE_mul,_FlxqE_rand,hash_GEN,zvV_equal,ell_is_inf, NULL};
     279             : 
     280             : const struct bb_group *
     281          34 : get_FlxqE_group(void ** pt_E, GEN a4, GEN a6, GEN T, ulong p)
     282             : {
     283          34 :   struct _FlxqE *e = (struct _FlxqE *) stack_malloc(sizeof(struct _FlxqE));
     284          34 :   e->a4 = a4; e->a6 = a6; e->T = Flx_get_red(T, p); e->p = p;
     285          34 :   *pt_E = (void *) e;
     286          34 :   return &FlxqE_group;
     287             : }
     288             : 
     289             : GEN
     290        2729 : FlxqE_order(GEN z, GEN o, GEN a4, GEN T, ulong p)
     291             : {
     292        2729 :   pari_sp av = avma;
     293             :   struct _FlxqE e;
     294        2729 :   e.a4=a4; e.T=T; e.p=p;
     295        2729 :   return gerepileuptoint(av, gen_order(z, o, (void*)&e, &FlxqE_group));
     296             : }
     297             : 
     298             : GEN
     299          49 : FlxqE_log(GEN a, GEN b, GEN o, GEN a4, GEN T, ulong p)
     300             : {
     301          49 :   pari_sp av = avma;
     302             :   struct _FlxqE e;
     303          49 :   e.a4=a4; e.T=T; e.p=p;
     304          49 :   return gerepileuptoint(av, gen_PH_log(a, b, o, (void*)&e, &FlxqE_group));
     305             : }
     306             : 
     307             : /***********************************************************************/
     308             : /**                                                                   **/
     309             : /**                            Pairings                               **/
     310             : /**                                                                   **/
     311             : /***********************************************************************/
     312             : 
     313             : /* Derived from APIP from and by Jerome Milan, 2012 */
     314             : 
     315             : static GEN
     316       67918 : FlxqE_vert(GEN P, GEN Q, GEN a4, GEN T, ulong p)
     317             : {
     318       67918 :   long vT = get_Flx_var(T);
     319             :   GEN df;
     320       67918 :   if (ell_is_inf(P))
     321       20627 :     return pol1_Flx(vT);
     322       47291 :   if (!Flx_equal(gel(Q, 1), gel(P, 1)))
     323       45480 :     return Flx_sub(gel(Q, 1), gel(P, 1), p);
     324        1811 :   if (lgpol(gel(P,2))!=0) return pol1_Flx(vT);
     325         672 :   df = typ(a4)==t_VEC ? Flxq_mul(gel(P,1), Flx_mulu(gel(a4, 1), 2, p), T, p)
     326        1142 :                       : a4;
     327        1142 :   return Flxq_inv(Flx_add(Flx_mulu(Flxq_sqr(gel(P,1), T, p), 3, p),
     328             :                           df, p), T, p);
     329             : }
     330             : 
     331             : static GEN
     332       30587 : FlxqE_Miller_line(GEN R, GEN Q, GEN slope, GEN a4, GEN T, ulong p)
     333             : {
     334       30587 :   long vT = get_Flx_var(T);
     335       30587 :   GEN x = gel(Q, 1), y = gel(Q, 2);
     336       30587 :   GEN tmp1 = Flx_sub(x, gel(R, 1), p);
     337       30587 :   GEN tmp2 = Flx_add(Flxq_mul(tmp1, slope, T, p), gel(R, 2), p);
     338       30587 :   if (!Flx_equal(y, tmp2))
     339       29452 :     return Flx_sub(y, tmp2, p);
     340        1135 :   if (lgpol(y) == 0)
     341         571 :     return pol1_Flx(vT);
     342             :   else
     343             :   {
     344         564 :     GEN s1, s2, a2 = typ(a4)==t_VEC ? gel(a4,1): NULL;
     345         564 :     GEN y2i = Flxq_inv(Flx_mulu(y, 2, p), T, p);
     346         564 :     GEN df = a2 ? Flxq_mul(x, Flx_mulu(a2, 2, p), T, p): a4;
     347             :     GEN x3, ddf;
     348         564 :     s1 = Flxq_mul(Flx_add(Flx_mulu(Flxq_sqr(x, T, p), 3, p), df, p), y2i, T, p);
     349         564 :     if (!Flx_equal(s1, slope))
     350         330 :       return Flx_sub(s1, slope, p);
     351         234 :     x3 = Flx_mulu(x, 3, p);
     352         234 :     ddf = a2 ? Flx_add(x3, a2, p): x3;
     353         234 :     s2 = Flxq_mul(Flx_sub(ddf, Flxq_sqr(s1, T, p), p), y2i, T, p);
     354         234 :     return lgpol(s2)!=0 ? s2: y2i;
     355             :   }
     356             : }
     357             : 
     358             : /* Computes the equation of the line tangent to R and returns its
     359             :  * evaluation at the point Q. Also doubles the point R. */
     360             : static GEN
     361       46863 : FlxqE_tangent_update(GEN R, GEN Q, GEN a4, GEN T, ulong p, GEN *pt_R)
     362             : {
     363       46863 :   if (ell_is_inf(R))
     364             :   {
     365        4077 :     *pt_R = ellinf();
     366        4077 :     return pol1_Flx(get_Flx_var(T));
     367             :   }
     368       42786 :   else if (!lgpol(gel(R,2)))
     369             :   {
     370       15226 :     *pt_R = ellinf();
     371       15226 :     return FlxqE_vert(R, Q, a4, T, p);
     372             :   } else {
     373             :     GEN slope;
     374       27560 :     *pt_R = FlxqE_dbl_slope(R, a4, T, p, &slope);
     375       27560 :     return FlxqE_Miller_line(R, Q, slope, a4, T, p);
     376             :   }
     377             : }
     378             : 
     379             : /* Computes the equation of the line through R and P, and returns its
     380             :  * evaluation at the point Q. Also adds P to the point R. */
     381             : static GEN
     382        4435 : FlxqE_chord_update(GEN R, GEN P, GEN Q, GEN a4, GEN T, ulong p, GEN *pt_R)
     383             : {
     384        4435 :   if (ell_is_inf(R))
     385             :   {
     386          77 :     *pt_R = gcopy(P);
     387          77 :     return FlxqE_vert(P, Q, a4, T, p);
     388             :   }
     389        4358 :   else if (ell_is_inf(P))
     390             :   {
     391           0 :     *pt_R = gcopy(R);
     392           0 :     return FlxqE_vert(R, Q, a4, T, p);
     393             :   }
     394        4358 :   else if (Flx_equal(gel(P, 1), gel(R, 1)))
     395             :   {
     396        1331 :     if (Flx_equal(gel(P, 2), gel(R, 2)))
     397           7 :       return FlxqE_tangent_update(R, Q, a4, T, p, pt_R);
     398             :     else
     399             :     {
     400        1324 :       *pt_R = ellinf();
     401        1324 :       return FlxqE_vert(R, Q, a4, T, p);
     402             :     }
     403             :   } else {
     404             :     GEN slope;
     405        3027 :     *pt_R = FlxqE_add_slope(P, R, a4, T, p, &slope);
     406        3027 :     return FlxqE_Miller_line(R, Q, slope, a4, T, p);
     407             :   }
     408             : }
     409             : 
     410             : struct _FlxqE_miller
     411             : {
     412             :   ulong p;
     413             :   GEN T, a4, P;
     414             : };
     415             : 
     416             : static GEN
     417       46856 : FlxqE_Miller_dbl(void* E, GEN d)
     418             : {
     419       46856 :   struct _FlxqE_miller *m = (struct _FlxqE_miller *)E;
     420       46856 :   ulong p  = m->p;
     421       46856 :   GEN T = m->T, a4 = m->a4, P = m->P;
     422       46856 :   GEN v, line, point = gel(d,3);
     423       46856 :   GEN N = Flxq_sqr(gel(d,1), T, p);
     424       46856 :   GEN D = Flxq_sqr(gel(d,2), T, p);
     425       46856 :   line = FlxqE_tangent_update(point, P, a4, T, p, &point);
     426       46856 :   N  = Flxq_mul(N, line, T, p);
     427       46856 :   v = FlxqE_vert(point, P, a4, T, p);
     428       46856 :   D = Flxq_mul(D, v, T, p); return mkvec3(N, D, point);
     429             : }
     430             : 
     431             : static GEN
     432        4435 : FlxqE_Miller_add(void* E, GEN va, GEN vb)
     433             : {
     434        4435 :   struct _FlxqE_miller *m = (struct _FlxqE_miller *)E;
     435        4435 :   ulong p = m->p;
     436        4435 :   GEN T = m->T, a4 = m->a4, P = m->P;
     437             :   GEN v, line, point;
     438        4435 :   GEN na = gel(va,1), da = gel(va,2), pa = gel(va,3);
     439        4435 :   GEN nb = gel(vb,1), db = gel(vb,2), pb = gel(vb,3);
     440        4435 :   GEN N = Flxq_mul(na, nb, T, p);
     441        4435 :   GEN D = Flxq_mul(da, db, T, p);
     442        4435 :   line = FlxqE_chord_update(pa, pb, P, a4, T, p, &point);
     443        4435 :   N  = Flxq_mul(N, line, T, p);
     444        4435 :   v = FlxqE_vert(point, P, a4, T, p);
     445        4435 :   D = Flxq_mul(D, v, T, p); return mkvec3(N, D, point);
     446             : }
     447             : 
     448             : /* Returns the Miller function f_{m, Q} evaluated at the point P using
     449             :  * the standard Miller algorithm. */
     450             : static GEN
     451       16473 : FlxqE_Miller(GEN Q, GEN P, GEN m, GEN a4, GEN T, ulong p)
     452             : {
     453       16473 :   pari_sp av = avma;
     454             :   struct _FlxqE_miller d;
     455             :   GEN v, N, D, g1;
     456             : 
     457       16473 :   d.a4 = a4; d.T = T; d.p = p; d.P = P;
     458       16473 :   g1 = pol1_Flx(get_Flx_var(T));
     459       16473 :   v = gen_pow_i(mkvec3(g1,g1,Q), m, (void*)&d,
     460             :                 FlxqE_Miller_dbl, FlxqE_Miller_add);
     461       16473 :   N = gel(v,1); D = gel(v,2);
     462       16473 :   return gerepileupto(av, Flxq_div(N, D, T, p));
     463             : }
     464             : 
     465             : GEN
     466       12792 : FlxqE_weilpairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, ulong p)
     467             : {
     468       12792 :   pari_sp av = avma;
     469             :   GEN N, D, result;
     470       12792 :   if (ell_is_inf(P) || ell_is_inf(Q)
     471       10085 :     || (Flx_equal(gel(P,1),gel(Q,1)) && Flx_equal(gel(P,2),gel(Q,2))))
     472        4587 :     return pol1_Flx(get_Flx_var(T));
     473        8205 :   N = FlxqE_Miller(P, Q, m, a4, T, p);
     474        8205 :   D = FlxqE_Miller(Q, P, m, a4, T, p);
     475        8205 :   result = Flxq_div(N, D, T, p);
     476        8205 :   if (mpodd(m)) result = Flx_neg(result, p);
     477        8205 :   return gerepileupto(av, result);
     478             : }
     479             : 
     480             : GEN
     481          63 : FlxqE_tatepairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, ulong p)
     482             : {
     483          63 :   if (ell_is_inf(P) || ell_is_inf(Q)) return pol1_Flx(get_Flx_var(T));
     484          63 :   return FlxqE_Miller(P, Q, m, a4, T, p);
     485             : }
     486             : 
     487             : static GEN
     488       12771 : _FlxqE_pairorder(void *E, GEN P, GEN Q, GEN m, GEN F)
     489             : {
     490       12771 :   struct _FlxqE *e = (struct _FlxqE *) E;
     491       12771 :   return  Flxq_order(FlxqE_weilpairing(P,Q,m,e->a4,e->T,e->p), F, e->T, e->p);
     492             : }
     493             : 
     494             : GEN
     495       15609 : Flxq_ellgroup(GEN a4, GEN a6, GEN N, GEN T, ulong p, GEN *pt_m)
     496             : {
     497             :   struct _FlxqE e;
     498       15609 :   GEN q = powuu(p, get_Flx_degree(T));
     499       15609 :   e.a4=a4; e.a6=a6; e.T=T; e.p=p;
     500       15609 :   return gen_ellgroup(N, subiu(q,1), pt_m, (void*)&e, &FlxqE_group, _FlxqE_pairorder);
     501             : }
     502             : 
     503             : GEN
     504       14363 : Flxq_ellgens(GEN a4, GEN a6, GEN ch, GEN D, GEN m, GEN T, ulong p)
     505             : {
     506             :   GEN P;
     507       14363 :   pari_sp av = avma;
     508             :   struct _FlxqE e;
     509       14363 :   e.a4=a4; e.a6=a6; e.T=T; e.p=p;
     510       14363 :   switch(lg(D)-1)
     511             :   {
     512          63 :   case 0:
     513          63 :     return cgetg(1,t_VEC);
     514       11794 :   case 1:
     515       11794 :     P = gen_gener(gel(D,1), (void*)&e, &FlxqE_group);
     516       11794 :     P = mkvec(FlxqE_changepoint(P, ch, T, p));
     517       11794 :     break;
     518        2506 :   default:
     519        2506 :     P = gen_ellgens(gel(D,1), gel(D,2), m, (void*)&e, &FlxqE_group, _FlxqE_pairorder);
     520        2506 :     gel(P,1) = FlxqE_changepoint(gel(P,1), ch, T, p);
     521        2506 :     gel(P,2) = FlxqE_changepoint(gel(P,2), ch, T, p);
     522        2506 :     break;
     523             :   }
     524       14300 :   return gerepilecopy(av, P);
     525             : }
     526             : /***********************************************************************/
     527             : /**                                                                   **/
     528             : /**                          Point counting                           **/
     529             : /**                                                                   **/
     530             : /***********************************************************************/
     531             : 
     532             : /* assume a and n  are coprime */
     533             : static GEN
     534       76251 : RgX_circular_shallow(GEN P, long a, long n)
     535             : {
     536       76251 :   long i, l = lgpol(P);
     537       76251 :   GEN Q = cgetg(2+n,t_POL);
     538       76251 :   Q[1] = P[1];
     539      512330 :   for(i=0; i<l; i++)
     540      436079 :     gel(Q,2+(i*a)%n) = gel(P,2+i);
     541      168693 :   for(   ; i<n; i++)
     542       92442 :     gel(Q,2+(i*a)%n) = gen_0;
     543       76251 :   return normalizepol_lg(Q,2+n);
     544             : }
     545             : 
     546             : static GEN
     547       76251 : ZpXQ_frob_cyc(GEN x, GEN T, GEN q, ulong p)
     548             : {
     549       76251 :   long n = get_FpX_degree(T);
     550       76251 :   return FpX_rem(RgX_circular_shallow(x,p,n+1), T, q);
     551             : }
     552             : 
     553             : static GEN
     554      113526 : ZpXQ_frob(GEN x, GEN Xm, GEN T, GEN q, ulong p)
     555             : {
     556      113526 :   if (lg(Xm)==1)
     557       43428 :     return ZpXQ_frob_cyc(x, T, q, p);
     558             :   else
     559             :   {
     560       70098 :     long n = get_FpX_degree(T);
     561       70098 :     GEN V = RgX_blocks(RgX_inflate(x, p), n, p);
     562       70098 :     GEN W = ZXV_dotproduct(V, Xm);
     563       70098 :     return FpX_rem(W, T, q);
     564             :   }
     565             : }
     566             : 
     567             : struct _lift_lin
     568             : {
     569             :   ulong p;
     570             :   GEN sqx, Tp;
     571             :   GEN ai, Xm;
     572             : };
     573             : 
     574       84035 : static GEN _lift_invl(void *E, GEN x)
     575             : {
     576       84035 :   struct _lift_lin *d = (struct _lift_lin *) E;
     577       84035 :   GEN T = d->Tp;
     578       84035 :   ulong p = d->p;
     579       84035 :   GEN xai = Flxq_mul(ZX_to_Flx(x, p), d->ai, T, p);
     580       84035 :   return Flx_to_ZX(Flxq_lroot_fast(xai, d->sqx, T, p));
     581             : }
     582             : 
     583       23744 : static GEN _lift_lin(void *E, GEN F, GEN x2, GEN q)
     584             : {
     585       23744 :   struct _lift_lin *d = (struct _lift_lin *) E;
     586       23744 :   pari_sp av = avma;
     587       23744 :   GEN T = gel(F,3), Xm = gel(F,4);
     588       23744 :   GEN y2  = ZpXQ_frob(x2, Xm, T, q, d->p);
     589       23744 :   GEN lin = FpX_add(ZX_mul(gel(F,1), y2), ZX_mul(gel(F,2), x2), q);
     590       23744 :   return gerepileupto(av, FpX_rem(lin, T, q));
     591             : }
     592             : 
     593             : static GEN
     594      180915 : FpM_FpXV_bilinear(GEN P, GEN X, GEN Y, GEN p)
     595             : {
     596      180915 :    pari_sp av = avma;
     597      180915 :    GEN s =  ZX_mul(FpXV_FpC_mul(X,gel(P,1),p),gel(Y,1));
     598      180915 :    long i, l = lg(P);
     599      850311 :    for(i=2; i<l; i++)
     600      669396 :      s = ZX_add(s, ZX_mul(FpXV_FpC_mul(X,gel(P,i),p),gel(Y,i)));
     601      180915 :    return gerepileupto(av, FpX_red(s, p));
     602             : }
     603             : 
     604             : static GEN
     605      180915 : FpM_FpXQV_bilinear(GEN P, GEN X, GEN Y, GEN T, GEN p)
     606             : {
     607      180915 :   return FpX_rem(FpM_FpXV_bilinear(P,X,Y,p),T,p);
     608             : }
     609             : 
     610             : static GEN
     611      120624 : FpXC_powderiv(GEN M, GEN p)
     612             : {
     613             :   long i, l;
     614      120624 :   long v = varn(gel(M,2));
     615      120624 :   GEN m = cgetg_copy(M, &l);
     616      120624 :   gel(m,1) = pol_0(v);
     617      120624 :   gel(m,2) = pol_1(v);
     618      446432 :   for(i=2; i<l-1; i++)
     619      325808 :     gel(m,i+1) = FpX_Fp_mul(gel(M,i),utoi(i), p);
     620      120624 :   return m;
     621             : }
     622             : 
     623             : struct _lift_iso
     624             : {
     625             :   GEN phi;
     626             :   GEN Xm,T;
     627             :   GEN sqx, Tp;
     628             :   ulong p;
     629             : };
     630             : 
     631             : static GEN
     632       60291 : _lift_iter(void *E, GEN x2, GEN q)
     633             : {
     634       60291 :   struct _lift_iso *d = (struct _lift_iso *) E;
     635       60291 :   ulong p = d->p;
     636       60291 :   long n = lg(d->phi)-2;
     637       60291 :   GEN TN = FpXT_red(d->T, q), XN = FpXV_red(d->Xm, q);
     638       60291 :   GEN y2 = ZpXQ_frob(x2, XN, TN, q, p);
     639       60291 :   GEN xp = FpXQ_powers(x2, n, TN, q);
     640       60291 :   GEN yp = FpXQ_powers(y2, n, TN, q);
     641       60291 :   GEN V  = FpM_FpXQV_bilinear(d->phi,xp,yp,TN,q);
     642       60291 :   return mkvec3(V,xp,yp);
     643             : }
     644             : 
     645             : static GEN
     646       60291 : _lift_invd(void *E, GEN V, GEN v, GEN qM, long M)
     647             : {
     648       60291 :   struct _lift_iso *d = (struct _lift_iso *) E;
     649             :   struct _lift_lin e;
     650       60291 :   ulong p = d->p;
     651       60291 :   GEN TM = FpXT_red(d->T, qM), XM = FpXV_red(d->Xm, qM);
     652       60291 :   GEN xp = FpXV_red(gel(v,2), qM);
     653       60291 :   GEN yp = FpXV_red(gel(v,3), qM);
     654       60291 :   GEN Dx = FpM_FpXQV_bilinear(d->phi, FpXC_powderiv(xp, qM), yp, TM, qM);
     655       60291 :   GEN Dy = FpM_FpXQV_bilinear(d->phi, xp, FpXC_powderiv(yp, qM), TM, qM);
     656       60291 :   GEN F = mkvec4(Dy, Dx, TM, XM);
     657       60291 :   e.ai = Flxq_inv(ZX_to_Flx(Dy,p),d->Tp,p);
     658       60291 :   e.sqx = d->sqx; e.Tp = d->Tp; e.p=p; e.Xm = XM;
     659       60291 :   return gen_ZpX_Dixon(F,V,qM,utoipos(p),M,(void*) &e, _lift_lin, _lift_invl);
     660             : }
     661             : 
     662             : static GEN
     663       25032 : lift_isogeny(GEN phi, GEN x0, long n, GEN Xm, GEN T, GEN sqx, GEN Tp, ulong p)
     664             : {
     665             :   struct _lift_iso d;
     666       25032 :   d.phi=phi;
     667       25032 :   d.Xm=Xm; d.T=T;
     668       25032 :   d.sqx=sqx; d.Tp=Tp; d.p=p;
     669       25032 :   return gen_ZpX_Newton(x0, utoipos(p), n,(void*)&d, _lift_iter, _lift_invd);
     670             : }
     671             : 
     672             : static GEN
     673       25011 : getc2(GEN act, GEN X, GEN T, GEN q, ulong p, long N)
     674             : {
     675       25011 :   GEN A1 = RgV_to_RgX(gel(act,1),0), A2 =  RgV_to_RgX(gel(act,2),0);
     676       25011 :   long n = brent_kung_optpow(maxss(degpol(A1),degpol(A2)),2,1);
     677       25011 :   GEN xp = FpXQ_powers(X,n,T,q);
     678       25011 :   GEN P  = FpX_FpXQV_eval(A1, xp, T, q);
     679       25011 :   GEN Q  = FpX_FpXQV_eval(A2, xp, T, q);
     680       25011 :   return ZpXQ_div(P, Q, T, q, utoipos(p), N);
     681             : }
     682             : 
     683             : struct _ZpXQ_norm
     684             : {
     685             :   long n;
     686             :   GEN T, p;
     687             : };
     688             : 
     689             : static GEN
     690       32823 : ZpXQ_norm_mul(void *E, GEN x, GEN y)
     691             : {
     692       32823 :   struct _ZpXQ_norm *D = (struct _ZpXQ_norm*)E;
     693       32823 :   GEN P = gel(x,1), Q = gel(y,1);
     694       32823 :   long a = mael(x,2,1), b = mael(y,2,1);
     695       32823 :   retmkvec2(FpXQ_mul(P,ZpXQ_frob_cyc(Q, D->T, D->p, a), D->T, D->p),
     696             :             mkvecsmall((a*b)%D->n));
     697             : }
     698             : 
     699             : static GEN
     700       22715 : ZpXQ_norm_sqr(void *E, GEN x)
     701             : {
     702       22715 :   return ZpXQ_norm_mul(E, x, x);
     703             : }
     704             : 
     705             : /* Assume T = Phi_(n) and n prime */
     706             : GEN
     707       11340 : ZpXQ_norm_pcyc(GEN x, GEN T, GEN q, GEN p)
     708             : {
     709             :   GEN z;
     710             :   struct _ZpXQ_norm D;
     711       11340 :   long d = get_FpX_degree(T);
     712       11340 :   D.T = T; D.p = q; D.n = d+1;
     713       11340 :   if (d==1) return ZX_copy(x);
     714       11340 :   z = mkvec2(x,mkvecsmall(p[2]));
     715       11340 :   z = gen_powu_i(z,d,(void*)&D,ZpXQ_norm_sqr,ZpXQ_norm_mul);
     716       11340 :   return gmael(z,1,2);
     717             : }
     718             : 
     719             : /* Assume T = Phi_(n) and n prime */
     720             : static GEN
     721       11102 : ZpXQ_sqrtnorm_pcyc(GEN x, GEN T, GEN q, GEN p, long e)
     722             : {
     723       11102 :   GEN z = ZpXQ_norm_pcyc(x, T, q, p);
     724       11102 :   return Zp_sqrtlift(z,Fp_sqrt(z,p),p,e);
     725             : }
     726             : 
     727             : /* Assume a = 1 [p], return the square root of the norm */
     728             : static GEN
     729       13930 : ZpXQ_sqrtnorm(GEN a, GEN T, GEN q, GEN p, long e)
     730             : {
     731       13930 :   GEN s = Fp_div(FpXQ_trace(ZpXQ_log(a, T, p, e), T, q), gen_2, q);
     732       13930 :   return modii(gel(Qp_exp(cvtop(s, p, e-1)),4), q);
     733             : }
     734             : 
     735             : struct _teich_lin
     736             : {
     737             :   ulong p;
     738             :   GEN sqx, Tp;
     739             :   long m;
     740             : };
     741             : 
     742             : static GEN
     743       29470 : _teich_invl(void *E, GEN x)
     744             : {
     745       29470 :   struct _teich_lin *d = (struct _teich_lin *) E;
     746       29470 :   ulong p = d->p;
     747       29470 :   GEN T = d->Tp;
     748       29470 :   return Flx_to_ZX(Flxq_lroot_fast(ZX_to_Flx(x, p), d->sqx, T, p));
     749             : }
     750             : 
     751             : static GEN
     752        8953 : _teich_lin(void *E, GEN F, GEN x2, GEN q)
     753             : {
     754        8953 :   struct _teich_lin *d = (struct _teich_lin *) E;
     755        8953 :   pari_sp av = avma;
     756        8953 :   GEN T = gel(F,2), Xm = gel(F,3);
     757        8953 :   GEN y2  = ZpXQ_frob(x2, Xm, T, q, d->p);
     758        8953 :   GEN lin = FpX_sub(y2, ZX_mulu(ZX_mul(gel(F,1), x2), d->p), q);
     759        8953 :   return gerepileupto(av, FpX_rem(lin, T, q));
     760             : }
     761             : 
     762             : struct _teich_iso
     763             : {
     764             :   GEN Xm, T;
     765             :   GEN sqx, Tp;
     766             :   ulong p;
     767             : };
     768             : 
     769             : static GEN
     770       20517 : _teich_iter(void *E, GEN x2, GEN q)
     771             : {
     772       20517 :   struct _teich_iso *d = (struct _teich_iso *) E;
     773       20517 :   ulong p = d->p;
     774       20517 :   GEN TN = FpXT_red(d->T, q), XN = FpXV_red(d->Xm, q);
     775       20517 :   GEN y2 = ZpXQ_frob(x2, XN, TN, q, d->p);
     776       20517 :   GEN x1 = FpXQ_powu(x2, p-1, TN, q);
     777       20517 :   GEN xp = FpXQ_mul(x2, x1, TN, q);
     778       20517 :   GEN V = FpX_sub(y2,xp,q);
     779       20517 :   return mkvec2(V,x1);
     780             : }
     781             : 
     782             : static GEN
     783       20517 : _teich_invd(void *E, GEN V, GEN v, GEN qM, long M)
     784             : {
     785       20517 :   struct _teich_iso *d = (struct _teich_iso *) E;
     786             :   struct _teich_lin e;
     787       20517 :   ulong p = d->p;
     788       20517 :   GEN TM = FpXT_red(d->T, qM), XM = FpXV_red(d->Xm, qM);
     789       20517 :   GEN x1 = FpX_red(gel(v,2), qM);
     790       20517 :   GEN F = mkvec3(x1, TM, XM);
     791       20517 :   e.sqx = d->sqx; e.Tp = d->Tp; e.p=p;
     792       20517 :   return gen_ZpX_Dixon(F,V,qM,utoipos(p),M,(void*) &e, _teich_lin, _teich_invl);
     793             : }
     794             : 
     795             : static GEN
     796       10213 : Teichmuller_lift(GEN x, GEN Xm, GEN T, GEN sqx, GEN Tp, ulong p, long N)
     797             : {
     798             :   struct _teich_iso d;
     799       10213 :   d.Xm = Xm; d.T = T; d.sqx = sqx; d.Tp = Tp; d.p = p;
     800       10213 :   return gen_ZpX_Newton(x,utoipos(p), N,(void*)&d, _teich_iter, _teich_invd);
     801             : }
     802             : 
     803             : static GEN
     804       25032 : get_norm(GEN a4, GEN a6, GEN T, ulong p, long N)
     805             : {
     806       25032 :   long sv=T[1];
     807             :   GEN a;
     808       25032 :   if (p==3) a = gel(a4,1);
     809             :   else
     810             :   {
     811       10227 :     GEN P = mkpoln(4, pol1_Flx(sv), pol0_Flx(sv), a4, a6);
     812       10227 :     a = gel(FlxqX_powu(P,p>>1,T,p),2+p-1);
     813             :   }
     814       25032 :   return Zp_sqrtnlift(gen_1,subss(p,1),utoi(Flxq_norm(a,T,p)),utoipos(p), N);
     815             : }
     816             : 
     817             : static GEN
     818       25011 : fill_pols(long n, const long *v, long m, const long *vn,
     819             :           const long *vd, GEN *act)
     820             : {
     821             :   long i, j;
     822       25011 :   long d = upowuu(n,12/(n-1));
     823       25011 :   GEN N, D, M = zeromatcopy(n+1,n+1);
     824       25011 :   gmael(M,1,n+1) = gen_1;
     825      120568 :   for(i=2;i<=n+1;i++)
     826      338373 :     for(j=i-1;j<=n;j++)
     827      242816 :       gmael(M,i,j) = mulis(powuu(d,i-2),v[j-i+1]);
     828       25011 :   N = cgetg(m+1,t_COL);
     829       25011 :   D = cgetg(m+1,t_COL);
     830      135359 :   for(i=1;i<=m;i++)
     831             :   {
     832      110348 :     gel(N,i) = stoi(*vn++);
     833      110348 :     gel(D,i) = stoi(*vd++);
     834             :   }
     835       25011 :   *act = mkmat2(N,D);
     836       25011 :   return M;
     837             : }
     838             : 
     839             : /*
     840             :   These polynomials were extracted from the ECHIDNA databases
     841             :   available at <http://echidna.maths.usyd.edu.au/echidna/>
     842             :   and computed by David R. Kohel.
     843             :   Return the matrix of the modular polynomial, set act to the parametrization,
     844             :   and set dj to the opposite of the supersingular j-invariant.
     845             : */
     846             : static GEN
     847       25011 : get_Kohel_polynomials(ulong p, GEN *act, long *dj)
     848             : {
     849       25011 :   const long mat3[] = {-1,-36,-270};
     850       25011 :   const long num3[] = {1,-483,-21141,-59049};
     851       25011 :   const long den3[] = {1,261, 4347, -6561};
     852       25011 :   const long mat5[] = {-1,-30,-315,-1300,-1575};
     853       25011 :   const long num5[] = {-1,490,20620,158750,78125};
     854       25011 :   const long den5[] = {-1,-254,-4124,-12250,3125};
     855       25011 :   const long mat7[] = {-1,-28,-322,-1904,-5915,-8624,-4018};
     856       25011 :   const long num7[] = {1,-485,-24058,-343833,-2021642,-4353013,-823543};
     857       25011 :   const long den7[] = {1,259,5894,49119,168406,166355,-16807};
     858       25011 :   const long mat13[]= {-1,-26,-325,-2548,-13832,-54340,-157118,-333580,-509366,
     859             :                        -534820,-354536,-124852,-15145};
     860       25011 :   const long num13[]= {1,-487,-24056,-391463,-3396483,-18047328,-61622301,
     861             :                        -133245853,-168395656,-95422301,-4826809};
     862       25011 :   const long den13[]= {1,257,5896,60649,364629,1388256,3396483,5089019,4065464,
     863             :                        1069939,-28561};
     864       25011 :   switch(p)
     865             :   {
     866       14805 :   case 3:
     867       14805 :     *dj = 0;
     868       14805 :     return fill_pols(3,mat3,4,num3,den3,act);
     869       10171 :   case 5:
     870       10171 :     *dj = 0;
     871       10171 :     return fill_pols(5,mat5,5,num5,den5,act);
     872          28 :   case 7:
     873          28 :     *dj = 1;
     874          28 :     return fill_pols(7,mat7,7,num7,den7,act);
     875           7 :   case 13:
     876           7 :     *dj = 8;
     877           7 :     return fill_pols(13,mat13,11,num13,den13,act);
     878             :   }
     879             :   *dj=0; *act = NULL; return NULL; /* LCOV_EXCL_LINE */
     880             : }
     881             : 
     882             : long
     883       32476 : zx_is_pcyc(GEN T)
     884             : {
     885       32476 :   long i, n = degpol(T);
     886       32476 :   if (!uisprime(n+1))
     887       11840 :     return 0;
     888       99148 :   for (i=0; i<=n; i++)
     889       87808 :     if (T[i+2]!=1UL)
     890        9296 :       return 0;
     891       11340 :   return 1;
     892             : }
     893             : 
     894             : static GEN
     895       25011 : Flxq_ellcard_Kohel(GEN a4, GEN a6, GEN T, ulong p)
     896             : {
     897       25011 :   pari_sp av = avma, av2;
     898             :   pari_timer ti;
     899       25011 :   long n = get_Flx_degree(T), N = (n+4)/2, dj;
     900       25011 :   GEN q = powuu(p, N);
     901             :   GEN T2, Xm, s1, c2, t, lr;
     902             :   GEN S1, sqx;
     903             :   GEN Nc2, Np;
     904       25011 :   GEN act, phi = get_Kohel_polynomials(p, &act, &dj);
     905       25011 :   long ispcyc = zx_is_pcyc(get_Flx_mod(T));
     906       25011 :   timer_start(&ti);
     907       25011 :   if (!ispcyc)
     908             :   {
     909       13916 :     T2 = Flx_Teichmuller(get_Flx_mod(T),p,N);
     910       13916 :     if (DEBUGLEVEL) timer_printf(&ti,"Teich");
     911             :   } else
     912       11095 :     T2 = Flx_to_ZX(get_Flx_mod(T));
     913       25011 :   T2 = FpX_get_red(T2, q); T = ZXT_to_FlxT(T2, p);
     914       25011 :   av2 = avma;
     915       25011 :   if (DEBUGLEVEL) timer_printf(&ti,"Barrett");
     916       25011 :   if (!ispcyc)
     917             :   {
     918       13916 :     Xm = FpXQ_powers(pol_xn(n,get_FpX_var(T2)),p-1,T2,q);
     919       13916 :     if (DEBUGLEVEL) timer_printf(&ti,"Xm");
     920             :   } else
     921       11095 :     Xm = cgetg(1,t_VEC);
     922       25011 :   s1 = Flxq_inv(Flx_Fl_add(Flxq_ellj(a4,a6,T,p),dj, p),T,p);
     923       25011 :   lr = Flxq_lroot(polx_Flx(get_Flx_var(T)), T, p);
     924       25011 :   sqx = Flxq_powers(lr, p-1, T, p);
     925       25011 :   S1 = lift_isogeny(phi, Flx_to_ZX(s1), N, Xm, T2, sqx, T ,p);
     926       25011 :   if (DEBUGLEVEL) timer_printf(&ti,"Lift isogeny");
     927       25011 :   c2 = getc2(act, S1, T2, q, p, N);
     928       25011 :   if (DEBUGLEVEL) timer_printf(&ti,"c^2");
     929       25011 :   if (p>3 && !ispcyc)
     930             :   {
     931       10199 :     GEN c2p = Flx_to_ZX(Flxq_inv(ZX_to_Flx(c2,p),T,p));
     932       10199 :     GEN tc2 = Teichmuller_lift(c2p,Xm, T2,sqx,T,p,N);
     933       10199 :     if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fq");
     934       10199 :     c2 = FpX_rem(FpX_mul(tc2,c2,q),T2,q);
     935             :   }
     936       25011 :   c2 = gerepileupto(av2, c2);
     937       25011 :   if (DEBUGLEVEL) timer_printf(&ti,"tc2");
     938       25011 :   Nc2 = (ispcyc? ZpXQ_sqrtnorm_pcyc: ZpXQ_sqrtnorm)(c2, T2, q, utoipos(p), N);
     939       25011 :   if (DEBUGLEVEL) timer_printf(&ti,"Norm");
     940       25011 :   Np = get_norm(a4,a6,T,p,N);
     941       25011 :   if (p>3 && ispcyc)
     942             :   {
     943           7 :     GEN Ncpi =  utoi(Fl_inv(umodiu(Nc2,p), p));
     944           7 :     GEN tNc2 = Zp_sqrtnlift(gen_1, subss(p,1), Ncpi, utoipos(p),N);
     945           7 :     if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fp");
     946           7 :     Nc2 = Fp_mul(Nc2,tNc2,q);
     947             :   }
     948       25011 :   t = Fp_center_i(Fp_mul(Nc2,Np,q),q,shifti(q,-1));
     949       25011 :   return gerepileupto(av, subii(addiu(powuu(p,n),1),t));
     950             : }
     951             : 
     952             : /* Use Damien Robert method */
     953             : 
     954             : static GEN
     955          21 : get_trace_Robert(GEN J, GEN phi, GEN Xm, GEN T, GEN q, ulong p, long e)
     956             : {
     957          21 :   long n = lg(phi)-2;
     958          21 :   GEN K = ZpXQ_frob(J, Xm, T, q, p);
     959          21 :   GEN Jp = FpXQ_powers(J, n, T, q);
     960          21 :   GEN Kp = FpXQ_powers(K, n, T, q);
     961          21 :   GEN Jd = FpXC_powderiv(Jp, q);
     962          21 :   GEN Kd = FpXC_powderiv(Kp, q);
     963          21 :   GEN Dx = FpM_FpXQV_bilinear(phi, Kd, Jp, T, q);
     964          21 :   GEN Dy = FpM_FpXQV_bilinear(phi, Kp, Jd, T, q);
     965          21 :   GEN C = ZpXQ_inv(ZX_divuexact(Dy, p), T, utoi(p), e);
     966          21 :   return FpX_neg(FpXQ_mul(Dx, C, T, q), q);
     967             : }
     968             : 
     969             : static GEN
     970          21 : Flxq_ellcard_Harley(GEN a4, GEN a6, GEN T, ulong p)
     971             : {
     972          21 :   pari_sp av = avma, av2;
     973             :   pari_timer ti;
     974          21 :   long n = get_Flx_degree(T), N = (n+5)/2;
     975          21 :   GEN pp = utoipos(p), q = powuu(p, N);
     976             :   GEN T2, j, t, phi;
     977             :   GEN J1,sqx,Xm;
     978             :   GEN c2, tc2, c2p, Nc2, Np;
     979          21 :   long ispcyc = zx_is_pcyc(get_Flx_mod(T));
     980          21 :   timer_start(&ti);
     981          21 :   if (!ispcyc)
     982             :   {
     983          14 :     T2 = Flx_Teichmuller(get_Flx_mod(T),p,N);
     984          14 :     if (DEBUGLEVEL) timer_printf(&ti,"Teich");
     985             :   } else
     986           7 :     T2 = Flx_to_ZX(get_Flx_mod(T));
     987          21 :   T2 = FpX_get_red(T2, q); T = ZXT_to_FlxT(T2, p);
     988          21 :   av2 = avma;
     989          21 :   if (DEBUGLEVEL) timer_printf(&ti,"Barrett");
     990          21 :   if (!ispcyc)
     991             :   {
     992          14 :     Xm = FpXQ_powers(pol_xn(n,get_FpX_var(T2)),p-1,T2,q);
     993          14 :     if (DEBUGLEVEL) timer_printf(&ti,"Xm");
     994             :   } else
     995           7 :     Xm = cgetg(1,t_VEC);
     996          21 :   j = Flxq_ellj(a4,a6,T,p);
     997          21 :   sqx = Flxq_powers(Flxq_lroot(polx_Flx(T[1]), T, p), p-1, T, p);
     998          21 :   phi = polmodular_ZM(p, 0);
     999          21 :   if (DEBUGLEVEL) timer_printf(&ti,"phi");
    1000          21 :   J1 = lift_isogeny(phi, Flx_to_ZX(j), N, Xm, T2,sqx,T,p);
    1001          21 :   if (DEBUGLEVEL) timer_printf(&ti,"Lift isogeny");
    1002          21 :   c2 = get_trace_Robert(J1, phi, Xm, T2, q, p, N);
    1003          21 :   q = diviuexact(q,p); N--;
    1004          21 :   if (DEBUGLEVEL) timer_printf(&ti,"c^2");
    1005          21 :   if (!ispcyc)
    1006             :   {
    1007          14 :     c2p = Flx_to_ZX(Flxq_inv(ZX_to_Flx(c2,p),T,p));
    1008          14 :     tc2 = Teichmuller_lift(c2p,Xm, T2,sqx,T,p,N);
    1009          14 :     if (DEBUGLEVEL) timer_printf(&ti,"teichmuller");
    1010          14 :     c2 = FpX_rem(FpX_mul(tc2,c2,q),T2,q);
    1011             :   }
    1012          21 :   c2 = gerepileupto(av2, c2);
    1013          21 :   q = powuu(p, N);
    1014          21 :   Nc2 = (ispcyc? ZpXQ_sqrtnorm_pcyc: ZpXQ_sqrtnorm)(c2, T2, q, pp, N);
    1015          21 :   if (DEBUGLEVEL) timer_printf(&ti,"Norm");
    1016          21 :   Np = get_norm(a4,a6,T,p,N);
    1017          21 :   if (ispcyc)
    1018             :   {
    1019           7 :     GEN Ncpi = utoi(Fl_inv(umodiu(Nc2,p), p));
    1020           7 :     GEN tNc2 = Zp_sqrtnlift(gen_1, subss(p,1), Ncpi, pp, N);
    1021           7 :     if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fp");
    1022           7 :     Nc2 = Fp_mul(Nc2,tNc2,q);
    1023             :   }
    1024          21 :   t = Fp_center_i(Fp_mul(Nc2,Np,q),q,shifti(q,-1));
    1025          21 :   return gerepileupto(av, subii(addiu(powuu(p,n),1),t));
    1026             : }
    1027             : 
    1028             : /***************************************************************************/
    1029             : /*                                                                         */
    1030             : /*                          Shanks Mestre                                  */
    1031             : /*                                                                         */
    1032             : /***************************************************************************/
    1033             : 
    1034             : /* Return the lift of a (mod b), which is closest to h */
    1035             : static GEN
    1036        2309 : closest_lift(GEN a, GEN b, GEN h)
    1037             : {
    1038        2309 :   return addii(a, mulii(b, diviiround(subii(h,a), b)));
    1039             : }
    1040             : 
    1041             : static GEN
    1042        1273 : FlxqE_find_order(GEN f, GEN h, GEN bound, GEN B, GEN a4, GEN T, ulong p)
    1043             : {
    1044        1273 :   pari_sp av = avma, av1;
    1045             :   pari_timer Ti;
    1046        1273 :   long s = itos( gceil(gsqrt(gdiv(bound,B),DEFAULTPREC)) ) >> 1;
    1047             :   GEN tx, ti;
    1048        1273 :   GEN fh = FlxqE_mul(f, h, a4, T, p);
    1049        1273 :   GEN F, P = fh, fg;
    1050             :   long i;
    1051        1273 :   if (DEBUGLEVEL >= 6) timer_start(&Ti);
    1052        1273 :   if (ell_is_inf(fh)) return h;
    1053        1149 :   F = FlxqE_mul(f, B, a4, T, p);
    1054        1149 :   if (s < 3)
    1055             :   { /* we're nearly done: naive search */
    1056         246 :     GEN Q = P;
    1057         246 :     for (i=1;; i++)
    1058             :     {
    1059         690 :       P = FlxqE_add(P, F, a4, T, p); /* h.f + i.F */
    1060         690 :       if (ell_is_inf(P)) return gerepileupto(av, addii(h, mului(i,B)));
    1061         620 :       Q = FlxqE_sub(Q, F, a4, T, p); /* h.f - i.F */
    1062         620 :       if (ell_is_inf(Q)) return gerepileupto(av, subii(h, mului(i,B)));
    1063             :     }
    1064             :   }
    1065         903 :   tx = cgetg(s+1,t_VECSMALL);
    1066             :   /* Baby Step/Giant Step */
    1067         903 :   av1 = avma;
    1068        5224 :   for (i=1; i<=s; i++)
    1069             :   { /* baby steps */
    1070        4453 :     tx[i] = hash_GEN(gel(P, 1));
    1071        4453 :     P = FlxqE_add(P, F, a4, T, p); /* h.f + i.F */
    1072        4453 :     if (ell_is_inf(P)) return gerepileupto(av, addii(h, mului(i,B)));
    1073        4321 :     if (gc_needed(av1,3))
    1074             :     {
    1075           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"[Flxq_ellcard] baby steps, i=%ld",i);
    1076           0 :       P = gerepileupto(av1,P);
    1077             :     }
    1078             :   }
    1079         771 :   if (DEBUGLEVEL >= 6) timer_printf(&Ti, "[Flxq_ellcard] baby steps, s = %ld",s);
    1080             :   /* giant steps: fg = s.F */
    1081         771 :   fg = gerepileupto(av1, FlxqE_sub(P, fh, a4, T, p));
    1082         771 :   if (ell_is_inf(fg)) return gerepileupto(av,mului(s,B));
    1083         771 :   ti = vecsmall_indexsort(tx); /* = permutation sorting tx */
    1084         771 :   tx = perm_mul(tx,ti);
    1085         771 :   if (DEBUGLEVEL >= 6) timer_printf(&Ti, "[Flxq_ellcard] sorting");
    1086         771 :   av1 = avma;
    1087         771 :   for (P=fg, i=1; ; i++)
    1088        2917 :   {
    1089        3688 :     long k = hash_GEN(gel(P,1));
    1090        3688 :     long r = zv_search(tx, k);
    1091        3688 :     if (r)
    1092             :     {
    1093        1543 :       while (r && tx[r] == k) r--;
    1094         771 :       for (r++; r <= s && tx[r] == k; r++)
    1095             :       {
    1096         771 :         long j = ti[r]-1;
    1097         771 :         GEN Q = FlxqE_add(FlxqE_mul(F, stoi(j), a4, T, p), fh, a4, T, p);
    1098         771 :         if (DEBUGLEVEL >= 6)
    1099           0 :           timer_printf(&Ti, "[Flxq_ellcard] giant steps, i = %ld",i);
    1100         771 :         if (Flx_equal(gel(P,1), gel(Q,1)))
    1101             :         {
    1102         771 :           if (Flx_equal(gel(P,2), gel(Q,2))) i = -i;
    1103         771 :           return gerepileupto(av,addii(h, mulii(addis(mulss(s,i), j), B)));
    1104             :         }
    1105             :       }
    1106             :     }
    1107        2917 :     P = FlxqE_add(P,fg,a4,T,p);
    1108        2917 :     if (gc_needed(av1,3))
    1109             :     {
    1110           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"[Flxq_ellcard] giants steps, i=%ld",i);
    1111           0 :       P = gerepileupto(av1,P);
    1112             :     }
    1113             :   }
    1114             : }
    1115             : 
    1116             : static void
    1117       33614 : Flx_next(GEN t, ulong p)
    1118             : {
    1119             :   long i;
    1120       33614 :   for(i=2;;i++)
    1121       41573 :     if (uel(t,i)==p-1)
    1122        7959 :       t[i]=0;
    1123             :     else
    1124             :     {
    1125       33614 :       t[i]++;
    1126       33614 :       break;
    1127             :     }
    1128       33614 : }
    1129             : 
    1130             : static void
    1131       33614 : Flx_renormalize_ip(GEN x, long lx)
    1132             : {
    1133             :   long i;
    1134       41573 :   for (i = lx-1; i>=2; i--)
    1135       38605 :     if (x[i]) break;
    1136       33614 :   setlg(x, i+1);
    1137       33614 : }
    1138             : 
    1139             : static ulong
    1140        2282 : F3xq_ellcard_naive(GEN a2, GEN a6, GEN T)
    1141             : {
    1142        2282 :   pari_sp av = avma;
    1143        2282 :   long i, d = get_Flx_degree(T), lx = d+2;
    1144        2282 :   long q = upowuu(3, d), a;
    1145        2282 :   GEN x = zero_zv(lx); x[1] = get_Flx_var(T);
    1146       11354 :   for(a=1, i=0; i<q; i++)
    1147             :   {
    1148             :     GEN rhs;
    1149        9072 :     Flx_renormalize_ip(x, lx);
    1150        9072 :     rhs = Flx_add(Flxq_mul(Flxq_sqr(x, T, 3), Flx_add(x, a2, 3), T, 3), a6, 3);
    1151        9072 :     if (!lgpol(rhs)) a++; else if (Flxq_issquare(rhs, T, 3)) a+=2;
    1152        9072 :     Flx_next(x, 3);
    1153             :   }
    1154        2282 :   set_avma(av);
    1155        2282 :   return a;
    1156             : }
    1157             : 
    1158             : static ulong
    1159         686 : Flxq_ellcard_naive(GEN a4, GEN a6, GEN T, ulong p)
    1160             : {
    1161         686 :   pari_sp av = avma;
    1162         686 :   long i, d = get_Flx_degree(T), lx = d+2;
    1163         686 :   long q = upowuu(p, d), a;
    1164         686 :   GEN x = zero_zv(lx); x[1] = get_Flx_var(T);
    1165       25228 :   for(a=1, i=0; i<q; i++)
    1166             :   {
    1167             :     GEN x2, rhs;
    1168       24542 :     Flx_renormalize_ip(x, lx);
    1169       24542 :     x2  = Flxq_sqr(x, T, p);
    1170       24542 :     rhs = Flx_add(Flxq_mul(x, Flx_add(x2, a4, p), T, p), a6, p);
    1171       24542 :     if (!lgpol(rhs)) a++; else if (Flxq_issquare(rhs,T,p)) a+=2;
    1172       24542 :     Flx_next(x,p);
    1173             :   }
    1174         686 :   set_avma(av);
    1175         686 :   return a;
    1176             : }
    1177             : 
    1178             : /* assume T irreducible mod p, m = (q-1)/(p-1) */
    1179             : static long
    1180        2611 : Flxq_kronecker(GEN x, GEN m, GEN T, ulong p)
    1181             : {
    1182             :   pari_sp av;
    1183             :   ulong z;
    1184        2611 :   if (lgpol(x) == 0) return 0;
    1185        2598 :   av = avma; z = Flxq_pow(x, m, T, p)[2];
    1186        2598 :   return gc_long(av, krouu(z, p));
    1187             : }
    1188             : 
    1189             : /* Find x such that kronecker(u = x^3+a4x+a6, p) is KRO.
    1190             :  * Return point [x*u,u^2] on E (KRO=1) / E^twist (KRO=-1) */
    1191             : static GEN
    1192        2611 : Flxq_ellpoint(long KRO, GEN a4, GEN a6, GEN m, long n, long vn, GEN T, ulong p)
    1193             : {
    1194             :   for(;;)
    1195        1338 :   {
    1196        2611 :     GEN x = random_Flx(n,vn,p);
    1197        2611 :     GEN u = Flx_add(a6, Flxq_mul(Flx_add(a4, Flxq_sqr(x,T,p), p), x, T,p), p);
    1198        2611 :     if (Flxq_kronecker(u, m,T,p) == KRO)
    1199        1273 :       return mkvec2(Flxq_mul(u,x, T,p), Flxq_sqr(u, T,p));
    1200             :   }
    1201             : }
    1202             : 
    1203             : static GEN
    1204        1036 : Flxq_ellcard_Shanks(GEN a4, GEN a6, GEN q, GEN T, ulong p)
    1205             : {
    1206        1036 :   pari_sp av = avma;
    1207        1036 :   long vn = get_Flx_var(T), n = get_Flx_degree(T), KRO = -1;
    1208             :   GEN h,f, ta4, A, B, m;
    1209        1036 :   GEN q1p = addiu(q,1), q2p = shifti(q1p, 1);
    1210        1036 :   GEN bound = addiu(sqrti(gmul2n(q,4)), 1); /* ceil( 4sqrt(q) ) */
    1211             :   /* once #E(Flxq) is know mod B >= bound, it is completely determined */
    1212             :   /* how many 2-torsion points ? */
    1213        1036 :   switch(FlxqX_nbroots(mkpoln(4, pol1_Flx(vn), pol0_Flx(vn), a4, a6), T, p))
    1214             :   {
    1215         385 :   case 3:  A = gen_0; B = utoipos(4); break;
    1216         245 :   case 1:  A = gen_0; B = gen_2; break;
    1217         406 :   default: A = gen_1; B = gen_2; break; /* 0 */
    1218             :   }
    1219        1036 :   m = diviuexact(subiu(powuu(p,n), 1), p-1);
    1220             :   for(;;)
    1221             :   {
    1222        1273 :     h = closest_lift(A, B, q1p);
    1223             :     /* [ux, u^2] is on E_u: y^2 = x^3 + c4 u^2 x + c6 u^3
    1224             :      * E_u isomorphic to E (resp. E') iff KRO = 1 (resp. -1)
    1225             :      * #E(F_p) = p+1 - a_p, #E'(F_p) = p+1 + a_p
    1226             :      *
    1227             :      * #E_u(Flxq) = A (mod B),  h is close to #E_u(Flxq) */
    1228        1273 :     KRO = -KRO;
    1229        1273 :     f = Flxq_ellpoint(KRO, a4,a6, m,n,vn, T,p);
    1230             : 
    1231        1273 :     ta4 = Flxq_mul(a4, gel(f,2), T, p); /* a4 for E_u */
    1232        1273 :     h = FlxqE_find_order(f, h, bound, B, ta4,T,p);
    1233        1273 :     h = FlxqE_order(f, h, ta4, T, p);
    1234             :     /* h | #E_u(Flxq) = A (mod B) */
    1235        1273 :     A = Z_chinese_all(A, gen_0, B, h, &B);
    1236        1273 :     if (cmpii(B, bound) >= 0) break;
    1237             :     /* not done, update A mod B for the _next_ curve, isomorphic to
    1238             :      * the quadratic twist of this one */
    1239         237 :     A = remii(subii(q2p,A), B); /* #E(Fq)+#E'(Fq) = 2q+2 */
    1240             :   }
    1241        1036 :   h = closest_lift(A, B, q1p);
    1242        1036 :   return gerepileuptoint(av, KRO == 1? h: subii(q2p,h));
    1243             : }
    1244             : 
    1245             : static GEN
    1246       17087 : F3xq_ellcard(GEN a2, GEN a6, GEN T)
    1247             : {
    1248       17087 :   long n = get_Flx_degree(T);
    1249       17087 :   if (n <= 2)
    1250        1981 :     return utoi(F3xq_ellcard_naive(a2, a6, T));
    1251             :   else
    1252             :   {
    1253       15106 :     GEN q1 = addiu(powuu(3, get_Flx_degree(T)), 1), t;
    1254       15106 :     GEN a = Flxq_div(a6,Flxq_powu(a2,3,T,3),T,3);
    1255       15106 :     if (Flx_equal1(Flxq_powu(a, 8, T, 3)))
    1256             :     {
    1257         301 :       GEN P = Flxq_minpoly(a,T,3);
    1258         301 :       long dP = degpol(P); /* dP <= 2 */
    1259         301 :       ulong q = upowuu(3,dP);
    1260         301 :       GEN A2 = pol1_Flx(P[1]), A6 = Flx_rem(polx_Flx(P[1]), P, 3);
    1261         301 :       long tP = q + 1 - F3xq_ellcard_naive(A2, A6, P);
    1262         301 :       t = elltrace_extension(stoi(tP), n/dP, utoi(q));
    1263         301 :       if (umodiu(t, 3)!=1) t = negi(t);
    1264         301 :       return Flx_equal1(a2) || Flxq_issquare(a2,T,3) ? subii(q1,t): addii(q1,t);
    1265             :     }
    1266       14805 :     else return Flxq_ellcard_Kohel(mkvec(a2), a6, T, 3);
    1267             :   }
    1268             : }
    1269             : 
    1270             : static GEN
    1271       10913 : Flxq_ellcard_Satoh(GEN a4, GEN a6, GEN j, GEN T, ulong p)
    1272             : {
    1273       10913 :   long n = get_Flx_degree(T);
    1274       10913 :   if (n <= 2)
    1275         406 :     return utoi(Flxq_ellcard_naive(a4, a6, T, p));
    1276             :   else
    1277             :   {
    1278       10507 :     GEN jp = Flxq_powu(j, p, T, p);
    1279       10507 :     GEN s = Flx_add(j, jp, p);
    1280       10507 :     if (degpol(s) <= 0)
    1281             :     { /* it is assumed j not in F_p */
    1282         280 :       GEN m = Flxq_mul(j, jp, T, p);
    1283         280 :       if (degpol(m) <= 0)
    1284             :       {
    1285         280 :         GEN q = sqru(p);
    1286         280 :         GEN q1 = addiu(powuu(p, get_Flx_degree(T)), 1);
    1287         280 :         GEN sk = Flx_Fl_add(Flx_neg(j, p), 1728%p, p);
    1288         280 :         GEN sA4 = Flx_triple(Flxq_mul(sk, j, T, p), p);
    1289         280 :         GEN u = Flxq_div(a4, sA4, T, p);
    1290         280 :         ulong ns = lgpol(s) ? Fl_neg(s[2], p): 0UL;
    1291         280 :         GEN P = mkvecsmall4(T[1], m[2], ns, 1L);
    1292             :         GEN A4, A6, t, tP;
    1293         280 :         Flxq_ellj_to_a4a6(polx_Flx(T[1]), P, p, &A4, &A6);
    1294         280 :         tP = addis(q, 1 - Flxq_ellcard_naive(A4, A6, P, p));
    1295         280 :         t = elltrace_extension(tP, n>>1, q);
    1296         280 :         return Flxq_is2npower(u, 2, T, p) ? subii(q1,t): addii(q1,t);
    1297             :       }
    1298             :     }
    1299       10227 :     if (p<=7 || p==13 ) return Flxq_ellcard_Kohel(a4, a6, T, p);
    1300          21 :     else return Flxq_ellcard_Harley(a4, a6, T, p);
    1301             :   }
    1302             : }
    1303             : 
    1304             : static GEN
    1305           0 : Flxq_ellcard_Kedlaya(GEN a4, GEN a6, GEN T, ulong p)
    1306             : {
    1307           0 :   pari_sp av = avma;
    1308           0 :   GEN H = mkpoln(4, gen_1, gen_0, Flx_to_ZX(a4), Flx_to_ZX(a6));
    1309           0 :   GEN Tp = Flx_to_ZX(get_Flx_mod(T));
    1310           0 :   long n = degpol(Tp), e = ((p < 16 ? n+1: n)>>1)+1;
    1311           0 :   GEN M = ZlXQX_hyperellpadicfrobenius(H, Tp, p, e);
    1312           0 :   GEN N = ZpXQM_prodFrobenius(M, Tp, utoipos(p), e);
    1313           0 :   GEN q = powuu(p, e);
    1314           0 :   GEN tp = Fq_add(gcoeff(N,1,1), gcoeff(N,2,2), Tp, q);
    1315           0 :   GEN t = Fp_center_i(typ(tp)==t_INT ? tp: leading_coeff(tp), q, shifti(q,-1));
    1316           0 :   return gerepileupto(av, subii(addiu(powuu(p, n), 1), t));
    1317             : }
    1318             : 
    1319             : GEN
    1320       51820 : Flxq_ellj(GEN a4, GEN a6, GEN T, ulong p)
    1321             : {
    1322       51820 :   pari_sp av=avma;
    1323       51820 :   if (p==3)
    1324             :   {
    1325             :     GEN J;
    1326       14805 :     if (typ(a4)!=t_VEC) return pol0_Flx(get_Flx_var(T));
    1327       14805 :     J = Flxq_div(Flxq_powu(gel(a4,1),3, T, p),Flx_neg(a6,p), T, p);
    1328       14805 :     return gerepileuptoleaf(av, J);
    1329             :   }
    1330             :   else
    1331             :   {
    1332       37015 :     pari_sp av=avma;
    1333       37015 :     GEN a43 = Flxq_mul(a4,Flxq_sqr(a4,T,p),T,p);
    1334       37015 :     GEN a62 = Flxq_sqr(a6,T,p);
    1335       37015 :     GEN num = Flx_mulu(a43,6912,p);
    1336       37015 :     GEN den = Flx_add(Flx_mulu(a43,4,p),Flx_mulu(a62,27,p),p);
    1337       37015 :     return gerepileuptoleaf(av, Flxq_div(num, den, T, p));
    1338             :   }
    1339             : }
    1340             : 
    1341             : void
    1342         280 : Flxq_ellj_to_a4a6(GEN j, GEN T, ulong p, GEN *pt_a4, GEN *pt_a6)
    1343             : {
    1344         280 :   ulong zagier = 1728 % p;
    1345         280 :   if (lgpol(j)==0)
    1346           0 :     { *pt_a4 = pol0_Flx(T[1]); *pt_a6 =pol1_Flx(T[1]); }
    1347         280 :   else if (lgpol(j)==1 && uel(j,2) == zagier)
    1348           0 :     { *pt_a4 = pol1_Flx(T[1]); *pt_a6 =pol0_Flx(T[1]); }
    1349             :   else
    1350             :   {
    1351         280 :     GEN k = Flx_Fl_add(Flx_neg(j, p), zagier, p);
    1352         280 :     GEN kj = Flxq_mul(k, j, T, p);
    1353         280 :     GEN k2j = Flxq_mul(kj, k, T, p);
    1354         280 :     *pt_a4 = Flx_triple(kj, p);
    1355         280 :     *pt_a6 = Flx_double(k2j, p);
    1356             :   }
    1357         280 : }
    1358             : 
    1359             : static GEN
    1360        6426 : F3xq_ellcardj(GEN a4, GEN a6, GEN T, GEN q, long n)
    1361             : {
    1362        6426 :   const ulong p = 3;
    1363             :   ulong t;
    1364        6426 :   GEN q1 = addiu(q,1);
    1365        6426 :   GEN na4 = Flx_neg(a4,p), ra4;
    1366        6426 :   if (!Flxq_issquare(na4,T,p))
    1367        3234 :     return q1;
    1368        3192 :   ra4 = Flxq_sqrt(na4,T,p);
    1369        3192 :   t = Flxq_trace(Flxq_div(a6,Flxq_mul(na4,ra4,T,p),T,p),T,p);
    1370        3192 :   if (n%2==1)
    1371             :   {
    1372             :     GEN q3;
    1373        1078 :     if (t==0) return q1;
    1374         301 :     q3 = powuu(p,(n+1)>>1);
    1375         301 :     return (t==1)^(n%4==1) ? subii(q1,q3): addii(q1,q3);
    1376             :   }
    1377             :   else
    1378             :   {
    1379        2114 :     GEN q22, q2 = powuu(p,n>>1);
    1380        2114 :     GEN W = Flxq_pow(a4,shifti(q,-2),T,p);
    1381        2114 :     long s = (W[2]==1)^(n%4==2);
    1382        2114 :     if (t!=0) return s ? addii(q1,q2): subii(q1, q2);
    1383        2114 :     q22 = shifti(q2,1);
    1384        2114 :     return s ? subii(q1,q22):  addii(q1, q22);
    1385             :   }
    1386             : }
    1387             : 
    1388             : static GEN
    1389       14812 : Flxq_ellcardj(GEN a4, GEN a6, ulong j, GEN T, GEN q, ulong p, long n)
    1390             : {
    1391       14812 :   GEN q1 = addiu(q,1);
    1392       14812 :   if (j==0)
    1393             :   {
    1394             :     ulong w;
    1395             :     GEN W, t, N;
    1396        5600 :     if (umodiu(q,6)!=1) return q1;
    1397        4200 :     N = Fp_ffellcard(gen_0,gen_1,q,n,utoipos(p));
    1398        4200 :     t = subii(q1, N);
    1399        4200 :     W = Flxq_pow(a6,diviuexact(shifti(q,-1), 3),T,p);
    1400        4200 :     if (degpol(W)>0) /*p=5 mod 6*/
    1401        1428 :       return Flx_equal1(Flxq_powu(W,3,T,p)) ? addii(q1,shifti(t,-1)):
    1402         476 :                                               subii(q1,shifti(t,-1));
    1403        3248 :     w = W[2];
    1404        3248 :     if (w==1)   return N;
    1405        2576 :     if (w==p-1) return addii(q1,t);
    1406             :     else /*p=1 mod 6*/
    1407             :     {
    1408        1904 :       GEN u = shifti(t,-1), v = sqrtint(diviuexact(subii(q,sqri(u)),3));
    1409        1904 :       GEN a = addii(u,v), b = shifti(v,1);
    1410        1904 :       if (Fl_powu(w,3,p)==1)
    1411             :       {
    1412         952 :         if (Fl_add(umodiu(a,p),Fl_mul(w,umodiu(b,p),p),p)==0)
    1413         511 :           return subii(q1,subii(shifti(b,1),a));
    1414             :         else
    1415         441 :           return addii(q1,addii(a,b));
    1416             :       }
    1417             :       else
    1418             :       {
    1419         952 :         if (Fl_sub(umodiu(a,p),Fl_mul(w,umodiu(b,p),p),p)==0)
    1420         511 :           return subii(q1,subii(a,shifti(b,1)));
    1421             :         else
    1422         441 :           return subii(q1,addii(a,b));
    1423             :       }
    1424             :     }
    1425        9212 :   } else if (j==1728%p)
    1426             :   {
    1427             :     ulong w;
    1428             :     GEN W, N, t;
    1429        5607 :     if (mod4(q)==3) return q1;
    1430        4207 :     W = Flxq_pow(a4,shifti(q,-2),T,p);
    1431        4207 :     if (degpol(W)>0) return q1; /*p=3 mod 4*/
    1432        3521 :     w = W[2];
    1433        3521 :     N = Fp_ffellcard(gen_1,gen_0,q,n,utoipos(p));
    1434        3521 :     if(w==1) return N;
    1435        2464 :     t = subii(q1, N);
    1436        2464 :     if(w==p-1) return addii(q1, t);
    1437             :     else /*p=1 mod 4*/
    1438             :     {
    1439        1400 :       GEN u = shifti(t,-1), v = sqrtint(subii(q,sqri(u)));
    1440        1400 :       if (Fl_add(umodiu(u,p),Fl_mul(w,umodiu(v,p),p),p)==0)
    1441         700 :         return subii(q1,shifti(v,1));
    1442             :       else
    1443         700 :         return addii(q1,shifti(v,1));
    1444             :     }
    1445             :   } else
    1446             :   {
    1447        3605 :     ulong g = Fl_div(j, Fl_sub(1728%p, j, p), p);
    1448        3605 :     GEN l = Flxq_div(Flx_triple(a6,p),Flx_double(a4,p),T,p);
    1449        3605 :     GEN N = Fp_ffellcard(utoi(Fl_triple(g,p)),utoi(Fl_double(g,p)),q,n,utoipos(p));
    1450        3605 :     if (Flxq_issquare(l,T,p)) return N;
    1451        2184 :     return subii(shifti(q1,1),N);
    1452             :   }
    1453             : }
    1454             : 
    1455             : GEN
    1456       50518 : Flxq_ellcard(GEN a4, GEN a6, GEN T, ulong p)
    1457             : {
    1458       50518 :   pari_sp av = avma;
    1459       50518 :   long n = get_Flx_degree(T);
    1460       50518 :   GEN J, r, q = powuu(p,  n);
    1461       50518 :   if (typ(a4)==t_VEC)
    1462       17087 :     r = F3xq_ellcard(gel(a4,1), a6, T);
    1463       33431 :   else if (p==3)
    1464        6426 :     r = F3xq_ellcardj(a4, a6, T, q, n);
    1465       27005 :   else if (degpol(a4)<=0 && degpol(a6)<=0)
    1466         217 :     r = Fp_ffellcard(utoi(Flx_eval(a4,0,p)),utoi(Flx_eval(a6,0,p)),q,n,utoipos(p));
    1467       26788 :   else if (degpol(J=Flxq_ellj(a4,a6,T,p))<=0)
    1468       14812 :     r = Flxq_ellcardj(a4,a6,lgpol(J)?J[2]:0,T,q,p,n);
    1469       11976 :   else if (p <= 7)
    1470       10843 :     r = Flxq_ellcard_Satoh(a4, a6, J, T, p);
    1471        1133 :   else if (cmpis(q,100)<0)
    1472           0 :     r = utoi(Flxq_ellcard_naive(a4, a6, T, p));
    1473        1133 :   else if (p == 13 || (7*p <= (ulong)10*n && (BITS_IN_LONG==64 || p <= 103)))
    1474          70 :     r = Flxq_ellcard_Satoh(a4, a6, J, T, p);
    1475        1063 :   else if (p <= (ulong)2*n)
    1476           0 :     r = Flxq_ellcard_Kedlaya(a4, a6, T, p);
    1477        1063 :   else if (expi(q)<=62)
    1478        1036 :     r = Flxq_ellcard_Shanks(a4, a6, q, T, p);
    1479             :   else
    1480          27 :     r = Fq_ellcard_SEA(Flx_to_ZX(a4),Flx_to_ZX(a6),q,Flx_to_ZX(T),utoipos(p),0);
    1481       50518 :   return gerepileuptoint(av, r);
    1482             : }

Generated by: LCOV version 1.13