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 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.8.0 lcov report (development 19226-b907b8d) Lines: 906 946 95.8 %
Date: 2016-07-29 07:10:27 Functions: 96 97 99.0 %
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. It is distributed in the hope that it will be useful, but WITHOUT
       8             : ANY WARRANTY WHATSOEVER.
       9             : 
      10             : Check the License for details. You should have received a copy of it, along
      11             : with the package; see the file 'COPYING'. If not, write to the Free Software
      12             : Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. */
      13             : 
      14             : #include "pari.h"
      15             : #include "paripriv.h"
      16             : 
      17             : /* Not so fast arithmetic with points over elliptic curves over Fq,
      18             : small characteristic. */
      19             : 
      20             : /***********************************************************************/
      21             : /**                                                                   **/
      22             : /**                              FlxqE                                **/
      23             : /**                                                                   **/
      24             : /***********************************************************************/
      25             : 
      26             : /* Theses functions deal with point over elliptic curves over Fq defined
      27             :  * by an equation of the form y^2=x^3+a4*x+a6.
      28             :  * Most of the time a6 is omitted since it can be recovered from any point
      29             :  * on the curve.
      30             :  */
      31             : 
      32             : GEN
      33       63854 : RgE_to_FlxqE(GEN x, GEN T, ulong p)
      34             : {
      35       63854 :   if (ell_is_inf(x)) return x;
      36       63854 :   retmkvec2(Rg_to_Flxq(gel(x,1),T,p),Rg_to_Flxq(gel(x,2),T,p));
      37             : }
      38             : 
      39             : GEN
      40      152915 : FlxqE_changepoint(GEN x, GEN ch, GEN T, ulong p)
      41             : {
      42      152915 :   pari_sp av = avma;
      43             :   GEN p1,z,u,r,s,t,v,v2,v3;
      44      152915 :   if (ell_is_inf(x)) return x;
      45       90762 :   u = gel(ch,1); r = gel(ch,2);
      46       90762 :   s = gel(ch,3); t = gel(ch,4);
      47       90762 :   v = Flxq_inv(u, T, p); v2 = Flxq_sqr(v, T, p); v3 = Flxq_mul(v,v2, T, p);
      48       90762 :   p1 = Flx_sub(gel(x,1),r, p);
      49       90762 :   z = cgetg(3,t_VEC);
      50       90762 :   gel(z,1) = Flxq_mul(v2, p1, T, p);
      51       90762 :   gel(z,2) = Flxq_mul(v3, Flx_sub(gel(x,2), Flx_add(Flxq_mul(s, p1, T, p),t, p), p), T, p);
      52       90762 :   return gerepileupto(av, z);
      53             : }
      54             : 
      55             : GEN
      56       63854 : FlxqE_changepointinv(GEN x, GEN ch, GEN T, ulong p)
      57             : {
      58             :   GEN u, r, s, t, X, Y, u2, u3, u2X, z;
      59       63854 :   if (ell_is_inf(x)) return x;
      60       63854 :   X = gel(x,1); Y = gel(x,2);
      61       63854 :   u = gel(ch,1); r = gel(ch,2);
      62       63854 :   s = gel(ch,3); t = gel(ch,4);
      63       63854 :   u2 = Flxq_sqr(u, T, p); u3 = Flxq_mul(u,u2, T, p);
      64       63854 :   u2X = Flxq_mul(u2,X, T, p);
      65       63854 :   z = cgetg(3, t_VEC);
      66       63854 :   gel(z,1) = Flx_add(u2X,r, p);
      67       63854 :   gel(z,2) = Flx_add(Flxq_mul(u3,Y, T, p), Flx_add(Flxq_mul(s,u2X, T, p), t, p), p);
      68       63854 :   return z;
      69             : }
      70             : 
      71             : static ulong
      72       20440 : nonsquare_Fl(ulong p)
      73             : {
      74             :   ulong a;
      75             :   do
      76       20440 :     a = random_Fl(p);
      77       20440 :   while (krouu(a, p) >= 0);
      78        7686 :   return a;
      79             : }
      80             : 
      81             : static GEN
      82       22834 : nonsquare_Flxq(GEN T, ulong p)
      83             : {
      84       22834 :   pari_sp av = avma;
      85       22834 :   long n = degpol(T), vs = T[1];
      86             :   GEN a;
      87       22834 :   if (odd(n))
      88        7686 :     return mkvecsmall2(vs, nonsquare_Fl(p));
      89             :   do
      90             :   {
      91       30891 :     avma = av;
      92       30891 :     a = random_Flx(n, vs, p);
      93       30891 :   } while (Flxq_issquare(a, T, p));
      94       15148 :   return a;
      95             : }
      96             : 
      97             : void
      98       22834 : Flxq_elltwist(GEN a, GEN a6, GEN T, ulong p, GEN *pt_a, GEN *pt_a6)
      99             : {
     100       22834 :   GEN d = nonsquare_Flxq(T, p);
     101       22834 :   GEN d2 = Flxq_sqr(d, T, p), d3 = Flxq_mul(d2, d, T, p);
     102       22834 :   if (typ(a)==t_VECSMALL)
     103             :   {
     104       15232 :     *pt_a  = Flxq_mul(a,  d2, T, p);
     105       15232 :     *pt_a6 = Flxq_mul(a6, d3, T, p);
     106             :   } else
     107             :   {
     108        7602 :     *pt_a  = mkvec(Flxq_mul(gel(a,1), d, T, p));
     109        7602 :     *pt_a6 = Flxq_mul(a6, d3, T, p);
     110             :   }
     111       22834 : }
     112             : 
     113             : static GEN
     114     1281730 : FlxqE_dbl_slope(GEN P, GEN a4, GEN T, ulong p, GEN *slope)
     115             : {
     116             :   GEN x, y, Q;
     117     1281730 :   if (ell_is_inf(P) || !lgpol(gel(P,2))) return ellinf();
     118     1182690 :   x = gel(P,1); y = gel(P,2);
     119     1182690 :   if (p==3UL)
     120     1570492 :     *slope = typ(a4)==t_VEC ? Flxq_div(Flxq_mul(x, gel(a4, 1), T, p), y, T, p)
     121     1040907 :                             : Flxq_div(a4, Flx_neg(y, p), T, p);
     122             :   else
     123             :   {
     124      653105 :     GEN sx = Flx_add(Flx_triple(Flxq_sqr(x, T, p), p), a4, p);
     125      653105 :     *slope = Flxq_div(sx, Flx_double(y, p), T, p);
     126             :   }
     127     1182690 :   Q = cgetg(3,t_VEC);
     128     1182690 :   gel(Q, 1) = Flx_sub(Flxq_sqr(*slope, T, p), Flx_double(x, p), p);
     129     1182690 :   if (typ(a4)==t_VEC) gel(Q, 1) = Flx_sub(gel(Q, 1), gel(a4, 1), p);
     130     1182690 :   gel(Q, 2) = Flx_sub(Flxq_mul(*slope, Flx_sub(x, gel(Q, 1), p), T, p), y, p);
     131     1182690 :   return Q;
     132             : }
     133             : 
     134             : GEN
     135     1255028 : FlxqE_dbl(GEN P, GEN a4, GEN T, ulong p)
     136             : {
     137     1255028 :   pari_sp av = avma;
     138             :   GEN slope;
     139     1255028 :   return gerepileupto(av, FlxqE_dbl_slope(P,a4, T, p,&slope));
     140             : }
     141             : 
     142             : static GEN
     143      532823 : FlxqE_add_slope(GEN P, GEN Q, GEN a4, GEN T, ulong p, GEN *slope)
     144             : {
     145             :   GEN Px, Py, Qx, Qy, R;
     146      532823 :   if (ell_is_inf(P)) return Q;
     147      529330 :   if (ell_is_inf(Q)) return P;
     148      529162 :   Px = gel(P,1); Py = gel(P,2);
     149      529162 :   Qx = gel(Q,1); Qy = gel(Q,2);
     150      529162 :   if (Flx_equal(Px, Qx))
     151             :   {
     152       47366 :     if (Flx_equal(Py, Qy))
     153        1329 :       return FlxqE_dbl_slope(P, a4, T, p, slope);
     154             :     else
     155       46037 :       return ellinf();
     156             :   }
     157      481796 :   *slope = Flxq_div(Flx_sub(Py, Qy, p), Flx_sub(Px, Qx, p), T, p);
     158      481796 :   R = cgetg(3,t_VEC);
     159      481796 :   gel(R, 1) = Flx_sub(Flx_sub(Flxq_sqr(*slope, T, p), Px, p), Qx, p);
     160      481796 :   if (typ(a4)==t_VEC) gel(R, 1) = Flx_sub(gel(R, 1),gel(a4, 1), p);
     161      481796 :   gel(R, 2) = Flx_sub(Flxq_mul(*slope, Flx_sub(Px, gel(R, 1), p), T, p), Py, p);
     162      481796 :   return R;
     163             : }
     164             : 
     165             : GEN
     166      529119 : FlxqE_add(GEN P, GEN Q, GEN a4, GEN T, ulong p)
     167             : {
     168      529119 :   pari_sp av = avma;
     169             :   GEN slope;
     170      529119 :   return gerepileupto(av, FlxqE_add_slope(P,Q,a4, T, p,&slope));
     171             : }
     172             : 
     173             : static GEN
     174         987 : FlxqE_neg_i(GEN P, ulong p)
     175             : {
     176         987 :   if (ell_is_inf(P)) return P;
     177         987 :   return mkvec2(gel(P,1), Flx_neg(gel(P,2), p));
     178             : }
     179             : 
     180             : GEN
     181         399 : FlxqE_neg(GEN P, GEN T, ulong p)
     182             : {
     183             :   (void) T;
     184         399 :   if (ell_is_inf(P)) return ellinf();
     185         399 :   return mkvec2(gcopy(gel(P,1)), Flx_neg(gel(P,2), p));
     186             : }
     187             : 
     188             : GEN
     189         987 : FlxqE_sub(GEN P, GEN Q, GEN a4, GEN T, ulong p)
     190             : {
     191         987 :   pari_sp av = avma;
     192             :   GEN slope;
     193         987 :   return gerepileupto(av, FlxqE_add_slope(P, FlxqE_neg_i(Q, p), a4, T, p, &slope));
     194             : }
     195             : 
     196             : struct _FlxqE
     197             : {
     198             :   GEN a4, a6;
     199             :   GEN T;
     200             :   ulong p;
     201             : };
     202             : 
     203             : static GEN
     204     1255028 : _FlxqE_dbl(void *E, GEN P)
     205             : {
     206     1255028 :   struct _FlxqE *ell = (struct _FlxqE *) E;
     207     1255028 :   return FlxqE_dbl(P, ell->a4, ell->T, ell->p);
     208             : }
     209             : 
     210             : static GEN
     211      521559 : _FlxqE_add(void *E, GEN P, GEN Q)
     212             : {
     213      521559 :   struct _FlxqE *ell=(struct _FlxqE *) E;
     214      521559 :   return FlxqE_add(P, Q, ell->a4, ell->T, ell->p);
     215             : }
     216             : 
     217             : static GEN
     218      215575 : _FlxqE_mul(void *E, GEN P, GEN n)
     219             : {
     220      215575 :   pari_sp av = avma;
     221      215575 :   struct _FlxqE *e=(struct _FlxqE *) E;
     222      215575 :   long s = signe(n);
     223      215575 :   if (!s || ell_is_inf(P)) return ellinf();
     224      215281 :   if (s<0) P = FlxqE_neg(P, e->T, e->p);
     225      215281 :   if (is_pm1(n)) return s>0? gcopy(P): P;
     226      210864 :   return gerepileupto(av, gen_pow(P, n, e, &_FlxqE_dbl, &_FlxqE_add));
     227             : }
     228             : 
     229             : GEN
     230       64848 : FlxqE_mul(GEN P, GEN n, GEN a4, GEN T, ulong p)
     231             : {
     232             :   struct _FlxqE E;
     233       64848 :   E.a4= a4; E.T = T; E.p = p;
     234       64848 :   return _FlxqE_mul(&E, P, n);
     235             : }
     236             : 
     237             : /* 3*x^2+2*a2*x = -a2*x, and a2!=0 */
     238             : 
     239             : /* Finds a random non-singular point on E */
     240             : static GEN
     241       76048 : random_F3xqE(GEN a2, GEN a6, GEN T)
     242             : {
     243       76048 :   pari_sp ltop = avma;
     244             :   GEN x, y, rhs;
     245       76048 :   const ulong p=3;
     246             :   do
     247             :   {
     248      152159 :     avma= ltop;
     249      152159 :     x   = random_Flx(get_Flx_degree(T),get_Flx_var(T),p);
     250      152159 :     rhs = Flx_add(Flxq_mul(Flxq_sqr(x, T, p), Flx_add(x, a2, p), T, p), a6, p);
     251      152159 :   } while ((!lgpol(rhs) && !lgpol(x)) || !Flxq_issquare(rhs, T, p));
     252       76048 :   y = Flxq_sqrt(rhs, T, p);
     253       76048 :   if (!y) pari_err_PRIME("random_F3xqE", T);
     254       76048 :   return gerepilecopy(ltop, mkvec2(x, y));
     255             : }
     256             : 
     257             : /* Finds a random non-singular point on E */
     258             : GEN
     259      141129 : random_FlxqE(GEN a4, GEN a6, GEN T, ulong p)
     260             : {
     261      141129 :   pari_sp ltop = avma;
     262             :   GEN x, x2, y, rhs;
     263      141129 :   if (typ(a4)==t_VEC)
     264       76048 :     return random_F3xqE(gel(a4,1), a6, T);
     265             :   do
     266             :   {
     267      130758 :     avma= ltop;
     268      130758 :     x   = random_Flx(get_Flx_degree(T),get_Flx_var(T),p);
     269      130758 :     x2  = Flxq_sqr(x, T, p); /*  x^3+a4*x+a6 = x*(x^2+a4)+a6  */
     270      130758 :     rhs = Flx_add(Flxq_mul(x, Flx_add(x2, a4, p), T, p), a6, p);
     271      131780 :   } while ((!lgpol(rhs) && !lgpol(Flx_add(Flx_triple(x2, p), a4, p)))
     272      261516 :           || !Flxq_issquare(rhs, T, p));
     273       65081 :   y = Flxq_sqrt(rhs, T, p);
     274       65081 :   if (!y) pari_err_PRIME("random_FlxqE", T);
     275       65081 :   return gerepilecopy(ltop, mkvec2(x, y));
     276             : }
     277             : 
     278             : static GEN
     279       66019 : _FlxqE_rand(void *E)
     280             : {
     281       66019 :   struct _FlxqE *ell=(struct _FlxqE *) E;
     282       66019 :   return random_FlxqE(ell->a4, ell->a6, ell->T, ell->p);
     283             : }
     284             : 
     285             : static const struct bb_group FlxqE_group={_FlxqE_add,_FlxqE_mul,_FlxqE_rand,hash_GEN,zvV_equal,ell_is_inf, NULL};
     286             : 
     287             : const struct bb_group *
     288          34 : get_FlxqE_group(void ** pt_E, GEN a4, GEN a6, GEN T, ulong p)
     289             : {
     290          34 :   struct _FlxqE *e = (struct _FlxqE *) stack_malloc(sizeof(struct _FlxqE));
     291          34 :   e->a4 = a4; e->a6 = a6; e->T = Flx_get_red(T, p); e->p = p;
     292          34 :   *pt_E = (void *) e;
     293          34 :   return &FlxqE_group;
     294             : }
     295             : 
     296             : GEN
     297        2471 : FlxqE_order(GEN z, GEN o, GEN a4, GEN T, ulong p)
     298             : {
     299        2471 :   pari_sp av = avma;
     300             :   struct _FlxqE e;
     301        2471 :   e.a4=a4; e.T=T; e.p=p;
     302        2471 :   return gerepileuptoint(av, gen_order(z, o, (void*)&e, &FlxqE_group));
     303             : }
     304             : 
     305             : GEN
     306          42 : FlxqE_log(GEN a, GEN b, GEN o, GEN a4, GEN T, ulong p)
     307             : {
     308          42 :   pari_sp av = avma;
     309             :   struct _FlxqE e;
     310          42 :   e.a4=a4; e.T=T; e.p=p;
     311          42 :   return gerepileuptoint(av, gen_PH_log(a, b, o, (void*)&e, &FlxqE_group));
     312             : }
     313             : 
     314             : /***********************************************************************/
     315             : /**                                                                   **/
     316             : /**                            Pairings                               **/
     317             : /**                                                                   **/
     318             : /***********************************************************************/
     319             : 
     320             : /* Derived from APIP from and by Jerome Milan, 2012 */
     321             : 
     322             : static GEN
     323       69400 : FlxqE_vert(GEN P, GEN Q, GEN a4, GEN T, ulong p)
     324             : {
     325       69400 :   long vT = get_Flx_var(T);
     326             :   GEN df;
     327       69400 :   if (ell_is_inf(P))
     328       22505 :     return pol1_Flx(vT);
     329       46895 :   if (!Flx_equal(gel(Q, 1), gel(P, 1)))
     330       42407 :     return Flx_sub(gel(Q, 1), gel(P, 1), p);
     331        4488 :   if (lgpol(gel(P,2))!=0) return pol1_Flx(vT);
     332        9828 :   df = typ(a4)==t_VEC ? Flxq_mul(gel(P,1), Flx_mulu(gel(a4, 1), 2, p), T, p)
     333        6006 :                       : a4;
     334        3822 :   return Flxq_inv(Flx_add(Flx_mulu(Flxq_sqr(gel(P,1), T, p), 3, p),
     335             :                           df, p), T, p);
     336             : }
     337             : 
     338             : static GEN
     339       28090 : FlxqE_Miller_line(GEN R, GEN Q, GEN slope, GEN a4, GEN T, ulong p)
     340             : {
     341       28090 :   long vT = get_Flx_var(T);
     342       28090 :   GEN x = gel(Q, 1), y = gel(Q, 2);
     343       28090 :   GEN tmp1 = Flx_sub(x, gel(R, 1), p);
     344       28090 :   GEN tmp2 = Flx_add(Flxq_mul(tmp1, slope, T, p), gel(R, 2), p);
     345       28090 :   if (!Flx_equal(y, tmp2))
     346       26332 :     return Flx_sub(y, tmp2, p);
     347        1758 :   if (lgpol(y) == 0)
     348         483 :     return pol1_Flx(vT);
     349             :   else
     350             :   {
     351        1275 :     GEN s1, s2, a2 = typ(a4)==t_VEC ? gel(a4,1): NULL;
     352        1275 :     GEN y2i = Flxq_inv(Flx_mulu(y, 2, p), T, p);
     353        1275 :     GEN df = a2 ? Flxq_mul(x, Flx_mulu(a2, 2, p), T, p): a4;
     354             :     GEN x3, ddf;
     355        1275 :     s1 = Flxq_mul(Flx_add(Flx_mulu(Flxq_sqr(x, T, p), 3, p), df, p), y2i, T, p);
     356        1275 :     if (!Flx_equal(s1, slope))
     357         344 :       return Flx_sub(s1, slope, p);
     358         931 :     x3 = Flx_mulu(x, 3, p);
     359         931 :     ddf = a2 ? Flx_add(x3, a2, p): x3;
     360         931 :     s2 = Flxq_mul(Flx_sub(ddf, Flxq_sqr(s1, T, p), p), y2i, T, p);
     361         931 :     return lgpol(s2)!=0 ? s2: y2i;
     362             :   }
     363             : }
     364             : 
     365             : /* Computes the equation of the line tangent to R and returns its
     366             :    evaluation at the point Q. Also doubles the point R.
     367             :  */
     368             : 
     369             : static GEN
     370       46534 : FlxqE_tangent_update(GEN R, GEN Q, GEN a4, GEN T, ulong p, GEN *pt_R)
     371             : {
     372       46534 :   if (ell_is_inf(R))
     373             :   {
     374        3798 :     *pt_R = ellinf();
     375        3798 :     return pol1_Flx(get_Flx_var(T));
     376             :   }
     377       42736 :   else if (!lgpol(gel(R,2)))
     378             :   {
     379       17363 :     *pt_R = ellinf();
     380       17363 :     return FlxqE_vert(R, Q, a4, T, p);
     381             :   } else {
     382             :     GEN slope;
     383       25373 :     *pt_R = FlxqE_dbl_slope(R, a4, T, p, &slope);
     384       25373 :     return FlxqE_Miller_line(R, Q, slope, a4, T, p);
     385             :   }
     386             : }
     387             : 
     388             : /* Computes the equation of the line through R and P, and returns its
     389             :    evaluation at the point Q. Also adds P to the point R.
     390             :  */
     391             : 
     392             : static GEN
     393        4110 : FlxqE_chord_update(GEN R, GEN P, GEN Q, GEN a4, GEN T, ulong p, GEN *pt_R)
     394             : {
     395        4110 :   if (ell_is_inf(R))
     396             :   {
     397          49 :     *pt_R = gcopy(P);
     398          49 :     return FlxqE_vert(P, Q, a4, T, p);
     399             :   }
     400        4061 :   else if (ell_is_inf(P))
     401             :   {
     402           0 :     *pt_R = gcopy(R);
     403           0 :     return FlxqE_vert(R, Q, a4, T, p);
     404             :   }
     405        4061 :   else if (Flx_equal(gel(P, 1), gel(R, 1)))
     406             :   {
     407        1344 :     if (Flx_equal(gel(P, 2), gel(R, 2)))
     408           0 :       return FlxqE_tangent_update(R, Q, a4, T, p, pt_R);
     409             :     else
     410             :     {
     411        1344 :       *pt_R = ellinf();
     412        1344 :       return FlxqE_vert(R, Q, a4, T, p);
     413             :     }
     414             :   } else {
     415             :     GEN slope;
     416        2717 :     *pt_R = FlxqE_add_slope(P, R, a4, T, p, &slope);
     417        2717 :     return FlxqE_Miller_line(R, Q, slope, a4, T, p);
     418             :   }
     419             : }
     420             : 
     421             : /* Returns the Miller function f_{m, Q} evaluated at the point P using
     422             :    the standard Miller algorithm.
     423             :  */
     424             : 
     425             : struct _FlxqE_miller
     426             : {
     427             :   ulong p;
     428             :   GEN T, a4, P;
     429             : };
     430             : 
     431             : static GEN
     432       46534 : FlxqE_Miller_dbl(void* E, GEN d)
     433             : {
     434       46534 :   struct _FlxqE_miller *m = (struct _FlxqE_miller *)E;
     435       46534 :   ulong p  = m->p;
     436       46534 :   GEN T = m->T, a4 = m->a4, P = m->P;
     437             :   GEN v, line;
     438       46534 :   GEN num = Flxq_sqr(gel(d,1), T, p);
     439       46534 :   GEN denom = Flxq_sqr(gel(d,2), T, p);
     440       46534 :   GEN point = gel(d,3);
     441       46534 :   line = FlxqE_tangent_update(point, P, a4, T, p, &point);
     442       46534 :   num  = Flxq_mul(num, line, T, p);
     443       46534 :   v = FlxqE_vert(point, P, a4, T, p);
     444       46534 :   denom = Flxq_mul(denom, v, T, p);
     445       46534 :   return mkvec3(num, denom, point);
     446             : }
     447             : 
     448             : static GEN
     449        4110 : FlxqE_Miller_add(void* E, GEN va, GEN vb)
     450             : {
     451        4110 :   struct _FlxqE_miller *m = (struct _FlxqE_miller *)E;
     452        4110 :   ulong p = m->p;
     453        4110 :   GEN T = m->T, a4 = m->a4, P = m->P;
     454             :   GEN v, line, point;
     455        4110 :   GEN na = gel(va,1), da = gel(va,2), pa = gel(va,3);
     456        4110 :   GEN nb = gel(vb,1), db = gel(vb,2), pb = gel(vb,3);
     457        4110 :   GEN num   = Flxq_mul(na, nb, T, p);
     458        4110 :   GEN denom = Flxq_mul(da, db, T, p);
     459        4110 :   line = FlxqE_chord_update(pa, pb, P, a4, T, p, &point);
     460        4110 :   num  = Flxq_mul(num, line, T, p);
     461        4110 :   v = FlxqE_vert(point, P, a4, T, p);
     462        4110 :   denom = Flxq_mul(denom, v, T, p);
     463        4110 :   return mkvec3(num, denom, point);
     464             : }
     465             : 
     466             : static GEN
     467       18658 : FlxqE_Miller(GEN Q, GEN P, GEN m, GEN a4, GEN T, ulong p)
     468             : {
     469       18658 :   pari_sp ltop = avma;
     470             :   struct _FlxqE_miller d;
     471             :   GEN v, num, denom, g1;
     472             : 
     473       18658 :   d.a4 = a4; d.T = T; d.p = p; d.P = P;
     474       18658 :   g1 = pol1_Flx(get_Flx_var(T));
     475       18658 :   v = gen_pow(mkvec3(g1,g1,Q), m, (void*)&d, FlxqE_Miller_dbl, FlxqE_Miller_add);
     476       18658 :   num = gel(v,1); denom = gel(v,2);
     477       18658 :   return gerepileupto(ltop, Flxq_div(num, denom, T, p));
     478             : }
     479             : 
     480             : GEN
     481       12228 : FlxqE_weilpairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, ulong p)
     482             : {
     483       12228 :   pari_sp ltop = avma;
     484             :   GEN num, denom, result;
     485       12228 :   if (ell_is_inf(P) || ell_is_inf(Q) || Flx_equal(P,Q))
     486        2927 :     return pol1_Flx(get_Flx_var(T));
     487        9301 :   num    = FlxqE_Miller(P, Q, m, a4, T, p);
     488        9301 :   denom  = FlxqE_Miller(Q, P, m, a4, T, p);
     489        9301 :   result = Flxq_div(num, denom, T, p);
     490        9301 :   if (mpodd(m))
     491         602 :     result  = Flx_neg(result, p);
     492        9301 :   return gerepileupto(ltop, result);
     493             : }
     494             : 
     495             : GEN
     496          56 : FlxqE_tatepairing(GEN P, GEN Q, GEN m, GEN a4, GEN T, ulong p)
     497             : {
     498          56 :   if (ell_is_inf(P) || ell_is_inf(Q))
     499           0 :     return pol1_Flx(get_Flx_var(T));
     500          56 :   return FlxqE_Miller(P, Q, m, a4, T, p);
     501             : }
     502             : 
     503             : static GEN
     504       12214 : _FlxqE_pairorder(void *E, GEN P, GEN Q, GEN m, GEN F)
     505             : {
     506       12214 :   struct _FlxqE *e = (struct _FlxqE *) E;
     507       12214 :   return  Flxq_order(FlxqE_weilpairing(P,Q,m,e->a4,e->T,e->p), F, e->T, e->p);
     508             : }
     509             : 
     510             : GEN
     511       14385 : Flxq_ellgroup(GEN a4, GEN a6, GEN N, GEN T, ulong p, GEN *pt_m)
     512             : {
     513             :   struct _FlxqE e;
     514       14385 :   GEN q = powuu(p, get_Flx_degree(T));
     515       14385 :   e.a4=a4; e.a6=a6; e.T=T; e.p=p;
     516       14385 :   return gen_ellgroup(N, subis(q,1), pt_m, (void*)&e, &FlxqE_group, _FlxqE_pairorder);
     517             : }
     518             : 
     519             : GEN
     520       13153 : Flxq_ellgens(GEN a4, GEN a6, GEN ch, GEN D, GEN m, GEN T, ulong p)
     521             : {
     522             :   GEN P;
     523       13153 :   pari_sp av = avma;
     524             :   struct _FlxqE e;
     525       13153 :   e.a4=a4; e.a6=a6; e.T=T; e.p=p;
     526       13153 :   switch(lg(D)-1)
     527             :   {
     528             :   case 1:
     529       10703 :     P = gen_gener(gel(D,1), (void*)&e, &FlxqE_group);
     530       10703 :     P = mkvec(FlxqE_changepoint(P, ch, T, p));
     531       10703 :     break;
     532             :   default:
     533        2450 :     P = gen_ellgens(gel(D,1), gel(D,2), m, (void*)&e, &FlxqE_group, _FlxqE_pairorder);
     534        2450 :     gel(P,1) = FlxqE_changepoint(gel(P,1), ch, T, p);
     535        2450 :     gel(P,2) = FlxqE_changepoint(gel(P,2), ch, T, p);
     536        2450 :     break;
     537             :   }
     538       13153 :   return gerepilecopy(av, P);
     539             : }
     540             : /***********************************************************************/
     541             : /**                                                                   **/
     542             : /**                          Point counting                           **/
     543             : /**                                                                   **/
     544             : /***********************************************************************/
     545             : 
     546       11130 : static GEN _can_invl(void *E, GEN V) {(void) E; return V; }
     547             : 
     548        3654 : static GEN _can_lin(void *E, GEN F, GEN V, GEN q)
     549             : {
     550        3654 :   GEN v = RgX_splitting(V, 3);
     551             :   (void) E;
     552        3654 :   return FpX_sub(V,ZXV_dotproduct(v, F), q);
     553             : }
     554             : 
     555             : static GEN
     556        7476 : _can_iter(void *E, GEN f, GEN q)
     557             : {
     558        7476 :   GEN h = RgX_splitting(f,3);
     559        7476 :   GEN h1s = ZX_sqr(gel(h,1)), h2s = ZX_sqr(gel(h,2)), h3s = ZX_sqr(gel(h,3));
     560        7476 :   GEN h12 = ZX_mul(gel(h,1), gel(h,2));
     561        7476 :   GEN h13 = ZX_mul(gel(h,1), gel(h,3));
     562        7476 :   GEN h23 = ZX_mul(gel(h,2), gel(h,3));
     563        7476 :   GEN h1c = ZX_mul(gel(h,1), h1s);
     564        7476 :   GEN h3c = ZX_mul(gel(h,3), h3s);
     565        7476 :   GEN th = ZX_mul(ZX_sub(h2s,ZX_mulu(h13,3)),gel(h,2));
     566        7476 :   GEN y = FpX_sub(f,ZX_add(RgX_shift_shallow(h3c,2),ZX_add(RgX_shift_shallow(th,1),h1c)),q);
     567             :   (void) E;
     568        7476 :   return mkvecn(7,y,h1s,h2s,h3s,h12,h13,h23);
     569             : }
     570             : 
     571             : static GEN
     572        7476 : _can_invd(void *E, GEN V, GEN v, GEN qM, long M)
     573             : {
     574        7476 :   GEN h1s=gel(v,2), h2s=gel(v,3), h3s=gel(v,4);
     575        7476 :   GEN h12=gel(v,5), h13=gel(v,6), h23=gel(v,7);
     576        7476 :   GEN F = mkvec3(ZX_sub(h1s,RgX_shift_shallow(h23,1)),RgX_shift_shallow(ZX_sub(h2s,h13),1),
     577             :                  ZX_sub(RgX_shift_shallow(h3s,2),RgX_shift_shallow(h12,1)));
     578             :   (void)E;
     579        7476 :   return gen_ZpX_Dixon(ZXV_Z_mul(F, utoi(3)), V, qM, utoi(3), M, NULL,
     580             :                                                  _can_lin, _can_invl);
     581             : }
     582             : 
     583             : static GEN
     584        3717 : F3x_canonlift(GEN P, long n)
     585        3717 : { return gen_ZpX_Newton(Flx_to_ZX(P),utoi(3), n, NULL, _can_iter, _can_invd); }
     586             : 
     587       29316 : static GEN _can5_invl(void *E, GEN V) {(void) E; return V; }
     588             : 
     589        8820 : static GEN _can5_lin(void *E, GEN F, GEN V, GEN q)
     590             : {
     591        8820 :   ulong p = *(ulong*)E;
     592        8820 :   GEN v = RgX_splitting(V, p);
     593        8820 :   return FpX_sub(V,ZXV_dotproduct(v, F), q);
     594             : }
     595             : 
     596             : /* P(X,t) -> P(X*t^n,t) mod (t^p-1) */
     597             : static GEN
     598       61775 : _shift(GEN P, long n, ulong p, long v)
     599             : {
     600       61775 :   long i, l=lg(P);
     601       61775 :   GEN r = cgetg(l,t_POL); r[1] = P[1];
     602      475020 :   for(i=2;i<l;i++)
     603             :   {
     604      413245 :     long s = n*(i-2)%p;
     605      413245 :     GEN ci = gel(P,i);
     606      413245 :     if (typ(ci)==t_INT)
     607      102942 :       gel(r,i) = monomial(ci, s, v);
     608             :     else
     609      310303 :       gel(r,i) = RgX_rotate_shallow(ci, s, p);
     610             :   }
     611       61775 :   return FpXX_renormalize(r, l);
     612             : }
     613             : 
     614             : struct _can_mul
     615             : {
     616             :   GEN T, q;
     617             :   ulong p;
     618             : };
     619             : 
     620             : static GEN
     621       41279 : _can5_mul(void *E, GEN A, GEN B)
     622             : {
     623       41279 :   struct _can_mul *d = (struct _can_mul *)E;
     624       41279 :   GEN a = gel(A,1), b = gel(B,1);
     625       41279 :   long n = itos(gel(A,2));
     626       41279 :   GEN bn = _shift(b, n, d->p, get_FpX_var(d->T));
     627       41279 :   GEN c = FpXQX_mul(a, bn, d->T, d->q);
     628       41279 :   return mkvec2(c, addii(gel(A,2), gel(B,2)));
     629             : }
     630             : 
     631             : static GEN
     632       41083 : _can5_sqr(void *E, GEN A)
     633             : {
     634       41083 :   return _can5_mul(E,A,A);
     635             : }
     636             : 
     637             : static GEN
     638       20496 : _can5_iter(void *E, GEN f, GEN q)
     639             : {
     640       20496 :   pari_sp av = avma;
     641             :   struct _can_mul D;
     642       20496 :   ulong p = *(ulong*)E;
     643       20496 :   long i, vT = fetch_var();
     644             :   GEN N, P, d, V, fs;
     645       20496 :   D.q = q; D.T = ZX_Z_sub(monomial(gen_1,p,vT),gen_1);
     646       20496 :   D.p = p;
     647       20496 :   fs = mkvec2(_shift(f, 1, p, vT), gen_1);
     648       20496 :   N = gel(gen_powu(fs,p-1,(void*)&D,_can5_sqr,_can5_mul),1);
     649       20496 :   N = simplify_shallow(FpXQX_red(N,polcyclo(p,vT),q));
     650       20496 :   P = FpX_mul(N,f,q);
     651       20496 :   P = RgX_deflate(P, p);
     652       20496 :   d = RgX_splitting(N, p);
     653       20496 :   V = cgetg(p+1,t_VEC);
     654       20496 :   gel(V,1) = ZX_mulu(gel(d,1), p);
     655      103306 :   for(i=2; i<= (long)p; i++)
     656       82810 :     gel(V,i) = ZX_mulu(RgX_shift_shallow(gel(d,p+2-i), 1), p);
     657       20496 :   (void)delete_var(); return gerepilecopy(av, mkvec2(ZX_sub(f,P),V));
     658             : }
     659             : 
     660             : static GEN
     661       20496 : _can5_invd(void *E, GEN H, GEN v, GEN qM, long M)
     662             : {
     663       20496 :   ulong p = *(long*)E;
     664       20496 :   return gen_ZpX_Dixon(gel(v,2), H, qM, utoi(p), M, E, _can5_lin, _can5_invl);
     665             : }
     666             : 
     667             : static GEN
     668       13930 : Flx_canonlift(GEN P, long n, ulong p)
     669             : {
     670       24143 :   return p==3 ? F3x_canonlift(P,n):
     671       10213 :          gen_ZpX_Newton(Flx_to_ZX(P),utoi(p), n, &p, _can5_iter, _can5_invd);
     672             : }
     673             : 
     674             : /* assume a and n  are coprime */
     675             : static GEN
     676       76118 : RgX_circular_shallow(GEN P, long a, long n)
     677             : {
     678       76118 :   long i, l = lgpol(P);
     679       76118 :   GEN Q = cgetg(2+n,t_POL);
     680       76118 :   Q[1] = P[1];
     681      510923 :   for(i=0; i<l; i++)
     682      434805 :     gel(Q,2+(i*a)%n) = gel(P,2+i);
     683      168399 :   for(   ; i<n; i++)
     684       92281 :     gel(Q,2+(i*a)%n) = gen_0;
     685       76118 :   return normalizepol_lg(Q,2+n);
     686             : }
     687             : 
     688             : static GEN
     689       76118 : ZpXQ_frob_cyc(GEN x, GEN T, GEN q, ulong p)
     690             : {
     691       76118 :   long n = get_FpX_degree(T);
     692       76118 :   return FpX_rem(RgX_circular_shallow(x,p,n+1), T, q);
     693             : }
     694             : 
     695             : static GEN
     696      113120 : ZpXQ_frob(GEN x, GEN Xm, GEN T, GEN q, ulong p)
     697             : {
     698      113120 :   if (lg(Xm)==1)
     699       43344 :     return ZpXQ_frob_cyc(x, T, q, p);
     700             :   else
     701             :   {
     702       69776 :     long n = get_FpX_degree(T);
     703       69776 :     GEN V = RgX_blocks(RgX_inflate(x, p), n, p);
     704       69776 :     GEN W = ZXV_dotproduct(V, Xm);
     705       69776 :     return FpX_rem(W, T, q);
     706             :   }
     707             : }
     708             : 
     709             : struct _lift_lin
     710             : {
     711             :   ulong p;
     712             :   GEN sqx, Tp;
     713             :   GEN ai, Xm;
     714             : };
     715             : 
     716       83790 : static GEN _lift_invl(void *E, GEN x)
     717             : {
     718       83790 :   struct _lift_lin *d = (struct _lift_lin *) E;
     719       83790 :   GEN T = d->Tp;
     720       83790 :   ulong p = d->p;
     721       83790 :   GEN xai = Flxq_mul(ZX_to_Flx(x, p), d->ai, T, p);
     722       83790 :   return Flx_to_ZX(Flxq_lroot_fast(xai, d->sqx, T, p));
     723             : }
     724             : 
     725       23562 : static GEN _lift_lin(void *E, GEN F, GEN x2, GEN q)
     726             : {
     727       23562 :   struct _lift_lin *d = (struct _lift_lin *) E;
     728       23562 :   pari_sp av = avma;
     729       23562 :   GEN T = gel(F,3), Xm = gel(F,4);
     730       23562 :   GEN y2  = ZpXQ_frob(x2, Xm, T, q, d->p);
     731       23562 :   GEN lin = FpX_add(ZX_mul(gel(F,1), y2), ZX_mul(gel(F,2), x2), q);
     732       23562 :   return gerepileupto(av, FpX_rem(lin, T, q));
     733             : }
     734             : 
     735             : static GEN
     736      180684 : FpM_FpXV_bilinear(GEN P, GEN X, GEN Y, GEN p)
     737             : {
     738      180684 :    pari_sp av = avma;
     739      180684 :    GEN s =  ZX_mul(FpXV_FpC_mul(X,gel(P,1),p),gel(Y,1));
     740      180684 :    long i, l = lg(P);
     741      848358 :    for(i=2; i<l; i++)
     742      667674 :      s = ZX_add(s, ZX_mul(FpXV_FpC_mul(X,gel(P,i),p),gel(Y,i)));
     743      180684 :    return gerepileupto(av, FpX_red(s, p));
     744             : }
     745             : 
     746             : static GEN
     747      180684 : FpM_FpXQV_bilinear(GEN P, GEN X, GEN Y, GEN T, GEN p)
     748             : {
     749      180684 :   return FpX_rem(FpM_FpXV_bilinear(P,X,Y,p),T,p);
     750             : }
     751             : 
     752             : static GEN
     753      120456 : FpXC_powderiv(GEN M, GEN p)
     754             : {
     755             :   long i, l;
     756      120456 :   long v = varn(gel(M,2));
     757      120456 :   GEN m = cgetg_copy(M, &l);
     758      120456 :   gel(m,1) = pol_0(v);
     759      120456 :   gel(m,2) = pol_1(v);
     760      445116 :   for(i=2; i<l-1; i++)
     761      324660 :     gel(m,i+1) = FpX_Fp_mul(gel(M,i),utoi(i), p);
     762      120456 :   return m;
     763             : }
     764             : 
     765             : struct _lift_iso
     766             : {
     767             :   GEN phi;
     768             :   GEN Xm,T;
     769             :   GEN sqx, Tp;
     770             :   ulong p;
     771             : };
     772             : 
     773             : static GEN
     774       60228 : _lift_iter(void *E, GEN x2, GEN q)
     775             : {
     776       60228 :   struct _lift_iso *d = (struct _lift_iso *) E;
     777       60228 :   ulong p = d->p;
     778       60228 :   long n = lg(d->phi)-2;
     779       60228 :   GEN TN = FpXT_red(d->T, q), XN = FpXV_red(d->Xm, q);
     780       60228 :   GEN y2 = ZpXQ_frob(x2, XN, TN, q, p);
     781       60228 :   GEN xp = FpXQ_powers(x2, n, TN, q);
     782       60228 :   GEN yp = FpXQ_powers(y2, n, TN, q);
     783       60228 :   GEN V  = FpM_FpXQV_bilinear(d->phi,xp,yp,TN,q);
     784       60228 :   return mkvec3(V,xp,yp);
     785             : }
     786             : 
     787             : static GEN
     788       60228 : _lift_invd(void *E, GEN V, GEN v, GEN qM, long M)
     789             : {
     790       60228 :   struct _lift_iso *d = (struct _lift_iso *) E;
     791             :   struct _lift_lin e;
     792       60228 :   ulong p = d->p;
     793       60228 :   GEN TM = FpXT_red(d->T, qM), XM = FpXV_red(d->Xm, qM);
     794       60228 :   GEN xp = FpXV_red(gel(v,2), qM);
     795       60228 :   GEN yp = FpXV_red(gel(v,3), qM);
     796       60228 :   GEN Dx = FpM_FpXQV_bilinear(d->phi, FpXC_powderiv(xp, qM), yp, TM, qM);
     797       60228 :   GEN Dy = FpM_FpXQV_bilinear(d->phi, xp, FpXC_powderiv(yp, qM), TM, qM);
     798       60228 :   GEN F = mkvec4(Dy, Dx, TM, XM);
     799       60228 :   e.ai = Flxq_inv(ZX_to_Flx(Dy,p),d->Tp,p);
     800       60228 :   e.sqx = d->sqx; e.Tp = d->Tp; e.p=p; e.Xm = XM;
     801       60228 :   return gen_ZpX_Dixon(F,V,qM,utoi(p),M,(void*) &e, _lift_lin, _lift_invl);
     802             : }
     803             : 
     804             : static GEN
     805       25018 : lift_isogeny(GEN phi, GEN x0, long n, GEN Xm, GEN T, GEN sqx, GEN Tp, ulong p)
     806             : {
     807             :   struct _lift_iso d;
     808       25018 :   d.phi=phi;
     809       25018 :   d.Xm=Xm; d.T=T;
     810       25018 :   d.sqx=sqx; d.Tp=Tp; d.p=p;
     811       25018 :   return gen_ZpX_Newton(x0, utoi(p), n,(void*)&d, _lift_iter, _lift_invd);
     812             : }
     813             : 
     814             : static GEN
     815       25004 : getc2(GEN act, GEN X, GEN T, GEN q, ulong p, long N)
     816             : {
     817       25004 :   GEN A1 = RgV_to_RgX(gel(act,1),0), A2 =  RgV_to_RgX(gel(act,2),0);
     818       25004 :   long n = brent_kung_optpow(maxss(degpol(A1),degpol(A2)),2,1);
     819       25004 :   GEN xp = FpXQ_powers(X,n,T,q);
     820       25004 :   GEN P  = FpX_FpXQV_eval(A1, xp, T, q);
     821       25004 :   GEN Q  = FpX_FpXQV_eval(A2, xp, T, q);
     822       25004 :   return ZpXQ_div(P, Q, T, q, utoi(p), N);
     823             : }
     824             : 
     825             : struct _ZpXQ_norm
     826             : {
     827             :   long n;
     828             :   GEN T, p;
     829             : };
     830             : 
     831             : static GEN
     832       32774 : ZpXQ_norm_mul(void *E, GEN x, GEN y)
     833             : {
     834       32774 :   struct _ZpXQ_norm *D = (struct _ZpXQ_norm*)E;
     835       32774 :   GEN P = gel(x,1), Q = gel(y,1);
     836       32774 :   long a = mael(x,2,1), b = mael(y,2,1);
     837       32774 :   retmkvec2(FpXQ_mul(P,ZpXQ_frob_cyc(Q, D->T, D->p, a), D->T, D->p),
     838             :             mkvecsmall((a*b)%D->n));
     839             : }
     840             : 
     841             : static GEN
     842       22680 : ZpXQ_norm_sqr(void *E, GEN x)
     843             : {
     844       22680 :   return ZpXQ_norm_mul(E, x, x);
     845             : }
     846             : 
     847             : /* Assume T = Phi_(n) and n prime */
     848             : GEN
     849       11326 : ZpXQ_norm_pcyc(GEN x, GEN T, GEN q, GEN p)
     850             : {
     851             :   GEN z;
     852             :   struct _ZpXQ_norm D;
     853       11326 :   long d = get_FpX_degree(T);
     854       11326 :   D.T = T; D.p = q; D.n = d+1;
     855       11326 :   if (d==1) return ZX_copy(x);
     856       11326 :   z = mkvec2(x,mkvecsmall(p[2]));
     857       11326 :   z = gen_powu(z,d,(void*)&D,ZpXQ_norm_sqr,ZpXQ_norm_mul);
     858       11326 :   return gmael(z,1,2);
     859             : }
     860             : 
     861             : /* Assume T = Phi_(n) and n prime */
     862             : static GEN
     863       11088 : ZpXQ_sqrtnorm_pcyc(GEN x, GEN T, GEN q, GEN p, long e)
     864             : {
     865       11088 :   GEN z = ZpXQ_norm_pcyc(x, T, q, p);
     866       11088 :   return Zp_sqrtlift(z,Fp_sqrt(z,p),p,e);
     867             : }
     868             : 
     869             : /* Assume a = 1 [p], return the square root of the norm */
     870             : static GEN
     871       13930 : ZpXQ_sqrtnorm(GEN a, GEN T, GEN q, GEN p, long e)
     872             : {
     873       13930 :   GEN s = Fp_div(FpXQ_trace(ZpXQ_log(a, T, p, e), T, q), gen_2, q);
     874       13930 :   return modii(gel(Qp_exp(cvtop(s, p, e-1)),4), q);
     875             : }
     876             : 
     877             : struct _teich_lin
     878             : {
     879             :   ulong p;
     880             :   GEN sqx, Tp;
     881             :   long m;
     882             : };
     883             : 
     884             : static GEN
     885       29302 : _teich_invl(void *E, GEN x)
     886             : {
     887       29302 :   struct _teich_lin *d = (struct _teich_lin *) E;
     888       29302 :   ulong p = d->p;
     889       29302 :   GEN T = d->Tp;
     890       29302 :   return Flx_to_ZX(Flxq_lroot_fast(ZX_to_Flx(x, p), d->sqx, T, p));
     891             : }
     892             : 
     893             : static GEN
     894        8806 : _teich_lin(void *E, GEN F, GEN x2, GEN q)
     895             : {
     896        8806 :   struct _teich_lin *d = (struct _teich_lin *) E;
     897        8806 :   pari_sp av = avma;
     898        8806 :   GEN T = gel(F,2), Xm = gel(F,3);
     899        8806 :   GEN y2  = ZpXQ_frob(x2, Xm, T, q, d->p);
     900        8806 :   GEN lin = FpX_sub(y2, ZX_mulu(ZX_mul(gel(F,1), x2), d->p), q);
     901        8806 :   return gerepileupto(av, FpX_rem(lin, T, q));
     902             : }
     903             : 
     904             : struct _teich_iso
     905             : {
     906             :   GEN Xm, T;
     907             :   GEN sqx, Tp;
     908             :   ulong p;
     909             : };
     910             : 
     911             : static GEN
     912       20496 : _teich_iter(void *E, GEN x2, GEN q)
     913             : {
     914       20496 :   struct _teich_iso *d = (struct _teich_iso *) E;
     915       20496 :   ulong p = d->p;
     916       20496 :   GEN TN = FpXT_red(d->T, q), XN = FpXV_red(d->Xm, q);
     917       20496 :   GEN y2 = ZpXQ_frob(x2, XN, TN, q, d->p);
     918       20496 :   GEN x1 = FpXQ_powu(x2, p-1, TN, q);
     919       20496 :   GEN xp = FpXQ_mul(x2, x1, TN, q);
     920       20496 :   GEN V = FpX_sub(y2,xp,q);
     921       20496 :   return mkvec2(V,x1);
     922             : }
     923             : 
     924             : static GEN
     925       20496 : _teich_invd(void *E, GEN V, GEN v, GEN qM, long M)
     926             : {
     927       20496 :   struct _teich_iso *d = (struct _teich_iso *) E;
     928             :   struct _teich_lin e;
     929       20496 :   ulong p = d->p;
     930       20496 :   GEN TM = FpXT_red(d->T, qM), XM = FpXV_red(d->Xm, qM);
     931       20496 :   GEN x1 = FpX_red(gel(v,2), qM);
     932       20496 :   GEN F = mkvec3(x1, TM, XM);
     933       20496 :   e.sqx = d->sqx; e.Tp = d->Tp; e.p=p;
     934       20496 :   return gen_ZpX_Dixon(F,V,qM,utoi(p),M,(void*) &e, _teich_lin, _teich_invl);
     935             : }
     936             : 
     937             : static GEN
     938       10213 : Teichmuller_lift(GEN x, GEN Xm, GEN T, GEN sqx, GEN Tp, ulong p, long N)
     939             : {
     940             :   struct _teich_iso d;
     941       10213 :   d.Xm = Xm; d.T = T; d.sqx = sqx; d.Tp = Tp; d.p = p;
     942       10213 :   return gen_ZpX_Newton(x,utoi(p), N,(void*)&d, _teich_iter, _teich_invd);
     943             : }
     944             : 
     945             : static GEN
     946       25018 : get_norm(GEN a4, GEN a6, GEN T, ulong p, long N)
     947             : {
     948       25018 :   long sv=T[1];
     949             :   GEN a;
     950       25018 :   if (p==3) a = gel(a4,1);
     951             :   else
     952             :   {
     953       10213 :     GEN P = mkpoln(4, pol1_Flx(sv), pol0_Flx(sv), a4, a6);
     954       10213 :     a = gel(FlxqX_powu(P,p>>1,T,p),2+p-1);
     955             :   }
     956       25018 :   return Zp_sqrtnlift(gen_1,subss(p,1),utoi(Flxq_norm(a,T,p)),utoi(p), N);
     957             : }
     958             : 
     959             : static GEN
     960       25004 : fill_pols(long n, const long *v, long m, const long *vn,
     961             :           const long *vd, GEN *act)
     962             : {
     963             :   long i, j;
     964       25004 :   long d = upowuu(n,12/(n-1));
     965       25004 :   GEN N, D, M = zeromatcopy(n+1,n+1);
     966       25004 :   gmael(M,1,n+1) = gen_1;
     967      120554 :   for(i=2;i<=n+1;i++)
     968      338443 :     for(j=i-1;j<=n;j++)
     969      242893 :       gmael(M,i,j) = mulis(powuu(d,i-2),v[j-i+1]);
     970       25004 :   N = cgetg(m+1,t_COL);
     971       25004 :   D = cgetg(m+1,t_COL);
     972      135345 :   for(i=1;i<=m;i++)
     973             :   {
     974      110341 :     gel(N,i) = stoi(*vn++);
     975      110341 :     gel(D,i) = stoi(*vd++);
     976             :   }
     977       25004 :   *act = mkmat2(N,D);
     978       25004 :   return M;
     979             : }
     980             : 
     981             : /*
     982             :   These polynomials were extracted from the ECHIDNA databases
     983             :   available at <http://echidna.maths.usyd.edu.au/echidna/>
     984             :   and computed by David R. Kohel.
     985             :   Return the matrix of the modular polynomial, set act to the parametrization,
     986             :   and set dj to the opposite of the supersingular j-invariant.
     987             : */
     988             : static GEN
     989       25004 : get_Kohel_polynomials(ulong p, GEN *act, long *dj)
     990             : {
     991       25004 :   const long mat3[] = {-1,-36,-270};
     992       25004 :   const long num3[] = {1,-483,-21141,-59049};
     993       25004 :   const long den3[] = {1,261, 4347, -6561};
     994       25004 :   const long mat5[] = {-1,-30,-315,-1300,-1575};
     995       25004 :   const long num5[] = {-1,490,20620,158750,78125};
     996       25004 :   const long den5[] = {-1,-254,-4124,-12250,3125};
     997       25004 :   const long mat7[] = {-1,-28,-322,-1904,-5915,-8624,-4018};
     998       25004 :   const long num7[] = {1,-485,-24058,-343833,-2021642,-4353013,-823543};
     999       25004 :   const long den7[] = {1,259,5894,49119,168406,166355,-16807};
    1000       25004 :   const long mat13[]= {-1,-26,-325,-2548,-13832,-54340,-157118,-333580,-509366,
    1001             :                        -534820,-354536,-124852,-15145};
    1002       25004 :   const long num13[]= {1,-487,-24056,-391463,-3396483,-18047328,-61622301,
    1003             :                        -133245853,-168395656,-95422301,-4826809};
    1004       25004 :   const long den13[]= {1,257,5896,60649,364629,1388256,3396483,5089019,4065464,
    1005             :                        1069939,-28561};
    1006       25004 :   switch(p)
    1007             :   {
    1008             :   case 3:
    1009       14805 :     *dj = 0;
    1010       14805 :     return fill_pols(3,mat3,4,num3,den3,act);
    1011             :   case 5:
    1012       10150 :     *dj = 0;
    1013       10150 :     return fill_pols(5,mat5,5,num5,den5,act);
    1014             :   case 7:
    1015          42 :     *dj = 1;
    1016          42 :     return fill_pols(7,mat7,7,num7,den7,act);
    1017             :   case 13:
    1018           7 :     *dj = 8;
    1019           7 :     return fill_pols(13,mat13,11,num13,den13,act);
    1020             :   }
    1021           0 :   *dj=0; *act = NULL;
    1022           0 :   return NULL;
    1023             : }
    1024             : 
    1025             : long
    1026       32161 : zx_is_pcyc(GEN T)
    1027             : {
    1028       32161 :   long i, n = degpol(T);
    1029       32161 :   if (!uisprime(n+1))
    1030       11539 :     return 0;
    1031       98994 :   for (i=0; i<=n; i++)
    1032       87668 :     if (T[i+2]!=1UL)
    1033        9296 :       return 0;
    1034       11326 :   return 1;
    1035             : }
    1036             : 
    1037             : static GEN
    1038       25004 : Flxq_ellcard_Kohel(GEN a4, GEN a6, GEN T, ulong p)
    1039             : {
    1040       25004 :   pari_sp av = avma, av2;
    1041             :   pari_timer ti;
    1042       25004 :   long n = get_Flx_degree(T), N = (n+4)/2, dj;
    1043       25004 :   GEN q = powuu(p, N);
    1044             :   GEN T2, Xm, s1, c2, t, lr;
    1045             :   GEN S1, sqx;
    1046             :   GEN Nc2, Np;
    1047       25004 :   GEN act, phi = get_Kohel_polynomials(p, &act, &dj);
    1048       25004 :   long ispcyc = zx_is_pcyc(get_Flx_mod(T));
    1049       25004 :   timer_start(&ti);
    1050       25004 :   if (!ispcyc)
    1051             :   {
    1052       13916 :     T2 = Flx_canonlift(get_Flx_mod(T),N,p);
    1053       13916 :     if (DEBUGLEVEL) timer_printf(&ti,"Teich");
    1054             :   } else
    1055       11088 :     T2 = Flx_to_ZX(get_Flx_mod(T));
    1056       25004 :   T2 = FpX_get_red(T2, q); T = ZXT_to_FlxT(T2, p);
    1057       25004 :   av2 = avma;
    1058       25004 :   if (DEBUGLEVEL) timer_printf(&ti,"Barrett");
    1059       25004 :   if (!ispcyc)
    1060             :   {
    1061       13916 :     Xm = FpXQ_powers(monomial(gen_1,n,get_FpX_var(T2)),p-1,T2,q);
    1062       13916 :     if (DEBUGLEVEL) timer_printf(&ti,"Xm");
    1063             :   } else
    1064       11088 :     Xm = cgetg(1,t_VEC);
    1065       25004 :   s1 = Flxq_inv(Flx_Fl_add(Flxq_ellj(a4,a6,T,p),dj, p),T,p);
    1066       25004 :   lr = Flxq_lroot(polx_Flx(get_Flx_var(T)), T, p);
    1067       25004 :   sqx = Flxq_powers(lr, p-1, T, p);
    1068       25004 :   S1 = lift_isogeny(phi, Flx_to_ZX(s1), N, Xm, T2, sqx, T ,p);
    1069       25004 :   if (DEBUGLEVEL) timer_printf(&ti,"Lift isogeny");
    1070       25004 :   c2 = getc2(act, S1, T2, q, p, N);
    1071       25004 :   if (DEBUGLEVEL) timer_printf(&ti,"c^2");
    1072       25004 :   if (p>3 && !ispcyc)
    1073             :   {
    1074       10199 :     GEN c2p = Flx_to_ZX(Flxq_inv(ZX_to_Flx(c2,p),T,p));
    1075       10199 :     GEN tc2 = Teichmuller_lift(c2p,Xm, T2,sqx,T,p,N);
    1076       10199 :     if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fq");
    1077       10199 :     c2 = FpX_rem(FpX_mul(tc2,c2,q),T2,q);
    1078             :   }
    1079       25004 :   c2 = gerepileupto(av2, c2);
    1080       25004 :   if (DEBUGLEVEL) timer_printf(&ti,"tc2");
    1081       25004 :   Nc2 = (ispcyc? ZpXQ_sqrtnorm_pcyc: ZpXQ_sqrtnorm)(c2, T2, q, utoi(p), N);
    1082       25004 :   if (DEBUGLEVEL) timer_printf(&ti,"Norm");
    1083       25004 :   Np = get_norm(a4,a6,T,p,N);
    1084       25004 :   if (p>3 && ispcyc)
    1085             :   {
    1086           0 :     GEN Ncpi =  utoi(Fl_inv(umodiu(Nc2,p), p));
    1087           0 :     GEN tNc2 = Zp_sqrtnlift(gen_1, subss(p,1), Ncpi, utoi(p),N);
    1088           0 :     if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fp");
    1089           0 :     Nc2 = Fp_mul(Nc2,tNc2,q);
    1090             :   }
    1091       25004 :   t = Fp_center(Fp_mul(Nc2,Np,q),q,shifti(q,-1));
    1092       25004 :   return gerepileupto(av, subii(addis(powuu(p,n),1),t));
    1093             : }
    1094             : 
    1095             : static void
    1096          14 : liftcurve(GEN J, GEN T, GEN q, ulong p, long N, GEN *A4, GEN *A6)
    1097             : {
    1098          14 :   pari_sp av = avma;
    1099          14 :   GEN r = ZpXQ_inv(Z_ZX_sub(utoi(1728),J),T,utoi(p),N);
    1100          14 :   GEN g = FpXQ_mul(J,r,T,q);
    1101          14 :   *A4 = FpX_mulu(g,3,q);
    1102          14 :   *A6 = FpX_mulu(g,2,q);
    1103          14 :   gerepileall(av,2,A4,A6);
    1104          14 : }
    1105             : 
    1106             : static GEN
    1107          14 : getc5(GEN H, GEN A40, GEN A60, GEN A41, GEN A61, GEN T, GEN q, ulong p, long N)
    1108             : {
    1109          14 :   long d = lg(H)-1;
    1110          14 :   GEN s1 = gel(H,d-1), s2 = gel(H,d-2), s3 = d<5 ? pol_0(varn(T)): gel(H,d-3);
    1111          14 :   GEN s12 = FpXQ_sqr(s1,T,q);
    1112          14 :   GEN h2 = ZX_sub(ZX_shifti(s2,1),s12); /*2*s2-s1^2*/
    1113          14 :   GEN h3 = ZX_sub(FpXQ_mul(ZX_add(h2,s2),s1,T,q),ZX_mulu(s3,3));
    1114             :                                         /*3*s2*s1-s1^3-3s3*/
    1115          14 :   GEN alpha= ZX_sub(ZX_mulu(h2,30), ZX_mulu(A40,5*p-6)); /* 30*h2+A40*(6-5*p)*/
    1116          14 :   GEN beta = ZX_sub(ZX_sub(ZX_mulu(FpXQ_mul(A40,s1,T,q),42),ZX_mulu(A60,14*p-15)),
    1117             :                     ZX_mulu(h3,70)); /* 42*A40*s1-A60*(14*p-15)-70*h3 */
    1118          14 :   GEN u2 = FpXQ_mul(FpXQ_mul(A41,beta,T,q),
    1119             :                     ZpXQ_inv(FpXQ_mul(A61,alpha,T,q),T,utoi(p),N),T,q);
    1120          14 :   return u2;
    1121             : }
    1122             : 
    1123             : static GEN
    1124          14 : ZpXQX_liftrootmod_vald(GEN f, GEN H, long v, GEN T, GEN p, long e)
    1125             : {
    1126          14 :   pari_sp av = avma, av2, lim;
    1127          14 :   GEN pv = p, q, qv, W, df, Tq, fr, dfr;
    1128             :   ulong mask;
    1129             :   pari_timer ti;
    1130          14 :   if (e <= v+1) return H;
    1131          14 :   df = RgX_deriv(f);
    1132          14 :   if (v) { pv = powiu(p,v); qv = mulii(pv,p); df = ZXX_Z_divexact(df, pv); }
    1133           0 :   else qv = p;
    1134          14 :   mask = quadratic_prec_mask(e-v);
    1135          14 :   Tq = FpXT_red(T, qv); dfr = FpXQX_red(df, Tq, p);
    1136          14 :   if (DEBUGLEVEL) timer_start(&ti);
    1137          14 :   W = FpXQXQ_inv(FpXQX_rem(dfr, H, Tq, p), H, Tq, p); /* 1/f'(a) mod (T,p) */
    1138          14 :   if (DEBUGLEVEL) timer_printf(&ti,"FpXQXQ_inv");
    1139          14 :   q = p;
    1140          14 :   av2 = avma; lim = stack_lim(av2, 2);
    1141             :   for (;;)
    1142             :   {
    1143             :     GEN u, fa, qv, q2v, Tq2, fadH;
    1144          56 :     GEN H2 = H, q2 = q;
    1145          56 :     q = sqri(q);
    1146          56 :     if (mask & 1) q = diviiexact(q,p);
    1147          56 :     mask >>= 1;
    1148          56 :     if (v) { qv = mulii(q, pv); q2v = mulii(q2, pv); }
    1149           0 :     else { qv = q; q2v = q2; }
    1150          56 :     Tq2 = FpXT_red(T, q2v); Tq = FpXT_red(T, qv);
    1151          56 :     fr = FpXQX_red(f, Tq, qv);
    1152          56 :     fa = FpXQX_rem(fr, H, Tq, qv);
    1153          56 :     fa = ZXX_Z_divexact(fa, q2v);
    1154          56 :     fadH = FpXQXQ_mul(RgX_deriv(H),fa,H,Tq2,q2);
    1155          56 :     H = FpXX_add(H, gmul(FpXQXQ_mul(W, fadH, H, Tq2, q2v), q2), qv);
    1156          56 :     if (mask == 1) return gerepileupto(av, H);
    1157          42 :     dfr = FpXQX_rem(FpXQX_red(df, Tq, q),H,Tq,q);
    1158          42 :     u = ZXX_Z_divexact(ZXX_Z_add_shallow(FpXQXQ_mul(W,dfr,H,Tq,q),gen_m1),q2);
    1159          42 :     W = gsub(W,gmul(FpXQXQ_mul(u,W,H2,Tq2,q2),q2));
    1160          42 :     if (low_stack(lim, stack_lim(av2,2)))
    1161             :     {
    1162           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"ZpXQX_liftroot, e = %ld", e);
    1163           0 :       gerepileall(av2, 3, &H, &W, &q);
    1164             :     }
    1165          42 :   }
    1166             : }
    1167             : 
    1168             : static GEN
    1169          14 : get_H1(GEN A41, GEN A61, GEN T2, ulong p)
    1170             : {
    1171          14 :   GEN q = utoi(p), T = FpXT_red(T2,q);
    1172          14 :   GEN pol = FpXQ_elldivpol(FpX_red(A41,q),FpX_red(A61,q),p,T,q);
    1173          14 :   return FpXQX_normalize(RgX_deflate(pol,p),T,q);
    1174             : }
    1175             : 
    1176             : static GEN
    1177          14 : Flxq_ellcard_Harley(GEN a4, GEN a6, GEN T, ulong p)
    1178             : {
    1179          14 :   pari_sp av = avma, av2;
    1180             :   pari_timer ti;
    1181          14 :   long n = get_Flx_degree(T), N = (n+5)/2;
    1182          14 :   GEN q = powuu(p, N);
    1183             :   GEN T2, j, t;
    1184             :   GEN J1,A40,A41,A60,A61, sqx,Xm;
    1185             :   GEN pol, h1, H;
    1186             :   GEN c2, tc2, c2p, Nc2, Np;
    1187          14 :   long ispcyc = zx_is_pcyc(get_Flx_mod(T));
    1188          14 :   timer_start(&ti);
    1189          14 :   if (!ispcyc)
    1190             :   {
    1191          14 :     T2 = Flx_canonlift(get_Flx_mod(T),N,p);
    1192          14 :     if (DEBUGLEVEL) timer_printf(&ti,"Teich");
    1193             :   } else
    1194           0 :     T2 = Flx_to_ZX(get_Flx_mod(T));
    1195          14 :   T2 = FpX_get_red(T2, q); T = ZXT_to_FlxT(T2, p);
    1196          14 :   av2 = avma;
    1197          14 :   if (DEBUGLEVEL) timer_printf(&ti,"Barrett");
    1198          14 :   if (!ispcyc)
    1199             :   {
    1200          14 :     Xm = FpXQ_powers(monomial(gen_1,n,get_FpX_var(T2)),p-1,T2,q);
    1201          14 :     if (DEBUGLEVEL) timer_printf(&ti,"Xm");
    1202             :   } else
    1203           0 :     Xm = cgetg(1,t_VEC);
    1204          14 :   if (DEBUGLEVEL) timer_printf(&ti,"Xm");
    1205          14 :   j = Flxq_ellj(a4,a6,T,p);
    1206          14 :   sqx = Flxq_powers(Flxq_lroot(polx_Flx(T[1]), T, p), p-1, T, p);
    1207          14 :   J1 = lift_isogeny(polmodular_ZM(p, 0), Flx_to_ZX(j), N, Xm, T2,sqx,T,p);
    1208          14 :   if (DEBUGLEVEL) timer_printf(&ti,"Lift isogeny");
    1209          14 :   liftcurve(J1,T2,q,p,N,&A41,&A61);
    1210          14 :   A40 = ZpXQ_frob(A41, Xm, T2, q, p);
    1211          14 :   A60 = ZpXQ_frob(A61, Xm, T2, q, p);
    1212          14 :   if (DEBUGLEVEL) timer_printf(&ti,"liftcurve");
    1213          14 :   pol = FpXQ_elldivpol(A40,A60,p,T2,q);
    1214          14 :   if (DEBUGLEVEL) timer_printf(&ti,"p-division");
    1215          14 :   h1 = get_H1(A41,A61,T2,p);
    1216          14 :   H = ZpXQX_liftrootmod_vald(pol,h1,1,T2,utoi(p),N);
    1217          14 :   q = diviuexact(q,p); N--;
    1218          14 :   if (DEBUGLEVEL) timer_printf(&ti,"kernel");
    1219          14 :   c2 = getc5(H,A40,A60,A41,A61,T2,q,p,N);
    1220          14 :   if (DEBUGLEVEL) timer_printf(&ti,"c^2");
    1221          14 :   if (!ispcyc)
    1222             :   {
    1223          14 :     c2p = Flx_to_ZX(Flxq_inv(ZX_to_Flx(c2,p),T,p));
    1224          14 :     tc2 = Teichmuller_lift(c2p,Xm, T2,sqx,T,p,N);
    1225          14 :     if (DEBUGLEVEL) timer_printf(&ti,"teichmuller");
    1226          14 :     c2 = FpX_rem(FpX_mul(tc2,c2,q),T2,q);
    1227             :   }
    1228          14 :   c2 = gerepileupto(av2, c2);
    1229          14 :   q = powuu(p, N);
    1230          14 :   Nc2 = (ispcyc? ZpXQ_sqrtnorm_pcyc: ZpXQ_sqrtnorm)(c2, T2, q, utoi(p), N);
    1231          14 :   if (DEBUGLEVEL) timer_printf(&ti,"Norm");
    1232          14 :   Np = get_norm(a4,a6,T,p,N);
    1233          14 :   if (ispcyc)
    1234             :   {
    1235           0 :     GEN Ncpi = utoi(Fl_inv(umodiu(Nc2,p), p));
    1236           0 :     GEN tNc2 = Zp_sqrtnlift(gen_1, subss(p,1), Ncpi, utoi(p), N);
    1237           0 :     if (DEBUGLEVEL) timer_printf(&ti,"Teichmuller/Fp");
    1238           0 :     Nc2 = Fp_mul(Nc2,tNc2,q);
    1239             :   }
    1240          14 :   t = Fp_center(Fp_mul(Nc2,Np,q),q,shifti(q,-1));
    1241          14 :   return gerepileupto(av, subii(addis(powuu(p,n),1),t));
    1242             : }
    1243             : 
    1244             : /***************************************************************************/
    1245             : /*                                                                         */
    1246             : /*                          Shanks Mestre                                  */
    1247             : /*                                                                         */
    1248             : /***************************************************************************/
    1249             : 
    1250             : /* Return the lift of a (mod b), which is closest to h */
    1251             : static GEN
    1252        1988 : closest_lift(GEN a, GEN b, GEN h)
    1253             : {
    1254        1988 :   return addii(a, mulii(b, diviiround(subii(h,a), b)));
    1255             : }
    1256             : 
    1257             : static GEN
    1258        1043 : FlxqE_find_order(GEN f, GEN h, GEN bound, GEN B, GEN a4, GEN T, ulong p)
    1259             : {
    1260        1043 :   pari_sp av = avma, av1;
    1261             :   pari_timer Ti;
    1262        1043 :   long s = itos( gceil(gsqrt(gdiv(bound,B),DEFAULTPREC)) ) >> 1;
    1263             :   GEN tx, ti;
    1264        1043 :   GEN fh = FlxqE_mul(f, h, a4, T, p);
    1265        1043 :   GEN F, P = fh, fg;
    1266             :   long i;
    1267        1043 :   if (DEBUGLEVEL >= 6) timer_start(&Ti);
    1268        1043 :   if (ell_is_inf(fh)) return h;
    1269         952 :   F = FlxqE_mul(f, B, a4, T, p);
    1270         952 :   if (s < 3)
    1271             :   { /* we're nearly done: naive search */
    1272         119 :     GEN Q = P;
    1273         385 :     for (i=1;; i++)
    1274             :     {
    1275         385 :       P = FlxqE_add(P, F, a4, T, p); /* h.f + i.F */
    1276         385 :       if (ell_is_inf(P)) return gerepileupto(av, addii(h, mului(i,B)));
    1277         336 :       Q = FlxqE_sub(Q, F, a4, T, p); /* h.f - i.F */
    1278         336 :       if (ell_is_inf(Q)) return gerepileupto(av, subii(h, mului(i,B)));
    1279         266 :     }
    1280             :   }
    1281         833 :   tx = cgetg(s+1,t_VECSMALL);
    1282             :   /* Baby Step/Giant Step */
    1283         833 :   av1 = avma;
    1284        4697 :   for (i=1; i<=s; i++)
    1285             :   { /* baby steps */
    1286        4046 :     tx[i] = hash_GEN(gel(P, 1));
    1287        4046 :     P = FlxqE_add(P, F, a4, T, p); /* h.f + i.F */
    1288        4046 :     if (ell_is_inf(P)) return gerepileupto(av, addii(h, mului(i,B)));
    1289        3864 :     if (gc_needed(av1,3))
    1290             :     {
    1291           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"[Flxq_ellcard] baby steps, i=%ld",i);
    1292           0 :       P = gerepileupto(av1,P);
    1293             :     }
    1294             :   }
    1295         651 :   if (DEBUGLEVEL >= 6) timer_printf(&Ti, "[Flxq_ellcard] baby steps, s = %ld",s);
    1296             :   /* giant steps: fg = s.F */
    1297         651 :   fg = gerepileupto(av1, FlxqE_sub(P, fh, a4, T, p));
    1298         651 :   if (ell_is_inf(fg)) return gerepileupto(av,mului(s,B));
    1299         651 :   ti = vecsmall_indexsort(tx); /* = permutation sorting tx */
    1300         651 :   tx = perm_mul(tx,ti);
    1301         651 :   if (DEBUGLEVEL >= 6) timer_printf(&Ti, "[Flxq_ellcard] sorting");
    1302         651 :   av1 = avma;
    1303        3129 :   for (P=fg, i=1; ; i++)
    1304             :   {
    1305        3129 :     long k = hash_GEN(gel(P,1));
    1306        3129 :     long r = zv_search(tx, k);
    1307        3129 :     if (r)
    1308             :     {
    1309         651 :       while (r && tx[r] == k) r--;
    1310         651 :       for (r++; r <= s && tx[r] == k; r++)
    1311             :       {
    1312         651 :         long j = ti[r]-1;
    1313         651 :         GEN Q = FlxqE_add(FlxqE_mul(F, stoi(j), a4, T, p), fh, a4, T, p);
    1314         651 :         if (DEBUGLEVEL >= 6)
    1315           0 :           timer_printf(&Ti, "[Flxq_ellcard] giant steps, i = %ld",i);
    1316         651 :         if (Flx_equal(gel(P,1), gel(Q,1)))
    1317             :         {
    1318         651 :           if (Flx_equal(gel(P,2), gel(Q,2))) i = -i;
    1319         651 :           return gerepileupto(av,addii(h, mulii(addis(mulss(s,i), j), B)));
    1320             :         }
    1321             :       }
    1322             :     }
    1323        2478 :     P = FlxqE_add(P,fg,a4,T,p);
    1324        2478 :     if (gc_needed(av1,3))
    1325             :     {
    1326           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"[Flxq_ellcard] giants steps, i=%ld",i);
    1327           0 :       P = gerepileupto(av1,P);
    1328             :     }
    1329        2478 :   }
    1330             : }
    1331             : 
    1332             : static void
    1333       33145 : Flx_next(GEN t, ulong p)
    1334             : {
    1335             :   long i;
    1336       40628 :   for(i=2;;i++)
    1337       40628 :     if (uel(t,i)==p-1)
    1338        7483 :       t[i]=0;
    1339             :     else
    1340             :     {
    1341       33145 :       t[i]++;
    1342       33145 :       break;
    1343        7483 :     }
    1344       33145 : }
    1345             : 
    1346             : static void
    1347       33145 : Flx_renormalize_ip(GEN x, long lx)
    1348             : {
    1349             :   long i;
    1350       40628 :   for (i = lx-1; i>=2; i--)
    1351       38381 :     if (x[i]) break;
    1352       33145 :   setlg(x, i+1);
    1353       33145 : }
    1354             : 
    1355             : static ulong
    1356        1547 : F3xq_ellcard_naive(GEN a2, GEN a6, GEN T)
    1357             : {
    1358        1547 :   pari_sp av = avma;
    1359        1547 :   long i, d = get_Flx_degree(T), lx = d+2;
    1360        1547 :   long q = upowuu(3, d), a;
    1361        1547 :   GEN x = zero_zv(lx); x[1] = get_Flx_var(T);
    1362        8624 :   for(a=1, i=0; i<q; i++)
    1363             :   {
    1364             :     GEN rhs;
    1365        7077 :     Flx_renormalize_ip(x, lx);
    1366        7077 :     rhs = Flx_add(Flxq_mul(Flxq_sqr(x, T, 3), Flx_add(x, a2, 3), T, 3), a6, 3);
    1367        7077 :     if (!lgpol(rhs)) a++; else if (Flxq_issquare(rhs, T, 3)) a+=2;
    1368        7077 :     Flx_next(x, 3);
    1369             :   }
    1370        1547 :   avma = av;
    1371        1547 :   return a;
    1372             : }
    1373             : 
    1374             : static ulong
    1375         700 : Flxq_ellcard_naive(GEN a4, GEN a6, GEN T, ulong p)
    1376             : {
    1377         700 :   pari_sp av = avma;
    1378         700 :   long i, d = get_Flx_degree(T), lx = d+2;
    1379         700 :   long q = upowuu(p, d), a;
    1380         700 :   GEN x = zero_zv(lx); x[1] = get_Flx_var(T);
    1381       26768 :   for(a=1, i=0; i<q; i++)
    1382             :   {
    1383             :     GEN x2, rhs;
    1384       26068 :     Flx_renormalize_ip(x, lx);
    1385       26068 :     x2  = Flxq_sqr(x, T, p);
    1386       26068 :     rhs = Flx_add(Flxq_mul(x, Flx_add(x2, a4, p), T, p), a6, p);
    1387       26068 :     if (!lgpol(rhs)) a++; else if (Flxq_issquare(rhs,T,p)) a+=2;
    1388       26068 :     Flx_next(x,p);
    1389             :   }
    1390         700 :   avma = av;
    1391         700 :   return a;
    1392             : }
    1393             : 
    1394             : /* assume T irreducible mod p, m = (q-1)/(p-1) */
    1395             : static int
    1396        1860 : Flxq_kronecker(GEN x, GEN m, GEN T, ulong p)
    1397             : {
    1398             :   pari_sp av;
    1399             :   ulong z;
    1400        1860 :   if (lgpol(x) == 0) return 0;
    1401        1860 :   av = avma;
    1402        1860 :   z = Flxq_pow(x, m, T, p)[2];
    1403        1860 :   avma = av; return krouu(z, p);
    1404             : }
    1405             : 
    1406             : /* Find x such that kronecker(u = x^3+a4x+a6, p) is KRO.
    1407             :  * Return point [x*u,u^2] on E (KRO=1) / E^twist (KRO=-1) */
    1408             : static GEN
    1409        1860 : Flxq_ellpoint(long KRO, GEN a4, GEN a6, GEN m, long n, long vn, GEN T, ulong p)
    1410             : {
    1411             :   for(;;)
    1412             :   {
    1413        1860 :     GEN x = random_Flx(n,vn,p);
    1414        1860 :     GEN u = Flx_add(a6, Flxq_mul(Flx_add(a4, Flxq_sqr(x,T,p), p), x, T,p), p);
    1415        1860 :     if (Flxq_kronecker(u, m,T,p) == KRO)
    1416        2086 :       return mkvec2(Flxq_mul(u,x, T,p), Flxq_sqr(u, T,p));
    1417         817 :   }
    1418             : }
    1419             : 
    1420             : static GEN
    1421         945 : Flxq_ellcard_Shanks(GEN a4, GEN a6, GEN q, GEN T, ulong p)
    1422             : {
    1423         945 :   pari_sp av = avma;
    1424         945 :   long vn = get_Flx_var(T), n = get_Flx_degree(T), KRO = -1;
    1425             :   GEN h,f, ta4, A, B, m;
    1426         945 :   GEN q1p = addsi(1, q), q2p = shifti(q1p, 1);
    1427         945 :   GEN bound = addis(sqrti(gmul2n(q,4)), 1); /* ceil( 4sqrt(q) ) */
    1428             :   /* once #E(Flxq) is know mod B >= bound, it is completely determined */
    1429             :   /* how many 2-torsion points ? */
    1430         945 :   switch(FlxqX_nbroots(mkpoln(4, pol1_Flx(vn), pol0_Flx(vn), a4, a6), T, p))
    1431             :   {
    1432         182 :   case 3:  A = gen_0; B = utoipos(4); break;
    1433         364 :   case 1:  A = gen_0; B = gen_2; break;
    1434         399 :   default: A = gen_1; B = gen_2; break; /* 0 */
    1435             :   }
    1436         945 :   m = diviuexact(subiu(powuu(p,n), 1), p-1);
    1437             :   for(;;)
    1438             :   {
    1439        1043 :     h = closest_lift(A, B, q1p);
    1440             :     /* [ux, u^2] is on E_u: y^2 = x^3 + c4 u^2 x + c6 u^3
    1441             :      * E_u isomorphic to E (resp. E') iff KRO = 1 (resp. -1)
    1442             :      * #E(F_p) = p+1 - a_p, #E'(F_p) = p+1 + a_p
    1443             :      *
    1444             :      * #E_u(Flxq) = A (mod B),  h is close to #E_u(Flxq) */
    1445        1043 :     KRO = -KRO;
    1446        1043 :     f = Flxq_ellpoint(KRO, a4,a6, m,n,vn, T,p);
    1447             : 
    1448        1043 :     ta4 = Flxq_mul(a4, gel(f,2), T, p); /* a4 for E_u */
    1449        1043 :     h = FlxqE_find_order(f, h, bound, B, ta4,T,p);
    1450        1043 :     h = FlxqE_order(f, h, ta4, T, p);
    1451             :     /* h | #E_u(Flxq) = A (mod B) */
    1452        1043 :     A = Z_chinese_all(A, gen_0, B, h, &B);
    1453        1043 :     if (cmpii(B, bound) >= 0) break;
    1454             :     /* not done, update A mod B for the _next_ curve, isomorphic to
    1455             :      * the quadratic twist of this one */
    1456          98 :     A = remii(subii(q2p,A), B); /* #E(Fq)+#E'(Fq) = 2q+2 */
    1457          98 :   }
    1458         945 :   h = closest_lift(A, B, q1p);
    1459         945 :   return gerepileuptoint(av, KRO == 1? h: subii(q2p,h));
    1460             : }
    1461             : 
    1462             : static GEN
    1463       16352 : F3xq_ellcard(GEN a2, GEN a6, GEN T)
    1464             : {
    1465       16352 :   long n = get_Flx_degree(T);
    1466       16352 :   if (n <= 2)
    1467        1232 :     return utoi(F3xq_ellcard_naive(a2, a6, T));
    1468             :   else
    1469             :   {
    1470       15120 :     GEN q1 = addis(powuu(3, get_Flx_degree(T)), 1), t;
    1471       15120 :     GEN a = Flxq_div(a6,Flxq_powu(a2,3,T,3),T,3);
    1472       15120 :     if (Flx_equal1(Flxq_powu(a, 8, T, 3)))
    1473             :     {
    1474         315 :       GEN P = Flxq_minpoly(a,T,3);
    1475         315 :       long dP = degpol(P); /* dP <= 2 */
    1476         315 :       ulong q = upowuu(3,dP);
    1477         315 :       GEN A2 = pol1_Flx(P[1]), A6 = Flx_rem(polx_Flx(P[1]), P, 3);
    1478         315 :       long tP = q + 1 - F3xq_ellcard_naive(A2, A6, P);
    1479         315 :       t = elltrace_extension(stoi(tP), n/dP, utoi(q));
    1480         315 :       if (umodiu(t, 3)!=1) t = negi(t);
    1481         315 :       return Flx_equal1(a2) || Flxq_issquare(a2,T,3) ? subii(q1,t): addii(q1,t);
    1482             :     }
    1483       14805 :     else return Flxq_ellcard_Kohel(mkvec(a2), a6, T, 3);
    1484             :   }
    1485             : }
    1486             : 
    1487             : static GEN
    1488       10913 : Flxq_ellcard_Satoh(GEN a4, GEN a6, GEN j, GEN T, ulong p)
    1489             : {
    1490       10913 :   long n = get_Flx_degree(T);
    1491       10913 :   if (n <= 2)
    1492         420 :     return utoi(Flxq_ellcard_naive(a4, a6, T, p));
    1493             :   else
    1494             :   {
    1495       10493 :     GEN jp = Flxq_powu(j, p, T, p);
    1496       10493 :     GEN s = Flx_add(j, jp, p);
    1497       10493 :     if (degpol(s) <= 0)
    1498             :     { /* it is assumed j not in F_p */
    1499         280 :       GEN m = Flxq_mul(j, jp, T, p);
    1500         280 :       if (degpol(m) <= 0)
    1501             :       {
    1502         280 :         GEN q = sqru(p);
    1503         280 :         GEN q1 = addis(powuu(p, get_Flx_degree(T)), 1);
    1504         280 :         GEN sk = Flx_Fl_add(Flx_neg(j, p), 1728%p, p);
    1505         280 :         GEN sA4 = Flx_triple(Flxq_mul(sk, j, T, p), p);
    1506         280 :         GEN u = Flxq_div(a4, sA4, T, p);
    1507         280 :         ulong ns = lgpol(s) ? Fl_neg(s[2], p): 0UL;
    1508         280 :         GEN P = mkvecsmall4(T[1], m[2], ns, 1L);
    1509             :         GEN A4, A6, t, tP;
    1510         280 :         Flxq_ellj_to_a4a6(polx_Flx(T[1]), P, p, &A4, &A6);
    1511         280 :         tP = addis(q, 1 - Flxq_ellcard_naive(A4, A6, P, p));
    1512         280 :         t = elltrace_extension(tP, n>>1, q);
    1513         280 :         return Flxq_is2npower(u, 2, T, p) ? subii(q1,t): addii(q1,t);
    1514             :       }
    1515             :     }
    1516       10213 :     if (p<=7 || p==13 ) return Flxq_ellcard_Kohel(a4, a6, T, p);
    1517          14 :     else return Flxq_ellcard_Harley(a4, a6, T, p);
    1518             :   }
    1519             : }
    1520             : 
    1521             : static GEN
    1522           0 : Flxq_ellcard_Kedlaya(GEN a4, GEN a6, GEN T, ulong p)
    1523             : {
    1524           0 :   pari_sp av = avma;
    1525           0 :   GEN H = mkpoln(4, gen_1, gen_0, Flx_to_ZX(a4), Flx_to_ZX(a6));
    1526           0 :   GEN Tp = Flx_to_ZX(get_Flx_mod(T));
    1527           0 :   long n = degpol(Tp), e = ((p < 16 ? n+1: n)>>1)+1;
    1528           0 :   GEN M = ZlXQX_hyperellpadicfrobenius(H, Tp, p, e);
    1529           0 :   GEN N = ZpXQM_prodFrobenius(M, Tp, utoi(p), e);
    1530           0 :   GEN q = powuu(p, e);
    1531           0 :   GEN tp = Fq_add(gcoeff(N,1,1), gcoeff(N,2,2), Tp, q);
    1532           0 :   GEN t = Fp_center(typ(tp)==t_INT ? tp: leading_coeff(tp), q, shifti(q,-1));
    1533           0 :   return gerepileupto(av, subii(addis(powuu(p, n), 1), t));
    1534             : }
    1535             : 
    1536             : GEN
    1537       51715 : Flxq_ellj(GEN a4, GEN a6, GEN T, ulong p)
    1538             : {
    1539       51715 :   pari_sp av=avma;
    1540       51715 :   if (p==3)
    1541             :   {
    1542             :     GEN J;
    1543       14805 :     if (typ(a4)!=t_VEC) return pol0_Flx(get_Flx_var(T));
    1544       14805 :     J = Flxq_div(Flxq_powu(gel(a4,1),3, T, p),Flx_neg(a6,p), T, p);
    1545       14805 :     return gerepileuptoleaf(av, J);
    1546             :   }
    1547             :   else
    1548             :   {
    1549       36910 :     pari_sp av=avma;
    1550       36910 :     GEN a43 = Flxq_mul(a4,Flxq_sqr(a4,T,p),T,p);
    1551       36910 :     GEN a62 = Flxq_sqr(a6,T,p);
    1552       36910 :     GEN num = Flx_mulu(a43,6912,p);
    1553       36910 :     GEN den = Flx_add(Flx_mulu(a43,4,p),Flx_mulu(a62,27,p),p);
    1554       36910 :     return gerepileuptoleaf(av, Flxq_div(num, den, T, p));
    1555             :   }
    1556             : }
    1557             : 
    1558             : void
    1559         280 : Flxq_ellj_to_a4a6(GEN j, GEN T, ulong p, GEN *pt_a4, GEN *pt_a6)
    1560             : {
    1561         280 :   ulong zagier = 1728 % p;
    1562         280 :   if (lgpol(j)==0)
    1563           0 :     { *pt_a4 = pol0_Flx(T[1]); *pt_a6 =pol1_Flx(T[1]); }
    1564         280 :   else if (lgpol(j)==1 && uel(j,2) == zagier)
    1565           0 :     { *pt_a4 = pol1_Flx(T[1]); *pt_a6 =pol0_Flx(T[1]); }
    1566             :   else
    1567             :   {
    1568         280 :     GEN k = Flx_Fl_add(Flx_neg(j, p), zagier, p);
    1569         280 :     GEN kj = Flxq_mul(k, j, T, p);
    1570         280 :     GEN k2j = Flxq_mul(kj, k, T, p);
    1571         280 :     *pt_a4 = Flx_triple(kj, p);
    1572         280 :     *pt_a6 = Flx_double(k2j, p);
    1573             :   }
    1574         280 : }
    1575             : 
    1576             : static GEN
    1577        6027 : F3xq_ellcardj(GEN a4, GEN a6, GEN T, GEN q, long n)
    1578             : {
    1579        6027 :   const ulong p = 3;
    1580             :   ulong t;
    1581        6027 :   GEN q1 = addis(q,1);
    1582        6027 :   GEN na4 = Flx_neg(a4,p), ra4;
    1583        6027 :   if (!Flxq_issquare(na4,T,p))
    1584        3101 :     return q1;
    1585        2926 :   ra4 = Flxq_sqrt(na4,T,p);
    1586        2926 :   t = Flxq_trace(Flxq_div(a6,Flxq_mul(na4,ra4,T,p),T,p),T,p);
    1587        2926 :   if (n%2==1)
    1588             :   {
    1589             :     GEN q3;
    1590         854 :     if (t==0) return q1;
    1591         168 :     q3 = powuu(p,(n+1)>>1);
    1592         168 :     return (t==1)^(n%4==1) ? subii(q1,q3): addii(q1,q3);
    1593             :   }
    1594             :   else
    1595             :   {
    1596        2072 :     GEN q22, q2 = powuu(p,n>>1);
    1597        2072 :     GEN W = Flxq_pow(a4,shifti(q,-2),T,p);
    1598        2072 :     long s = (W[2]==1)^(n%4==2);
    1599        2072 :     if (t!=0) return s ? addii(q1,q2): subii(q1, q2);
    1600        2072 :     q22 = shifti(q2,1);
    1601        2072 :     return s ? subii(q1,q22):  addii(q1, q22);
    1602             :   }
    1603             : }
    1604             : 
    1605             : static GEN
    1606       14812 : Flxq_ellcardj(GEN a4, GEN a6, ulong j, GEN T, GEN q, ulong p, long n)
    1607             : {
    1608       14812 :   GEN q1 = addis(q,1);
    1609       14812 :   if (j==0)
    1610             :   {
    1611             :     ulong w;
    1612             :     GEN W, t, N;
    1613        5614 :     if (umodiu(q,6)!=1) return q1;
    1614        4207 :     N = Fp_ffellcard(gen_0,gen_1,q,n,utoi(p));
    1615        4207 :     t = subii(q1, N);
    1616        4207 :     W = Flxq_pow(a6,diviuexact(shifti(q,-1), 3),T,p);
    1617        4207 :     if (degpol(W)>0) /*p=5 mod 6*/
    1618        1372 :       return Flx_equal1(Flxq_powu(W,3,T,p)) ? addii(q1,shifti(t,-1)):
    1619         455 :                                               subii(q1,shifti(t,-1));
    1620        3290 :     w = W[2];
    1621        3290 :     if (w==1)   return N;
    1622        2590 :     if (w==p-1) return addii(q1,t);
    1623             :     else /*p=1 mod 6*/
    1624             :     {
    1625        1890 :       GEN u = shifti(t,-1), v = sqrtint(diviuexact(subii(q,sqri(u)),3));
    1626        1890 :       GEN a = addii(u,v), b = shifti(v,1);
    1627        1890 :       if (Fl_powu(w,3,p)==1)
    1628             :       {
    1629         945 :         if (Fl_add(umodiu(a,p),Fl_mul(w,umodiu(b,p),p),p)==0)
    1630         574 :           return subii(q1,subii(shifti(b,1),a));
    1631             :         else
    1632         371 :           return addii(q1,addii(a,b));
    1633             :       }
    1634             :       else
    1635             :       {
    1636         945 :         if (Fl_sub(umodiu(a,p),Fl_mul(w,umodiu(b,p),p),p)==0)
    1637         574 :           return subii(q1,subii(a,shifti(b,1)));
    1638             :         else
    1639         371 :           return subii(q1,addii(a,b));
    1640             :       }
    1641             :     }
    1642        9198 :   } else if (j==1728%p)
    1643             :   {
    1644             :     ulong w;
    1645             :     GEN W, N, t;
    1646        5614 :     if (mod4(q)==3) return q1;
    1647        4214 :     W = Flxq_pow(a4,shifti(q,-2),T,p);
    1648        4214 :     if (degpol(W)>0) return q1; /*p=3 mod 4*/
    1649        3500 :     w = W[2];
    1650        3500 :     N = Fp_ffellcard(gen_1,gen_0,q,n,utoi(p));
    1651        3500 :     if(w==1) return N;
    1652        2457 :     t = subii(q1, N);
    1653        2457 :     if(w==p-1) return addii(q1, t);
    1654             :     else /*p=1 mod 4*/
    1655             :     {
    1656        1400 :       GEN u = shifti(t,-1), v = sqrtint(subii(q,sqri(u)));
    1657        1400 :       if (Fl_add(umodiu(u,p),Fl_mul(w,umodiu(v,p),p),p)==0)
    1658         700 :         return subii(q1,shifti(v,1));
    1659             :       else
    1660         700 :         return addii(q1,shifti(v,1));
    1661             :     }
    1662             :   } else
    1663             :   {
    1664        3584 :     ulong g = Fl_div(j, Fl_sub(1728%p, j, p), p);
    1665        3584 :     GEN l = Flxq_div(Flx_triple(a6,p),Flx_double(a4,p),T,p);
    1666        3584 :     GEN N = Fp_ffellcard(utoi(Fl_triple(g,p)),utoi(Fl_double(g,p)),q,n,utoi(p));
    1667        3584 :     if (Flxq_issquare(l,T,p)) return N;
    1668        2184 :     return subii(shifti(q1,1),N);
    1669             :   }
    1670             : }
    1671             : 
    1672             : GEN
    1673       49286 : Flxq_ellcard(GEN a4, GEN a6, GEN T, ulong p)
    1674             : {
    1675       49286 :   pari_sp av = avma;
    1676       49286 :   long n = get_Flx_degree(T);
    1677       49286 :   GEN J, r, q = powuu(p,  n);
    1678       49286 :   if (typ(a4)==t_VEC)
    1679       16352 :     r = F3xq_ellcard(gel(a4,1), a6, T);
    1680       32934 :   else if (p==3)
    1681        6027 :     r = F3xq_ellcardj(a4, a6, T, q, n);
    1682       26907 :   else if (degpol(a4)<=0 && degpol(a6)<=0)
    1683         210 :     r = Fp_ffellcard(utoi(Flx_eval(a4,0,p)),utoi(Flx_eval(a6,0,p)),q,n,utoi(p));
    1684       26697 :   else if (degpol(J=Flxq_ellj(a4,a6,T,p))<=0)
    1685       14812 :     r = Flxq_ellcardj(a4,a6,lgpol(J)?J[2]:0,T,q,p,n);
    1686       11885 :   else if (p <= 7)
    1687       10843 :     r = Flxq_ellcard_Satoh(a4, a6, J, T, p);
    1688        1042 :   else if (cmpis(q,100)<0)
    1689           0 :     r = utoi(Flxq_ellcard_naive(a4, a6, T, p));
    1690        1042 :   else if (p == 13 || (7*p <= (ulong)10*n && (BITS_IN_LONG==64 || p <= 103)))
    1691          70 :     r = Flxq_ellcard_Satoh(a4, a6, J, T, p);
    1692         972 :   else if (p <= (ulong)2*n)
    1693           0 :     r = Flxq_ellcard_Kedlaya(a4, a6, T, p);
    1694         972 :   else if (expi(q)<=62)
    1695         945 :     r = Flxq_ellcard_Shanks(a4, a6, q, T, p);
    1696             :   else
    1697          27 :     r = Fq_ellcard_SEA(Flx_to_ZX(a4),Flx_to_ZX(a6),q,Flx_to_ZX(T),utoi(p),0);
    1698       49286 :   return gerepileuptoint(av, r);
    1699             : }

Generated by: LCOV version 1.11