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

Generated by: LCOV version 1.13