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 - ellsea.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.18.0 lcov report (development 29785-cb4530e674) Lines: 1208 1260 95.9 %
Date: 2024-12-02 09:09:34 Functions: 94 97 96.9 %
Legend: Lines: hit not hit

          Line data    Source code
       1             : /* Copyright (C) 2008  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             : /* This file is a C version by Bill Allombert of the 'ellsea' GP package
      16             :  * whose copyright statement is as follows:
      17             : Authors:
      18             :   Christophe Doche   <cdoche@math.u-bordeaux.fr>
      19             :   Sylvain Duquesne <duquesne@math.u-bordeaux.fr>
      20             : 
      21             : Universite Bordeaux I, Laboratoire A2X
      22             : For the AREHCC project, see http://www.arehcc.com/
      23             : 
      24             : Contributors:
      25             :   Karim Belabas (code cleanup and package release, faster polynomial arithmetic)
      26             : 
      27             : 'ellsea' is free software; you can redistribute it and/or modify it under the
      28             : terms of the GNU General Public License as published by the Free Software
      29             : Foundation. It is distributed in the hope that it will be useful, but WITHOUT
      30             : ANY WARRANTY WHATSOEVER. */
      31             : 
      32             : /* Extension to non prime finite fields by Bill Allombert 2012 */
      33             : 
      34             : #include "pari.h"
      35             : #include "paripriv.h"
      36             : 
      37             : #define DEBUGLEVEL DEBUGLEVEL_ellsea
      38             : 
      39             : static THREAD GEN modular_eqn;
      40             : 
      41             : void
      42      326178 : pari_set_seadata(GEN mod)  { modular_eqn = mod; }
      43             : GEN
      44      324524 : pari_get_seadata(void)  { return modular_eqn; }
      45             : 
      46             : static char *
      47          91 : seadata_filename(ulong ell)
      48          91 : { return stack_sprintf("%s/seadata/sea%ld", pari_datadir, ell); }
      49             : 
      50             : static GEN
      51          91 : get_seadata(ulong ell)
      52             : {
      53          91 :   pari_sp av = avma;
      54             :   GEN eqn;
      55          91 :   char *s = seadata_filename(ell);
      56          91 :   pariFILE *F = pari_fopengz(s);
      57          91 :   if (!F) return NULL;
      58          35 :   if (ell) /* large single polynomial */
      59           7 :     eqn = gp_read_stream(F->file);
      60             :   else
      61             :   { /* table of polynomials of small level */
      62          28 :     eqn = gp_readvec_stream(F->file);
      63          28 :     modular_eqn = eqn = gclone(eqn);
      64          28 :     set_avma(av);
      65             :   }
      66          35 :   pari_fclose(F);
      67          35 :   return eqn;
      68             : }
      69             : 
      70             : /*Builds the modular equation corresponding to the vector list. Shallow */
      71             : static GEN
      72        9968 : list_to_pol(GEN list, long vx, long vy)
      73             : {
      74        9968 :   long i, l = lg(list);
      75        9968 :   GEN P = cgetg(l, t_VEC);
      76      205121 :   for (i = 1; i < l; i++)
      77             :   {
      78      195153 :     GEN L = gel(list,i);
      79      195153 :     if (typ(L) == t_VEC) L = RgV_to_RgX_reverse(L, vy);
      80      195153 :     gel(P, i) = L;
      81             :   }
      82        9968 :   return RgV_to_RgX_reverse(P, vx);
      83             : }
      84             : 
      85             : struct meqn {
      86             :   char type;
      87             :   GEN eq, eval;
      88             :   long vx,vy;
      89             : };
      90             : 
      91             : static GEN
      92       10024 : seadata_cache(ulong ell)
      93             : {
      94       10024 :   long n = uprimepi(ell)-1;
      95             :   GEN C;
      96       10024 :   if (!modular_eqn && !get_seadata(0))
      97          56 :     C = NULL;
      98        9968 :   else if (n && n < lg(modular_eqn))
      99        9961 :     C = gel(modular_eqn, n);
     100             :   else
     101           7 :     C = get_seadata(ell);
     102       10024 :   return C;
     103             : }
     104             : /* C = [prime level, type "A" or "C", pol. coeffs] */
     105             : static void
     106        9968 : seadata_parse(struct meqn *M, GEN C, long vx, long vy)
     107             : {
     108        9968 :   M->type = *GSTR(gel(C,2));
     109        9968 :   M->eq = list_to_pol(gel(C,3), vx, vy);
     110        9968 : }
     111             : static void
     112       10003 : get_modular_eqn(struct meqn *M, ulong ell, long vx, long vy)
     113             : {
     114       10003 :   GEN C = seadata_cache(ell);
     115       10003 :   M->vx = vx;
     116       10003 :   M->vy = vy;
     117       10003 :   M->eval = gen_0;
     118       10003 :   if (C) seadata_parse(M, C, vx, vy);
     119             :   else
     120             :   {
     121          56 :     M->type = 'J'; /* j^(1/3) for ell != 3, j for 3 */
     122          56 :     M->eq = polmodular_ZXX(ell, ell==3? 0: 5, vx, vy);
     123             :   }
     124       10003 : }
     125             : 
     126             : GEN
     127          35 : ellmodulareqn(long ell, long vx, long vy)
     128             : {
     129          35 :   pari_sp av = avma;
     130             :   struct meqn meqn;
     131             :   GEN C;
     132          35 :   if (vx < 0) vx = 0;
     133          35 :   if (vy < 0) vy = 1;
     134          35 :   if (varncmp(vx,vy) >= 0)
     135           7 :     pari_err_PRIORITY("ellmodulareqn", pol_x(vx), ">=", vy);
     136          28 :   if (ell < 2 || !uisprime(ell))
     137           7 :     pari_err_PRIME("ellmodulareqn (level)", stoi(ell));
     138          21 :   C = seadata_cache(ell);
     139          21 :   if (!C) pari_err_FILE("seadata file", seadata_filename(ell));
     140          21 :   seadata_parse(&meqn, C, vx, vy);
     141          21 :   return gerepilecopy(av, mkvec2(meqn.eq, meqn.type=='A'? gen_1: gen_0));
     142             : }
     143             : 
     144             : /***********************************************************************/
     145             : /**                                                                   **/
     146             : /**                           FqE_group                               **/
     147             : /**                                                                   **/
     148             : /***********************************************************************/
     149             : 
     150             : static GEN
     151         122 : Fq_to_Flx(GEN a4, GEN T, ulong p)
     152         122 : { return typ(a4)==t_INT ? Z_to_Flx(a4, p, get_Flx_var(T)): ZX_to_Flx(a4, p); }
     153             : 
     154             : /*FIXME: the name of the function does not quite match what it does*/
     155             : static const struct bb_group *
     156         980 : get_FqE_group(void ** pt_E, GEN a4, GEN a6, GEN T, GEN p)
     157             : {
     158         980 :   if (!T) return get_FpE_group(pt_E,a4,a6,p);
     159          77 :   else if (lgefint(p)==3)
     160             :   {
     161          61 :     ulong pp = uel(p,2);
     162          61 :     GEN Tp = ZXT_to_FlxT(T,pp);
     163          61 :     return get_FlxqE_group(pt_E, Fq_to_Flx(a4, Tp, pp), Fq_to_Flx(a6, Tp, pp),
     164             :                            Tp, pp);
     165             :   }
     166          16 :   return get_FpXQE_group(pt_E,a4,a6,T,p);
     167             : }
     168             : 
     169             : /***********************************************************************/
     170             : /**                                                                   **/
     171             : /**   Handle curves with CM by small order                            **/
     172             : /**                                                                   **/
     173             : /***********************************************************************/
     174             : 
     175             : /* l odd prime. Return the list of discriminants D such that
     176             :  *   polclass(D) | poldisc(polmodular(l)) */
     177             : static GEN
     178          14 : list_singular_discs(long l)
     179             : {
     180          14 :   const long _4l2 = 4*l*l;
     181             :   long v;
     182          14 :   GEN V = zero_F2v(_4l2);
     183             :   /* special cased for efficiency + remove factor l^2 from conductor */
     184          14 :   F2v_set(V, 4); /* v = 0 */
     185          14 :   F2v_set(V, 3); /* v = l */
     186        1232 :   for (v = 1; v < 2*l; v++)
     187        1218 :     if (v != l)
     188             :     { /* l does not divide _4l2 - v*v */
     189        1204 :       GEN F = factoru(_4l2 - v*v), P, E, c;
     190        1204 :       ulong d = coredisc2u_fact(F, -1, &P, &E);
     191             :       long i, lc;
     192        1204 :       c = divisorsu_fact(mkvec2(P,E));
     193        1204 :       lc = lg(c);
     194        3528 :       for (i = 1; i < lc; i++)
     195        2324 :         F2v_set(V, d * uel(c,i)*uel(c,i));
     196             :     }
     197          14 :   return V;
     198             : }
     199             : 
     200             : /* l odd prime. Find D such that j has CM by D, assuming
     201             :  * subst(polmodular(l),x,j) has a double root */
     202             : static long
     203          14 : find_CM(long l, GEN j, GEN T, GEN p)
     204             : {
     205          14 :   const long inv = 0;
     206          14 :   GEN v = list_singular_discs(l);
     207          14 :   long i, n = v[1];
     208          14 :   GEN db = polmodular_db_init(inv);
     209         861 :   for (i = 1; i < n; i++)
     210         861 :     if (F2v_coeff(v,i))
     211             :     {
     212         161 :       GEN C = polclass0(-i, inv, 0, &db);
     213         161 :       GEN F = FqX_eval(C, j, T, p);
     214         161 :       if (signe(F)==0) break;
     215             :     }
     216          14 :   gunclone_deep(db); return i < n ? -i: 0;
     217             : }
     218             : 
     219             : static GEN
     220          14 : vecpoints_to_vecx(GEN x, GEN q1)
     221             : {
     222          42 :   pari_APPLY_type(t_COL, gadd(q1, signe(gmael(x,i,2)) > 0 ? gmael(x,i,1)
     223             :                                                           : negi(gmael(x,i,1))));
     224             : }
     225             : 
     226             : static GEN
     227          14 : Fq_ellcard_CM(long disc, GEN a4, GEN a6, GEN T, GEN p)
     228             : {
     229             :   const struct bb_group *grp;
     230             :   void *E;
     231          14 :   long d = T ? degpol(T): 1;
     232          14 :   GEN q = powiu(p, d), q1 = addiu(q, 1), Q, S;
     233          14 :   Q = qfbsolve(Qfb0(gen_1,gen_0,stoi(-disc)), mkmat22(gen_2, gen_2, p, utoi(d)), 3);
     234          14 :   if (lg(Q)==1) return q1;
     235          14 :   S = vecpoints_to_vecx(Q, q1);
     236          14 :   grp = get_FqE_group(&E, a4, a6, T, p);
     237          14 :   return gen_select_order(S, E, grp);
     238             : }
     239             : 
     240             : /***********************************************************************/
     241             : /**                                                                   **/
     242             : /**                      n-division polynomial                        **/
     243             : /**                                                                   **/
     244             : /***********************************************************************/
     245             : 
     246             : static GEN divpol(GEN t, GEN r2, long n, void *E, const struct bb_algebra *ff);
     247             : 
     248             : /* f_n^2, return ff->(zero|one) or a clone */
     249             : static GEN
     250      145208 : divpol_f2(GEN t, GEN r2, long n, void *E, const struct bb_algebra *ff)
     251             : {
     252      145208 :   if (n==0) return ff->zero(E);
     253      145208 :   if (n<=2) return ff->one(E);
     254      120428 :   if (gmael(t,2,n)) return gmael(t,2,n);
     255       44408 :   gmael(t,2,n) = gclone(ff->sqr(E,divpol(t,r2,n,E,ff)));
     256       44408 :   return gmael(t,2,n);
     257             : }
     258             : 
     259             : /* f_n f_{n-2}, return ff->zero or a clone */
     260             : static GEN
     261       88214 : divpol_ff(GEN t, GEN r2, long n, void *E, const struct bb_algebra *ff)
     262             : {
     263       88214 :   if (n<=2) return ff->zero(E);
     264       88214 :   if (gmael(t,3,n)) return gmael(t,3,n);
     265       56784 :   if (n<=4) return divpol(t,r2,n,E,ff);
     266       25011 :   gmael(t,3,n) = gclone(ff->mul(E,divpol(t,r2,n,E,ff), divpol(t,r2,n-2,E,ff)));
     267       25011 :   return gmael(t,3,n);
     268             : }
     269             : 
     270             : /* f_n, return ff->zero or a clone */
     271             : static GEN
     272      188048 : divpol(GEN t, GEN r2, long n, void *E, const struct bb_algebra *ff)
     273             : {
     274      188048 :   long m = n/2;
     275      188048 :   pari_sp av = avma;
     276             :   GEN f;
     277      188048 :   if (n==0) return ff->zero(E);
     278      184240 :   if (gmael(t,1,n)) return gmael(t,1,n);
     279       51331 :   switch(n)
     280             :   {
     281        7224 :   case 1:
     282             :   case 2:
     283        7224 :     f = ff->one(E);
     284        7224 :     break;
     285       44107 :   default:
     286       44107 :     if (odd(n))
     287       25627 :       if (odd(m))
     288       11382 :         f = ff->sub(E, ff->mul(E, divpol_ff(t,r2,m+2,E,ff),
     289             :                                   divpol_f2(t,r2,m,E,ff)),
     290       11382 :                        ff->mul(E, r2,
     291       11382 :                                   ff->mul(E,divpol_ff(t,r2,m+1,E,ff),
     292             :                                             divpol_f2(t,r2,m+1,E,ff))));
     293             :       else
     294       28490 :         f = ff->sub(E, ff->mul(E, r2,
     295       14245 :                                   ff->mul(E, divpol_ff(t,r2,m+2,E,ff),
     296             :                                              divpol_f2(t,r2,m,E,ff))),
     297       14245 :                        ff->mul(E, divpol_ff(t,r2,m+1,E,ff),
     298             :                                   divpol_f2(t,r2,m+1,E,ff)));
     299             :     else
     300       18480 :       f = ff->sub(E, ff->mul(E, divpol_ff(t,r2,m+2,E,ff),
     301             :                                 divpol_f2(t,r2,m-1,E,ff)),
     302       18480 :                      ff->mul(E, divpol_ff(t,r2,m,E,ff),
     303             :                                 divpol_f2(t,r2,m+1,E,ff)));
     304             :   }
     305       51331 :   gmael(t,1,n) = f = gclone( ff->red(E, f) );
     306       51331 :   set_avma(av); return f;
     307             : }
     308             : 
     309             : static void
     310       11312 : divpol_free(GEN t)
     311             : {
     312       11312 :   long i, l = lg(gel(t,1));
     313      209685 :   for (i=1; i<l; i++)
     314             :   {
     315      198373 :     guncloneNULL(gmael(t,1,i));
     316      198373 :     guncloneNULL(gmael(t,2,i));
     317      198373 :     guncloneNULL(gmael(t,3,i));
     318             :   }
     319       11312 : }
     320             : 
     321             : static GEN
     322        1522 : Flxq_elldivpol34(long n, GEN a4, GEN a6, GEN S, GEN T, ulong p)
     323             : {
     324             :   GEN res;
     325        1522 :   long vs = T[1];
     326        1522 :   switch(n)
     327             :   {
     328         761 :   case 3:
     329         761 :     res = mkpoln(5, Fl_to_Flx(3%p,vs), pol0_Flx(vs), Flx_mulu(a4, 6, p),
     330             :                     Flx_mulu(a6, 12, p), Flx_neg(Flxq_sqr(a4, T, p), p));
     331         761 :     break;
     332         761 :   case 4:
     333             :     {
     334         761 :       GEN a42 = Flxq_sqr(a4, T, p);
     335        1522 :       res = mkpoln(7, pol1_Flx(vs), pol0_Flx(vs), Flx_mulu(a4, 5, p),
     336             :           Flx_mulu(a6, 20, p), Flx_mulu(a42,p-5, p),
     337             :           Flx_mulu(Flxq_mul(a4, a6, T, p), p-4, p),
     338         761 :           Flx_sub(Flx_mulu(Flxq_sqr(a6, T, p), p-8%p, p),
     339             :             Flxq_mul(a4, a42, T, p), p));
     340         761 :       res = FlxX_double(res, p);
     341             :     }
     342         761 :     break;
     343           0 :     default:
     344           0 :       pari_err_BUG("Flxq_elldivpol34");
     345             :       return NULL;/*LCOV_EXCL_LINE*/
     346             :   }
     347        1522 :   setvarn(res, get_FlxqX_var(S));
     348        1522 :   return FlxqX_rem(res, S, T, p);
     349             : }
     350             : 
     351             : static GEN
     352       21102 : Fq_elldivpol34(long n, GEN a4, GEN a6, GEN S, GEN T, GEN p)
     353             : {
     354             :   GEN res;
     355       21102 :   switch(n)
     356             :   {
     357       10551 :   case 3:
     358       10551 :     res = mkpoln(5, utoi(3), gen_0, Fq_mulu(a4, 6, T, p),
     359             :         Fq_mulu(a6, 12, T, p), Fq_neg(Fq_sqr(a4, T, p), T, p));
     360       10551 :     break;
     361       10551 :   case 4:
     362             :     {
     363       10551 :       GEN a42 = Fq_sqr(a4, T, p);
     364       10551 :       res = mkpoln(7, gen_1, gen_0, Fq_mulu(a4, 5, T, p),
     365             :           Fq_mulu(a6, 20, T, p), Fq_Fp_mul(a42,stoi(-5), T, p),
     366             :           Fq_Fp_mul(Fq_mul(a4, a6, T, p), stoi(-4), T, p),
     367             :           Fq_sub(Fq_Fp_mul(Fq_sqr(a6, T, p), stoi(-8), T, p),
     368             :             Fq_mul(a4,a42, T, p), T, p));
     369       10551 :       res = FqX_mulu(res, 2, T, p);
     370             :     }
     371       10551 :     break;
     372           0 :     default:
     373           0 :       pari_err_BUG("Fq_elldivpol34");
     374             :       return NULL;/*LCOV_EXCL_LINE*/
     375             :   }
     376       21102 :   if (S)
     377             :   {
     378       21102 :     setvarn(res, get_FpXQX_var(S));
     379       21102 :     res = FqX_rem(res, S, T, p);
     380             :   }
     381       21102 :   return res;
     382             : }
     383             : 
     384             : static GEN
     385       17670 : rhs(GEN a4, GEN a6, long v)
     386             : {
     387       17670 :   GEN RHS = mkpoln(4, gen_1, gen_0, a4, a6);
     388       17670 :   setvarn(RHS, v); return RHS;
     389             : }
     390             : 
     391             : static GEN
     392        1132 : Flxq_rhs(GEN a4, GEN a6, long v, long vs)
     393             : {
     394        1132 :   GEN RHS = mkpoln(4, pol1_Flx(vs),  pol0_Flx(vs), a4, a6);
     395        1132 :   setvarn(RHS, v); return RHS;
     396             : }
     397             : 
     398             : struct divpolmod_red
     399             : {
     400             :   const struct bb_algebra *ff;
     401             :   void *E;
     402             :   GEN t, r2;
     403             : };
     404             : 
     405             : static void
     406       11312 : divpolmod_init(struct divpolmod_red *d, GEN D3, GEN D4, GEN RHS, long n,
     407             :                void *E, const struct bb_algebra *ff)
     408             : {
     409       11312 :   long k = n+2;
     410       11312 :   d->ff = ff; d->E = E;
     411       11312 :   d->t  = mkvec3(const_vec(k, NULL),const_vec(k, NULL),const_vec(k, NULL));
     412       11312 :   if (k>=3) gmael(d->t,1,3) = gclone(D3);
     413       11312 :   if (k>=4) gmael(d->t,1,4) = gclone(D4);
     414       11312 :   d->r2 = ff->sqr(E, RHS);
     415       11312 : }
     416             : 
     417             : static void
     418       10551 : Fq_elldivpolmod_init(struct divpolmod_red *d, GEN a4, GEN a6, long n, GEN h, GEN T, GEN p)
     419             : {
     420             :   void *E;
     421             :   const struct bb_algebra *ff;
     422       10551 :   GEN RHS, D3 = NULL, D4 = NULL;
     423       10551 :   long v = h ? get_FpXQX_var(h): 0;
     424       10551 :   D3 = n>=0 ? Fq_elldivpol34(3, a4, a6, h, T, p): NULL;
     425       10551 :   D4 = n>=1 ? Fq_elldivpol34(4, a4, a6, h, T, p): NULL;
     426       10551 :   RHS = rhs(a4, a6, v);
     427       10551 :   RHS = h ? FqX_rem(RHS, h, T, p): RHS;
     428       10551 :   RHS = FqX_mulu(RHS, 4, T, p);
     429       10551 :   ff = h ? T ? get_FpXQXQ_algebra(&E, h, T, p): get_FpXQ_algebra(&E, h, p):
     430           0 :            T ? get_FpXQX_algebra(&E, T, p, v): get_FpX_algebra(&E, p, v);
     431       10551 :   divpolmod_init(d, D3, D4, RHS, n, E, ff);
     432       10551 : }
     433             : 
     434             : static void
     435         761 : Flxq_elldivpolmod_init(struct divpolmod_red *d, GEN a4, GEN a6, long n, GEN h, GEN T, ulong p)
     436             : {
     437             :   void *E;
     438             :   const struct bb_algebra *ff;
     439         761 :   GEN RHS, D3 = NULL, D4 = NULL;
     440         761 :   long v = get_FlxqX_var(h), vT = get_Flx_var(T);
     441         761 :   D3 = n>=0 ? Flxq_elldivpol34(3, a4, a6, h, T, p): NULL;
     442         761 :   D4 = n>=1 ? Flxq_elldivpol34(4, a4, a6, h, T, p): NULL;
     443         761 :   RHS = FlxX_Fl_mul(FlxqX_rem(Flxq_rhs(a4, a6, v, vT), h, T, p), 4, p);
     444         761 :   ff = get_FlxqXQ_algebra(&E, h, T, p);
     445         761 :   divpolmod_init(d, D3, D4, RHS, n, E, ff);
     446         761 : }
     447             : 
     448             : /*Computes the n-division polynomial modulo the polynomial h \in Fq[x] */
     449             : GEN
     450         390 : Flxq_elldivpolmod(GEN a4, GEN a6, long n, GEN h, GEN T, ulong p)
     451             : {
     452             :   struct divpolmod_red d;
     453         390 :   pari_sp ltop = avma;
     454             :   GEN res;
     455         390 :   Flxq_elldivpolmod_init(&d, a4, a6, n, h, T, p);
     456         390 :   res = gcopy(divpol(d.t,d.r2,n,d.E,d.ff));
     457         390 :   divpol_free(d.t);
     458         390 :   return gerepileupto(ltop, res);
     459             : }
     460             : 
     461             : /*Computes the n-division polynomial modulo the polynomial h \in Fq[x] */
     462             : GEN
     463        4851 : Fq_elldivpolmod(GEN a4, GEN a6, long n, GEN h, GEN T, GEN p)
     464             : {
     465             :   struct divpolmod_red d;
     466        4851 :   pari_sp ltop = avma;
     467             :   GEN res;
     468        4851 :   if (lgefint(p)==3 && T)
     469             :   {
     470         390 :     ulong pp = p[2];
     471         390 :     GEN a4p = ZX_to_Flx(a4,pp), a6p = ZX_to_Flx(a6,pp);
     472         390 :     GEN hp = ZXX_to_FlxX(h, pp, get_FpX_var(T)), Tp = ZXT_to_FlxT(T , pp);
     473         390 :     res = Flxq_elldivpolmod(a4p, a6p, n, hp, Tp, pp);
     474         390 :     return gerepileupto(ltop, FlxX_to_ZXX(res));
     475             :   }
     476        4461 :   Fq_elldivpolmod_init(&d, a4, a6, n, h, T, p);
     477        4461 :   res = gcopy(divpol(d.t,d.r2,n,d.E,d.ff));
     478        4461 :   divpol_free(d.t);
     479        4461 :   return gerepileupto(ltop, res);
     480             : }
     481             : 
     482             : GEN
     483           0 : FpXQ_elldivpol(GEN a4, GEN a6, long n, GEN T, GEN p)
     484           0 : { return Fq_elldivpolmod(a4,a6,n,NULL,T,p); }
     485             : 
     486             : GEN
     487           0 : Fp_elldivpol(GEN a4, GEN a6, long n, GEN p)
     488           0 : { return Fq_elldivpolmod(a4,a6,n,NULL,NULL,p); }
     489             : 
     490             : static GEN
     491       24451 : Fq_ellyn(struct divpolmod_red *d, long k)
     492             : {
     493       24451 :   pari_sp av = avma;
     494       24451 :   void *E = d->E;
     495       24451 :   const struct bb_algebra *ff = d->ff;
     496       24451 :   if (k==1) return mkvec2(ff->one(E), ff->one(E));
     497             :   else
     498             :   {
     499       18998 :     GEN t = d->t, r2 = d->r2;
     500       18998 :     GEN pn2 = divpol(t,r2,k-2,E,ff);
     501       18998 :     GEN pp2 = divpol(t,r2,k+2,E,ff);
     502       18998 :     GEN pn12 = divpol_f2(t,r2,k-1,E,ff);
     503       18998 :     GEN pp12 = divpol_f2(t,r2,k+1,E,ff);
     504       18998 :     GEN on = ff->red(E,ff->sub(E, ff->mul(E,pp2,pn12), ff->mul(E,pn2,pp12)));
     505       18998 :     GEN f  = divpol(t,r2,k,E,ff);
     506       18998 :     GEN f2 = divpol_f2(t,r2,k,E,ff);
     507       18998 :     GEN f3 = ff->mul(E,f,f2);
     508       18998 :     if (!odd(k)) f3 = ff->mul(E,f3,r2);
     509       18998 :     return gerepilecopy(av,mkvec2(on, f3));
     510             :   }
     511             : }
     512             : 
     513             : static void
     514        6461 : Fq_elldivpolmod_close(struct divpolmod_red *d)
     515        6461 : { divpol_free(d->t); }
     516             : static GEN
     517        1540 : Fq_elldivpol2(GEN a4, GEN a6, GEN T, GEN p)
     518        1540 : { return mkpoln(4, utoi(4), gen_0, Fq_mulu(a4, 4, T, p), Fq_mulu(a6, 4, T, p)); }
     519             : 
     520             : static GEN
     521        1540 : Fq_elldivpol2d(GEN a4, GEN T, GEN p)
     522        1540 : { return mkpoln(3, utoi(6), gen_0, Fq_mulu(a4, 2, T, p)); }
     523             : 
     524             : static GEN
     525        1540 : FqX_numer_isog_abscissa(GEN h, GEN a4, GEN a6, GEN T, GEN p, long vx)
     526             : {
     527             :   GEN mp1, dh, ddh, t, u, t1, t2, t3, t4, f0;
     528        1540 :   long m = degpol(h);
     529        1540 :   mp1 = gel(h, m + 1); /* negative of first power sum */
     530        1540 :   dh = FqX_deriv(h, T, p);
     531        1540 :   ddh = FqX_deriv(dh, T, p);
     532        1540 :   t  = Fq_elldivpol2(a4, a6, T, p);
     533        1540 :   u  = Fq_elldivpol2d(a4, T, p);
     534        1540 :   t1 = FqX_sub(FqX_sqr(dh, T, p), FqX_mul(ddh, h, T, p), T, p);
     535        1540 :   t2 = FqX_mul(u, FqX_mul(h, dh, T, p), T, p);
     536        1540 :   t3 = FqX_mul(FqX_sqr(h, T, p),
     537             :                deg1pol_shallow(stoi(2*m), Fq_mulu(mp1, 2, T, p), vx), T, p);
     538        1540 :   f0 = FqX_add(FqX_sub(FqX_mul(t, t1, T, p), t2, T, p), t3, T, p);
     539        1540 :   t4 = FqX_mul(pol_x(vx),  FqX_sqr(h, T, p), T, p);
     540        1540 :   return FqX_add(t4, f0, T, p);
     541             : }
     542             : 
     543             : static GEN
     544        1036 : Zq_inv(GEN b, GEN T, GEN p, long e)
     545             : {
     546        2023 :   return e==1 ? Fq_inv(b, T, p):
     547         987 :          typ(b)==t_INT ? Zp_inv(b, p, e):  ZpXQ_inv(b, T, p, e);
     548             : }
     549             : 
     550             : static GEN
     551       98441 : Zq_div(GEN a, GEN b, GEN T, GEN q, GEN p, long e)
     552             : {
     553       98441 :   if (e==1) return Fq_div(a, b, T, p);
     554         987 :   return Fq_mul(a, Zq_inv(b, T, p, e), T, q);
     555             : }
     556             : 
     557             : static GEN
     558           0 : Zq_sqrt(GEN b, GEN T, GEN p, long e)
     559             : {
     560           0 :   return e==1 ? Fq_sqrt(b, T, p):
     561           0 :          typ(b)==t_INT ? Zp_sqrt(b, p, e):  ZpXQ_sqrt(b, T, p, e);
     562             : }
     563             : 
     564             : static GEN
     565          14 : Zq_divexact(GEN a, GEN b)
     566          14 : { return typ(a)==t_INT ? diviiexact(a, b): ZX_Z_divexact(a, b); }
     567             : 
     568             : static long
     569          14 : Zq_pval(GEN a, GEN p)
     570          14 : { return typ(a)==t_INT ? Z_pval(a, p): ZX_pval(a, p); }
     571             : 
     572             : static GEN
     573      120204 : Zq_divu_safe(GEN a, ulong b, GEN T, GEN q, GEN p, long e)
     574             : {
     575             :   long v, w;
     576      120204 :   if (e==1) return Fq_div(a, utoi(b), T, q);
     577        2611 :   v = u_pvalrem(b, p, &b);
     578        2611 :   if (v > 0)
     579             :   {
     580          14 :     if (signe(a)==0) return gen_0;
     581          14 :     w = Zq_pval(a, p);
     582          14 :     if (v > w) return NULL;
     583          14 :     a = Zq_divexact(a, powiu(p,v));
     584             :   }
     585        2611 :   return Fq_Fp_mul(a, Fp_inv(utoi(b), q), T, q);
     586             : }
     587             : 
     588             : static GEN
     589      164381 : FqX_shift(GEN P,long n)
     590      164381 : { return RgX_shift_shallow(P, n); }
     591             : 
     592             : static GEN
     593       38822 : FqX_mulhigh_i(GEN f, GEN g, long n, GEN T, GEN p)
     594       38822 : { return FqX_shift(FqX_mul(f,g,T, p),-n); }
     595             : 
     596             : static GEN
     597       38822 : FqX_mulhigh(GEN f, GEN g, long n2, long n, GEN T, GEN p)
     598             : {
     599       38822 :   GEN F = RgX_blocks(f, n2, 2), fl = gel(F,1), fh = gel(F,2);
     600       38822 :   return FqX_add(FqX_mulhigh_i(fl, g, n2, T, p), FqXn_mul(fh, g, n - n2, T, p), T, p);
     601             : }
     602             : 
     603             : static GEN
     604       19411 : FqX_invlift1(GEN Q, GEN P, long t1, long t2, GEN T, GEN p)
     605             : {
     606       19411 :   GEN H = FqXn_mul(FqX_mulhigh(Q, P, t1, t2, T, p), Q, t2-t1, T, p);
     607       19411 :   return FqX_sub(Q, FqX_shift(H, t1), T, p);
     608             : }
     609             : 
     610             : static GEN
     611       19411 : FqX_invsqrtlift1(GEN Q, GEN P, long t1, long t2, GEN T, GEN p)
     612             : {
     613       19411 :   GEN D = FqX_mulhigh(P, FqX_sqr(Q, T, p), t1, t2, T, p);
     614       19411 :   GEN H = FqXn_mul(Q, FqX_halve(D, T, p), t2-t1, T, p);
     615       19411 :   return FqX_sub(Q, FqX_shift(H, t1), T, p);
     616             : }
     617             : 
     618             : /* Q(x^2) = intformal(subst(x^N*P,x,x^2)) */
     619             : static GEN
     620       26537 : ZqX_integ2Xn(GEN P, long N, GEN T, GEN p, GEN pp, long e)
     621             : {
     622       26537 :   long d = degpol(P), v = varn(P);
     623             :   long k;
     624             :   GEN Q;
     625       26537 :   if(d==-1) return pol_0(v);
     626       19411 :   Q = cgetg(d+3,t_POL);
     627       19411 :   Q[1] = evalsigne(1) | evalvarn(v);
     628       83076 :   for (k = 0; k <= d; k++)
     629             :   {
     630       63665 :     GEN q = Zq_divu_safe(gel(P,2+k), 2*(k+N)+1, T, p, pp, e);
     631       63665 :     if (!q) return NULL;
     632       63665 :     gel(Q, 2+k) = q;
     633             :   }
     634       19411 :   return ZXX_renormalize(Q,d+3);
     635             : }
     636             : 
     637             : /* solution of G*(S'^2)=(S/x)*(HoS) mod x^m */
     638             : static GEN
     639        7126 : Zq_Weierstrass(GEN a4, GEN a6, GEN b4, GEN b6, long m, GEN T, GEN p, GEN pp, long n)
     640             : {
     641        7126 :   pari_sp av = avma;
     642        7126 :   long v = 0;
     643        7126 :   ulong mask = quadratic_prec_mask(m);
     644        7126 :   GEN iGdS2 = pol_1(v);
     645        7126 :   GEN G = mkpoln(4, a6, a4, gen_0, gen_1);
     646        7126 :   GEN GdS2 = G, S = pol_x(v), sG = pol_1(v), isG = sG, dS = sG;
     647        7126 :   long N = 1;
     648       26537 :   for (;mask>1;)
     649             :   {
     650             :     GEN S2, HS, K, dK, E;
     651       26537 :     long N2 = N, d;
     652       26537 :     N<<=1; if (mask & 1) N--;
     653       26537 :     mask >>= 1;
     654       26537 :     d = N-N2;
     655       26537 :     S2 = FqX_sqr(S, T, p);
     656       26537 :     HS = FqX_Fq_add(FqX_Fq_mul(S, b6, T, p), b4, T, p);
     657       26537 :     HS = FqX_Fq_add(FqXn_mul(S2, HS, N, T, p), gen_1, T, p);
     658       26537 :     HS = FqXn_mul(HS, FqX_shift(S,-1), N, T, p);
     659       26537 :     sG  = FqXn_mul(G, isG, N2, T, p);
     660             :     /* (HS-Gds2)/(Gds2*sG) */
     661       26537 :     dK = FqXn_mul(FqX_shift(FqX_sub(HS, GdS2, T, p), -N2),
     662             :                   FqXn_mul(iGdS2, isG, d, T, p), d, T, p);
     663       26537 :     K = ZqX_integ2Xn(dK, N2, T, p, pp, n);
     664       26537 :     if (!K) return gc_NULL(av);
     665       26537 :     E = FqXn_mul(FqXn_mul(K, sG, d, T, p), dS, d, T, p);
     666       26537 :     S = FqX_add(S, FqX_shift(E, N2+1), T, p);
     667       26537 :     if (mask <= 1) break;
     668       19411 :     isG = FqX_invsqrtlift1(isG, G, N2, N, T, p);
     669       19411 :     dS = FqX_deriv(S, T, p);
     670       19411 :     GdS2 = FqX_mul(G, FqX_sqr(dS, T, p), T, p);
     671       19411 :     iGdS2 = FqX_invlift1(iGdS2, GdS2, N2, N, T, p);
     672             :   }
     673        7126 :   return gerepileupto(av, S);
     674             : }
     675             : 
     676             : static GEN
     677        7126 : ZqXn_WNewton(GEN S, long l, GEN a4, GEN a6, GEN pp1, GEN T, GEN p, GEN pp, long e)
     678             : {
     679        7126 :   long d = degpol(S);
     680             :   long k;
     681        7126 :   GEN Ge = cgetg(2+d,t_POL);
     682        7126 :   Ge[1] = evalsigne(1);
     683        7126 :   gel(Ge,2) = pp1;
     684        7126 :   if (d >= 2)
     685             :   {
     686        7126 :     GEN g = Zq_divu_safe(Fq_sub(gel(S,4), Fq_mulu(a4,(l-1),T,p),T,p), 6,T,p,pp,e);
     687        7126 :     if (!g) return NULL;
     688        7126 :     gel(Ge, 3) = g;
     689             :   }
     690        7126 :   if (d >= 3)
     691             :   {
     692        7126 :     GEN g = Zq_divu_safe(Fq_sub(Fq_sub(gel(S,5),
     693             :             Fq_mul(a4,Fq_mulu(pp1,6,T,p),T,p),T,p),
     694        7126 :             Fq_mulu(a6,(l-1)*2,T,p),T,p),10,T,p,pp,e);
     695        7126 :     if (!g) return NULL;
     696        7126 :     gel(Ge, 4) = g;
     697             :   }
     698       49413 :   for (k = 4; k <= d; k++)
     699             :   {
     700       84574 :     GEN g = Zq_divu_safe(Fq_sub(Fq_sub(gel(S,4+k-2),
     701       42287 :             Fq_mul(a4,Fq_mulu(gel(Ge,k-1),4*k-6,T,p),T,p),T,p),
     702       42287 :             Fq_mul(a6,Fq_mulu(gel(Ge,k-2),4*k-8,T,p),T,p),T,p),
     703       42287 :             4*k-2, T, p, pp, e);
     704       42287 :     if (!g) return NULL;
     705       42287 :     gel(Ge, k+1) = g;
     706             :   }
     707        7126 :   return ZXX_renormalize(Ge, 2+d);
     708             : }
     709             : 
     710             : /****************************************************************************/
     711             : /*               SIMPLE ELLIPTIC CURVE OVER Fq                              */
     712             : /****************************************************************************/
     713             : 
     714             : static GEN
     715        2604 : Fq_ellj(GEN a4, GEN a6, GEN T, GEN p)
     716             : {
     717        2604 :   pari_sp ltop=avma;
     718        2604 :   GEN a43 = Fq_mulu(Fq_powu(a4, 3, T, p), 4, T, p);
     719        2604 :   GEN j   = Fq_div(Fq_mulu(a43, 1728, T, p),
     720             :                    Fq_add(a43, Fq_mulu(Fq_sqr(a6, T, p), 27, T, p), T, p), T, p);
     721        2604 :   return gerepileupto(ltop, j);
     722             : }
     723             : 
     724             : static GEN
     725        2688 : Zq_ellj(GEN a4, GEN a6, GEN T, GEN p, GEN pp, long e)
     726             : {
     727        2688 :   pari_sp ltop=avma;
     728        2688 :   GEN a43 = Fq_mulu(Fq_powu(a4, 3, T, p), 4, T, p);
     729        2688 :   GEN j   = Zq_div(Fq_mulu(a43, 1728, T, p),
     730             :                    Fq_add(a43, Fq_mulu(Fq_sqr(a6, T, p), 27, T, p), T, p), T, p, pp, e);
     731        2688 :   return gerepileupto(ltop, j);
     732             : }
     733             : /****************************************************************************/
     734             : /*                              EIGENVALUE                                  */
     735             : /****************************************************************************/
     736             : 
     737             : static GEN
     738         371 : Flxq_find_eigen_Frobenius(GEN a4, GEN a6, GEN h, GEN T, ulong p)
     739             : {
     740         371 :   long v = get_FlxqX_var(h), vT = get_Flx_var(T);
     741         371 :   GEN RHS = FlxqX_rem(Flxq_rhs(a4, a6, v, vT), h, T, p);
     742         371 :   return FlxqXQ_halfFrobenius(RHS, h, T, p);
     743             : }
     744             : 
     745             : static GEN
     746        6090 : Fq_find_eigen_Frobenius(GEN a4, GEN a6, GEN h, GEN T, GEN p)
     747             : {
     748        6090 :   long v = T ? get_FpXQX_var(h): get_FpX_var(h);
     749        6090 :   GEN RHS  = FqX_rem(rhs(a4, a6, v), h, T, p);
     750       11942 :   return T ? FpXQXQ_halfFrobenius(RHS, h, T, p):
     751        5852 :              FpXQ_pow(RHS, shifti(p, -1), h, p);
     752             : }
     753             : /*Finds the eigenvalue of the Frobenius given E, ell odd prime, h factor of the
     754             :  *ell-division polynomial, p and tr the possible values for the trace
     755             :  *(useful for primes with one root)*/
     756             : static ulong
     757         504 : find_eigen_value_oneroot(GEN a4, GEN a6, ulong ell, GEN tr, GEN h, GEN T, GEN p)
     758             : {
     759         504 :   pari_sp ltop = avma;
     760             :   ulong t;
     761             :   struct divpolmod_red d;
     762             :   GEN f, Dy, Gy;
     763         504 :   h = FqX_get_red(h, T, p);
     764         504 :   Gy = Fq_find_eigen_Frobenius(a4, a6, h, T, p);
     765         504 :   t = Fl_div(tr[1], 2, ell);
     766         504 :   if (t < (ell>>1)) t = ell - t;
     767         504 :   Fq_elldivpolmod_init(&d, a4, a6, t, h, T, p);
     768         504 :   f = Fq_ellyn(&d, t);
     769         504 :   Dy = FqXQ_mul(Gy, gel(f,2), h, T, p);
     770         504 :   if (!gequal(gel(f,1), Dy)) t = ell-t;
     771         504 :   Fq_elldivpolmod_close(&d);
     772         504 :   return gc_ulong(ltop, t);
     773             : }
     774             : 
     775             : static ulong
     776         371 : Flxq_find_eigen_value_power(GEN a4, GEN a6, ulong ell, long k, ulong lambda,
     777             :                             GEN h, GEN T, ulong p)
     778             : {
     779         371 :   pari_sp ltop = avma;
     780         371 :   ulong t, ellk1 = upowuu(ell, k-1), ellk = ell*ellk1;
     781             :   pari_timer ti;
     782             :   struct divpolmod_red d;
     783             :   GEN Gy;
     784         371 :   timer_start(&ti);
     785         371 :   h = FlxqX_get_red(h, T, p);
     786         371 :   Gy = Flxq_find_eigen_Frobenius(a4, a6, h, T, p);
     787         371 :   if (DEBUGLEVEL>2) err_printf(" (%ld ms)",timer_delay(&ti));
     788         371 :   Flxq_elldivpolmod_init(&d, a4, a6, ellk, h, T, p);
     789        1685 :   for (t = lambda; t < ellk; t += ellk1)
     790             :   {
     791        1685 :     GEN f = Fq_ellyn(&d, t);
     792        1685 :     GEN Dr = FlxqXQ_mul(Gy, gel(f,2), h, T, p);
     793        1685 :     if (varn(gel(f,1))!=varn(Dr)) pari_err_BUG("find_eigen_value_power");
     794        1685 :     if (gequal(gel(f,1), Dr)) break;
     795        1441 :     if (gequal(gel(f,1), FlxX_neg(Dr,p))) { t = ellk-t; break; }
     796             :   }
     797         371 :   if (DEBUGLEVEL>2) err_printf(" (%ld ms)",timer_delay(&ti));
     798         371 :   Fq_elldivpolmod_close(&d);
     799         371 :   return gc_ulong(ltop, t);
     800             : }
     801             : 
     802             : /*Finds the eigenvalue of the Frobenius modulo ell^k given E, ell, k, h factor
     803             :  *of the ell-division polynomial, lambda the previous eigen value and p */
     804             : static ulong
     805        5586 : Fq_find_eigen_value_power(GEN a4, GEN a6, ulong ell, long k, ulong lambda, GEN h, GEN T, GEN p)
     806             : {
     807        5586 :   pari_sp ltop = avma;
     808        5586 :   ulong t, ellk1 = upowuu(ell, k-1), ellk = ell*ellk1;
     809             :   pari_timer ti;
     810             :   struct divpolmod_red d;
     811             :   GEN Gy;
     812        5586 :   timer_start(&ti);
     813        5586 :   h = FqX_get_red(h, T, p);
     814        5586 :   Gy = Fq_find_eigen_Frobenius(a4, a6, h, T, p);
     815        5586 :   if (DEBUGLEVEL>2) err_printf(" (%ld ms)",timer_delay(&ti));
     816        5586 :   Fq_elldivpolmod_init(&d, a4, a6, ellk, h, T, p);
     817       22262 :   for (t = lambda; t < ellk; t += ellk1)
     818             :   {
     819       22262 :     GEN f = Fq_ellyn(&d, t);
     820       22262 :     GEN Dr = FqXQ_mul(Gy, gel(f,2), h, T, p);
     821       22262 :     if (varn(gel(f,1))!=varn(Dr)) pari_err_BUG("find_eigen_value_power");
     822       22262 :     if (gequal(gel(f,1), Dr)) break;
     823       17830 :     if (gequal(gel(f,1), FqX_neg(Dr,T,p))) { t = ellk-t; break; }
     824             :   }
     825        5586 :   if (DEBUGLEVEL>2) err_printf(" (%ld ms)",timer_delay(&ti));
     826        5586 :   Fq_elldivpolmod_close(&d);
     827        5586 :   return gc_ulong(ltop, t);
     828             : }
     829             : 
     830             : static ulong
     831        5957 : find_eigen_value_power(GEN a4, GEN a6, ulong ell, long k, ulong lambda, GEN hq, GEN T, GEN p)
     832             : {
     833        5957 :   ulong pp = itou_or_0(p);
     834        5957 :   if (pp && T)
     835             :   {
     836         371 :     GEN a4p = ZX_to_Flx(a4, pp);
     837         371 :     GEN a6p = ZX_to_Flx(a6, pp);
     838         371 :     GEN hp = ZXXT_to_FlxXT(hq, pp,varn(a4));
     839         371 :     GEN Tp = ZXT_to_FlxT(T, pp);
     840         371 :     return Flxq_find_eigen_value_power(a4p, a6p, ell, k, lambda, hp, Tp, pp);
     841             :   }
     842        5586 :   return Fq_find_eigen_value_power(a4, a6, ell, k, lambda, hq, T, p);
     843             : }
     844             : 
     845             : static GEN
     846        8939 : find_kernel(GEN a4, GEN a6, long l, GEN b4, GEN b6, GEN pp1, GEN T, GEN p, GEN pp, long e)
     847             : {
     848             :   GEN Ge, S, Sd;
     849        8939 :   long d = ((l+1)>>1)+1;
     850        8939 :   if(l==3)
     851        1813 :     return deg1pol(gen_1, Fq_neg(pp1, T, p), 0);
     852        7126 :   S = Zq_Weierstrass(a4, a6, b4, b6, d + 1, T, p, pp, e);
     853        7126 :   if (S==NULL) return NULL;
     854        7126 :   S  = FqX_shift(S, -1);
     855        7126 :   Sd = FqXn_inv(S, d, T, p);
     856        7126 :   Ge = ZqXn_WNewton(Sd, l, a4, a6, pp1, T, p, pp, e);
     857        7126 :   if (!Ge) return NULL;
     858        7126 :   Ge = FqX_neg(Ge, T, p);
     859         714 :   Ge = T && lgefint(pp)==3 ? ZlXQXn_expint(Ge, d, T, p, pp[2])
     860        7304 :                            : FqXn_expint(Ge, d, T, p);
     861        7126 :   Ge = RgX_recip(FqX_red(Ge, T, pp));
     862        7126 :   if (degpol(Ge)==(l-1)>>1) return Ge;
     863        1463 :   return NULL;
     864             : }
     865             : 
     866             : static GEN
     867        6573 : compute_u(GEN gprime, GEN Dxxg, GEN DxJg, GEN DJJg, GEN j, GEN pJ, GEN px, ulong q, GEN E4, GEN E6, GEN T, GEN p, GEN pp, long e)
     868             : {
     869        6573 :   pari_sp ltop = avma;
     870        6573 :   GEN dxxgj = FqX_eval(Dxxg, j, T, p);
     871        6573 :   GEN dxJgj = FqX_eval(DxJg, j, T, p);
     872        6573 :   GEN dJJgj = FqX_eval(DJJg, j, T, p);
     873        6573 :   GEN E42 = Fq_sqr(E4, T, p), E6ovE4 = Zq_div(E6, E4, T, p, pp, e);
     874        6573 :   GEN a = Fq_mul(gprime, dxxgj, T, p);
     875        6573 :   GEN b = Fq_mul(Fq_mul(Fq_mulu(j,2*q, T, p), dxJgj, T, p), E6ovE4, T, p);
     876        6573 :   GEN c = Fq_mul(Zq_div(Fq_sqr(E6ovE4, T, p), gprime, T, p, pp, e), j, T, p);
     877        6573 :   GEN d = Fq_mul(Fq_mul(c,sqru(q), T, p), Fq_add(pJ, Fq_mul(j, dJJgj, T, p), T, p), T, p);
     878        6573 :   GEN f = Fq_sub(Fq_div(E6ovE4,utoi(3), T, p),
     879             :                  Zq_div(E42, Fq_mulu(E6,2,T, p), T, p, pp, e), T, p);
     880        6573 :   GEN g = Fq_sub(Fq_sub(b,a,T,p), d, T, p);
     881        6573 :   return gerepileupto(ltop, Fq_add(Zq_div(g,px,T,p,pp,e), Fq_mulu(f,q,T,p), T, p));
     882             : }
     883             : 
     884             : static void
     885        8890 : a4a6t(GEN *a4t, GEN *a6t, ulong l, GEN E4t, GEN E6t, GEN T, GEN p)
     886             : {
     887        8890 :   GEN l2 = modii(sqru(l), p), l4 = Fp_sqr(l2, p), l6 = Fp_mul(l4, l2, p);
     888        8890 :   *a4t = Fq_mul(E4t, Fp_muls(l4, -3, p), T, p);
     889        8890 :   *a6t = Fq_mul(E6t, Fp_muls(l6, -2, p), T, p);
     890        8890 : }
     891             : static void
     892          49 : a4a6t_from_J(GEN *a4t, GEN *a6t, ulong l, GEN C4t, GEN C6t, GEN T, GEN p)
     893             : {
     894          49 :   GEN l2 = modii(sqru(l), p), l4 = Fp_sqr(l2, p), l6 = Fp_mul(l4, l2, p);
     895          49 :   GEN v = Fp_inv(stoi(-864), p), u = Fp_mulu(v, 18, p);
     896          49 :   *a4t = Fq_mul(C4t, Fp_mul(u, l4, p), T, p);
     897          49 :   *a6t = Fq_mul(C6t, Fp_mul(v, l6, p), T, p);
     898          49 : }
     899             : /* Finds the isogenous EC, and the sum of the x-coordinates of the points in
     900             :  * the kernel of the isogeny E -> Eb
     901             :  * E: elliptic curve, ell: a prime, meqn: Atkin modular equation
     902             :  * g: root of meqn defining isogenous curve Eb. */
     903             : static GEN
     904        2576 : find_isogenous_from_Atkin(GEN a4, GEN a6, ulong ell, struct meqn *MEQN, GEN g, GEN T, GEN pp, long e)
     905             : {
     906        2576 :   pari_sp ltop = avma, btop;
     907        2576 :   GEN meqn = MEQN->eq, meqnx, Dmeqnx, Roots, gprime, u1;
     908        2576 :   long k, vJ = MEQN->vy;
     909        2576 :   GEN p = e==1 ? pp: powiu(pp, e);
     910        2576 :   GEN j = Zq_ellj(a4, a6, T, p, pp, e);
     911        2576 :   GEN E4 = Fq_div(a4, stoi(-3), T, p);
     912        2576 :   GEN E6 = Fq_neg(Fq_halve(a6, T, p), T, p);
     913        2576 :   GEN Dx = RgX_deriv(meqn);
     914        2576 :   GEN DJ = deriv(meqn, vJ);
     915        2576 :   GEN Dxg = FpXY_Fq_evaly(Dx, g, T, p, vJ);
     916        2576 :   GEN px = FqX_eval(Dxg, j, T, p), dx = Fq_mul(px, g, T, p);
     917        2576 :   GEN DJg = FpXY_Fq_evaly(DJ, g, T, p, vJ);
     918        2576 :   GEN pJ = FqX_eval(DJg, j, T, p), dJ = Fq_mul(pJ, j, T, p);
     919        2576 :   GEN Dxx = RgX_deriv(Dx);
     920        2576 :   GEN DxJg = FqX_deriv(Dxg, T, p);
     921             : 
     922        2576 :   GEN Dxxg = FpXY_Fq_evaly(Dxx, g, T, p, vJ);
     923        2576 :   GEN DJJg = FqX_deriv(DJg, T, p);
     924             :   GEN a, b;
     925        2576 :   if (!signe(Fq_red(dJ,T,pp)) || !signe(Fq_red(dx,T,pp)))
     926             :   {
     927          21 :     if (DEBUGLEVEL>0) err_printf("[A: d%c=0]",signe(dJ)? 'x': 'J');
     928          21 :     return gc_NULL(ltop);
     929             :   }
     930        2555 :   a = Fq_mul(dJ, Fq_mul(g, E6, T, p), T, p);
     931        2555 :   b = Fq_mul(E4, dx, T, p);
     932        2555 :   gprime = Zq_div(a, b, T, p, pp, e);
     933             : 
     934        2555 :   u1 = compute_u(gprime, Dxxg, DxJg, DJJg, j, pJ, px, 1, E4, E6, T, p, pp, e);
     935        2555 :   meqnx = FpXY_Fq_evaly(meqn, g, T, p, vJ);
     936        2555 :   Dmeqnx = FqX_deriv(meqnx, T, pp);
     937        2555 :   Roots = FqX_roots(meqnx, T, pp);
     938             : 
     939        2555 :   btop = avma;
     940        4032 :   for (k = lg(Roots)-1; k >= 1; k--, set_avma(btop))
     941             :   {
     942        4032 :     GEN jt = gel(Roots, k);
     943        4032 :     if (signe(FqX_eval(Dmeqnx, jt, T, pp))==0)
     944           0 :       continue;
     945        4032 :     if (e > 1)
     946          91 :       jt = ZqX_liftroot(meqnx, gel(Roots, k), T, pp, e);
     947        4032 :     if (signe(Fq_red(jt, T, pp)) == 0 || signe(Fq_sub(jt, utoi(1728), T, pp)) == 0)
     948             :     {
     949          14 :       if (DEBUGLEVEL>0) err_printf("[A: jt=%ld]",signe(Fq_red(jt,T,p))? 1728: 0);
     950          14 :       return gc_NULL(ltop);
     951             :     }
     952             :     else
     953             :     {
     954        4018 :       GEN pxstar = FqX_eval(Dxg, jt, T, p);
     955        4018 :       GEN dxstar = Fq_mul(pxstar, g, T, p);
     956        4018 :       GEN pJstar = FqX_eval(DJg, jt, T, p);
     957        4018 :       GEN dJstar = Fq_mul(Fq_mulu(jt, ell, T, p), pJstar, T, p);
     958        4018 :       GEN u = Fq_mul(Fq_mul(dxstar, dJ, T, p), E6, T, p);
     959        4018 :       GEN v = Fq_mul(Fq_mul(dJstar, dx, T, p), E4, T, p);
     960        4018 :       GEN E4t = Zq_div(Fq_mul(Fq_sqr(u, T, p), jt, T, p), Fq_mul(Fq_sqr(v, T, p), Fq_sub(jt, utoi(1728), T, p), T, p), T, p, pp, e);
     961        4018 :       GEN E6t = Zq_div(Fq_mul(u, E4t, T, p), v, T, p, pp, e);
     962        4018 :       GEN u2 = compute_u(gprime, Dxxg, DxJg, DJJg, jt, pJstar, pxstar, ell, E4t, E6t, T, p, pp, e);
     963        4018 :       GEN pp1 = Fq_mulu(Fq_sub(u1, u2, T, p), 3*ell, T, p);
     964             :       GEN a4t, a6t, h;
     965        4018 :       a4a6t(&a4t, &a6t, ell, E4t, E6t, T, p);
     966        4018 :       h = find_kernel(a4, a6, ell, a4t, a6t, pp1, T, p, pp, e);
     967        4018 :       if (h && signe(Fq_elldivpolmod(a4, a6, ell, h, T, pp))==0)
     968        2541 :         return gerepilecopy(ltop, mkvec3(a4t, a6t, h));
     969             :     }
     970             :   }
     971           0 :   pari_err_BUG("find_isogenous_from_Atkin, kernel not found");
     972             :   return NULL;/*LCOV_EXCL_LINE*/
     973             : }
     974             : 
     975             : /* Finds E' ell-isogenous to E and the trace term p1 from canonical modular
     976             :  *   equation meqn
     977             :  * E: elliptic curve, ell: a prime, meqn: canonical modular equation
     978             :  * g: root of meqn defining isogenous curve Eb. */
     979             : static GEN
     980        4879 : find_isogenous_from_canonical(GEN a4, GEN a6, ulong ell, struct meqn *MEQN, GEN g, GEN T, GEN pp, long e)
     981             : {
     982        4879 :   pari_sp ltop = avma;
     983        4879 :   GEN meqn = MEQN->eq;
     984        4879 :   long vJ = MEQN->vy;
     985        4879 :   GEN p = e==1 ? pp: powiu(pp, e);
     986             :   GEN h;
     987        4879 :   GEN E4 = Fq_div(a4, stoi(-3), T, p);
     988        4879 :   GEN E6 = Fq_neg(Fq_halve(a6, T, p), T, p);
     989        4879 :   GEN E42 = Fq_sqr(E4, T, p);
     990        4879 :   GEN E43 = Fq_mul(E4, E42, T, p);
     991        4879 :   GEN E62 = Fq_sqr(E6, T, p);
     992        4879 :   GEN delta = Fq_div(Fq_sub(E43, E62, T, p), utoi(1728), T, p);
     993        4879 :   GEN j = Zq_div(E43, delta, T, p, pp, e);
     994        4879 :   GEN Dx = RgX_deriv(meqn);
     995        4879 :   GEN DJ = deriv(meqn, vJ);
     996        4879 :   GEN Dxg = FpXY_Fq_evaly(Dx, g, T, p, vJ);
     997        4879 :   GEN px  = FqX_eval(Dxg, j, T, p), dx  = Fq_mul(px, g, T, p);
     998        4879 :   GEN DJg = FpXY_Fq_evaly(DJ, g, T, p, vJ);
     999        4879 :   GEN pJ = FqX_eval(DJg, j, T, p), dJ = Fq_mul(j, pJ, T, p);
    1000        4879 :   GEN Dxx = RgX_deriv(Dx);
    1001        4879 :   GEN DxJg = FqX_deriv(Dxg, T, p);
    1002             : 
    1003        4879 :   GEN ExJ = FqX_eval(DxJg, j, T, p);
    1004        4879 :   ulong tis = ugcd(12, ell-1), is = 12 / tis;
    1005        4879 :   GEN itis = Fq_inv(stoi(-tis), T, p);
    1006        4879 :   GEN deltal = Fq_div(Fq_mul(delta, Fq_powu(g, tis, T, p), T, p), powuu(ell, 12), T, p);
    1007             :   GEN E4l, E6l, a4t, a6t, p_1;
    1008        4879 :   if (signe(Fq_red(dx,T, pp))==0)
    1009             :   {
    1010           0 :     if (DEBUGLEVEL>0) err_printf("[C: dx=0]");
    1011           0 :     return gc_NULL(ltop);
    1012             :   }
    1013        4879 :   if (signe(Fq_red(dJ, T, pp))==0)
    1014             :   {
    1015             :     GEN jl;
    1016           0 :     if (DEBUGLEVEL>0) err_printf("[C: dJ=0]");
    1017           0 :     E4l = Fq_div(E4, sqru(ell), T, p);
    1018           0 :     jl  = Zq_div(Fq_powu(E4l, 3, T, p), deltal, T, p, pp, e);
    1019           0 :     E6l = Zq_sqrt(Fq_mul(Fq_sub(jl, utoi(1728), T, p),
    1020             :                          deltal, T, p), T, pp, e);
    1021           0 :     p_1 = gen_0;
    1022             :   }
    1023             :   else
    1024             :   {
    1025             :     GEN jl, f, fd, Dgs, Djs, jld;
    1026        4879 :     GEN E2s = Zq_div(Fq_mul(Fq_neg(Fq_mulu(E6, 12, T, p), T, p), dJ, T, p),
    1027             :                      Fq_mul(Fq_mulu(E4, is, T, p), dx, T, p), T, p, pp, e);
    1028        4879 :     GEN gd = Fq_mul(Fq_mul(E2s, itis, T, p), g, T, p);
    1029        4879 :     GEN jd = Zq_div(Fq_mul(Fq_neg(E42, T, p), E6, T, p), delta, T, p, pp, e);
    1030        4879 :     GEN E0b = Zq_div(E6, Fq_mul(E4, E2s, T, p), T, p, pp, e);
    1031        4879 :     GEN Dxxgj = FqXY_eval(Dxx, g, j, T, p);
    1032        4879 :     GEN Dgd = Fq_add(Fq_mul(gd, px, T, p), Fq_mul(g, Fq_add(Fq_mul(gd, Dxxgj, T, p), Fq_mul(jd, ExJ, T, p), T, p), T, p), T, p);
    1033        4879 :     GEN DJgJj = FqX_eval(FqX_deriv(DJg, T, p), j, T, p);
    1034        4879 :     GEN Djd = Fq_add(Fq_mul(jd, pJ, T, p), Fq_mul(j, Fq_add(Fq_mul(jd, DJgJj, T, p), Fq_mul(gd, ExJ, T, p), T, p), T, p), T, p);
    1035        4879 :     GEN E0bd = Zq_div(Fq_sub(Fq_mul(Dgd, itis, T, p), Fq_mul(E0b, Djd, T, p), T, p), dJ, T, p, pp, e);
    1036        4879 :     E4l = Fq_div(Fq_sub(E4, Fq_mul(E2s, Fq_sub(Fq_sub(Fq_add(Zq_div(Fq_mulu(E0bd, 12, T, p), E0b, T, p, pp, e), Zq_div(Fq_mulu(E42, 6, T, p), E6, T, p, pp, e), T, p), Zq_div(Fq_mulu(E6, 4, T, p), E4, T, p, pp, e), T, p), E2s, T, p), T, p), T, p), sqru(ell), T, p);
    1037        4879 :     jl = Zq_div(Fq_powu(E4l, 3, T, p), deltal, T, p, pp, e);
    1038        4879 :     if (signe(Fq_red(jl,T,pp))==0)
    1039             :     {
    1040           7 :       if (DEBUGLEVEL>0) err_printf("[C: jl=0]");
    1041           7 :       return gc_NULL(ltop);
    1042             :     }
    1043        4872 :     f =  Zq_div(powuu(ell, is), g, T, p, pp, e);
    1044        4872 :     fd = Fq_neg(Fq_mul(Fq_mul(E2s, f, T, p), itis, T, p), T, p);
    1045        4872 :     Dgs = FqXY_eval(Dx, f, jl, T, p);
    1046        4872 :     Djs = FqXY_eval(DJ, f, jl, T, p);
    1047        4872 :     jld = Zq_div(Fq_mul(Fq_neg(fd, T, p), Dgs, T, p),
    1048             :                  Fq_mulu(Djs, ell, T, p), T, p, pp, e);
    1049        4872 :     E6l = Zq_div(Fq_mul(Fq_neg(E4l, T, p), jld, T, p), jl, T, p, pp, e);
    1050        4872 :     p_1 = Fq_neg(Fq_halve(Fq_mulu(E2s, ell, T, p), T, p),T,p);
    1051             :   }
    1052        4872 :   a4a6t(&a4t, &a6t, ell, E4l, E6l, T, p);
    1053        4872 :   h = find_kernel(a4, a6, ell, a4t, a6t, p_1, T, p, pp, e);
    1054        4872 :   if (!h) return NULL;
    1055        4872 :   return gerepilecopy(ltop, mkvec3(a4t, a6t, h));
    1056             : }
    1057             : 
    1058             : static GEN
    1059          98 : corr(GEN c4, GEN c6, GEN T, GEN p, GEN pp, long e)
    1060             : {
    1061          98 :   GEN c46 = Zq_div(Fq_sqr(c4, T, p), c6, T, p, pp, e);
    1062          98 :   GEN c64 = Zq_div(c6, c4, T, p, pp, e);
    1063          98 :   GEN a = Fp_divu(gen_2, 3, p);
    1064          98 :   return Fq_add(Fq_halve(c46, T, p), Fq_mul(a, c64, T, p), T, p);
    1065             : }
    1066             : 
    1067             : static GEN
    1068         168 : RgXY_deflatex(GEN H, long n, long d)
    1069             : {
    1070         168 :   long i, l = lg(H);
    1071         168 :   GEN R = cgetg(l, t_POL);
    1072         168 :   R[1] = H[1];
    1073         980 :   for(i = 2; i < l; i++)
    1074             :   {
    1075         812 :     GEN Hi = gel(H, i);
    1076         812 :     gel(R,i) = typ(Hi)==t_POL? RgX_deflate(RgX_shift_shallow(Hi, d), n): Hi;
    1077             :   }
    1078         168 :   return RgX_renormalize_lg(R, l);
    1079             : }
    1080             : 
    1081             : static GEN
    1082          70 : Fq_polmodular_eval(GEN meqn, GEN j, long N, GEN T, GEN p, long vJ)
    1083             : {
    1084          70 :   pari_sp av = avma;
    1085             :   GEN R, dR, ddR;
    1086          70 :   long t0 = N%3 == 1 ? 2: 0;
    1087          70 :   long t2 = N%3 == 1 ? 0: 2;
    1088          70 :   if (N == 3)
    1089             :   {
    1090          14 :     GEN P = FpXX_red(meqn, p);
    1091          14 :     GEN dP = deriv(P, -1), ddP = deriv(dP, -1);
    1092          14 :     R = FpXY_Fq_evaly(P, j, T, p, vJ);
    1093          14 :     dR = FpXY_Fq_evaly(dP, j, T, p, vJ);
    1094          14 :     ddR = FpXY_Fq_evaly(ddP, j, T, p, vJ);
    1095          14 :     return gerepilecopy(av, mkvec3(R,dR,ddR));
    1096             :   }
    1097             :   else
    1098             :   {
    1099          56 :     GEN P5 = FpXX_red(meqn, p);
    1100          56 :     GEN H = RgX_splitting(P5, 3);
    1101          56 :     GEN H0 = RgXY_deflatex(gel(H,1), 3, -t0);
    1102          56 :     GEN H1 = RgXY_deflatex(gel(H,2), 3, -1);
    1103          56 :     GEN H2 = RgXY_deflatex(gel(H,3), 3, -t2);
    1104          56 :     GEN h0 = FpXY_Fq_evaly(H0, j, T, p, vJ);
    1105          56 :     GEN h1 = FpXY_Fq_evaly(H1, j, T, p, vJ);
    1106          56 :     GEN h2 = FpXY_Fq_evaly(H2, j, T, p, vJ);
    1107          56 :     GEN dH0 = RgX_deriv(H0);
    1108          56 :     GEN dH1 = RgX_deriv(H1);
    1109          56 :     GEN dH2 = RgX_deriv(H2);
    1110          56 :     GEN ddH0 = RgX_deriv(dH0);
    1111          56 :     GEN ddH1 = RgX_deriv(dH1);
    1112          56 :     GEN ddH2 = RgX_deriv(dH2);
    1113          56 :     GEN d0 = FpXY_Fq_evaly(dH0, j, T, p, vJ);
    1114          56 :     GEN d1 = FpXY_Fq_evaly(dH1, j, T, p, vJ);
    1115          56 :     GEN d2 = FpXY_Fq_evaly(dH2, j, T, p, vJ);
    1116          56 :     GEN dd0 = FpXY_Fq_evaly(ddH0, j, T, p, vJ);
    1117          56 :     GEN dd1 = FpXY_Fq_evaly(ddH1, j, T, p, vJ);
    1118          56 :     GEN dd2 = FpXY_Fq_evaly(ddH2, j, T, p, vJ);
    1119             :     GEN h02, h12, h22, h03, h13, h23, h012, dh03, dh13, dh23, dh012;
    1120             :     GEN ddh03, ddh13, ddh23, ddh012;
    1121             :     GEN R1, dR1, ddR1, ddR2;
    1122          56 :     h02 = FqX_sqr(h0, T, p);
    1123          56 :     h12 = FqX_sqr(h1, T, p);
    1124          56 :     h22 = FqX_sqr(h2, T, p);
    1125          56 :     h03 = FqX_mul(h0, h02, T, p);
    1126          56 :     h13 = FqX_mul(h1, h12, T, p);
    1127          56 :     h23 = FqX_mul(h2, h22, T, p);
    1128          56 :     h012 = FqX_mul(FqX_mul(h0, h1, T, p), h2, T, p);
    1129          56 :     dh03 = FqX_mul(FqX_mulu(d0, 3, T, p), h02, T, p);
    1130          56 :     dh13 = FqX_mul(FqX_mulu(d1, 3, T, p), h12, T, p);
    1131          56 :     dh23 = FqX_mul(FqX_mulu(d2, 3, T, p), h22, T, p);
    1132          56 :     dh012 = FqX_add(FqX_add(FqX_mul(FqX_mul(d0, h1, T, p), h2, T, p), FqX_mul(FqX_mul(h0, d1, T, p), h2, T, p), T, p), FqX_mul(FqX_mul(h0, h1, T, p), d2, T, p), T, p);
    1133          56 :     R1 = FqX_sub(h13, FqX_mulu(h012, 3, T, p), T, p);
    1134          56 :     R = FqX_add(FqX_add(FqX_Fq_mul(RgX_shift_shallow(h23, t2), Fq_sqr(j, T, p), T, p), FqX_Fq_mul(RgX_shift_shallow(R1, 1), j, T, p), T, p), RgX_shift_shallow(h03, t0), T, p);
    1135          56 :     dR1 = FqX_sub(dh13, FqX_mulu(dh012, 3, T, p), T, p);
    1136          56 :     dR = FqX_add(FqX_add(RgX_shift_shallow(FqX_add(FqX_Fq_mul(dh23, Fq_sqr(j, T, p), T, p), FqX_Fq_mul(h23, Fq_mulu(j, 2, T, p), T, p), T, p), t2), RgX_shift_shallow(FqX_add(FqX_Fq_mul(dR1, j, T, p), R1, T, p), 1), T, p), RgX_shift_shallow(dh03, t0), T, p);
    1137          56 :     ddh03 = FqX_mulu(FqX_add(FqX_mul(dd0, h02, T, p), FqX_mul(FqX_mulu(FqX_sqr(d0, T, p), 2, T, p), h0, T, p), T, p), 3, T, p);
    1138          56 :     ddh13 = FqX_mulu(FqX_add(FqX_mul(dd1, h12, T, p), FqX_mul(FqX_mulu(FqX_sqr(d1, T, p), 2, T, p), h1, T, p), T, p), 3, T, p);
    1139          56 :     ddh23 = FqX_mulu(FqX_add(FqX_mul(dd2, h22, T, p), FqX_mul(FqX_mulu(FqX_sqr(d2, T, p), 2, T, p), h2, T, p), T, p), 3, T, p);
    1140          56 :     ddh012 = FqX_add(FqX_add(FqX_add(FqX_mul(FqX_mul(dd0, h1, T, p), h2, T, p), FqX_mul(FqX_mul(h0, dd1, T, p), h2, T, p), T, p), FqX_mul(FqX_mul(h0, h1, T, p), dd2, T, p), T, p), FqX_mulu(FqX_add(FqX_add(FqX_mul(FqX_mul(d0, d1, T, p), h2, T, p), FqX_mul(FqX_mul(d0, h1, T, p), d2, T, p), T, p), FqX_mul(FqX_mul(h0, d1, T, p), d2, T, p), T, p), 2, T, p), T, p);
    1141          56 :     ddR1 = FqX_sub(ddh13, FqX_mulu(ddh012, 3, T, p), T, p);
    1142          56 :     ddR2 = FqX_add(FqX_add(FqX_Fq_mul(ddh23, Fq_sqr(j, T, p), T, p), FqX_Fq_mul(dh23, Fq_mulu(j, 4, T, p), T, p), T, p), FqX_mulu(h23, 2, T, p), T, p);
    1143          56 :     ddR = FqX_add(FqX_add(RgX_shift_shallow(ddR2, t2), RgX_shift_shallow(FqX_add(FqX_mulu(dR1, 2, T, p), FqX_Fq_mul(ddR1, j, T, p), T, p), 1), T, p), RgX_shift_shallow(ddh03, t0), T, p);
    1144          56 :     return gerepilecopy(av, mkvec3(R ,dR ,ddR));
    1145             :   }
    1146             : }
    1147             : 
    1148             : static GEN
    1149       11606 : meqn_j(struct meqn *MEQN, GEN j, long ell, GEN T, GEN p)
    1150             : {
    1151       11606 :   if (MEQN->type=='J')
    1152             :   {
    1153          70 :     MEQN->eval = Fq_polmodular_eval(MEQN->eq, j, ell, T, p, MEQN->vy);
    1154          70 :     return gel(MEQN->eval, 1);
    1155             :   }
    1156             :   else
    1157       11536 :     return FqXY_evalx(MEQN->eq, j, T, p);
    1158             : }
    1159             : 
    1160             : static GEN
    1161          49 : find_isogenous_from_J(GEN a4, GEN a6, ulong ell, struct meqn *MEQN, GEN g, GEN T, GEN pp, long e)
    1162             : {
    1163          49 :   pari_sp ltop = avma;
    1164          49 :   GEN meqn = MEQN->eval;
    1165          49 :   GEN p = e==1 ? pp: powiu(pp, e);
    1166             :   GEN h, a4t, a6t;
    1167             :   GEN C4, C6, C4t, C6t;
    1168             :   GEN j, jp, jtp, jtp2, jtp3;
    1169             :   GEN Py, Pxy, Pyy, Pxj, Pyj, Pxxj, Pxyj, Pyyj;
    1170             :   GEN s0, s1, s2, s3;
    1171             :   GEN den, D, co, cot, c0, p_1;
    1172          49 :   if (signe(g) == 0 || signe(Fq_sub(g, utoi(1728), T, p)) == 0)
    1173             :   {
    1174           0 :     if (DEBUGLEVEL>0) err_printf("[J: g=%ld]",signe(g)==0 ?0: 1728);
    1175           0 :     return gc_NULL(ltop);
    1176             :   }
    1177          49 :   C4 = Fq_mul(a4, stoi(-48), T, p);
    1178          49 :   C6 = Fq_mul(a6, stoi(-864), T, p);
    1179          49 :   if (signe(C4)==0 || signe(C6)==0)
    1180             :   {
    1181           0 :     if (DEBUGLEVEL>0) err_printf("[J: C%ld=0]",signe(C4)==0 ?4: 6);
    1182           0 :     return gc_NULL(ltop);
    1183             :   }
    1184          49 :   j = Zq_ellj(a4, a6, T, p, pp, e);
    1185          49 :   jp = Fq_mul(j, Zq_div(C6, C4, T, p, pp, e), T, p);
    1186          49 :   co = corr(C4, C6, T, p, pp, e);
    1187          49 :   Py = RgX_deriv(gel(meqn, 1));
    1188          49 :   Pxy = RgX_deriv(gel(meqn,2));
    1189          49 :   Pyy = RgX_deriv(Py);
    1190          49 :   Pxj = FqX_eval(gel(meqn, 2), g, T, p);
    1191          49 :   if (signe(Pxj)==0)
    1192             :   {
    1193           0 :     if (DEBUGLEVEL>0) err_printf("[J: Pxj=0]");
    1194           0 :     return gc_NULL(ltop);
    1195             :   }
    1196          49 :   Pyj = FqX_eval(Py, g, T, p);
    1197          49 :   Pxxj = FqX_eval(gel(meqn, 3), g, T, p);
    1198          49 :   Pxyj = FqX_eval(Pxy, g, T, p);
    1199          49 :   Pyyj = FqX_eval(Pyy, g, T, p);
    1200          49 :   jtp = Fq_div(Fq_mul(jp, Zq_div(Pxj, Pyj, T, p, pp, e), T, p),
    1201             :                utoineg(ell), T, p);
    1202          49 :   jtp2 = Fq_sqr(jtp,T,p);
    1203          49 :   jtp3 = Fq_mul(jtp,jtp2,T,p);
    1204          49 :   den = Fq_mul(Fq_sqr(g,T,p),Fq_sub(g,utoi(1728),T,p),T, p);
    1205          49 :   D  =  Zq_inv(den, T, pp, e);
    1206          49 :   C4t = Fq_mul(jtp2,Fq_mul(g, D, T, p), T, p);
    1207          49 :   C6t = Fq_mul(jtp3, D, T, p);
    1208          49 :   s0 = Fq_mul(Fq_sqr(jp, T, p), Pxxj, T, p);
    1209          49 :   s1 = Fq_mul(Fq_mulu(Fq_mul(jp,jtp,T,p),2*ell,T,p), Pxyj, T, p);
    1210          49 :   s2 = Fq_mul(Fq_mulu(jtp2,ell*ell,T,p), Pyyj, T, p);
    1211          49 :   s3 = Zq_div(Fq_add(s0, Fq_add(s1, s2, T, p), T, p),Fq_mul(jp, Pxj, T, p),T,p,pp,e);
    1212          49 :   cot = corr(C4t, C6t, T, p, pp, e);
    1213          49 :   c0 = Fq_sub(co,Fq_mulu(cot,ell,T,p),T,p);
    1214          49 :   p_1 = Fq_div(Fq_mulu(Fq_add(s3, c0, T, p),ell,T,p),stoi(-4),T,p);
    1215          49 :   a4a6t_from_J(&a4t, &a6t, ell, C4t, C6t, T, p);
    1216          49 :   h = find_kernel(a4, a6, ell, a4t, a6t, p_1, T, p, pp, e);
    1217          49 :   if (!h) return NULL;
    1218          49 :   return gerepilecopy(ltop, mkvec3(a4t, a6t, h));
    1219             : }
    1220             : 
    1221             : static GEN
    1222        7511 : find_isogenous(GEN a4,GEN a6, ulong ell, struct meqn *MEQN, GEN g, GEN T,GEN p)
    1223             : {
    1224        7511 :   ulong pp = itou_or_0(p);
    1225        7511 :   long e = pp ? ulogint(((ell+1)>>1)+1, pp) + ulogint(2*ell+4, pp) + 1: 1;
    1226        7511 :   if (signe(a4)==0 || signe(a6)==0)
    1227             :   {
    1228           7 :     if (DEBUGLEVEL>0) err_printf("[%c: j=%ld]",MEQN->type,signe(a4)==0 ?0: 1728);
    1229           7 :     return NULL;
    1230             :   }
    1231        7504 :   if (e > 1)
    1232             :   {
    1233          63 :     GEN pe = powiu(p, e);
    1234          63 :     GEN meqnj = meqn_j(MEQN, Zq_ellj(a4, a6, T, pe, p, e), ell, T, pe);
    1235          63 :     g = ZqX_liftroot(meqnj, g, T, p, e);
    1236             :   }
    1237        7504 :   switch(MEQN->type)
    1238             :   {
    1239        4879 :     case 'C': return find_isogenous_from_canonical(a4,a6,ell, MEQN, g, T,p,e);
    1240        2576 :     case 'A': return find_isogenous_from_Atkin(a4,a6,ell, MEQN, g, T,p,e);
    1241          49 :     default:  return find_isogenous_from_J(a4,a6,ell, MEQN, g, T,p,e);
    1242             :   }
    1243             : }
    1244             : 
    1245             : static GEN
    1246        6181 : FqX_homogenous_eval(GEN P, GEN A, GEN B, GEN T, GEN p)
    1247             : {
    1248        6181 :   long d = degpol(P), i, v = varn(A);
    1249        6181 :   GEN s =  scalar_ZX_shallow(gel(P, d+2), v), Bn = pol_1(v);
    1250       20454 :   for (i = d-1; i >= 0; i--)
    1251             :   {
    1252       14273 :     Bn = FqX_mul(Bn, B, T, p);
    1253       14273 :     s = FqX_add(FqX_mul(s, A, T, p), FqX_Fq_mul(Bn, gel(P,i+2), T, p), T, p);
    1254             :   }
    1255        6181 :   return s;
    1256             : }
    1257             : 
    1258             : static GEN
    1259        1295 : FqX_homogenous_div(GEN P, GEN Q, GEN A, GEN B, GEN T, GEN p)
    1260             : {
    1261        1295 :   GEN z = cgetg(3, t_RFRAC);
    1262        1295 :   long d = degpol(Q)-degpol(P);
    1263        1295 :   gel(z, 1) = FqX_homogenous_eval(P, A, B, T, p);
    1264        1295 :   gel(z, 2) = FqX_homogenous_eval(Q, A, B, T, p);
    1265        1295 :   if (d > 0)
    1266           0 :     gel(z, 1) = FqX_mul(gel(z, 1), FqX_powu(B, d, T, p), T, p);
    1267        1295 :   else if (d < 0)
    1268        1295 :     gel(z, 2) = FqX_mul(gel(z, 2), FqX_powu(B, -d, T, p), T, p);
    1269        1295 :   return z;
    1270             : }
    1271             : 
    1272             : static GEN
    1273        1540 : find_kernel_power(GEN Eba4, GEN Eba6, GEN Eca4, GEN Eca6, ulong ell, struct meqn *MEQN, GEN kpoly, GEN Ib, GEN T, GEN p)
    1274             : {
    1275        1540 :   pari_sp ltop = avma, btop;
    1276             :   GEN a4t, a6t, gtmp;
    1277        1540 :   GEN num_iso = FqX_numer_isog_abscissa(kpoly, Eba4, Eba6, T, p, 0);
    1278        1540 :   GEN mpoly = meqn_j(MEQN, Fq_ellj(Eca4, Eca6, T, p), ell, T, p);
    1279        1540 :   GEN mroots = FqX_roots(mpoly, T, p);
    1280        1540 :   GEN kpoly2 = FqX_sqr(kpoly, T, p);
    1281        1540 :   long i, l1 = lg(mroots);
    1282        1540 :   btop = avma;
    1283        2541 :   for (i = 1; i < l1; i++)
    1284             :   {
    1285             :     GEN h;
    1286        2303 :     GEN tmp = find_isogenous(Eca4, Eca6, ell, MEQN, gel(mroots, i), T, p);
    1287        2303 :     if (!tmp) return gc_NULL(ltop);
    1288        2296 :     a4t =  gel(tmp, 1);
    1289        2296 :     a6t =  gel(tmp, 2);
    1290        2296 :     gtmp = gel(tmp, 3);
    1291             : 
    1292             :     /*check that the kernel kpoly is the good one */
    1293        2296 :     h = FqX_homogenous_eval(gtmp, num_iso, kpoly2, T, p);
    1294        2296 :     if (signe(Fq_elldivpolmod(Eba4, Eba6, ell, h, T, p)))
    1295             :     {
    1296        1295 :       GEN Ic = FqX_homogenous_div(num_iso,kpoly2, numer_i(Ib),denom_i(Ib), T,p);
    1297        1295 :       GEN kpoly_new = FqX_homogenous_eval(gtmp,   numer_i(Ic),denom_i(Ic), T,p);
    1298        1295 :       return gerepilecopy(ltop, mkvecn(5, a4t, a6t, kpoly_new, gtmp, Ic));
    1299             :     }
    1300        1001 :     set_avma(btop);
    1301             :   }
    1302         238 :   return gc_NULL(ltop);
    1303             : }
    1304             : 
    1305             : /****************************************************************************/
    1306             : /*                                  TRACE                                   */
    1307             : /****************************************************************************/
    1308             : enum mod_type {MTcm, MTpathological, MTAtkin, MTElkies, MTone_root, MTroots};
    1309             : 
    1310             : static GEN
    1311         678 : Flxq_study_eqn(GEN mpoly, GEN T, ulong p, long *pt_dG, long *pt_r)
    1312             : {
    1313         678 :   GEN Xq = FlxqX_Frobenius(mpoly, T, p);
    1314         678 :   GEN G  = FlxqX_gcd(FlxX_sub(Xq, pol_x(0), p), mpoly, T, p);
    1315         678 :   *pt_dG = degpol(G);
    1316         678 :   if (!*pt_dG) { *pt_r = FlxqX_ddf_degree(mpoly, Xq, T, p); return NULL; }
    1317         410 :   return gel(FlxqX_roots(G, T, p), 1);
    1318             : }
    1319             : 
    1320             : static GEN
    1321        8988 : Fp_study_eqn(GEN mpoly, GEN p, long *pt_dG, long *pt_r)
    1322             : {
    1323        8988 :   GEN T  = FpX_get_red(mpoly, p);
    1324        8988 :   GEN XP = FpX_Frobenius(T, p);
    1325        8988 :   GEN G  = FpX_gcd(FpX_sub(XP, pol_x(0), p), mpoly, p);
    1326        8988 :   *pt_dG = degpol(G);
    1327        8988 :   if (!*pt_dG) { *pt_r = FpX_ddf_degree(T, XP, p); return NULL; }
    1328        4732 :   return FpX_oneroot(G, p);
    1329             : }
    1330             : 
    1331             : static GEN
    1332        9989 : Fq_study_eqn(GEN mpoly, GEN T, GEN p, long *pt_dG, long *pt_r)
    1333             : {
    1334             :   GEN G;
    1335        9989 :   if (!T) return Fp_study_eqn(mpoly, p, pt_dG, pt_r);
    1336        1001 :   if (lgefint(p)==3)
    1337             :   {
    1338         678 :     ulong pp = p[2];
    1339         678 :     GEN Tp = ZXT_to_FlxT(T,pp);
    1340         678 :     GEN mpolyp = ZXX_to_FlxX(mpoly,pp,get_FpX_var(T));
    1341         678 :     G = Flxq_study_eqn(mpolyp, Tp, pp, pt_dG, pt_r);
    1342         678 :     return G ? Flx_to_ZX(G): NULL;
    1343             :   }
    1344             :   else
    1345             :   {
    1346         323 :     GEN Xq = FpXQX_Frobenius(mpoly, T, p);
    1347         323 :     G  = FpXQX_gcd(FpXX_sub(Xq, pol_x(0), p), mpoly, T, p);
    1348         323 :     *pt_dG = degpol(G);
    1349         323 :     if (!*pt_dG) { *pt_r = FpXQX_ddf_degree(mpoly, Xq, T, p); return NULL; }
    1350         136 :     return gel(FpXQX_roots(G, T, p), 1);
    1351             :   }
    1352             : }
    1353             : 
    1354             : /* Berlekamp variant */
    1355             : static GEN
    1356       10003 : study_modular_eqn(long ell, GEN mpoly, GEN T, GEN p, enum mod_type *mt, long *ptr_r)
    1357             : {
    1358       10003 :   pari_sp ltop = avma;
    1359       10003 :   GEN g = gen_0;
    1360       10003 :   *ptr_r = 0; /*gcc -Wall*/
    1361       10003 :   if (!FqX_is_squarefree(mpoly, T, p)) *mt = MTcm;
    1362             :   else
    1363             :   {
    1364             :     long dG;
    1365        9989 :     g = Fq_study_eqn(mpoly, T, p, &dG, ptr_r);
    1366        9989 :     switch(dG)
    1367             :     {
    1368        4711 :       case 0:  *mt = MTAtkin; break;
    1369         539 :       case 1:  *mt = MTone_root; break;
    1370        4669 :       case 2:  *mt = MTElkies;   break;
    1371          70 :       default: *mt = (dG == ell + 1)? MTroots: MTpathological;
    1372             :     }
    1373             :   }
    1374       10003 :   if (DEBUGLEVEL) switch(*mt)
    1375             :   {
    1376           0 :     case MTone_root: err_printf("One root\t"); break;
    1377           0 :     case MTElkies: err_printf("Elkies\t"); break;
    1378           0 :     case MTroots: err_printf("l+1 roots\t"); break;
    1379           0 :     case MTAtkin: err_printf("Atkin\t"); break;
    1380           0 :     case MTpathological: err_printf("Pathological\n"); break;
    1381           0 :     case MTcm: err_printf("CM\t"); break;
    1382             :   }
    1383       10003 :   return g ? gerepilecopy(ltop, g): NULL;
    1384             : }
    1385             : 
    1386             : /*Returns the trace modulo ell^k when ell is an Elkies prime */
    1387             : static GEN
    1388        5208 : find_trace_Elkies_power(GEN a4, GEN a6, ulong ell, long *pt_k, struct meqn *MEQN, GEN g, GEN tr, GEN q, GEN T, GEN p, long smallfact, pari_timer *ti)
    1389             : {
    1390        5208 :   pari_sp ltop = avma, btop;
    1391             :   GEN tmp, Eba4, Eba6, Eca4, Eca6, Ib, kpoly;
    1392        5208 :   long k = *pt_k;
    1393        5208 :   ulong lambda, ellk = upowuu(ell, k), pellk = umodiu(q, ellk);
    1394             :   long cnt;
    1395             : 
    1396        5208 :   if (DEBUGLEVEL) { err_printf("mod %ld", ell); }
    1397        5208 :   Eba4 = a4;
    1398        5208 :   Eba6 = a6;
    1399        5208 :   tmp = find_isogenous(a4,a6, ell, MEQN, g, T, p);
    1400        5208 :   if (!tmp) return gc_NULL(ltop);
    1401        5166 :   Eca4 =  gel(tmp, 1);
    1402        5166 :   Eca6 =  gel(tmp, 2);
    1403        5166 :   kpoly = gel(tmp, 3);
    1404        5166 :   Ib = pol_x(0);
    1405        5166 :   lambda = tr ? find_eigen_value_oneroot(a4, a6, ell, tr, kpoly, T, p):
    1406        4662 :                 find_eigen_value_power(a4, a6, ell, 1, 1, kpoly, T, p);
    1407        5166 :   if (DEBUGLEVEL>1) err_printf(" [%ld ms]", timer_delay(ti));
    1408        5166 :   if (smallfact && smallfact%(long)ell!=0)
    1409             :   {
    1410         378 :     ulong pell = pellk%ell;
    1411         378 :     ulong ap = Fl_add(lambda, Fl_div(pell, lambda, ell), ell);
    1412         378 :     if (Fl_sub(pell, ap, ell)==ell-1) { set_avma(ltop); return mkvecsmall(ap); }
    1413         364 :     if (smallfact < 0 && Fl_add(pell, ap, ell)==ell-1) { set_avma(ltop); return mkvecsmall(ap); }
    1414             :   }
    1415        5138 :   btop = avma;
    1416        6433 :   for (cnt = 2; cnt <= k; cnt++)
    1417             :   {
    1418        1540 :     GEN tmp = find_kernel_power(Eba4, Eba6, Eca4, Eca6, ell, MEQN, kpoly, Ib, T, p);
    1419        1540 :     if (!tmp) { k = cnt-1; break; }
    1420        1295 :     if (DEBUGLEVEL) err_printf(", %Ps", powuu(ell, cnt));
    1421        1295 :     lambda = find_eigen_value_power(a4, a6, ell, cnt, lambda, gel(tmp,3), T, p);
    1422        1295 :     Eba4 = Eca4;
    1423        1295 :     Eba6 = Eca6;
    1424        1295 :     Eca4 = gel(tmp,1);
    1425        1295 :     Eca6 = gel(tmp,2);
    1426        1295 :     kpoly = gel(tmp,4);
    1427        1295 :     Ib = gel(tmp, 5);
    1428        1295 :     if (gc_needed(btop, 1))
    1429             :     {
    1430           0 :       if(DEBUGMEM>1) pari_warn(warnmem,"find_trace_Elkies_power");
    1431           0 :       gerepileall(btop, 6, &Eba4, &Eba6, &Eca4, &Eca6, &kpoly, &Ib);
    1432             :     }
    1433        1295 :     if (DEBUGLEVEL>1) err_printf(" [%ld ms]", timer_delay(ti));
    1434             :   }
    1435        5138 :   set_avma(ltop);
    1436        5138 :   ellk = upowuu(ell, k);
    1437        5138 :   pellk = umodiu(q, ellk);
    1438        5138 :   *pt_k = k;
    1439        5138 :   return mkvecsmall(Fl_add(lambda, Fl_div(pellk, lambda, ellk), ellk));
    1440             : }
    1441             : 
    1442             : /*Returns the possible values of the trace when ell is an Atkin prime, */
    1443             : /*given r the splitting degree of the modular equation at J = E.j */
    1444             : static GEN
    1445        4711 : find_trace_Atkin(ulong ell, long r, GEN q)
    1446             : {
    1447        4711 :   pari_sp ltop = avma;
    1448        4711 :   long nval = 0;
    1449        4711 :   ulong teta, pell = umodiu(q, ell), invp = Fl_inv(pell, ell);
    1450        4711 :   GEN val_pos = cgetg(1+ell, t_VECSMALL), P = gel(factoru(r), 1);
    1451        4711 :   GEN S = mkvecsmall4(0, pell, 0, 1);
    1452        4711 :   GEN U = mkvecsmall3(0, ell-1, 0);
    1453        4711 :   pari_sp btop = avma;
    1454        4711 :   if (r==2 && krouu(ell-pell, ell) < 0)
    1455         707 :     val_pos[++nval] = 0;
    1456       92099 :   for (teta = 1; teta < ell; teta++, set_avma(btop))
    1457             :   {
    1458       87388 :     ulong disc = Fl_sub(Fl_sqr(teta,ell), Fl_mul(4UL,pell,ell), ell);
    1459             :     GEN a;
    1460       87388 :     if (krouu(disc, ell) >= 0) continue;
    1461       43162 :     S[3] = Fl_neg(teta, ell);
    1462       43162 :     U[3] = Fl_mul(invp, teta, ell);
    1463       43162 :     a = Flxq_powu(U, r/P[1], S, ell);
    1464       43162 :     if (!Flx_equal1(a) && Flx_equal1(Flxq_powu(a, P[1], S, ell)))
    1465             :     {
    1466       29260 :       pari_sp av = avma;
    1467       29260 :       long i, l=lg(P);
    1468       49924 :       for (i = 2; i < l; i++, set_avma(av))
    1469       26250 :         if (Flx_equal1(Flxq_powu(U, r/P[i], S, ell))) break;
    1470       29260 :       if (i==l) val_pos[++nval] = teta;
    1471             :     }
    1472             :   }
    1473        4711 :   return gerepileupto(ltop, vecsmall_shorten(val_pos, nval));
    1474             : }
    1475             : 
    1476             : /*Returns the possible traces when there is only one root */
    1477             : static GEN
    1478         539 : find_trace_one_root(ulong ell, GEN q)
    1479             : {
    1480         539 :   ulong a = Fl_double(Fl_sqrt(umodiu(q,ell), ell), ell);
    1481         539 :   return mkvecsmall2(a, ell - a);
    1482             : }
    1483             : 
    1484             : static GEN
    1485          70 : find_trace_lp1_roots(long ell, GEN q)
    1486             : {
    1487          70 :   ulong ell2 = ell * ell, pell = umodiu(q, ell2);
    1488          70 :   ulong a  = Fl_sqrt(pell%ell, ell);
    1489          70 :   ulong pa = Fl_add(Fl_div(pell, a, ell2), a, ell2);
    1490          70 :   return mkvecsmall2(pa, ell2 - pa);
    1491             : }
    1492             : 
    1493             : /*ell odd prime; trace modulo ell^k: [], [t] or [t1,...,td] */
    1494             : static GEN
    1495       10003 : find_trace(GEN a4, GEN a6, GEN j, ulong ell, GEN q, GEN T, GEN p, long *ptr_kt,
    1496             :   long smallfact, long vx, long vy)
    1497             : {
    1498       10003 :   pari_sp ltop = avma;
    1499             :   GEN g, meqnj, tr, tr2;
    1500             :   long kt, r;
    1501             :   enum mod_type mt;
    1502             :   struct meqn MEQN;
    1503             :   pari_timer ti;
    1504             : 
    1505       10003 :   kt = maxss((long)(log(expi(q)*M_LN2)/log((double)ell)), 1);
    1506       10003 :   if (DEBUGLEVEL)
    1507           0 :   { err_printf("SEA: Prime %5ld ", ell); timer_start(&ti); }
    1508       10003 :   get_modular_eqn(&MEQN, ell, vx, vy);
    1509       10003 :   meqnj = meqn_j(&MEQN, j, ell, T, p);
    1510       10003 :   g = study_modular_eqn(ell, meqnj, T, p, &mt, &r);
    1511             :   /* If l is an Elkies prime, search for a factor of the l-division polynomial.
    1512             :   * Then deduce the trace by looking for eigenvalues of the Frobenius by
    1513             :   * computing modulo this factor */
    1514       10003 :   switch (mt)
    1515             :   {
    1516         539 :   case MTone_root:
    1517         539 :     tr2 = find_trace_one_root(ell, q);
    1518         539 :     tr = find_trace_Elkies_power(a4,a6,ell, &kt, &MEQN, g, tr2, q, T, p, smallfact, &ti);
    1519         539 :     if (!tr) { tr = tr2; kt = 1; }
    1520         539 :     break;
    1521        4669 :   case MTElkies:
    1522             :     /* Contrary to MTone_root, may look mod higher powers of ell */
    1523        4669 :     if (abscmpiu(p, 2*ell+3) <= 0)
    1524          49 :       kt = 1; /* Not implemented in this case */
    1525        4669 :     tr = find_trace_Elkies_power(a4,a6,ell, &kt, &MEQN, g, NULL, q, T, p, smallfact, &ti);
    1526        4669 :     if (!tr)
    1527             :     {
    1528           7 :       if (DEBUGLEVEL) err_printf("[fail]\n");
    1529           7 :       return gc_NULL(ltop);
    1530             :     }
    1531        4662 :     break;
    1532          70 :   case MTroots:
    1533          70 :     tr = find_trace_lp1_roots(ell, q);
    1534          70 :     kt = 2;
    1535          70 :     break;
    1536        4711 :   case MTAtkin:
    1537        4711 :     tr = find_trace_Atkin(ell, r, q);
    1538        4711 :     if (lg(tr)==1) pari_err_PRIME("ellap",p);
    1539        4711 :     kt = 1;
    1540        4711 :     break;
    1541          14 :   case MTcm:
    1542             :     {
    1543          14 :       long D = find_CM(ell, j, T, p);
    1544          14 :       GEN C = Fq_ellcard_CM(D, a4, a6, T, p);
    1545          14 :       if (DEBUGLEVEL>1) err_printf(" D=%ld [%ld ms]\n", D, timer_delay(&ti));
    1546          14 :       return gc_const(ltop, C);
    1547             :     }
    1548           0 :   default: /* case MTpathological: */
    1549           0 :     return gc_NULL(ltop);
    1550             :   }
    1551        9982 :   if (DEBUGLEVEL) {
    1552           0 :     long n = lg(tr)-1;
    1553           0 :     if (n > 1 || mt == MTAtkin)
    1554             :     {
    1555           0 :       err_printf("%3ld trace(s)",n);
    1556           0 :       if (DEBUGLEVEL>1) err_printf(" [%ld ms]", timer_delay(&ti));
    1557             :     }
    1558           0 :     if (n > 1) err_printf("\n");
    1559             :   }
    1560        9982 :   *ptr_kt = kt;
    1561        9982 :   return gerepileupto(ltop, tr);
    1562             : }
    1563             : 
    1564             : /* A partition of compile_atkin in baby and giant is represented as the binary
    1565             :    developpement of an integer; if the i-th bit is 1, the i-th prime in
    1566             :    compile-atkin is a baby. The optimum is obtained when the ratio between
    1567             :    the number of possibilities for traces modulo giants (p_g) and babies (p_b)
    1568             :    is near 3/4. */
    1569             : static long
    1570         910 : separation(GEN cnt)
    1571             : {
    1572             :   pari_sp btop;
    1573         910 :   long k = lg(cnt)-1, l = (1L<<k)-1, best_i, i, j;
    1574             :   GEN best_r, P, P3, r;
    1575             : 
    1576         910 :   P = gen_1;
    1577        4550 :   for (j = 1; j <= k; ++j) P = mulis(P, cnt[j]);
    1578             :   /* p_b * p_g = P is constant */
    1579         910 :   P3 = mulsi(3, P);
    1580         910 :   btop = avma;
    1581         910 :   best_i = 0;
    1582         910 :   best_r = P3;
    1583       44282 :   for (i = 1; i < l; i++)
    1584             :   {
    1585             :     /* scan all possibilities */
    1586       43463 :     GEN p_b = gen_1;
    1587      415947 :     for (j = 0; j < k; j++)
    1588      372484 :       if (i & (1L<<j)) p_b = mulis(p_b, cnt[1+j]);
    1589       43463 :     r = subii(shifti(sqri(p_b), 2), P3); /* (p_b/p_g - 3/4)*4*P */
    1590       43463 :     if (!signe(r)) { best_i = i; break; }
    1591       43372 :     if (abscmpii(r, best_r) < 0) { best_i = i; best_r = r; }
    1592       43372 :     if (gc_needed(btop, 1))
    1593           0 :       best_r = gerepileuptoint(btop, best_r);
    1594             :   }
    1595         910 :   return best_i;
    1596             : }
    1597             : 
    1598             : /* x VEC defined modulo P (= *P), y VECSMALL modulo q, (q,P) = 1. */
    1599             : /* Update in place:
    1600             :  *   x to vector mod q P congruent to x mod P (resp. y mod q). */
    1601             : /*   P ( <-- qP ) */
    1602             : static void
    1603        1820 : multiple_crt(GEN x, GEN y, GEN q, GEN P)
    1604             : {
    1605        1820 :   pari_sp ltop = avma, av;
    1606        1820 :   long i, j, k, lx = lg(x)-1, ly = lg(y)-1;
    1607             :   GEN  a1, a2, u, v, A2X;
    1608        1820 :   (void)bezout(P,q,&u,&v);
    1609        1820 :   a1 = mulii(P,u);
    1610        1820 :   a2 = mulii(q,v); A2X = ZC_Z_mul(x, a2);
    1611        1820 :   av = avma; affii(mulii(P,q), P);
    1612       73010 :   for (i = 1, k = 1; i <= lx; i++, set_avma(av))
    1613             :   {
    1614       71190 :     GEN a2x = gel(A2X,i);
    1615     1194718 :     for (j = 1; j <= ly; ++j)
    1616             :     {
    1617     1123528 :       GEN t = Fp_add(Fp_mulu(a1, y[j], P), a2x, P);
    1618     1123528 :       affii(t, gel(x, k++));
    1619             :     }
    1620             :   }
    1621        1820 :   setlg(x, k); set_avma(ltop);
    1622        1820 : }
    1623             : 
    1624             : /****************************************************************************/
    1625             : /*                              MATCH AND SORT                              */
    1626             : /****************************************************************************/
    1627             : 
    1628             : static GEN
    1629        1820 : possible_traces(GEN compile, GEN mask, GEN *P, int larger)
    1630             : {
    1631        1820 :   GEN V, Pfinal = gen_1, C = shallowextract(compile, mask);
    1632        1820 :   long i, lfinal = 1, lC = lg(C), lP;
    1633        1820 :   pari_sp av = avma;
    1634             : 
    1635        5460 :   for (i = 1; i < lC; i++)
    1636             :   {
    1637        3640 :     GEN c = gel(C,i), t;
    1638        3640 :     Pfinal = mulii(Pfinal, gel(c,1));
    1639        3640 :     t = muluu(lfinal, lg(gel(c,2))-1);
    1640        3640 :     lfinal = itou(t);
    1641             :   }
    1642        1820 :   Pfinal = gerepileuptoint(av, Pfinal);
    1643        1820 :   if (larger)
    1644         910 :     lP = lgefint(shifti(Pfinal,1));
    1645             :   else
    1646         910 :     lP = lgefint(Pfinal);
    1647        1820 :   lfinal++;
    1648             :   /* allocate room for final result */
    1649        1820 :   V = cgetg(lfinal, t_VEC);
    1650     1061382 :   for (i = 1; i < lfinal; i++) gel(V,i) = cgeti(lP);
    1651             : 
    1652             :   {
    1653        1820 :     GEN c = gel(C,1), v = gel(c,2);
    1654        1820 :     long l = lg(v);
    1655        9044 :     for (i = 1; i < l; i++) affsi(v[i], gel(V,i));
    1656        1820 :     setlg(V, l); affii(gel(c,1), Pfinal); /* reset Pfinal */
    1657             :   }
    1658        3640 :   for (i = 2; i < lC; i++)
    1659             :   {
    1660        1820 :     GEN c = gel(C,i);
    1661        1820 :     multiple_crt(V, gel(c,2), gel(c,1), Pfinal); /* Pfinal updated! */
    1662             :   }
    1663        1820 :   *P = Pfinal; return V;
    1664             : }
    1665             : 
    1666             : static GEN
    1667      459375 : cost(long mask, GEN cost_vec)
    1668             : {
    1669      459375 :   pari_sp ltop = avma;
    1670             :   long i;
    1671      459375 :   GEN c = gen_1;
    1672     7173831 :   for (i = 1; i < lg(cost_vec); i++)
    1673     6714456 :     if (mask&(1L<<(i-1)))
    1674     2976967 :       c = mulis(c, cost_vec[i]);
    1675      459375 :   return gerepileuptoint(ltop, c);
    1676             : }
    1677             : 
    1678             : static GEN
    1679      369894 : value(long mask, GEN atkin, long k)
    1680             : {
    1681      369894 :   pari_sp ltop = avma;
    1682             :   long i;
    1683      369894 :   GEN c = gen_1;
    1684     5777625 :   for (i = 1; i <= k; i++)
    1685     5407731 :     if (mask&(1L<<(i-1)))
    1686     2386237 :       c = mulii(c, gmael(atkin, i, 1));
    1687      369894 :   return gerepileuptoint(ltop, c);
    1688             : }
    1689             : 
    1690             : static void
    1691      182616 : set_cost(GEN B, long b, GEN cost_vec, long *pi)
    1692             : {
    1693      182616 :   pari_sp av = avma;
    1694      182616 :   GEN costb = cost(b, cost_vec);
    1695      182616 :   long i = *pi;
    1696      250474 :   while (cmpii(costb, cost(B[i], cost_vec)) < 0) --i;
    1697      182616 :   B[++i] = b;
    1698      182616 :   *pi = i; set_avma(av);
    1699      182616 : }
    1700             : 
    1701             : static GEN
    1702        1925 : get_lgatkin(GEN compile_atkin, long k)
    1703             : {
    1704        1925 :   GEN v = cgetg(k+1, t_VECSMALL);
    1705             :   long j;
    1706       10248 :   for (j = 1; j <= k; ++j) v[j] = lg(gmael(compile_atkin, j, 2))-1;
    1707        1925 :   return v;
    1708             : }
    1709             : 
    1710             : static GEN
    1711        1015 : champion(GEN atkin, long k, GEN bound_champ)
    1712             : {
    1713        1015 :   const long two_k = 1L<<k;
    1714        1015 :   pari_sp ltop = avma;
    1715             :   long i, j, n, i1, i2;
    1716        1015 :   GEN B, Bp, cost_vec, res = NULL;
    1717             : 
    1718        1015 :   cost_vec = get_lgatkin(atkin, k);
    1719        1015 :   if (k == 1) return mkvec2(gen_1, utoipos(cost_vec[1]));
    1720             : 
    1721        1001 :   B  = zero_zv(two_k);
    1722        1001 :   Bp = zero_zv(two_k);
    1723        1001 :   Bp[2] = 1;
    1724        4669 :   for (n = 2, j = 2; j <= k; j++)
    1725             :   {
    1726             :     long b;
    1727        3668 :     i = 1;
    1728      173418 :     for (i1 = 2, i2 = 1; i1 <= n; )
    1729             :     {
    1730      169750 :       pari_sp av = avma;
    1731      169750 :       long b1 = Bp[i1], b2 = Bp[i2]|(1L<<(j-1));
    1732      169750 :       if (cmpii(value(b1, atkin, k), value(b2, atkin, k)) < 0)
    1733      169750 :         { b = b1; i1++; } else { b = b2; i2++; }
    1734      169750 :       set_avma(av);
    1735      169750 :       set_cost(B, b, cost_vec, &i);
    1736             :     }
    1737       16534 :     for ( ; i2 <= n; i2++)
    1738             :     {
    1739       12866 :       b = Bp[i2]|(1L<<(j-1));
    1740       12866 :       set_cost(B, b, cost_vec, &i);
    1741             :     }
    1742        3668 :     n = i;
    1743      122094 :     for (i = 1; i <= n; i++)
    1744      118426 :       Bp[i] = B[i];
    1745             :   }
    1746     9631069 :   for (i = 1; i <= two_k; i++)
    1747     9630068 :     if (B[i])
    1748             :     {
    1749       26285 :       GEN b = cost (B[i], cost_vec);
    1750       26285 :       GEN v = value(B[i], atkin, k);
    1751       26285 :       if (cmpii(v, bound_champ) <=0) continue;
    1752        5005 :       if (res && gcmp(b, gel(res, 2)) >=0) continue;
    1753        1001 :       res = mkvec2(utoi(B[i]), b);
    1754             :     }
    1755        1001 :   return gerepilecopy(ltop, res);
    1756             : }
    1757             : 
    1758             : static GEN
    1759        1820 : compute_diff(GEN v)
    1760             : {
    1761        1820 :   long i, l = lg(v) - 1;
    1762        1820 :   GEN diff = cgetg(l, t_VEC);
    1763     1059562 :   for (i = 1; i < l; i++) gel(diff, i) = subii(gel(v, i+1), gel(v, i));
    1764        1820 :   return ZV_sort_uniq_shallow(diff);
    1765             : }
    1766             : 
    1767             : static int
    1768       17276 : cmp_atkin(void*E, GEN a, GEN b)
    1769             : {
    1770       17276 :   long ta=typ(a)==t_INT, tb=typ(b)==t_INT, c;
    1771             :   (void) E;
    1772       17276 :   if (ta || tb) return ta-tb;
    1773        5670 :   c = lg(gel(a,2)) - lg(gel(b,2));
    1774        5670 :   if (c) return c;
    1775         847 :   return cmpii(gel(b,1), gel(a,1));
    1776             : }
    1777             : 
    1778             : static void
    1779        4109 : add_atkin(GEN atkin, GEN trace, long *nb)
    1780             : {
    1781        4109 :   long l = lg(atkin)-1;
    1782        4109 :   long i, k = gen_search(atkin, trace, NULL, cmp_atkin);
    1783        4109 :   if (k > 0 || (k = -k) > l) return;
    1784       79926 :   for (i = l; i > k; i--) gel(atkin,i) = gel(atkin,i-1);
    1785        4109 :   if (typ(gel(atkin,l))==t_INT) (*nb)++;
    1786        4109 :   gel(atkin,k) = trace;
    1787             : }
    1788             : 
    1789             : /* V = baby / giant, P = Pb / Pg */
    1790             : static GEN
    1791        1820 : BSGS_pre(GEN *pdiff, GEN V, GEN P, void *E, const struct bb_group *grp)
    1792             : {
    1793        1820 :   GEN diff = compute_diff(V);
    1794        1820 :   GEN pre = cgetg(lg(diff), t_VEC);
    1795        1820 :   long i, l = lg(diff);
    1796        1820 :   gel(pre, 1) = grp->pow(E, P, gel(diff, 1));
    1797             :   /* what we'd _really_ want here is a hashtable diff[i] -> pre[i]  */
    1798       39018 :   for (i = 2; i < l; i++)
    1799             :   {
    1800       37198 :     pari_sp av = avma;
    1801       37198 :     GEN d = subii(gel(diff, i), gel(diff, i-1));
    1802       37198 :     GEN Q = grp->mul(E, gel(pre, i-1), grp->pow(E, P, d));
    1803       37198 :     gel(pre, i) = gerepilecopy(av, Q);
    1804             :   }
    1805        1820 :   *pdiff = diff; return pre;
    1806             : }
    1807             : 
    1808             : /* u = trace_elkies, Mu = prod_elkies. Let caller collect garbage */
    1809             : /* Match & sort: variant from Lercier's thesis, section 11.2.3 */
    1810             : /* baby/giant/table updated in place: this routines uses
    1811             :  *   size(baby)+size(giant)+size(table)+size(table_ind) + O(log p)
    1812             :  * bits of stack */
    1813             : static GEN
    1814         966 : match_and_sort(GEN compile_atkin, GEN Mu, GEN u, GEN q, void *E, const struct bb_group *grp)
    1815             : {
    1816             :   pari_sp av1, av2;
    1817         966 :   GEN baby, giant, SgMb, Mb, Mg, den, Sg, dec_inf, div, pp1 = addiu(q,1);
    1818             :   GEN P, Pb, Pg, point, diff, pre, table, table_ind;
    1819         966 :   long best_i, i, lbaby, lgiant, k = lg(compile_atkin)-1;
    1820         966 :   GEN bound = sqrti(shifti(q, 2)), card;
    1821         966 :   const long lcard = 100;
    1822         966 :   long lq = lgefint(q), nbcard;
    1823             :   pari_timer ti;
    1824             : 
    1825         966 :   if (k == 1)
    1826             :   { /*only one Atkin prime, check the cardinality with random points */
    1827          56 :     GEN r = gel(compile_atkin, 1), r1 = gel(r,1), r2 = gel(r,2);
    1828          56 :     long l = lg(r2), j;
    1829          56 :     GEN card = cgetg(l, t_VEC), Cs2, C, U;
    1830          56 :     Z_chinese_pre(Mu, r1, &C,&U, NULL);
    1831          56 :     Cs2 = shifti(C, -1);
    1832         378 :     for (j = 1, i = 1; i < l; i++)
    1833             :     {
    1834         322 :       GEN t = Z_chinese_post(u, stoi(r2[i]), C, U, NULL);
    1835         322 :       t = Fp_center_i(t, C, Cs2);
    1836         322 :       if (abscmpii(t, bound) <= 0) gel(card, j++) = subii(pp1, t);
    1837             :     }
    1838          56 :     setlg(card, j);
    1839          56 :     return gen_select_order(card, E, grp);
    1840             :   }
    1841         910 :   if (DEBUGLEVEL>=2) timer_start(&ti);
    1842         910 :   av1 = avma;
    1843         910 :   best_i = separation( get_lgatkin(compile_atkin, k) );
    1844         910 :   set_avma(av1);
    1845             : 
    1846         910 :   baby  = possible_traces(compile_atkin, utoi(best_i), &Mb, 1);
    1847         910 :   giant = possible_traces(compile_atkin, subiu(int2n(k), best_i+1), &Mg, 0);
    1848         910 :   lbaby = lg(baby);
    1849         910 :   lgiant = lg(giant);
    1850         910 :   den = Fp_inv(Fp_mul(Mu, Mb, Mg), Mg);
    1851         910 :   av2 = avma;
    1852      622790 :   for (i = 1; i < lgiant; i++, set_avma(av2))
    1853      621880 :     affii(Fp_mul(gel(giant,i), den, Mg), gel(giant,i));
    1854         910 :   ZV_sort_inplace(giant);
    1855         910 :   Sg = Fp_mul(negi(u), den, Mg);
    1856         910 :   den = Fp_inv(Fp_mul(Mu, Mg, Mb), Mb);
    1857         910 :   dec_inf = divii(mulii(Mb,addii(Mg,shifti(Sg,1))), shifti(Mg,1));
    1858         910 :   togglesign(dec_inf); /* now, dec_inf = ceil(- (Mb/2 + Sg Mb/Mg) ) */
    1859         910 :   div = mulii(truedivii(dec_inf, Mb), Mb);
    1860         910 :   av2 = avma;
    1861      438592 :   for (i = 1; i < lbaby; i++, set_avma(av2))
    1862             :   {
    1863      437682 :     GEN b = addii(Fp_mul(Fp_sub(gel(baby,i), u, Mb), den, Mb), div);
    1864      437682 :     if (cmpii(b, dec_inf) < 0) b = addii(b, Mb);
    1865      437682 :     affii(b, gel(baby,i));
    1866             :   }
    1867         910 :   ZV_sort_inplace(baby);
    1868             : 
    1869         910 :   SgMb = mulii(Sg, Mb);
    1870         910 :   card = cgetg(lcard+1,t_VEC);
    1871       91910 :   for (i = 1; i <= lcard; i++) gel(card,i) = cgetipos(lq+1);
    1872             : 
    1873         910 :   av2 = avma;
    1874         910 : MATCH_RESTART:
    1875         910 :   set_avma(av2);
    1876         910 :   nbcard = 0;
    1877         910 :   P = grp->rand(E);
    1878         910 :   point = grp->pow(E,P, Mu);
    1879         910 :   Pb = grp->pow(E,point, Mg);
    1880         910 :   Pg = grp->pow(E,point, Mb);
    1881             :   /* Precomputation for babies */
    1882         910 :   pre = BSGS_pre(&diff, baby, Pb, E, grp);
    1883             : 
    1884             :   /*Now we compute the table of babies, this table contains only the */
    1885             :   /*lifted x-coordinate of the points in order to use less memory */
    1886         910 :   table = cgetg(lbaby, t_VECSMALL);
    1887         910 :   av1 = avma;
    1888             :   /* (p+1 - u - Mu*Mb*Sg) P - (baby[1]) Pb */
    1889         910 :   point = grp->pow(E,P, subii(subii(pp1, u), mulii(Mu, addii(SgMb, mulii(Mg, gel(baby,1))))));
    1890         910 :   table[1] = grp->hash(gel(point,1));
    1891      437682 :   for (i = 2; i < lbaby; i++)
    1892             :   {
    1893      436772 :     GEN d = subii(gel(baby, i), gel(baby, i-1));
    1894      436772 :     point =  grp->mul(E, point, grp->pow(E, gel(pre, ZV_search(diff, d)), gen_m1));
    1895      436772 :     table[i] = grp->hash(gel(point,1));
    1896      436772 :     if (gc_needed(av1,3))
    1897             :     {
    1898          19 :       if(DEBUGMEM>1) pari_warn(warnmem,"match_and_sort, baby = %ld", i);
    1899          19 :       point = gerepileupto(av1, point);
    1900             :     }
    1901             :   }
    1902         910 :   set_avma(av1);
    1903             :   /* Precomputations for giants */
    1904         910 :   pre = BSGS_pre(&diff, giant, Pg, E, grp);
    1905             : 
    1906             :   /* Look for a collision among the x-coordinates */
    1907         910 :   table_ind = vecsmall_indexsort(table);
    1908         910 :   table = perm_mul(table,table_ind);
    1909             : 
    1910         910 :   av1 = avma;
    1911         910 :   point = grp->pow(E, Pg, gel(giant, 1));
    1912         910 :   for (i = 1; ; i++)
    1913      620970 :   {
    1914             :     GEN d;
    1915      621880 :     long h = grp->hash(gel(point, 1));
    1916      621880 :     long s = zv_search(table, h);
    1917      621880 :     if (s) {
    1918        1820 :       while (table[s] == h && s) s--;
    1919        1820 :       for (s++; s < lbaby && table[s] == h; s++)
    1920             :       {
    1921         910 :         GEN B = gel(baby,table_ind[s]), G = gel(giant,i);
    1922         910 :         GEN GMb = mulii(G, Mb), BMg = mulii(B, Mg);
    1923         910 :         GEN Be = subii(subii(pp1, u), mulii(Mu, addii(SgMb, BMg)));
    1924         910 :         GEN Bp = grp->pow(E,P, Be);
    1925             :         /* p+1 - u - Mu (Sg Mb + GIANT Mb + BABY Mg) */
    1926         910 :         if (gequal(gel(Bp,1),gel(point,1)))
    1927             :         {
    1928         910 :           GEN card1 = subii(Be, mulii(Mu, GMb));
    1929         910 :           GEN card2 = addii(card1, mulii(mulsi(2,Mu), GMb));
    1930         910 :           if (abscmpii(subii(pp1, card1), bound) <= 0)
    1931         798 :             affii(card1, gel(card, ++nbcard));
    1932         910 :           if (nbcard >= lcard) goto MATCH_RESTART;
    1933         910 :           if (abscmpii(subii(pp1, card2), bound) <= 0)
    1934         490 :             affii(card2, gel(card, ++nbcard));
    1935         910 :           if (nbcard >= lcard) goto MATCH_RESTART;
    1936             :         }
    1937             :       }
    1938             :     }
    1939      621880 :     if (i==lgiant-1) break;
    1940      620970 :     d = subii(gel(giant, i+1), gel(giant, i));
    1941      620970 :     point = grp->mul(E,point, gel(pre, ZV_search(diff, d)));
    1942      620970 :     if (gc_needed(av1,3))
    1943             :     {
    1944          26 :       if(DEBUGMEM>1) pari_warn(warnmem,"match_and_sort, giant = %ld", i);
    1945          26 :       point = gerepileupto(av1, point);
    1946             :     }
    1947             :   }
    1948         910 :   setlg(card, nbcard+1);
    1949         910 :   if (DEBUGLEVEL>=2) timer_printf(&ti,"match_and_sort");
    1950         910 :   return gen_select_order(card, E, grp);
    1951             : }
    1952             : 
    1953             : static GEN
    1954        1015 : get_bound_bsgs(long lp)
    1955             : {
    1956             :   GEN B;
    1957        1015 :   if (lp <= 160)
    1958         980 :     B = divru(powru(dbltor(1.048), lp), 9);
    1959          35 :   else if (lp <= 192)
    1960          28 :     B = divrr(powru(dbltor(1.052), lp), dbltor(16.65));
    1961             :   else
    1962           7 :     B = mulrr(powru(dbltor(1.035), minss(lp,307)), dbltor(1.35));
    1963        1015 :   return mulru(B, 1000000);
    1964             : }
    1965             : 
    1966             : /* E is an elliptic curve defined over Z or over Fp in ellinit format, defined
    1967             :  * by the equation E: y^2 + a1*x*y + a2*y = x^3 + a2*x^2 + a4*x + a6
    1968             :  * p is a prime number
    1969             :  * set smallfact to stop whenever a small factor of the order, not dividing smallfact,
    1970             :  * is detected. Useful when searching for a good curve for cryptographic
    1971             :  * applications */
    1972             : GEN
    1973        1064 : Fq_ellcard_SEA(GEN a4, GEN a6, GEN q, GEN T, GEN p, long smallfact)
    1974             : {
    1975        1064 :   const long MAX_ATKIN = 21;
    1976        1064 :   pari_sp ltop = avma, btop;
    1977             :   long ell, i, nb_atkin, vx,vy;
    1978             :   GEN TR, TR_mod, compile_atkin, bound, bound_bsgs, champ;
    1979        1064 :   GEN prod_atkin = gen_1, max_traces = gen_0;
    1980             :   GEN j;
    1981        1064 :   double bound_gr = 1.;
    1982        1064 :   const double growth_factor = 1.26;
    1983             :   forprime_t TT;
    1984             : 
    1985        1064 :   j = Fq_ellj(a4, a6, T, p);
    1986        1064 :   if (signe(j) == 0 || signe(Fq_sub(j, utoi(1728), T, p)) == 0)
    1987           0 :     return T ? FpXQ_ellcard(Fq_to_FpXQ(a4, T, p), Fq_to_FpXQ(a6, T, p), T, p)
    1988          14 :              : Fp_ellcard(a4, a6, p);
    1989        1050 :   if (Fq_elljissupersingular(j, T, p))
    1990          21 :     return Fq_ellcard_supersingular(a4, a6, T, p);
    1991             :   /*First compute the trace modulo 2 */
    1992        1029 :   switch(FqX_nbroots(rhs(a4, a6, 0), T, p))
    1993             :   {
    1994          70 :   case 3: /* bonus time: 4 | #E(Fq) = q+1 - t */
    1995          70 :     i = mod4(q)+1; if (i > 2) i -= 4;
    1996          70 :     TR_mod = utoipos(4);
    1997          70 :     TR = stoi(i); break;
    1998         511 :   case 1:
    1999         511 :     TR_mod = gen_2;
    2000         511 :     TR = gen_0; break;
    2001         448 :   default : /* 0 */
    2002         448 :     TR_mod = gen_2;
    2003         448 :     TR = gen_1; break;
    2004             :   }
    2005        1029 :   if (odd(smallfact) && !mpodd(TR))
    2006             :   {
    2007          14 :     if (DEBUGLEVEL) err_printf("Aborting: #E(Fq) divisible by 2\n");
    2008          14 :     set_avma(ltop); return gen_0;
    2009             :   }
    2010        1015 :   vy = fetch_var();
    2011        1015 :   vx = fetch_var_higher();
    2012             : 
    2013             :   /* compile_atkin is a vector containing informations about Atkin primes,
    2014             :    * informations about Elkies primes lie in Mod(TR, TR_mod). */
    2015        1015 :   u_forprime_init(&TT, 3, ULONG_MAX);
    2016        1015 :   bound = sqrti(shifti(q, 4));
    2017        1015 :   bound_bsgs = get_bound_bsgs(expi(q));
    2018        1015 :   compile_atkin = zerovec(MAX_ATKIN); nb_atkin = 0;
    2019        1015 :   btop = avma;
    2020       10010 :   while ( (ell = u_forprime_next(&TT)) )
    2021             :   {
    2022       10010 :     long ellkt, kt = 1, nbtrace;
    2023             :     GEN trace_mod;
    2024       10017 :     if (absequalui(ell, p)) continue;
    2025       10003 :     trace_mod = find_trace(a4, a6, j, ell, q, T, p, &kt, smallfact, vx,vy);
    2026       10003 :     if (!trace_mod) continue;
    2027        9996 :     if (typ(trace_mod)==t_INT)
    2028             :     {
    2029          14 :       delete_var(); delete_var();
    2030        1015 :       return gerepileuptoint(ltop, trace_mod);
    2031             :     }
    2032        9982 :     nbtrace = lg(trace_mod) - 1;
    2033        9982 :     ellkt = (long)upowuu(ell, kt);
    2034        9982 :     if (nbtrace == 1)
    2035             :     {
    2036        5873 :       long t_mod_ellkt = trace_mod[1];
    2037        5873 :       if (smallfact && smallfact%ell!=0)
    2038             :       { /* does ell divide q + 1 - t ? */
    2039         385 :         long q_mod_ell_plus_one = umodiu(q,ell) + 1;
    2040         385 :         ulong  card_mod_ell = umodsu(q_mod_ell_plus_one - t_mod_ellkt, ell);
    2041         385 :         ulong tcard_mod_ell = 1;
    2042         385 :         if (card_mod_ell && smallfact < 0)
    2043         133 :           tcard_mod_ell = umodsu(q_mod_ell_plus_one + t_mod_ellkt, ell);
    2044         385 :         if (!card_mod_ell || !tcard_mod_ell)
    2045             :         {
    2046          28 :           if (DEBUGLEVEL)
    2047           0 :             err_printf("\nAborting: #E%s(Fq) divisible by %ld\n",
    2048             :                        tcard_mod_ell ? "" : "_twist", ell);
    2049          28 :           delete_var(); delete_var();
    2050          28 :           return gc_const(ltop, gen_0);
    2051             :         }
    2052             :       }
    2053        5845 :       (void)Z_incremental_CRT(&TR, t_mod_ellkt, &TR_mod, ellkt);
    2054        5845 :       if (DEBUGLEVEL)
    2055           0 :         err_printf(", missing %ld bits\n",expi(bound)-expi(TR_mod));
    2056             :     }
    2057             :     else
    2058             :     {
    2059        4109 :       add_atkin(compile_atkin, mkvec2(utoipos(ellkt), trace_mod), &nb_atkin);
    2060        4109 :       prod_atkin = value(-1, compile_atkin, nb_atkin);
    2061             :     }
    2062        9954 :     if (cmpii(mulii(TR_mod, prod_atkin), bound) > 0)
    2063             :     {
    2064             :       GEN bound_tr;
    2065        1057 :       if (!nb_atkin)
    2066             :       {
    2067           7 :         delete_var(); delete_var();
    2068           7 :         return gerepileuptoint(ltop, subii(addiu(q, 1), TR));
    2069             :       }
    2070        1050 :       bound_tr = mulrr(bound_bsgs, dbltor(bound_gr));
    2071        1050 :       bound_gr *= growth_factor;
    2072        1050 :       if (signe(max_traces))
    2073             :       {
    2074          84 :         max_traces = divis(muliu(max_traces,nbtrace), ellkt);
    2075          84 :         if (DEBUGLEVEL>=3)
    2076           0 :           err_printf("At least %Ps remaining possibilities.\n",max_traces);
    2077             :       }
    2078        1050 :       if (cmpir(max_traces, bound_tr) < 0)
    2079             :       {
    2080        1015 :         GEN bound_atkin = truedivii(bound, TR_mod);
    2081        1015 :         champ = champion(compile_atkin, nb_atkin, bound_atkin);
    2082        1015 :         max_traces = gel(champ,2);
    2083        1015 :         if (DEBUGLEVEL>=2)
    2084           0 :           err_printf("%Ps remaining possibilities.\n", max_traces);
    2085        1015 :         if (cmpir(max_traces, bound_tr) < 0)
    2086             :         {
    2087         966 :           GEN res, cat = shallowextract(compile_atkin, gel(champ,1));
    2088             :           const struct bb_group *grp;
    2089             :           void *E;
    2090         966 :           if (DEBUGLEVEL)
    2091           0 :             err_printf("Match and sort for %Ps possibilities.\n", max_traces);
    2092         966 :           delete_var();
    2093         966 :           delete_var();
    2094         966 :           grp = get_FqE_group(&E,a4,a6,T,p);
    2095         966 :           res = match_and_sort(cat, TR_mod, TR, q, E, grp);
    2096         966 :           return gerepileuptoint(ltop, res);
    2097             :         }
    2098             :       }
    2099             :     }
    2100        8981 :     if (gc_needed(btop, 1))
    2101           0 :       gerepileall(btop,5, &TR,&TR_mod, &compile_atkin, &max_traces, &prod_atkin);
    2102             :   }
    2103             :   return NULL;/*LCOV_EXCL_LINE*/
    2104             : }
    2105             : 
    2106             : GEN
    2107         973 : Fp_ellcard_SEA(GEN a4, GEN a6, GEN p, long smallfact)
    2108         973 : { return Fq_ellcard_SEA(a4, a6, p, NULL, p, smallfact); }

Generated by: LCOV version 1.16