Code coverage tests

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

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

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

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

Generated by: LCOV version 1.11