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.18.0 lcov report (development 29644-f1a46412a2) Lines: 865 926 93.4 %
Date: 2024-10-12 09:06:38 Functions: 86 91 94.5 %
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             : /* These functions deal with point over elliptic curves over Fq defined
      29             :  * by an equation of the form y^2=x^3+a4*x+a6. Most of the time a6 is omitted
      30             :  * since it can be recovered from any point on the curve. */
      31             : 
      32             : GEN
      33       65804 : RgE_to_FlxqE(GEN x, GEN T, ulong p)
      34             : {
      35       65804 :   if (ell_is_inf(x)) return x;
      36       65804 :   retmkvec2(Rg_to_Flxq(gel(x,1),T,p), Rg_to_Flxq(gel(x,2),T,p));
      37             : }
      38             : 
      39             : GEN
      40      157822 : FlxqE_changepoint(GEN x, GEN ch, GEN T, ulong p)
      41             : {
      42      157822 :   pari_sp av = avma;
      43             :   GEN p1, p2, z, u, r, s, t, v, v2, v3;
      44             :   ulong pi;
      45      157822 :   if (ell_is_inf(x)) return x;
      46       93810 :   pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
      47       93810 :   u = gel(ch,1); r = gel(ch,2);
      48       93810 :   s = gel(ch,3); t = gel(ch,4);
      49       93810 :   v = Flxq_inv_pre(u, T, p, pi);
      50       93810 :   v2 = Flxq_sqr_pre(v, T, p, pi);
      51       93810 :   v3 = Flxq_mul_pre(v,v2, T, p, pi);
      52       93810 :   p1 = Flx_sub(gel(x,1), r, p);
      53       93810 :   p2 = Flx_sub(gel(x,2), Flx_add(Flxq_mul_pre(s, p1, T, p, pi),t, p), p);
      54       93810 :   z = cgetg(3,t_VEC);
      55       93810 :   gel(z,1) = Flxq_mul_pre(v2, p1, T, p, pi);
      56       93810 :   gel(z,2) = Flxq_mul_pre(v3, p2, T, p, pi);
      57       93810 :   return gerepileupto(av, z);
      58             : }
      59             : 
      60             : GEN
      61       65804 : FlxqE_changepointinv(GEN x, GEN ch, GEN T, ulong p)
      62             : {
      63       65804 :   pari_sp av = avma;
      64             :   GEN p1, p2, u, r, s, t, X, Y, u2, u3, u2X, z;
      65             :   ulong pi;
      66       65804 :   if (ell_is_inf(x)) return x;
      67       65804 :   pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
      68       65804 :   X = gel(x,1); Y = gel(x,2);
      69       65804 :   u = gel(ch,1); r = gel(ch,2);
      70       65804 :   s = gel(ch,3); t = gel(ch,4);
      71       65804 :   u2 = Flxq_sqr_pre(u, T, p, pi);
      72       65804 :   u3 = Flxq_mul_pre(u,u2, T, p, pi);
      73       65804 :   u2X = Flxq_mul_pre(u2,X, T, p, pi);
      74       65804 :   p1 = Flxq_mul_pre(u3,Y, T, p, pi);
      75       65804 :   p2 = Flx_add(Flxq_mul_pre(s,u2X, T, p, pi), t, p);
      76       65804 :   z = cgetg(3, t_VEC);
      77       65804 :   gel(z,1) = Flx_add(u2X, r, p);
      78       65804 :   gel(z,2) = Flx_add(p1, p2, p);
      79       65804 :   return gerepileupto(av, z);
      80             : }
      81             : 
      82             : static GEN
      83       22834 : nonsquare_Flxq(GEN T, ulong p)
      84             : {
      85       22834 :   pari_sp av = avma;
      86       22834 :   long n = degpol(T), vs = T[1];
      87             :   GEN a;
      88       22834 :   if (odd(n))
      89        7686 :     return mkvecsmall2(vs, nonsquare_Fl(p));
      90             :   do
      91             :   {
      92       30485 :     set_avma(av);
      93       30485 :     a = random_Flx(n, vs, p);
      94       30485 :   } while (Flxq_issquare(a, T, p));
      95       15148 :   return a;
      96             : }
      97             : 
      98             : void
      99       22834 : Flxq_elltwist(GEN a, GEN a6, GEN T, ulong p, GEN *pt_a, GEN *pt_a6)
     100             : {
     101       22834 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
     102       22834 :   GEN d = nonsquare_Flxq(T, p);
     103       22834 :   GEN d2 = Flxq_sqr_pre(d, T, p, pi), d3 = Flxq_mul_pre(d2, d, T, p, pi);
     104       22834 :   if (typ(a)==t_VECSMALL)
     105             :   {
     106       15232 :     *pt_a  = Flxq_mul_pre(a,  d2, T, p, pi);
     107       15232 :     *pt_a6 = Flxq_mul_pre(a6, d3, T, p, pi);
     108             :   } else
     109             :   {
     110        7602 :     *pt_a  = mkvec(Flxq_mul_pre(gel(a,1), d, T, p, pi));
     111        7602 :     *pt_a6 = Flxq_mul_pre(a6, d3, T, p, pi);
     112             :   }
     113       22834 : }
     114             : 
     115             : static GEN
     116     1594555 : FlxqE_dbl_slope(GEN P, GEN a4, GEN T, ulong p, ulong pi, GEN *ps)
     117             : {
     118             :   GEN x, y, Q, s;
     119     1594555 :   if (ell_is_inf(P) || !lgpol(gel(P,2))) return ellinf();
     120     1478012 :   x = gel(P,1); y = gel(P,2);
     121     1478012 :   if (p==3UL)
     122      531265 :     s = typ(a4)==t_VEC? Flxq_div_pre(Flxq_mul_pre(x, gel(a4,1), T,p,pi), y, T,p,pi)
     123      531265 :                       : Flxq_div_pre(a4, Flx_neg(y, p), T,p,pi);
     124             :   else
     125             :   {
     126      946747 :     GEN sx = Flx_add(Flx_triple(Flxq_sqr_pre(x, T, p, pi), p), a4, p);
     127      946747 :     s = Flxq_div_pre(sx, Flx_double(y, p), T, p, pi);
     128             :   }
     129     1478012 :   Q = cgetg(3,t_VEC);
     130     1478012 :   gel(Q,1) = Flx_sub(Flxq_sqr_pre(s, T, p, pi), Flx_double(x, p), p);
     131     1478012 :   if (typ(a4)==t_VEC) gel(Q, 1) = Flx_sub(gel(Q,1), gel(a4,1), p);
     132     1478012 :   gel(Q,2) = Flx_sub(Flxq_mul_pre(s, Flx_sub(x, gel(Q,1), p), T, p, pi),
     133             :                      y, p);
     134     1478012 :   if (ps) *ps = s;
     135     1478012 :   return Q;
     136             : }
     137             : 
     138             : GEN
     139           0 : FlxqE_dbl(GEN P, GEN a4, GEN T, ulong p)
     140             : {
     141           0 :   pari_sp av = avma;
     142           0 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
     143           0 :   return gerepileupto(av, FlxqE_dbl_slope(P,a4, T, p, pi, NULL));
     144             : }
     145             : 
     146             : static GEN
     147      713458 : FlxqE_add_slope(GEN P, GEN Q, GEN a4, GEN T, ulong p, ulong pi, GEN *ps)
     148             : {
     149             :   GEN Px, Py, Qx, Qy, R, s;
     150      713458 :   if (ell_is_inf(P)) return Q;
     151      708835 :   if (ell_is_inf(Q)) return P;
     152      708649 :   Px = gel(P,1); Py = gel(P,2);
     153      708649 :   Qx = gel(Q,1); Qy = gel(Q,2);
     154      708649 :   if (Flx_equal(Px, Qx))
     155             :   {
     156       53667 :     if (Flx_equal(Py, Qy))
     157        4325 :       return FlxqE_dbl_slope(P, a4, T, p, pi, ps);
     158             :     else
     159       49342 :       return ellinf();
     160             :   }
     161      654982 :   s = Flxq_div_pre(Flx_sub(Py, Qy, p), Flx_sub(Px, Qx, p), T, p, pi);
     162      654982 :   R = cgetg(3,t_VEC);
     163      654982 :   gel(R,1) = Flx_sub(Flx_sub(Flxq_sqr_pre(s, T, p, pi), Px, p), Qx, p);
     164      654982 :   if (typ(a4)==t_VEC) gel(R,1) = Flx_sub(gel(R,1), gel(a4,1), p);
     165      654982 :   gel(R,2) = Flx_sub(Flxq_mul_pre(s, Flx_sub(Px, gel(R,1), p), T, p, pi),
     166             :                      Py, p);
     167      654982 :   if (ps) *ps = s;
     168      654982 :   return R;
     169             : }
     170             : 
     171             : GEN
     172           0 : FlxqE_add(GEN P, GEN Q, GEN a4, GEN T, ulong p)
     173             : {
     174           0 :   pari_sp av = avma;
     175           0 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
     176           0 :   return gerepileupto(av, FlxqE_add_slope(P,Q,a4, T,p,pi, NULL));
     177             : }
     178             : 
     179             : static GEN
     180        6819 : FlxqE_neg_i(GEN P, ulong p)
     181             : {
     182        6819 :   if (ell_is_inf(P)) return P;
     183        6819 :   return mkvec2(gel(P,1), Flx_neg(gel(P,2), p));
     184             : }
     185             : 
     186             : GEN
     187        1193 : FlxqE_neg(GEN P, GEN T, ulong p)
     188             : {
     189             :   (void) T;
     190        1193 :   if (ell_is_inf(P)) return ellinf();
     191        1193 :   return mkvec2(gcopy(gel(P,1)), Flx_neg(gel(P,2), p));
     192             : }
     193             : 
     194             : GEN
     195           0 : FlxqE_sub(GEN P, GEN Q, GEN a4, GEN T, ulong p)
     196             : {
     197           0 :   pari_sp av = avma;
     198           0 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
     199           0 :   return gerepileupto(av, FlxqE_add_slope(P, FlxqE_neg_i(Q, p), a4, T,p,pi, NULL));
     200             : }
     201             : 
     202             : struct _FlxqE
     203             : {
     204             :   GEN a4, a6, T;
     205             :   ulong p, pi;
     206             : };
     207             : 
     208             : static GEN
     209     1566373 : _FlxqE_dbl(void *E, GEN P)
     210             : {
     211     1566373 :   struct _FlxqE *e = (struct _FlxqE *) E;
     212     1566373 :   return FlxqE_dbl_slope(P, e->a4, e->T, e->p, e->pi, NULL);
     213             : }
     214             : 
     215             : static GEN
     216      703790 : _FlxqE_add(void *E, GEN P, GEN Q)
     217             : {
     218      703790 :   struct _FlxqE *e = (struct _FlxqE *) E;
     219      703790 :   return FlxqE_add_slope(P, Q, e->a4, e->T, e->p, e->pi, NULL);
     220             : }
     221             : 
     222             : static GEN
     223        6819 : _FlxqE_sub(void *E, GEN P, GEN Q)
     224             : {
     225        6819 :   struct _FlxqE *e = (struct _FlxqE *) E;
     226        6819 :   return FlxqE_add_slope(P, FlxqE_neg_i(Q,e->p), e->a4, e->T,e->p,e->pi, NULL);
     227             : }
     228             : 
     229             : static GEN
     230      264696 : _FlxqE_mul(void *E, GEN P, GEN n)
     231             : {
     232      264696 :   pari_sp av = avma;
     233      264696 :   struct _FlxqE *e=(struct _FlxqE *) E;
     234      264696 :   long s = signe(n);
     235      264696 :   if (!s || ell_is_inf(P)) return ellinf();
     236      263985 :   if (s < 0) P = FlxqE_neg(P, e->T, e->p);
     237      263985 :   if (is_pm1(n)) return s>0? gcopy(P): P;
     238      251134 :   return gerepilecopy(av, gen_pow_i(P, n, e, &_FlxqE_dbl, &_FlxqE_add));
     239             : }
     240             : 
     241             : GEN
     242       64068 : FlxqE_mul(GEN P, GEN n, GEN a4, GEN T, ulong p)
     243             : {
     244             :   struct _FlxqE E;
     245       64068 :   E.a4= a4; E.T = T; E.p = p; E.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
     246       64068 :   return _FlxqE_mul(&E, P, n);
     247             : }
     248             : 
     249             : /* 3*x^2+2*a2*x = -a2*x, and a2!=0 */
     250             : 
     251             : /* Finds a random nonsingular point on E */
     252             : static GEN
     253       76839 : random_F3xqE(GEN a2, GEN a6, GEN T)
     254             : {
     255       76839 :   pari_sp ltop = avma;
     256             :   GEN x, y, rhs;
     257       76839 :   const ulong p = 3;
     258             :   do
     259             :   {
     260      151221 :     set_avma(ltop);
     261      151221 :     x   = random_Flx(get_Flx_degree(T),get_Flx_var(T),p);
     262      151221 :     rhs = Flx_add(Flxq_mul(Flxq_sqr(x, T, p), Flx_add(x, a2, p), T, p), a6, p);
     263      151221 :   } while ((!lgpol(rhs) && !lgpol(x)) || !Flxq_issquare(rhs, T, p));
     264       76839 :   y = Flxq_sqrt(rhs, T, p);
     265       76839 :   if (!y) pari_err_PRIME("random_F3xqE", T);
     266       76839 :   return gerepilecopy(ltop, mkvec2(x, y));
     267             : }
     268             : 
     269             : /* Finds a random nonsingular point on E */
     270             : static GEN
     271      150676 : random_FlxqE_pre(GEN a4, GEN a6, GEN T, ulong p, ulong pi)
     272             : {
     273      150676 :   pari_sp ltop = avma;
     274             :   GEN x, x2, y, rhs;
     275      150676 :   if (typ(a4)==t_VEC) return random_F3xqE(gel(a4,1), a6, T);
     276             :   do
     277             :   {
     278      142709 :     set_avma(ltop);
     279      142709 :     x   = random_Flx(get_Flx_degree(T),get_Flx_var(T),p);
     280      142709 :     x2  = Flxq_sqr_pre(x, T, p, pi); /*  x^3+a4*x+a6 = x*(x^2+a4)+a6  */
     281      142709 :     rhs = Flx_add(Flxq_mul_pre(x, Flx_add(x2, a4, p), T, p, pi), a6, p);
     282      148064 :   } while ((!lgpol(rhs) && !lgpol(Flx_add(Flx_triple(x2, p), a4, p)))
     283      148022 :           || !Flxq_issquare(rhs, T, p));
     284       73837 :   y = Flxq_sqrt(rhs, T, p);
     285       73837 :   if (!y) pari_err_PRIME("random_FlxqE", T);
     286       73837 :   return gerepilecopy(ltop, mkvec2(x, y));
     287             : }
     288             : GEN
     289       76941 : random_FlxqE(GEN a4, GEN a6, GEN T, ulong p)
     290       76941 : { return random_FlxqE_pre(a4, a6, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p)); }
     291             : 
     292             : static GEN
     293       73735 : _FlxqE_rand(void *E)
     294             : {
     295       73735 :   struct _FlxqE *e=(struct _FlxqE *) E;
     296       73735 :   return random_FlxqE_pre(e->a4, e->a6, e->T, e->p, e->pi);
     297             : }
     298             : 
     299             : static const struct bb_group FlxqE_group={_FlxqE_add,_FlxqE_mul,_FlxqE_rand,hash_GEN,zvV_equal,ell_is_inf, NULL};
     300             : 
     301             : const struct bb_group *
     302          61 : get_FlxqE_group(void ** pt_E, GEN a4, GEN a6, GEN T, ulong p)
     303             : {
     304          61 :   struct _FlxqE *e = (struct _FlxqE *) stack_malloc(sizeof(struct _FlxqE));
     305          61 :   e->a4 = a4; e->a6 = a6;
     306          61 :   e->pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
     307          61 :   e->p = p;
     308          61 :   e->T = Flx_get_red_pre(T, p, e->pi);
     309          61 :   *pt_E = (void *)e; return &FlxqE_group;
     310             : }
     311             : 
     312             : GEN
     313        1470 : FlxqE_order(GEN z, GEN o, GEN a4, GEN T, ulong p)
     314             : {
     315        1470 :   pari_sp av = avma;
     316             :   struct _FlxqE e;
     317        1470 :   e.a4 = a4; e.T = T; e.p = p; e.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
     318        1470 :   return gerepileuptoint(av, gen_order(z, o, (void*)&e, &FlxqE_group));
     319             : }
     320             : 
     321             : GEN
     322          49 : FlxqE_log(GEN a, GEN b, GEN o, GEN a4, GEN T, ulong p)
     323             : {
     324          49 :   pari_sp av = avma;
     325             :   struct _FlxqE e;
     326          49 :   e.a4 = a4; e.T = T; e.p = p; e.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
     327          49 :   return gerepileuptoint(av, gen_PH_log(a, b, o, (void*)&e, &FlxqE_group));
     328             : }
     329             : 
     330             : /***********************************************************************/
     331             : /**                            Pairings                               **/
     332             : /***********************************************************************/
     333             : /* Derived from APIP by Jerome Milan, 2012 */
     334             : static GEN
     335       65208 : FlxqE_vert(GEN P, GEN Q, GEN a4, GEN T, ulong p, ulong pi)
     336             : {
     337       65208 :   long vT = get_Flx_var(T);
     338             :   GEN df;
     339       65208 :   if (ell_is_inf(P)) return pol1_Flx(vT);
     340       44064 :   if (!Flx_equal(gel(Q,1), gel(P,1))) return Flx_sub(gel(Q,1), gel(P,1), p);
     341        1844 :   if (lgpol(gel(P,2))!=0) return pol1_Flx(vT);
     342         798 :   df = typ(a4)==t_VEC ? Flxq_mul_pre(gel(P,1), Flx_double(gel(a4,1), p), T,p,pi)
     343        1246 :                       : a4;
     344        1246 :   return Flxq_inv_pre(Flx_add(Flx_triple(Flxq_sqr_pre(gel(P,1), T,p, pi), p),
     345             :                               df, p), T, p, pi);
     346             : }
     347             : 
     348             : static GEN
     349       26706 : FlxqE_Miller_line(GEN R, GEN Q, GEN slope, GEN a4, GEN T, ulong p, ulong pi)
     350             : {
     351       26706 :   long vT = get_Flx_var(T);
     352       26706 :   GEN x = gel(Q,1), y = gel(Q,2);
     353       26706 :   GEN tmp1 = Flx_sub(x, gel(R,1), p);
     354       26706 :   GEN tmp2 = Flx_add(Flxq_mul_pre(tmp1, slope, T, p, pi), gel(R,2), p);
     355       26706 :   if (!Flx_equal(y, tmp2)) return Flx_sub(y, tmp2, p);
     356        1228 :   if (lgpol(y) == 0) return pol1_Flx(vT);
     357             :   else
     358             :   {
     359         605 :     GEN s1, s2, a2 = typ(a4)==t_VEC ? gel(a4,1): NULL;
     360         605 :     GEN y2i = Flxq_inv_pre(Flx_mulu(y, 2, p), T, p, pi);
     361         605 :     GEN df = a2 ? Flxq_mul_pre(x, Flx_mulu(a2, 2, p), T, p, pi): a4;
     362             :     GEN x3, ddf;
     363         605 :     s1 = Flxq_mul_pre(Flx_add(Flx_triple(Flxq_sqr_pre(x, T, p, pi), p), df, p), y2i, T, p, pi);
     364         605 :     if (!Flx_equal(s1, slope)) return Flx_sub(s1, slope, p);
     365         244 :     x3 = Flx_triple(x, p);
     366         244 :     ddf = a2 ? Flx_add(x3, a2, p): x3;
     367         244 :     s2 = Flxq_mul_pre(Flx_sub(ddf, Flxq_sqr_pre(s1, T,p,pi), p), y2i, T,p,pi);
     368         244 :     return lgpol(s2)!=0 ? s2: y2i;
     369             :   }
     370             : }
     371             : 
     372             : /* Computes the equation of the line tangent to R and returns its
     373             :  * evaluation at the point Q. Also doubles the point R. */
     374             : static GEN
     375       43748 : FlxqE_tangent_update(GEN R, GEN Q, GEN a4, GEN T, ulong p, ulong pi, GEN *pt_R)
     376             : {
     377       43748 :   if (ell_is_inf(R))
     378             :   {
     379        3940 :     *pt_R = ellinf();
     380        3940 :     return pol1_Flx(get_Flx_var(T));
     381             :   }
     382       39808 :   else if (!lgpol(gel(R,2)))
     383             :   {
     384       15951 :     *pt_R = ellinf();
     385       15951 :     return FlxqE_vert(R, Q, a4, T, p, pi);
     386             :   } else {
     387             :     GEN slope;
     388       23857 :     *pt_R = FlxqE_dbl_slope(R, a4, T, p, pi, &slope);
     389       23857 :     return FlxqE_Miller_line(R, Q, slope, a4, T, p, pi);
     390             :   }
     391             : }
     392             : 
     393             : /* Computes the equation of the line through R and P, and returns its
     394             :  * evaluation at the point Q. Also adds P to the point R. */
     395             : static GEN
     396        4179 : FlxqE_chord_update(GEN R, GEN P, GEN Q, GEN a4, GEN T, ulong p, ulong pi, GEN *pt_R)
     397             : {
     398        4179 :   if (ell_is_inf(R))
     399             :   {
     400          77 :     *pt_R = gcopy(P);
     401          77 :     return FlxqE_vert(P, Q, a4, T, p, pi);
     402             :   }
     403        4102 :   else if (ell_is_inf(P))
     404             :   {
     405           0 :     *pt_R = gcopy(R);
     406           0 :     return FlxqE_vert(R, Q, a4, T, p, pi);
     407             :   }
     408        4102 :   else if (Flx_equal(gel(P, 1), gel(R, 1)))
     409             :   {
     410        1253 :     if (Flx_equal(gel(P, 2), gel(R, 2)))
     411           0 :       return FlxqE_tangent_update(R, Q, a4, T, p, pi, pt_R);
     412             :     else
     413             :     {
     414        1253 :       *pt_R = ellinf();
     415        1253 :       return FlxqE_vert(R, Q, a4, T, p, pi);
     416             :     }
     417             :   } else {
     418             :     GEN slope;
     419        2849 :     *pt_R = FlxqE_add_slope(P, R, a4, T, p, pi, &slope);
     420        2849 :     return FlxqE_Miller_line(R, Q, slope, a4, T, p, pi);
     421             :   }
     422             : }
     423             : 
     424             : struct _FlxqE_miller
     425             : {
     426             :   ulong p, pi;
     427             :   GEN T, a4, P;
     428             : };
     429             : 
     430             : static GEN
     431       43748 : FlxqE_Miller_dbl(void* E, GEN d)
     432             : {
     433       43748 :   struct _FlxqE_miller *m = (struct _FlxqE_miller *)E;
     434       43748 :   ulong p = m->p, pi = m->pi;
     435       43748 :   GEN T = m->T, a4 = m->a4, P = m->P;
     436       43748 :   GEN v, line, point = gel(d,3);
     437       43748 :   GEN N = Flxq_sqr_pre(gel(d,1), T, p, pi);
     438       43748 :   GEN D = Flxq_sqr_pre(gel(d,2), T, p, pi);
     439       43748 :   line = FlxqE_tangent_update(point, P, a4, T, p, pi, &point);
     440       43748 :   N  = Flxq_mul_pre(N, line, T, p, pi);
     441       43748 :   v = FlxqE_vert(point, P, a4, T, p, pi);
     442       43748 :   D = Flxq_mul_pre(D, v, T, p, pi); return mkvec3(N, D, point);
     443             : }
     444             : 
     445             : static GEN
     446        4179 : FlxqE_Miller_add(void* E, GEN va, GEN vb)
     447             : {
     448        4179 :   struct _FlxqE_miller *m = (struct _FlxqE_miller *)E;
     449        4179 :   ulong p = m->p, pi = m->pi;
     450        4179 :   GEN T = m->T, a4 = m->a4, P = m->P;
     451             :   GEN v, line, point;
     452        4179 :   GEN na = gel(va,1), da = gel(va,2), pa = gel(va,3);
     453        4179 :   GEN nb = gel(vb,1), db = gel(vb,2), pb = gel(vb,3);
     454        4179 :   GEN N = Flxq_mul_pre(na, nb, T, p, pi);
     455        4179 :   GEN D = Flxq_mul_pre(da, db, T, p, pi);
     456        4179 :   line = FlxqE_chord_update(pa, pb, P, a4, T, p, pi, &point);
     457        4179 :   N  = Flxq_mul_pre(N, line, T, p, pi);
     458        4179 :   v = FlxqE_vert(point, P, a4, T, p, pi);
     459        4179 :   D = Flxq_mul_pre(D, v, T, p, pi); return mkvec3(N, D, point);
     460             : }
     461             : 
     462             : /* Returns the Miller function f_{m, Q} evaluated at the point P using
     463             :  * the standard Miller algorithm. */
     464             : static GEN
     465       17127 : FlxqE_Miller(GEN Q, GEN P, GEN m, GEN a4, GEN T, ulong p, ulong pi)
     466             : {
     467       17127 :   pari_sp av = avma;
     468             :   struct _FlxqE_miller d;
     469             :   GEN v, N, D, g1;
     470             : 
     471       17127 :   d.a4 = a4; d.T = T; d.p = p; d.P = P; d.pi = pi;
     472       17127 :   g1 = pol1_Flx(get_Flx_var(T));
     473       17127 :   v = gen_pow_i(mkvec3(g1,g1,Q), m, (void*)&d,
     474             :                 FlxqE_Miller_dbl, FlxqE_Miller_add);
     475       17127 :   N = gel(v,1); D = gel(v,2);
     476       17127 :   return gerepileupto(av, Flxq_div_pre(N, D, T, p, pi));
     477             : }
     478             : 
     479             : GEN
     480       13840 : FlxqE_weilpairing_pre(GEN P, GEN Q, GEN m, GEN a4, GEN T, ulong p, ulong pi)
     481             : {
     482       13840 :   pari_sp av = avma;
     483             :   GEN N, D, w;
     484       13840 :   if (ell_is_inf(P) || ell_is_inf(Q)
     485       11073 :     || (Flx_equal(gel(P,1),gel(Q,1)) && Flx_equal(gel(P,2),gel(Q,2))))
     486        5308 :     return pol1_Flx(get_Flx_var(T));
     487        8532 :   N = FlxqE_Miller(P, Q, m, a4, T, p, pi);
     488        8532 :   D = FlxqE_Miller(Q, P, m, a4, T, p, pi);
     489        8532 :   w = Flxq_div_pre(N, D, T, p, pi); if (mpodd(m)) w = Flx_neg(w, p);
     490        8532 :   return gerepileupto(av, w);
     491             : }
     492             : GEN
     493          21 : FlxqE_weilpairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, ulong p)
     494          21 : { return FlxqE_weilpairing_pre(P,Q,m,a4,T,p, SMALL_ULONG(p)?0:get_Fl_red(p)); }
     495             : 
     496             : GEN
     497          63 : FlxqE_tatepairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, ulong p)
     498             : {
     499          63 :   if (ell_is_inf(P) || ell_is_inf(Q)) return pol1_Flx(get_Flx_var(T));
     500          63 :   return FlxqE_Miller(P, Q, m, a4, T, p, SMALL_ULONG(p)? 0: get_Fl_red(p));
     501             : }
     502             : 
     503             : static GEN
     504       13819 : _FlxqE_pairorder(void *E, GEN P, GEN Q, GEN m, GEN F)
     505             : {
     506       13819 :   struct _FlxqE *e = (struct _FlxqE *) E;
     507       13819 :   return  Flxq_order(FlxqE_weilpairing_pre(P,Q,m,e->a4,e->T,e->p,e->pi), F, e->T, e->p);
     508             : }
     509             : 
     510             : GEN
     511       22294 : Flxq_ellgroup(GEN a4, GEN a6, GEN N, GEN T, ulong p, GEN *pt_m)
     512             : {
     513             :   struct _FlxqE e;
     514       22294 :   GEN q = powuu(p, get_Flx_degree(T));
     515       22294 :   e.a4=a4; e.a6=a6; e.T=T; e.p=p; e.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
     516       22294 :   return gen_ellgroup(N, subiu(q,1), pt_m, (void*)&e, &FlxqE_group, _FlxqE_pairorder);
     517             : }
     518             : 
     519             : GEN
     520       14370 : Flxq_ellgens(GEN a4, GEN a6, GEN ch, GEN D, GEN m, GEN T, ulong p)
     521             : {
     522             :   GEN P;
     523       14370 :   pari_sp av = avma;
     524             :   struct _FlxqE e;
     525       14370 :   e.a4=a4; e.a6=a6; e.T=T; e.p=p; e.pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
     526       14370 :   switch(lg(D)-1)
     527             :   {
     528          63 :   case 0:
     529          63 :     return cgetg(1,t_VEC);
     530       11801 :   case 1:
     531       11801 :     P = gen_gener(gel(D,1), (void*)&e, &FlxqE_group);
     532       11801 :     P = mkvec(FlxqE_changepoint(P, ch, T, p));
     533       11801 :     break;
     534        2506 :   default:
     535        2506 :     P = gen_ellgens(gel(D,1), gel(D,2), m, (void*)&e, &FlxqE_group, _FlxqE_pairorder);
     536        2506 :     gel(P,1) = FlxqE_changepoint(gel(P,1), ch, T, p);
     537        2506 :     gel(P,2) = FlxqE_changepoint(gel(P,2), ch, T, p);
     538        2506 :     break;
     539             :   }
     540       14307 :   return gerepilecopy(av, P);
     541             : }
     542             : /***********************************************************************/
     543             : /**                          Point counting                           **/
     544             : /***********************************************************************/
     545             : 
     546             : /* assume a and n are coprime */
     547             : static GEN
     548       76251 : RgX_circular_shallow(GEN P, long a, long n)
     549             : {
     550       76251 :   long i, l = lgpol(P);
     551       76251 :   GEN Q = cgetg(2+n,t_POL);
     552       76251 :   Q[1] = P[1];
     553      512330 :   for(i=0; i<l; i++)
     554      436079 :     gel(Q,2+(i*a)%n) = gel(P,2+i);
     555      168693 :   for(   ; i<n; i++)
     556       92442 :     gel(Q,2+(i*a)%n) = gen_0;
     557       76251 :   return normalizepol_lg(Q,2+n);
     558             : }
     559             : 
     560             : static GEN
     561       76251 : ZpXQ_frob_cyc(GEN x, GEN T, GEN q, ulong p)
     562             : {
     563       76251 :   long n = get_FpX_degree(T);
     564       76251 :   return FpX_rem(RgX_circular_shallow(x,p,n+1), T, q);
     565             : }
     566             : 
     567             : static GEN
     568      113610 : ZpXQ_frob(GEN x, GEN Xm, GEN T, GEN q, ulong p)
     569             : {
     570      113610 :   if (lg(Xm)==1)
     571       43428 :     return ZpXQ_frob_cyc(x, T, q, p);
     572             :   else
     573             :   {
     574       70182 :     long n = get_FpX_degree(T);
     575       70182 :     GEN V = RgX_blocks(RgX_inflate(x, p), n, p);
     576       70182 :     GEN W = ZXV_dotproduct(V, Xm);
     577       70182 :     return FpX_rem(W, T, q);
     578             :   }
     579             : }
     580             : 
     581             : struct _lift_lin
     582             : {
     583             :   ulong p, pi;
     584             :   GEN sqx, Tp, ai, Xm;
     585             : };
     586             : 
     587             : static GEN
     588       84077 : _lift_invl(void *E, GEN x)
     589             : {
     590       84077 :   struct _lift_lin *d = (struct _lift_lin *) E;
     591       84077 :   GEN T = d->Tp;
     592       84077 :   ulong p = d->p, pi = d->pi;
     593       84077 :   GEN xai = Flxq_mul_pre(ZX_to_Flx(x, p), d->ai, T, p, pi);
     594       84077 :   return Flx_to_ZX(Flxq_lroot_fast_pre(xai, d->sqx, T, p, pi));
     595             : }
     596             : static GEN
     597       23744 : _lift_lin(void *E, GEN F, GEN x2, GEN q)
     598             : {
     599       23744 :   struct _lift_lin *d = (struct _lift_lin *) E;
     600       23744 :   pari_sp av = avma;
     601       23744 :   GEN T = gel(F,3), Xm = gel(F,4);
     602       23744 :   GEN y2  = ZpXQ_frob(x2, Xm, T, q, d->p);
     603       23744 :   GEN lin = FpX_add(ZX_mul(gel(F,1), y2), ZX_mul(gel(F,2), x2), q);
     604       23744 :   return gerepileupto(av, FpX_rem(lin, T, q));
     605             : }
     606             : 
     607             : static GEN
     608      181041 : FpM_FpXV_bilinear(GEN P, GEN X, GEN Y, GEN p)
     609             : {
     610      181041 :    pari_sp av = avma;
     611      181041 :    GEN s =  ZX_mul(FpXV_FpC_mul(X,gel(P,1),p),gel(Y,1));
     612      181041 :    long i, l = lg(P);
     613      851487 :    for(i=2; i<l; i++)
     614      670446 :      s = ZX_add(s, ZX_mul(FpXV_FpC_mul(X,gel(P,i),p),gel(Y,i)));
     615      181041 :    return gerepileupto(av, FpX_red(s, p));
     616             : }
     617             : 
     618             : static GEN
     619      181041 : FpM_FpXQV_bilinear(GEN P, GEN X, GEN Y, GEN T, GEN p)
     620      181041 : { return FpX_rem(FpM_FpXV_bilinear(P,X,Y,p),T,p); }
     621             : 
     622             : static GEN
     623      120708 : FpXC_powderiv(GEN M, GEN p)
     624             : {
     625             :   long i, l;
     626      120708 :   long v = varn(gel(M,2));
     627      120708 :   GEN m = cgetg_copy(M, &l);
     628      120708 :   gel(m,1) = pol_0(v);
     629      120708 :   gel(m,2) = pol_1(v);
     630      447132 :   for(i=2; i<l-1; i++)
     631      326424 :     gel(m,i+1) = FpX_Fp_mul(gel(M,i),utoi(i), p);
     632      120708 :   return m;
     633             : }
     634             : 
     635             : struct _lift_iso
     636             : {
     637             :   GEN phi, Xm, T, sqx, Tp;
     638             :   ulong p, pi;
     639             : };
     640             : 
     641             : static GEN
     642       60333 : _lift_iter(void *E, GEN x2, GEN q)
     643             : {
     644       60333 :   struct _lift_iso *d = (struct _lift_iso *) E;
     645       60333 :   ulong p = d->p;
     646       60333 :   long n = lg(d->phi)-2;
     647       60333 :   GEN TN = FpXT_red(d->T, q), XN = FpXV_red(d->Xm, q);
     648       60333 :   GEN y2 = ZpXQ_frob(x2, XN, TN, q, p);
     649       60333 :   GEN xp = FpXQ_powers(x2, n, TN, q);
     650       60333 :   GEN yp = FpXQ_powers(y2, n, TN, q);
     651       60333 :   GEN V  = FpM_FpXQV_bilinear(d->phi,xp,yp,TN,q);
     652       60333 :   return mkvec3(V,xp,yp);
     653             : }
     654             : 
     655             : static GEN
     656       60333 : _lift_invd(void *E, GEN V, GEN v, GEN qM, long M)
     657             : {
     658       60333 :   struct _lift_iso *d = (struct _lift_iso *) E;
     659             :   struct _lift_lin e;
     660       60333 :   ulong p = d->p, pi = d->pi;
     661       60333 :   GEN TM = FpXT_red(d->T, qM), XM = FpXV_red(d->Xm, qM);
     662       60333 :   GEN xp = FpXV_red(gel(v,2), qM);
     663       60333 :   GEN yp = FpXV_red(gel(v,3), qM);
     664       60333 :   GEN Dx = FpM_FpXQV_bilinear(d->phi, FpXC_powderiv(xp, qM), yp, TM, qM);
     665       60333 :   GEN Dy = FpM_FpXQV_bilinear(d->phi, xp, FpXC_powderiv(yp, qM), TM, qM);
     666       60333 :   GEN F = mkvec4(Dy, Dx, TM, XM);
     667       60333 :   e.ai = Flxq_inv_pre(ZX_to_Flx(Dy,p),d->Tp, p, pi);
     668       60333 :   e.sqx = d->sqx; e.Tp = d->Tp; e.p=p; e.pi=pi; e.Xm = XM;
     669       60333 :   return gen_ZpX_Dixon(F,V,qM,utoipos(p),M,(void*) &e, _lift_lin, _lift_invl);
     670             : }
     671             : 
     672             : static GEN
     673       25053 : lift_isogeny(GEN phi, GEN x0, long n, GEN Xm, GEN T, GEN sqx, GEN Tp,
     674             :   ulong p, ulong pi)
     675             : {
     676             :   struct _lift_iso d;
     677       25053 :   d.phi = phi; d.Xm = Xm; d.T = T;
     678       25053 :   d.sqx = sqx; d.Tp = Tp; d.p = p; d.pi = pi;
     679       25053 :   return gen_ZpX_Newton(x0, utoipos(p), n,(void*)&d, _lift_iter, _lift_invd);
     680             : }
     681             : 
     682             : static GEN
     683       25032 : getc2(GEN act, GEN X, GEN T, GEN q, ulong p, long N)
     684             : {
     685       25032 :   GEN A1 = RgV_to_RgX(gel(act,1),0), A2 =  RgV_to_RgX(gel(act,2),0);
     686       25032 :   long n = brent_kung_optpow(maxss(degpol(A1),degpol(A2)),2,1);
     687       25032 :   GEN xp = FpXQ_powers(X,n,T,q);
     688       25032 :   GEN P  = FpX_FpXQV_eval(A1, xp, T, q);
     689       25032 :   GEN Q  = FpX_FpXQV_eval(A2, xp, T, q);
     690       25032 :   return ZpXQ_div(P, Q, T, q, utoipos(p), N);
     691             : }
     692             : 
     693             : struct _ZpXQ_norm
     694             : {
     695             :   long n;
     696             :   GEN T, p;
     697             : };
     698             : 
     699             : static GEN
     700       32823 : ZpXQ_norm_mul(void *E, GEN x, GEN y)
     701             : {
     702       32823 :   struct _ZpXQ_norm *D = (struct _ZpXQ_norm*)E;
     703       32823 :   GEN P = gel(x,1), Q = gel(y,1);
     704       32823 :   long a = mael(x,2,1), b = mael(y,2,1);
     705       32823 :   retmkvec2(FpXQ_mul(P,ZpXQ_frob_cyc(Q, D->T, D->p, a), D->T, D->p),
     706             :             mkvecsmall((a*b)%D->n));
     707             : }
     708             : static GEN
     709       22715 : ZpXQ_norm_sqr(void *E, GEN x) { return ZpXQ_norm_mul(E, x, x); }
     710             : 
     711             : /* Assume T = Phi_(n) and n prime */
     712             : GEN
     713       11340 : ZpXQ_norm_pcyc(GEN x, GEN T, GEN q, GEN p)
     714             : {
     715             :   GEN z;
     716             :   struct _ZpXQ_norm D;
     717       11340 :   long d = get_FpX_degree(T);
     718       11340 :   D.T = T; D.p = q; D.n = d+1;
     719       11340 :   if (d==1) return ZX_copy(x);
     720       11340 :   z = mkvec2(x,mkvecsmall(p[2]));
     721       11340 :   z = gen_powu_i(z,d,(void*)&D,ZpXQ_norm_sqr,ZpXQ_norm_mul);
     722       11340 :   return gmael(z,1,2);
     723             : }
     724             : 
     725             : /* Assume T = Phi_(n) and n prime */
     726             : static GEN
     727       11102 : ZpXQ_sqrtnorm_pcyc(GEN x, GEN T, GEN q, GEN p, long e)
     728             : {
     729       11102 :   GEN z = ZpXQ_norm_pcyc(x, T, q, p);
     730       11102 :   return Zp_sqrtlift(z,Fp_sqrt(z,p),p,e);
     731             : }
     732             : 
     733             : /* Assume a = 1 [p], return the square root of the norm */
     734             : static GEN
     735       13951 : ZpXQ_sqrtnorm(GEN a, GEN T, GEN q, GEN p, long e)
     736             : {
     737       13951 :   GEN s = Fp_div(FpXQ_trace(ZpXQ_log(a, T, p, e), T, q), gen_2, q);
     738       13951 :   return modii(gel(Qp_exp(cvtop(s, p, e-1)),4), q);
     739             : }
     740             : 
     741             : struct _teich_lin
     742             : {
     743             :   ulong p, pi;
     744             :   GEN sqx, Tp;
     745             :   long m;
     746             : };
     747             : 
     748             : static GEN
     749       29512 : _teich_invl(void *E, GEN x)
     750             : {
     751       29512 :   struct _teich_lin *d = (struct _teich_lin *) E;
     752       29512 :   ulong p = d->p, pi = d->pi;
     753       29512 :   return Flx_to_ZX(Flxq_lroot_fast_pre(ZX_to_Flx(x,p), d->sqx, d->Tp, p, pi));
     754             : }
     755             : 
     756             : static GEN
     757        8953 : _teich_lin(void *E, GEN F, GEN x2, GEN q)
     758             : {
     759        8953 :   struct _teich_lin *d = (struct _teich_lin *) E;
     760        8953 :   pari_sp av = avma;
     761        8953 :   GEN T = gel(F,2), Xm = gel(F,3);
     762        8953 :   GEN y2  = ZpXQ_frob(x2, Xm, T, q, d->p);
     763        8953 :   GEN lin = FpX_sub(y2, ZX_mulu(ZX_mul(gel(F,1), x2), d->p), q);
     764        8953 :   return gerepileupto(av, FpX_rem(lin, T, q));
     765             : }
     766             : 
     767             : struct _teich_iso
     768             : {
     769             :   GEN Xm, T, sqx, Tp;
     770             :   ulong p, pi;
     771             : };
     772             : 
     773             : static GEN
     774       20559 : _teich_iter(void *E, GEN x2, GEN q)
     775             : {
     776       20559 :   struct _teich_iso *d = (struct _teich_iso *) E;
     777       20559 :   ulong p = d->p;
     778       20559 :   GEN TN = FpXT_red(d->T, q), XN = FpXV_red(d->Xm, q);
     779       20559 :   GEN y2 = ZpXQ_frob(x2, XN, TN, q, d->p);
     780       20559 :   GEN x1 = FpXQ_powu(x2, p-1, TN, q);
     781       20559 :   GEN xp = FpXQ_mul(x2, x1, TN, q);
     782       20559 :   GEN V = FpX_sub(y2,xp,q);
     783       20559 :   return mkvec2(V,x1);
     784             : }
     785             : 
     786             : static GEN
     787       20559 : _teich_invd(void *E, GEN V, GEN v, GEN qM, long M)
     788             : {
     789       20559 :   struct _teich_iso *d = (struct _teich_iso *) E;
     790             :   struct _teich_lin e;
     791       20559 :   ulong p = d->p;
     792       20559 :   GEN TM = FpXT_red(d->T, qM), XM = FpXV_red(d->Xm, qM);
     793       20559 :   GEN x1 = FpX_red(gel(v,2), qM);
     794       20559 :   GEN F = mkvec3(x1, TM, XM);
     795       20559 :   e.sqx = d->sqx; e.Tp = d->Tp; e.p = p; e.pi = d->pi;
     796       20559 :   return gen_ZpX_Dixon(F,V,qM,utoipos(p),M,(void*) &e, _teich_lin, _teich_invl);
     797             : }
     798             : 
     799             : static GEN
     800       10234 : Teichmuller_lift(GEN x, GEN Xm, GEN T, GEN sqx, GEN Tp, ulong p, ulong pi,
     801             :   long N)
     802             : {
     803             :   struct _teich_iso d;
     804       10234 :   d.Xm = Xm; d.T = T; d.sqx = sqx; d.Tp = Tp; d.p = p; d.pi = pi;
     805       10234 :   return gen_ZpX_Newton(x,utoipos(p), N,(void*)&d, _teich_iter, _teich_invd);
     806             : }
     807             : 
     808             : static GEN
     809       25053 : get_norm(GEN a4, GEN a6, GEN T, ulong p, ulong pi, long N)
     810             : {
     811       25053 :   long sv=T[1];
     812             :   GEN a;
     813       25053 :   if (p==3) a = gel(a4,1);
     814             :   else
     815             :   {
     816       10248 :     GEN P = mkpoln(4, pol1_Flx(sv), pol0_Flx(sv), a4, a6);
     817       10248 :     a = gel(FlxqX_powu_pre(P, p>>1, T,p,pi), 2+p-1);
     818             :   }
     819       25053 :   return Zp_sqrtnlift(gen_1,subss(p,1),utoi(Flxq_norm(a,T,p)),utoipos(p), N);
     820             : }
     821             : 
     822             : static GEN
     823       25032 : fill_pols(long n, const long *v, long m, const long *vn,
     824             :           const long *vd, GEN *act)
     825             : {
     826             :   long i, j;
     827       25032 :   long d = upowuu(n,12/(n-1));
     828       25032 :   GEN N, D, M = zeromatcopy(n+1,n+1);
     829       25032 :   gmael(M,1,n+1) = gen_1;
     830      120764 :   for (i = 2; i <= n+1; i++)
     831      339486 :     for (j = i-1; j <= n; j++) gmael(M,i,j) = mulis(powuu(d,i-2), v[j-i+1]);
     832       25032 :   N = cgetg(m+1,t_COL);
     833       25032 :   D = cgetg(m+1,t_COL);
     834      135541 :   for(i = 1; i <= m; i++)
     835             :   {
     836      110509 :     gel(N,i) = stoi(*vn++);
     837      110509 :     gel(D,i) = stoi(*vd++);
     838             :   }
     839       25032 :   *act = mkmat2(N,D); return M;
     840             : }
     841             : 
     842             : /*
     843             :   These polynomials were extracted from the ECHIDNA databases
     844             :   available at <http://echidna.maths.usyd.edu.au/echidna/>
     845             :   and computed by David R. Kohel.
     846             :   Return the matrix of the modular polynomial, set act to the parametrization,
     847             :   and set dj to the opposite of the supersingular j-invariant.
     848             : */
     849             : static GEN
     850       25032 : get_Kohel_polynomials(ulong p, GEN *act, long *dj)
     851             : {
     852       25032 :   const long mat3[] = {-1,-36,-270};
     853       25032 :   const long num3[] = {1,-483,-21141,-59049};
     854       25032 :   const long den3[] = {1,261, 4347, -6561};
     855       25032 :   const long mat5[] = {-1,-30,-315,-1300,-1575};
     856       25032 :   const long num5[] = {-1,490,20620,158750,78125};
     857       25032 :   const long den5[] = {-1,-254,-4124,-12250,3125};
     858       25032 :   const long mat7[] = {-1,-28,-322,-1904,-5915,-8624,-4018};
     859       25032 :   const long num7[] = {1,-485,-24058,-343833,-2021642,-4353013,-823543};
     860       25032 :   const long den7[] = {1,259,5894,49119,168406,166355,-16807};
     861       25032 :   const long mat13[]= {-1,-26,-325,-2548,-13832,-54340,-157118,-333580,-509366,
     862             :                        -534820,-354536,-124852,-15145};
     863       25032 :   const long num13[]= {1,-487,-24056,-391463,-3396483,-18047328,-61622301,
     864             :                        -133245853,-168395656,-95422301,-4826809};
     865       25032 :   const long den13[]= {1,257,5896,60649,364629,1388256,3396483,5089019,4065464,
     866             :                        1069939,-28561};
     867       25032 :   switch(p)
     868             :   {
     869       14805 :   case 3:
     870       14805 :     *dj = 0;
     871       14805 :     return fill_pols(3,mat3,4,num3,den3,act);
     872       10178 :   case 5:
     873       10178 :     *dj = 0;
     874       10178 :     return fill_pols(5,mat5,5,num5,den5,act);
     875          35 :   case 7:
     876          35 :     *dj = 1;
     877          35 :     return fill_pols(7,mat7,7,num7,den7,act);
     878          14 :   case 13:
     879          14 :     *dj = 8;
     880          14 :     return fill_pols(13,mat13,11,num13,den13,act);
     881             :   }
     882             :   *dj=0; *act = NULL; return NULL; /* LCOV_EXCL_LINE */
     883             : }
     884             : 
     885             : long
     886       32465 : zx_is_pcyc(GEN T)
     887             : {
     888       32465 :   long i, n = degpol(T);
     889       32465 :   if (!uisprime(n+1)) return 0;
     890       99148 :   for (i = 0; i <= n; i++)
     891       87808 :     if (T[i+2] != 1UL) return 0;
     892       11340 :   return 1;
     893             : }
     894             : 
     895             : static GEN
     896       25032 : Flxq_ellcard_Kohel(GEN a4, GEN a6, GEN T, ulong p)
     897             : {
     898       25032 :   pari_sp av = avma, av2;
     899             :   pari_timer ti;
     900       25032 :   long n = get_Flx_degree(T), N = (n+4)/2, dj;
     901       25032 :   GEN q = powuu(p, N);
     902             :   GEN T2, Xm, s1, c2, t, lr, S1, sqx, Nc2, Np;
     903       25032 :   GEN act, phi = get_Kohel_polynomials(p, &act, &dj);
     904       25032 :   long ispcyc = zx_is_pcyc(get_Flx_mod(T));
     905       25032 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
     906       25032 :   timer_start(&ti);
     907       25032 :   if (!ispcyc)
     908             :   {
     909       13937 :     T2 = Flx_Teichmuller(get_Flx_mod(T),p,N);
     910       13937 :     if (DEBUGLEVEL) timer_printf(&ti,"Teich");
     911             :   } else
     912       11095 :     T2 = Flx_to_ZX(get_Flx_mod(T));
     913             : 
     914       25032 :   T2 = FpX_get_red(T2, q); T = ZXT_to_FlxT(T2, p);
     915       25032 :   av2 = avma;
     916       25032 :   if (DEBUGLEVEL) timer_printf(&ti,"Barrett");
     917       25032 :   if (!ispcyc)
     918             :   {
     919       13937 :     Xm = FpXQ_powers(pol_xn(n,get_FpX_var(T2)),p-1,T2,q);
     920       13937 :     if (DEBUGLEVEL) timer_printf(&ti,"Xm");
     921             :   } else
     922       11095 :     Xm = cgetg(1,t_VEC);
     923       25032 :   s1 = Flxq_inv_pre(Flx_Fl_add(Flxq_ellj(a4,a6,T,p),dj, p),T,p,pi);
     924       25032 :   lr = Flxq_lroot_pre(polx_Flx(get_Flx_var(T)), T,p,pi);
     925       25032 :   sqx = Flxq_powers_pre(lr, p-1, T, p, pi);
     926       25032 :   S1 = lift_isogeny(phi, Flx_to_ZX(s1), N, Xm, T2, sqx, T,p,pi);
     927       25032 :   if (DEBUGLEVEL) timer_printf(&ti,"Lift isogeny");
     928       25032 :   c2 = getc2(act, S1, T2, q, p, N);
     929       25032 :   if (DEBUGLEVEL) timer_printf(&ti,"c^2");
     930       25032 :   if (p>3 && !ispcyc)
     931             :   {
     932       10220 :     GEN c2p = Flx_to_ZX(Flxq_inv_pre(ZX_to_Flx(c2,p),T,p,pi));
     933       10220 :     GEN tc2 = Teichmuller_lift(c2p,Xm, T2,sqx,T,p,pi,N);
     934       10220 :     if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fq");
     935       10220 :     c2 = FpX_rem(FpX_mul(tc2,c2,q),T2,q);
     936             :   }
     937       25032 :   c2 = gerepileupto(av2, c2);
     938       25032 :   if (DEBUGLEVEL) timer_printf(&ti,"tc2");
     939       25032 :   Nc2 = (ispcyc? ZpXQ_sqrtnorm_pcyc: ZpXQ_sqrtnorm)(c2, T2, q, utoipos(p), N);
     940       25032 :   if (DEBUGLEVEL) timer_printf(&ti,"Norm");
     941       25032 :   Np = get_norm(a4,a6,T,p,pi,N);
     942       25032 :   if (p>3 && ispcyc)
     943             :   {
     944           7 :     GEN Ncpi =  utoi(Fl_inv(umodiu(Nc2,p), p));
     945           7 :     GEN tNc2 = Zp_sqrtnlift(gen_1, subss(p,1), Ncpi, utoipos(p),N);
     946           7 :     if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fp");
     947           7 :     Nc2 = Fp_mul(Nc2,tNc2,q);
     948             :   }
     949       25032 :   t = Fp_center_i(Fp_mul(Nc2,Np,q),q,shifti(q,-1));
     950       25032 :   return gerepileupto(av, subii(addiu(powuu(p,n),1),t));
     951             : }
     952             : 
     953             : /* Use Damien Robert's method */
     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             : /* in p^2, so p is tiny */
     970             : static GEN
     971          21 : Flxq_ellcard_Harley(GEN a4, GEN a6, GEN T, ulong p)
     972             : {
     973          21 :   pari_sp av = avma, av2;
     974             :   pari_timer ti;
     975          21 :   long n = get_Flx_degree(T), N = (n+5)/2;
     976          21 :   GEN pp = utoipos(p), q = powuu(p, N);
     977             :   GEN T2, j, t, phi, J1, sqx, Xm, c2, tc2, c2p, Nc2, Np;
     978          21 :   long ispcyc = zx_is_pcyc(get_Flx_mod(T));
     979          21 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p); /* = 0 here */
     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_pre(Flxq_lroot_pre(polx_Flx(T[1]), T,p,pi), p-1, T,p,pi);
     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,pi);
    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_pre(ZX_to_Flx(c2,p),T,p,pi));
    1008          14 :     tc2 = Teichmuller_lift(c2p,Xm, T2,sqx,T,p,pi,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,pi,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             : /*                          Shanks-Mestre                                  */
    1030             : /***************************************************************************/
    1031             : 
    1032             : /* Return the lift of a (mod b), which is closest to h */
    1033             : static GEN
    1034       10317 : closest_lift(GEN a, GEN b, GEN h)
    1035       10317 : { return addii(a, mulii(b, diviiround(subii(h,a), b))); }
    1036             : 
    1037             : /* find multiple of order of f using Baby Step/Giant Step, f^h close to 1,
    1038             :  * order lies in an interval of size <= 'bound' and known mod B */
    1039             : static GEN
    1040        5543 : _FlxqE_order_multiple(void *E, GEN f, GEN h, GEN bound, GEN B)
    1041             : {
    1042        5543 :   pari_sp av = avma, av1;
    1043             :   pari_timer Ti;
    1044        5543 :   long i, s = ceilsqrtdiv(bound, B) >> 1;
    1045             :   GEN P, F, tx, ti, fg, fh;
    1046             : 
    1047        5543 :   P = fh = _FlxqE_mul(E, f, h);
    1048        5543 :   if (DEBUGLEVEL >= 6) timer_start(&Ti);
    1049        5543 :   if (ell_is_inf(fh)) return h;
    1050        5151 :   F = _FlxqE_mul(E, f, B);
    1051        5151 :   if (s < 3)
    1052             :   { /* we're nearly done: naive search */
    1053        1098 :     GEN Q = P;
    1054        1098 :     for (i=1;; i++)
    1055             :     {
    1056        3641 :       P = _FlxqE_add(E, P, F); /* h.f + i.F */
    1057        3641 :       if (ell_is_inf(P)) return gerepileupto(av, addii(h, mului(i,B)));
    1058        3094 :       Q = _FlxqE_sub(E, Q, F); /* h.f - i.F */
    1059        3094 :       if (ell_is_inf(Q)) return gerepileupto(av, subii(h, mului(i,B)));
    1060             :     }
    1061             :   }
    1062        4053 :   tx = cgetg(s+1,t_VECSMALL); av1 = avma;
    1063       43300 :   for (i=1; i<=s; i++)
    1064             :   { /* baby steps */
    1065       39575 :     tx[i] = hash_GEN(gel(P, 1));
    1066       39575 :     P = _FlxqE_add(E, P, F); /* h.f + i.F */
    1067       39575 :     if (ell_is_inf(P)) return gerepileupto(av, addii(h, mului(i,B)));
    1068       39247 :     if (gc_needed(av1,3))
    1069             :     {
    1070           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"[Flxq_ellcard] baby steps, i=%ld",i);
    1071           0 :       P = gerepileupto(av1,P);
    1072             :     }
    1073             :   }
    1074        3725 :   if (DEBUGLEVEL >= 6) timer_printf(&Ti,"[Flxq_ellcard] baby steps, s = %ld",s);
    1075             :   /* giant steps: fg = s.F */
    1076        3725 :   fg = gerepileupto(av1, _FlxqE_sub(E, P, fh));
    1077        3725 :   if (ell_is_inf(fg)) return gerepileupto(av, mului(s,B));
    1078        3725 :   ti = vecsmall_indexsort(tx); /* = permutation sorting tx */
    1079        3725 :   tx = perm_mul(tx,ti);
    1080        3725 :   if (DEBUGLEVEL >= 6) timer_printf(&Ti, "[Flxq_ellcard] sorting");
    1081        3725 :   av1 = avma;
    1082        3725 :   for (P=fg, i=1; ; i++)
    1083       35479 :   {
    1084       39204 :     long k = hash_GEN(gel(P,1)), r = zv_search(tx, k);
    1085       39204 :     if (r)
    1086             :     {
    1087        7453 :       while (r && tx[r] == k) r--;
    1088        3725 :       for (r++; r <= s && tx[r] == k; r++)
    1089             :       {
    1090        3725 :         long j = ti[r]-1;
    1091        3725 :         GEN Q = _FlxqE_add(E, _FlxqE_mul(E, F, stoi(j)), fh);
    1092        3725 :         if (DEBUGLEVEL >= 6)
    1093           0 :           timer_printf(&Ti, "[Flxq_ellcard] giant steps, i = %ld",i);
    1094        3725 :         if (Flx_equal(gel(P,1), gel(Q,1)))
    1095             :         {
    1096        3725 :           if (Flx_equal(gel(P,2), gel(Q,2))) i = -i;
    1097        3725 :           return gerepileupto(av,addii(h, mulii(addis(mulss(s,i), j), B)));
    1098             :         }
    1099             :       }
    1100             :     }
    1101       35479 :     P = _FlxqE_add(E, P, fg);
    1102       35479 :     if (gc_needed(av1,3))
    1103             :     {
    1104           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"[Flxq_ellcard] giants steps, i=%ld",i);
    1105           0 :       P = gerepileupto(av1,P);
    1106             :     }
    1107             :   }
    1108             : }
    1109             : static GEN
    1110        5543 : _FlxqE_order(void *E, GEN f, GEN h, GEN bound, GEN B)
    1111             : {
    1112        5543 :   GEN o = _FlxqE_order_multiple(E, f, h, bound, B);
    1113        5543 :   return gen_order(f, o, E, &FlxqE_group);
    1114             : }
    1115             : 
    1116             : static void
    1117       64190 : Flx_next(GEN t, ulong p)
    1118             : {
    1119             :   long i;
    1120       78141 :   for(i=2;;i++)
    1121       78141 :     if (uel(t,i)==p-1) t[i]=0; else { t[i]++; break; }
    1122       64190 : }
    1123             : 
    1124             : static void
    1125       64190 : Flx_renormalize_ip(GEN x, long lx)
    1126             : {
    1127             :   long i;
    1128       78141 :   for (i = lx-1; i>=2; i--)
    1129       71057 :     if (x[i]) break;
    1130       64190 :   setlg(x, i+1);
    1131       64190 : }
    1132             : 
    1133             : static ulong
    1134        6188 : F3xq_ellcard_naive(GEN a2, GEN a6, GEN T)
    1135             : {
    1136        6188 :   pari_sp av = avma;
    1137        6188 :   long i, d = get_Flx_degree(T), lx = d+2;
    1138        6188 :   long q = upowuu(3, d), a;
    1139        6188 :   GEN x = zero_zv(lx); x[1] = get_Flx_var(T);
    1140       26978 :   for(a=1, i=0; i<q; i++)
    1141             :   {
    1142             :     GEN rhs;
    1143       20790 :     Flx_renormalize_ip(x, lx);
    1144       20790 :     rhs = Flx_add(Flxq_mul(Flxq_sqr(x, T, 3), Flx_add(x, a2, 3), T, 3), a6, 3);
    1145       20790 :     if (!lgpol(rhs)) a++; else if (Flxq_issquare(rhs, T, 3)) a+=2;
    1146       20790 :     Flx_next(x, 3);
    1147             :   }
    1148        6188 :   set_avma(av); return a;
    1149             : }
    1150             : 
    1151             : /* p^deg(T) is tiny */
    1152             : static ulong
    1153         896 : Flxq_ellcard_naive(GEN a4, GEN a6, GEN T, ulong p)
    1154             : {
    1155         896 :   pari_sp av = avma;
    1156         896 :   long i, d = get_Flx_degree(T), lx = d+2;
    1157         896 :   long q = upowuu(p, d), a;
    1158         896 :   GEN x = zero_zv(lx); x[1] = get_Flx_var(T);
    1159       44296 :   for(a = 1, i = 0; i < q; i++)
    1160             :   {
    1161             :     GEN x2, rhs;
    1162       43400 :     Flx_renormalize_ip(x, lx);
    1163       43400 :     x2  = Flxq_sqr_pre(x, T, p, 0);
    1164       43400 :     rhs = Flx_add(Flxq_mul_pre(x, Flx_add(x2, a4, p), T, p, 0), a6, p);
    1165       43400 :     if (!lgpol(rhs)) a++; else if (Flxq_issquare(rhs,T,p)) a += 2;
    1166       43400 :     Flx_next(x,p);
    1167             :   }
    1168         896 :   set_avma(av); return a;
    1169             : }
    1170             : 
    1171             : static long
    1172       11395 : Flxq_kronecker(GEN x, GEN T, ulong p)
    1173             : {
    1174             :   pari_sp av;
    1175       11395 :   if (lgpol(x) == 0) return 0;
    1176       11370 :   av = avma; return gc_long(av, krouu(Flxq_norm(x, T, p), p));
    1177             : }
    1178             : 
    1179             : /* Find x such that kronecker(u = x^3+a4x+a6, p) is KRO.
    1180             :  * Return point [x*u,u^2] on E (KRO=1) / E^twist (KRO=-1) */
    1181             : static GEN
    1182        5543 : Flxq_ellpoint(long KRO, GEN a4, GEN a6, GEN T, ulong p, ulong pi)
    1183             : {
    1184        5543 :   long v = get_Flx_var(T), n = get_Flx_degree(T);
    1185             :   for(;;)
    1186        5852 :   {
    1187       11395 :     GEN x = random_Flx(n, v, p), x2 = Flxq_sqr_pre(x,T,p,pi);
    1188       11395 :     GEN u = Flx_add(a6, Flxq_mul_pre(Flx_add(a4, x2, p), x, T,p, pi), p);
    1189       11395 :     if (Flxq_kronecker(u,T,p) == KRO)
    1190        5543 :       return mkvec2(Flxq_mul_pre(u,x, T,p,pi), Flxq_sqr_pre(u, T,p,pi));
    1191             :   }
    1192             : }
    1193             : 
    1194             : static GEN
    1195        4774 : Flxq_ellcard_Shanks(GEN a4, GEN a6, GEN q, GEN T, ulong p)
    1196             : {
    1197        4774 :   pari_sp av = avma;
    1198        4774 :   ulong pi = SMALL_ULONG(p)? 0: get_Fl_red(p);
    1199        4774 :   long v = get_Flx_var(T), KRO = -1;
    1200             :   GEN h,f, A, B;
    1201        4774 :   GEN q1p = addiu(q,1), q2p = shifti(q1p, 1);
    1202        4774 :   GEN bound = addiu(sqrti(gmul2n(q,4)), 1); /* ceil( 4sqrt(q) ) */
    1203             :   struct _FlxqE e;
    1204        4774 :   e.p = p; e.pi = pi; e.T = Flx_get_red_pre(T, p, pi);
    1205             :   /* once #E(Flxq) is known mod B >= bound, it is determined */
    1206        4774 :   switch(FlxqX_nbroots(mkpoln(4, pol1_Flx(v), pol0_Flx(v), a4, a6), T, p))
    1207             :   { /* how many 2-torsion points ? */
    1208        2891 :   case 3:  A = gen_0; B = utoipos(4); break;
    1209        1414 :   case 1:  A = gen_0; B = gen_2; break;
    1210         469 :   default: A = gen_1; B = gen_2; break; /* 0 */
    1211             :   }
    1212             :   for(;;)
    1213             :   {
    1214        5543 :     h = closest_lift(A, B, q1p);
    1215             :     /* [ux, u^2] is on E_u: y^2 = x^3 + c4 u^2 x + c6 u^3
    1216             :      * E_u isomorphic to E (resp. E') iff KRO = 1 (resp. -1)
    1217             :      * #E(F_p) = p+1 - a_p, #E'(F_p) = p+1 + a_p
    1218             :      *
    1219             :      * #E_u(Flxq) = A (mod B),  h is close to #E_u(Flxq) */
    1220        5543 :     KRO = -KRO;
    1221        5543 :     f = Flxq_ellpoint(KRO, a4,a6, T,p,pi);
    1222        5543 :     e.a4 = Flxq_mul_pre(a4, gel(f,2), T,p,pi); /* a4 for E_u */
    1223        5543 :     h = _FlxqE_order((void*)&e, f, h, bound, B);
    1224             :     /* h | #E_u(Flxq) = A (mod B) */
    1225        5543 :     A = Z_chinese_all(A, gen_0, B, h, &B);
    1226        5543 :     if (cmpii(B, bound) >= 0) break;
    1227             :     /* not done, update A mod B for the _next_ curve, isomorphic to
    1228             :      * the quadratic twist of this one */
    1229         769 :     A = remii(subii(q2p,A), B); /* #E(Fq)+#E'(Fq) = 2q+2 */
    1230             :   }
    1231        4774 :   h = closest_lift(A, B, q1p);
    1232        4774 :   return gerepileuptoint(av, KRO == 1? h: subii(q2p,h));
    1233             : }
    1234             : 
    1235             : static GEN
    1236       20993 : F3xq_ellcard(GEN a2, GEN a6, GEN T)
    1237             : {
    1238       20993 :   long n = get_Flx_degree(T);
    1239       20993 :   if (n <= 2)
    1240        5887 :     return utoi(F3xq_ellcard_naive(a2, a6, T));
    1241             :   else
    1242             :   {
    1243       15106 :     GEN q1 = addiu(powuu(3, get_Flx_degree(T)), 1), t;
    1244       15106 :     GEN a = Flxq_div(a6,Flxq_powu(a2,3,T,3),T,3);
    1245       15106 :     if (Flx_equal1(Flxq_powu(a, 8, T, 3)))
    1246             :     {
    1247         301 :       GEN P = Flxq_minpoly(a,T,3);
    1248         301 :       long dP = degpol(P); /* dP <= 2 */
    1249         301 :       ulong q = upowuu(3,dP);
    1250         301 :       GEN A2 = pol1_Flx(P[1]), A6 = Flx_rem(polx_Flx(P[1]), P, 3);
    1251         301 :       long tP = q + 1 - F3xq_ellcard_naive(A2, A6, P);
    1252         301 :       t = elltrace_extension(stoi(tP), n/dP, utoi(q));
    1253         301 :       if (umodiu(t, 3)!=1) t = negi(t);
    1254         301 :       return Flx_equal1(a2) || Flxq_issquare(a2,T,3) ? subii(q1,t): addii(q1,t);
    1255             :     }
    1256       14805 :     else return Flxq_ellcard_Kohel(mkvec(a2), a6, T, 3);
    1257             :   }
    1258             : }
    1259             : 
    1260             : static GEN
    1261       11144 : Flxq_ellcard_Satoh(GEN a4, GEN a6, GEN j, GEN T, ulong p)
    1262             : {
    1263       11144 :   long n = get_Flx_degree(T);
    1264       11144 :   if (n <= 2)
    1265         896 :     return utoi(Flxq_ellcard_naive(a4, a6, T, p));
    1266             :   else
    1267             :   {
    1268       10248 :     GEN jp = Flxq_powu(j, p, T, p);
    1269       10248 :     GEN s = Flx_add(j, jp, p);
    1270       10248 :     if (degpol(s) <= 0)
    1271             :     { /* it is assumed j not in F_p */
    1272           0 :       GEN m = Flxq_mul(j, jp, T, p);
    1273           0 :       if (degpol(m) <= 0)
    1274             :       {
    1275           0 :         GEN q = sqru(p);
    1276           0 :         GEN q1 = addiu(powuu(p, get_Flx_degree(T)), 1);
    1277           0 :         GEN sk = Flx_Fl_add(Flx_neg(j, p), 1728%p, p);
    1278           0 :         GEN sA4 = Flx_triple(Flxq_mul(sk, j, T, p), p);
    1279           0 :         GEN u = Flxq_div(a4, sA4, T, p);
    1280           0 :         ulong ns = lgpol(s) ? Fl_neg(s[2], p): 0UL;
    1281           0 :         GEN P = mkvecsmall4(T[1], m[2], ns, 1L);
    1282             :         GEN A4, A6, t, tP;
    1283           0 :         Flxq_ellj_to_a4a6(polx_Flx(T[1]), P, p, &A4, &A6);
    1284           0 :         tP = addis(q, 1 - Flxq_ellcard_naive(A4, A6, P, p));
    1285           0 :         t = elltrace_extension(tP, n>>1, q);
    1286           0 :         return Flxq_is2npower(u, 2, T, p) ? subii(q1,t): addii(q1,t);
    1287             :       }
    1288             :     }
    1289       10248 :     if (p<=7 || p==13) return Flxq_ellcard_Kohel(a4, a6, T, p);
    1290          21 :     else return Flxq_ellcard_Harley(a4, a6, T, p);
    1291             :   }
    1292             : }
    1293             : 
    1294             : static GEN
    1295           0 : Flxq_ellcard_Kedlaya(GEN a4, GEN a6, GEN T, ulong p)
    1296             : {
    1297           0 :   pari_sp av = avma;
    1298           0 :   GEN H = mkpoln(4, gen_1, gen_0, Flx_to_ZX(a4), Flx_to_ZX(a6));
    1299           0 :   GEN Tp = Flx_to_ZX(get_Flx_mod(T));
    1300           0 :   long n = degpol(Tp), e = ((p < 16 ? n+1: n)>>1)+1;
    1301           0 :   GEN M = ZlXQX_hyperellpadicfrobenius(H, Tp, p, e);
    1302           0 :   GEN N = ZpXQM_prodFrobenius(M, Tp, utoipos(p), e);
    1303           0 :   GEN q = powuu(p, e);
    1304           0 :   GEN tp = Fq_add(gcoeff(N,1,1), gcoeff(N,2,2), Tp, q);
    1305           0 :   GEN t = Fp_center_i(typ(tp)==t_INT ? tp: leading_coeff(tp), q, shifti(q,-1));
    1306           0 :   return gerepileupto(av, subii(addiu(powuu(p, n), 1), t));
    1307             : }
    1308             : 
    1309             : GEN
    1310       57299 : Flxq_ellj(GEN a4, GEN a6, GEN T, ulong p)
    1311             : {
    1312       57299 :   pari_sp av=avma;
    1313       57299 :   if (p==3)
    1314             :   {
    1315             :     GEN J;
    1316       14805 :     if (typ(a4)!=t_VEC) return pol0_Flx(get_Flx_var(T));
    1317       14805 :     J = Flxq_div(Flxq_powu(gel(a4,1),3, T, p),Flx_neg(a6,p), T, p);
    1318       14805 :     return gerepileuptoleaf(av, J);
    1319             :   }
    1320             :   else
    1321             :   {
    1322       42494 :     pari_sp av=avma;
    1323       42494 :     GEN a43 = Flxq_mul(a4,Flxq_sqr(a4,T,p),T,p);
    1324       42494 :     GEN a62 = Flxq_sqr(a6,T,p);
    1325       42494 :     GEN num = Flx_mulu(a43,6912,p);
    1326       42494 :     GEN den = Flx_add(Flx_mulu(a43,4,p),Flx_mulu(a62,27,p),p);
    1327       42494 :     return gerepileuptoleaf(av, Flxq_div(num, den, T, p));
    1328             :   }
    1329             : }
    1330             : 
    1331             : void
    1332           0 : Flxq_ellj_to_a4a6(GEN j, GEN T, ulong p, GEN *pt_a4, GEN *pt_a6)
    1333             : {
    1334           0 :   ulong zagier = 1728 % p;
    1335           0 :   if (lgpol(j)==0)
    1336           0 :     { *pt_a4 = pol0_Flx(T[1]); *pt_a6 =pol1_Flx(T[1]); }
    1337           0 :   else if (lgpol(j)==1 && uel(j,2) == zagier)
    1338           0 :     { *pt_a4 = pol1_Flx(T[1]); *pt_a6 =pol0_Flx(T[1]); }
    1339             :   else
    1340             :   {
    1341           0 :     GEN k = Flx_Fl_add(Flx_neg(j, p), zagier, p);
    1342           0 :     GEN kj = Flxq_mul(k, j, T, p);
    1343           0 :     GEN k2j = Flxq_mul(kj, k, T, p);
    1344           0 :     *pt_a4 = Flx_triple(kj, p);
    1345           0 :     *pt_a6 = Flx_double(k2j, p);
    1346             :   }
    1347           0 : }
    1348             : 
    1349             : static GEN
    1350        9205 : F3xq_ellcardj(GEN a4, GEN a6, GEN T, GEN q, long n)
    1351             : {
    1352        9205 :   const ulong p = 3;
    1353             :   ulong t;
    1354        9205 :   GEN q1 = addiu(q,1);
    1355        9205 :   GEN na4 = Flx_neg(a4,p), ra4;
    1356        9205 :   if (!Flxq_issquare(na4,T,p))
    1357        4823 :     return q1;
    1358        4382 :   ra4 = Flxq_sqrt(na4,T,p);
    1359        4382 :   t = Flxq_trace(Flxq_div(a6,Flxq_mul(na4,ra4,T,p),T,p),T,p);
    1360        4382 :   if (n%2==1)
    1361             :   {
    1362             :     GEN q3;
    1363        2261 :     if (t==0) return q1;
    1364         812 :     q3 = powuu(p,(n+1)>>1);
    1365         812 :     return (t==1)^(n%4==1) ? subii(q1,q3): addii(q1,q3);
    1366             :   }
    1367             :   else
    1368             :   {
    1369        2121 :     GEN q22, q2 = powuu(p,n>>1);
    1370        2121 :     GEN W = Flxq_pow(a4,shifti(q,-2),T,p);
    1371        2121 :     long s = (W[2]==1)^(n%4==2);
    1372        2121 :     if (t!=0) return s ? addii(q1,q2): subii(q1, q2);
    1373        2121 :     q22 = shifti(q2,1);
    1374        2121 :     return s ? subii(q1,q22):  addii(q1, q22);
    1375             :   }
    1376             : }
    1377             : 
    1378             : static GEN
    1379       15820 : Flxq_ellcardj(GEN a4, GEN a6, ulong j, GEN T, GEN q, ulong p, long n)
    1380             : {
    1381       15820 :   GEN q1 = addiu(q,1);
    1382       15820 :   if (j==0)
    1383             :   {
    1384             :     ulong w;
    1385             :     GEN W, t, N;
    1386        5677 :     if (umodiu(q,6)!=1) return q1;
    1387        4277 :     N = Fp_ffellcard(gen_0,gen_1,q,n,utoipos(p));
    1388        4277 :     t = subii(q1, N);
    1389        4277 :     W = Flxq_pow(a6,diviuexact(shifti(q,-1), 3),T,p);
    1390        4277 :     if (degpol(W)>0) /*p=5 mod 6*/
    1391        1323 :       return Flx_equal1(Flxq_powu(W,3,T,p)) ? addii(q1,shifti(t,-1)):
    1392         441 :                                               subii(q1,shifti(t,-1));
    1393        3395 :     w = W[2];
    1394        3395 :     if (w==1)   return N;
    1395        2653 :     if (w==p-1) return addii(q1,t);
    1396             :     else /*p=1 mod 6*/
    1397             :     {
    1398        1862 :       GEN u = shifti(t,-1), v = sqrtint(diviuexact(subii(q,sqri(u)),3));
    1399        1862 :       GEN a = addii(u,v), b = shifti(v,1);
    1400        1862 :       if (Fl_powu(w,3,p)==1)
    1401             :       {
    1402         931 :         if (Fl_add(umodiu(a,p),Fl_mul(w,umodiu(b,p),p),p)==0)
    1403         469 :           return subii(q1,subii(shifti(b,1),a));
    1404             :         else
    1405         462 :           return addii(q1,addii(a,b));
    1406             :       }
    1407             :       else
    1408             :       {
    1409         931 :         if (Fl_sub(umodiu(a,p),Fl_mul(w,umodiu(b,p),p),p)==0)
    1410         469 :           return subii(q1,subii(a,shifti(b,1)));
    1411             :         else
    1412         462 :           return subii(q1,addii(a,b));
    1413             :       }
    1414             :     }
    1415       10143 :   } else if (j==1728%p)
    1416             :   {
    1417             :     ulong w;
    1418             :     GEN W, N, t;
    1419        5656 :     if (mod4(q)==3) return q1;
    1420        4256 :     W = Flxq_pow(a4,shifti(q,-2),T,p);
    1421        4256 :     if (degpol(W)>0) return q1; /*p=3 mod 4*/
    1422        3612 :     w = W[2];
    1423        3612 :     N = Fp_ffellcard(gen_1,gen_0,q,n,utoipos(p));
    1424        3612 :     if(w==1) return N;
    1425        2534 :     t = subii(q1, N);
    1426        2534 :     if(w==p-1) return addii(q1, t);
    1427             :     else /*p=1 mod 4*/
    1428             :     {
    1429        1400 :       GEN u = shifti(t,-1), v = sqrtint(subii(q,sqri(u)));
    1430        1400 :       if (Fl_add(umodiu(u,p),Fl_mul(w,umodiu(v,p),p),p)==0)
    1431         700 :         return subii(q1,shifti(v,1));
    1432             :       else
    1433         700 :         return addii(q1,shifti(v,1));
    1434             :     }
    1435             :   } else
    1436             :   {
    1437        4487 :     ulong g = Fl_div(j, Fl_sub(1728%p, j, p), p);
    1438        4487 :     GEN N = Fp_ffellcard(utoi(Fl_triple(g,p)),utoi(Fl_double(g,p)),q,n,utoipos(p));
    1439        4487 :     GEN l = Flxq_mul(Flx_triple(a6,p),Flx_double(a4,p),T,p);
    1440        4487 :     if (Flxq_issquare(l,T,p)) return N;
    1441        2940 :     return subii(shifti(q1,1),N);
    1442             :   }
    1443             : }
    1444             : 
    1445             : static GEN
    1446         454 : Flxq_ffellcard(GEN a4, GEN a6, GEN M, GEN q, GEN T, ulong p, long n)
    1447             : {
    1448         454 :   long m = degpol(M);
    1449         454 :   GEN j = polx_Flx(M[1]);
    1450         454 :   GEN g = Flxq_div(j, mkvecsmall3(M[1],1728%p,p-1), M, p);
    1451         454 :   GEN N = Flxq_ellcard(Flx_triple(g, p), Flx_double(g, p), M, p);
    1452         454 :   GEN qm =  powuu(p, m), q1 = addiu(q, 1), qm1 = addiu(qm, 1);
    1453         454 :   GEN l = Flxq_mul(Flx_triple(a6,p), Flx_double(a4,p), T, p);
    1454         454 :   GEN te = elltrace_extension(subii(qm1, N), n/m, qm);
    1455         454 :   return Flxq_issquare(l,T,p) ? subii(q1, te): addii(q1, te);
    1456             : }
    1457             : 
    1458             : static GEN
    1459       62661 : Flxq_ellcard_i(GEN a4, GEN a6, GEN T, ulong p)
    1460             : {
    1461       62661 :   long n = get_Flx_degree(T);
    1462       62661 :   GEN J, M, q = powuu(p,  n);
    1463       62661 :   if (typ(a4)==t_VEC)
    1464       20993 :     return F3xq_ellcard(gel(a4,1), a6, T);
    1465       41668 :   if (p==3)
    1466        9205 :     return F3xq_ellcardj(a4, a6, T, q, n);
    1467       32463 :   if (degpol(a4)<=0 && degpol(a6)<=0)
    1468         217 :     return Fp_ffellcard(utoi(Flx_eval(a4,0,p)),utoi(Flx_eval(a6,0,p)),q,n,utoipos(p));
    1469       32246 :   J = Flxq_ellj(a4,a6,T,p);
    1470       32246 :   if (degpol(J)<=0)
    1471       15820 :     return Flxq_ellcardj(a4,a6,lgpol(J)?J[2]:0,T,q,p,n);
    1472       16426 :   M = Flxq_minpoly(J, T, p);
    1473       16426 :   if (degpol(M) < n)
    1474         454 :     return Flxq_ffellcard(a4, a6, M, q, T, p, n);
    1475       15972 :   if (p <= 7)
    1476       10990 :     return Flxq_ellcard_Satoh(a4, a6, J, T, p);
    1477        4982 :   if (cmpis(q,100)<0)
    1478           0 :     return utoi(Flxq_ellcard_naive(a4, a6, T, p));
    1479        4982 :   if (p == 13 || (7*p <= (ulong)10*n && (BITS_IN_LONG==64 || p <= 103)))
    1480         154 :     return Flxq_ellcard_Satoh(a4, a6, J, T, p);
    1481        4828 :   if (p <= (ulong)2*n)
    1482           0 :     return Flxq_ellcard_Kedlaya(a4, a6, T, p);
    1483        4828 :   if (expi(q)<=62)
    1484        4774 :     return Flxq_ellcard_Shanks(a4, a6, q, T, p);
    1485             :   else
    1486          54 :     return Fq_ellcard_SEA(Flx_to_ZX(a4),Flx_to_ZX(a6),q,Flx_to_ZX(T),utoipos(p),0);
    1487             : }
    1488             : 
    1489             : GEN
    1490       62661 : Flxq_ellcard(GEN a4, GEN a6, GEN T, ulong p)
    1491             : {
    1492       62661 :   pari_sp av = avma;
    1493       62661 :   return gerepileuptoint(av, Flxq_ellcard_i(a4, a6, T, p));
    1494             : }
    1495             : 
    1496             : static long
    1497         350 : Fl_ellj_trace(ulong j, ulong p)
    1498             : {
    1499             :   ulong a4, a6;
    1500         350 :   Fl_ellj_to_a4a6(j, p, &a4, &a6);
    1501         350 :   return Fl_elltrace(a4, a6, p);
    1502             : }
    1503             : 
    1504             : /* Given ordinary E/Fq, a prime ell, and the height of the ell-volcano
    1505             :  * containing j(E) (= v_ell(conductor of Z[pi_E]) returns the height of j(E)
    1506             :  * on its ell-volcano (= v_ell(conductor of the order End(E)). */
    1507             : static long
    1508         147 : Fl_ellheightabovefloor(ulong j, long ell, long e, ulong p)
    1509             : {
    1510         147 :   pari_sp av = avma;
    1511             :   GEN Xp, G, phi, phix, j0, j1;
    1512             :   long h, i, nj1;
    1513         147 :   if (e==0) return 0;
    1514         147 :   if (j==0 || j==1728%p) return e;
    1515         119 :   phi = ZXX_to_FlxX(polmodular_ZXX(ell, 0, 0, 1), p, 1);
    1516         119 :   phix = FlxY_evalx(phi, j, p);
    1517         119 :   Xp = Flx_Frobenius(phix, p);
    1518         119 :   G  = Flx_gcd(Flx_sub(Xp, polx_Flx(0), p), phix, p);
    1519         119 :   nj1 = degpol(G);
    1520         119 :   if (nj1 < ell) return 0;
    1521         112 :   if (e==1 || nj1 != ell+1) return e;
    1522          21 :   j1 = Flx_roots(G, p);
    1523          21 :   nj1 = lg(j1)-1;
    1524          21 :   if (nj1 < 3) return 0;
    1525          21 :   j0 = mkvecsmall3(j,j,j);
    1526          42 :   for (h = 1; ; h++)
    1527         133 :     for(i = 1; i <= 3; i++)
    1528             :     {
    1529         112 :       GEN P = Flx_div_by_X_x(FlxY_evalx(phi, uel(j1,i), p), uel(j0,i), p, NULL);
    1530         112 :       GEN r = Flx_roots(P, p);
    1531         112 :       if (lg(r) == 1) return gc_long(av, h);
    1532          91 :       j0[i] = j1[i];
    1533          91 :       j1[i] = r[1];
    1534             :     }
    1535             : }
    1536             : 
    1537             : /* Given an ordinary elliptic curve E/Fp and an integer h, returns
    1538             :  * D = disc(End(E)) assuming h(D) = h, using the approach sketched in
    1539             :  * Remark 13. If the algorithm returns 0 it has proved that h(D) != h, but it
    1540             :  * is under no obligation to do so and is allowed to return any value when the
    1541             :  * assumption h(d) = h is false. */
    1542             : static long
    1543         350 : Fl_end13(ulong j, ulong h, ulong p)
    1544             : {
    1545             :   ulong D0, v, h0;
    1546             :   long i, lL, lc, lD, nc;
    1547             :   GEN D, DF, cs, L, vP, vE;
    1548         350 :   ulong t = Fl_ellj_trace(j, p);
    1549             : 
    1550         350 :   D0 = coredisc2u_fact(factoru(4*p-t*t), -1, &vP, &vE);
    1551         350 :   h0 = itou(classno(stoi(-D0)));
    1552         350 :   if (h % h0 != 0) return 0;
    1553         350 :   h /= h0;
    1554         350 :   D = divisorsu_fact_factored(mkmat2(vP,vE));
    1555         350 :   DF = gel(D,2); D = gel(D,1);
    1556         350 :   lD = lg(D); v = D[lD-1];
    1557         350 :   cs = cgetg(lD,t_VECSMALL); nc = 0;
    1558        1848 :   for (i = 1; i < lD; i++)
    1559             :   {
    1560        1498 :     GEN F = gel(DF,i);
    1561        1498 :     ulong w = uquadclassnoF_fact(D0, -1, gel(F,1), gel(F,2));
    1562        1498 :     if (w == h) uel(cs,++nc) = v / uel(D,i);
    1563             :   }
    1564         350 :   if (nc==0) return 0;
    1565         350 :   if (nc==1) { v /= uel(cs,1); return -D0*v*v; }
    1566         147 :   L = cgetg(nc+1, t_VEC);
    1567         448 :   for (i = 1; i <= nc; i++) gel(L,i) = gel(factoru(uel(cs,i)), 1);
    1568         147 :   L = vecsmall_uniq(shallowconcat1(L));
    1569         147 :   lL = lg(L); lc = nc+1;
    1570         147 :   for (i = 1; i < lL; i++)
    1571             :   {
    1572         147 :     ulong ell = L[i];
    1573         147 :     long k, e = Fl_ellheightabovefloor(j, ell, z_lval(v,ell), p);
    1574         448 :     for (k = 1; k < lc; k++)
    1575         301 :       if(cs[k] && z_lval(cs[k], ell) != e) { cs[k] = 0; nc--; }
    1576         147 :     if (nc==0) return 0;
    1577         147 :     if (nc==1)
    1578             :     {
    1579         161 :       for (k = 1; k < lc; k++)
    1580         161 :         if (cs[k]) { v /= uel(cs,k); return -D0*v*v; }
    1581             :     }
    1582             :   }
    1583           0 :   return 0;
    1584             : }
    1585             : 
    1586             : INLINE int
    1587         350 : RgX_is_monic_ZX(GEN pol)
    1588         350 : { return RgX_is_ZX(pol) && ZX_is_monic(pol); }
    1589             : 
    1590             : long
    1591         357 : polisclass(GEN H)
    1592             : {
    1593         357 :   pari_sp av = avma, btop;
    1594         357 :   long h = degpol(H), hl, i, pmin, vH = varn(H), vh;
    1595             :   double lmin;
    1596             :   ulong p;
    1597             :   GEN h2list;
    1598             :   forprime_t T;
    1599             : 
    1600         357 :   if (typ(H)!= t_POL) pari_err_TYPE("polsisclass",H);
    1601         357 :   if (h <= 0 || !RgX_is_monic_ZX(H)) return 0;
    1602         350 :   vh = vals(h);
    1603         350 :   h2list = cgetg(vh+2, t_VECSMALL); hl = 1;
    1604        1071 :   for (i = 0; i <= vh; i++)
    1605             :   {
    1606         721 :     ulong d = 1UL<<i;
    1607         721 :     if (((d-h)&1)==0) h2list[hl++] = d;
    1608             :   }
    1609         350 :   setlg(h2list, hl);
    1610         350 :   lmin = h * (log(log(h+2))+2);
    1611         350 :   pmin = 33 * ceil(lmin*lmin);
    1612         350 :   u_forprime_init(&T, pmin, ULONG_MAX);
    1613         350 :   btop = avma;
    1614        3339 :   while((p = u_forprime_next(&T)))
    1615             :   {
    1616             :     ulong r;
    1617             :     long D, nroots;
    1618        3339 :     GEN Xp, G, Hp = ZX_to_Flx(H,p);
    1619        3339 :     if (!Flx_is_squarefree(Hp, p)) { set_avma(btop); continue; }
    1620        3339 :     Xp = Flx_Frobenius(Hp, p);
    1621        3339 :     G  = Flx_gcd(Flx_sub(Xp, polx_Flx(evalvarn(vH)), p), Hp, p);
    1622        3339 :     nroots = degpol(G);
    1623        3339 :     if (nroots==0) { set_avma(btop); continue; }
    1624        1274 :     if (nroots < h && !zv_search(h2list,nroots)) return gc_long(av, 0);
    1625        1274 :     r = Flx_oneroot(G, p);
    1626        1274 :     if (Fp_elljissupersingular(utoi(r), utoi(p))) { set_avma(btop); continue; }
    1627         350 :     D = Fl_end13(r, h, p);
    1628         350 :     if (D && gequal(H, polclass(stoi(D), 0, vH))) return gc_long(av, D);
    1629           0 :     return gc_long(av, 0);
    1630             :   }
    1631           0 :   pari_err_BUG("polisclass");
    1632             :   return 0; /* LCOV_EXCL_LINE */
    1633             : }

Generated by: LCOV version 1.16