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-bordeaux1.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 - modules - ellsea.c (source / functions) Hit Total Coverage
Test: PARI/GP v2.8.0 lcov report (development 17240-4611fa9) Lines: 954 1015 94.0 %
Date: 2014-12-21 Functions: 64 67 95.5 %
Legend: Lines: hit not hit | Branches: + taken - not taken # not executed Branches: 403 538 74.9 %

           Branch data     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 modular_eqn;
      37                 :            : 
      38                 :            : void
      39                 :       1243 : pari_init_seadata(void)  { modular_eqn = NULL; }
      40                 :            : void
      41         [ +  + ]:       1231 : pari_close_seadata(void) { if (modular_eqn) gunclone(modular_eqn); }
      42                 :            : 
      43                 :            : static int
      44                 :      24334 : FqX_equal(GEN x, GEN y) { return gequal(x,y); }
      45                 :            : 
      46                 :            : static int
      47                 :       1601 : FlxX_equal(GEN x, GEN y) { return gequal(x,y); }
      48                 :            : 
      49                 :            : static char *
      50                 :         21 : seadata_filename(ulong ell)
      51                 :         21 : { return stack_sprintf("%s/seadata/sea%ld", pari_datadir, ell); }
      52                 :            : 
      53                 :            : static GEN
      54                 :         21 : get_seadata(ulong ell)
      55                 :            : {
      56                 :         21 :   pari_sp av=avma;
      57                 :            :   GEN eqn;
      58                 :         21 :   char *s = seadata_filename(ell);
      59                 :         21 :   pariFILE *F = pari_fopengz(s);
      60         [ -  + ]:         21 :   if (!F) return NULL;
      61         [ +  - ]:         21 :   if (ell==0)
      62                 :            :   {
      63                 :         21 :     eqn = gp_readvec_stream(F->file);
      64                 :         21 :     pari_fclose(F);
      65                 :         21 :     modular_eqn = gclone(eqn);
      66                 :         21 :     avma=av;
      67                 :         21 :     return gen_0;
      68                 :            :   } else {
      69                 :          0 :     eqn = gp_read_stream(F->file);
      70                 :          0 :     pari_fclose(F);
      71                 :         21 :     return eqn;
      72                 :            :   }
      73                 :            : }
      74                 :            : 
      75                 :            : /*Builds the modular equation corresponding to the vector list. Shallow */
      76                 :            : static GEN
      77                 :       6209 : list_to_pol(GEN list, long vx, long vy)
      78                 :            : {
      79                 :       6209 :   long i, l = lg(list);
      80                 :       6209 :   GEN P = cgetg(l, t_VEC);
      81         [ +  + ]:     133784 :   for (i = 1; i < l; i++)
      82                 :            :   {
      83                 :     127575 :     GEN L = gel(list,i);
      84         [ +  + ]:     127575 :     if (typ(L) == t_VEC) L = RgV_to_RgX_reverse(L, vy);
      85                 :     127575 :     gel(P, i) = L;
      86                 :            :   }
      87                 :       6209 :   return RgV_to_RgX_reverse(P, vx);
      88                 :            : }
      89                 :            : 
      90                 :            : struct meqn {
      91                 :            :   char type;
      92                 :            :   GEN eq;
      93                 :            :   long vx,vy;
      94                 :            : };
      95                 :            : 
      96                 :            : static int
      97                 :       6209 : get_modular_eqn(struct meqn *M, ulong ell, long vx, long vy)
      98                 :            : {
      99                 :            :   GEN eqn;
     100                 :       6209 :   long idx = uprimepi(ell)-1;
     101 [ +  + ][ -  + ]:       6209 :   if (!modular_eqn && !get_seadata(0)) pari_err_PACKAGE("seadata");
     102 [ +  - ][ +  - ]:       6209 :   if (idx && idx<lg(modular_eqn))
     103                 :       6209 :     eqn = gel(modular_eqn, idx);
     104                 :            :   else
     105                 :          0 :     eqn = get_seadata(ell);
     106         [ -  + ]:       6209 :   if (!eqn) { M->type = 0; M->eq = NULL; return 0; }
     107                 :       6209 :   M->type = *GSTR(gel(eqn, 2));
     108                 :       6209 :   M->eq = list_to_pol(gel(eqn, 3), vx, vy);
     109                 :       6209 :   M->vx = vx;
     110                 :       6209 :   M->vy = vy; return 1;
     111                 :            : }
     112                 :            : 
     113                 :            : static void
     114                 :          0 : err_modular_eqn(long ell)
     115                 :          0 : { pari_err_FILE("seadata file", seadata_filename(ell)); }
     116                 :            : 
     117                 :            : GEN
     118                 :          7 : ellmodulareqn(long ell, long vx, long vy)
     119                 :            : {
     120                 :          7 :   pari_sp av = avma;
     121                 :            :   struct meqn meqn;
     122         [ +  - ]:          7 :   if (vx<0) vx=0;
     123         [ +  - ]:          7 :   if (vy<0) vy=1;
     124         [ -  + ]:          7 :   if (varncmp(vx,vy)>=0)
     125                 :          0 :     pari_err_PRIORITY("ellmodulareqn", pol_x(vx), ">=", vy);
     126 [ +  - ][ -  + ]:          7 :   if (ell < 0 || !uisprime(ell))
     127                 :          0 :     pari_err_PRIME("ellmodulareqn (level)", stoi(ell));
     128                 :            : 
     129         [ -  + ]:          7 :   if (!get_modular_eqn(&meqn, ell, vx, vy))
     130                 :          0 :     err_modular_eqn(ell);
     131                 :          7 :   return gerepilecopy(av,mkvec2(meqn.eq, stoi(meqn.type=='A')));
     132                 :            : }
     133                 :            : 
     134                 :            : static GEN
     135                 :       6705 : Fq_elldivpol2(GEN a4, GEN a6, GEN T, GEN p)
     136                 :            : {
     137                 :       6705 :   return mkpoln(4, utoi(4), gen_0, Fq_mulu(a4, 4, T, p), Fq_mulu(a6, 4, T, p));
     138                 :            : }
     139                 :            : 
     140                 :            : static GEN
     141                 :       6705 : Fq_elldivpol2d(GEN a4, GEN T, GEN p)
     142                 :            : {
     143                 :       6705 :   return mkpoln(3, utoi(6), gen_0, Fq_mulu(a4, 2, T, p));
     144                 :            : }
     145                 :            : 
     146                 :            : static GEN
     147                 :        980 : FqX_numer_isog_abscissa(GEN h, GEN a4, GEN a6, GEN T, GEN p, long vx)
     148                 :            : {
     149                 :            :   GEN mp1, dh, ddh, t, u, t1, t2, t3, t4, f0;
     150                 :        980 :   long m = degpol(h);
     151                 :        980 :   mp1 = gel(h, m + 1); /* negative of first power sum */
     152                 :        980 :   dh = FqX_deriv(h, T, p);
     153                 :        980 :   ddh = FqX_deriv(dh, T, p);
     154                 :        980 :   t  = Fq_elldivpol2(a4, a6, T, p);
     155                 :        980 :   u  = Fq_elldivpol2d(a4, T, p);
     156                 :        980 :   t1 = FqX_sub(FqX_sqr(dh, T, p), FqX_mul(ddh, h, T, p), T, p);
     157                 :        980 :   t2 = FqX_mul(u, FqX_mul(h, dh, T, p), T, p);
     158                 :        980 :   t3 = FqX_mul(FqX_sqr(h, T, p),
     159                 :            :                deg1pol_shallow(stoi(2*m), Fq_mulu(mp1, 2, T, p), vx), T, p);
     160                 :        980 :   f0 = FqX_add(FqX_sub(FqX_mul(t, t1, T, p), t2, T, p), t3, T, p);
     161                 :        980 :   t4 = FqX_mul(pol_x(vx),  FqX_sqr(h, T, p), T, p);
     162                 :        980 :   return FqX_add(t4, f0, T, p);
     163                 :            : }
     164                 :            : 
     165                 :            : static GEN
     166                 :       1463 : Zq_inv(GEN b, GEN T, GEN q, GEN p, long e)
     167                 :            : {
     168         [ -  + ]:       2926 :   return e==1 ? Fq_inv(b, T, p):
     169         [ +  + ]:       1463 :          typ(b)==t_INT ? Fp_inv(b, q):  ZpXQ_inv(b, T, p, e);
     170                 :            : }
     171                 :            : 
     172                 :            : static GEN
     173                 :     170119 : Zq_div(GEN a, GEN b, GEN T, GEN q, GEN p, long e)
     174                 :            : {
     175         [ +  + ]:     170119 :   if (e==1) return Fq_div(a, b, T, q);
     176                 :     170119 :   return Fq_mul(a, Zq_inv(b, T, q, p, e), T, q);
     177                 :            : }
     178                 :            : 
     179                 :            : static GEN
     180                 :          0 : Zq_sqrt(GEN b, GEN T, GEN q, GEN p, long e)
     181                 :            : {
     182         [ #  # ]:          0 :   return e==1 ? Fq_sqrt(b, T, q):
     183         [ #  # ]:          0 :          typ(b)==t_INT ? Zp_sqrt(b, p, e):  ZpXQ_sqrt(b, T, p, e);
     184                 :            : }
     185                 :            : 
     186                 :            : static GEN
     187                 :      63637 : Zq_divexact(GEN a, GEN b)
     188                 :            : {
     189         [ +  + ]:      63637 :   return typ(a)==t_INT ? diviiexact(a, b): ZX_Z_divexact(a, b);
     190                 :            : }
     191                 :            : 
     192                 :            : static long
     193                 :      63602 : Zq_pval(GEN a, GEN p)
     194                 :            : {
     195         [ +  + ]:      63602 :   return typ(a)==t_INT ? Z_pval(a, p): ZX_pval(a, p);
     196                 :            : }
     197                 :            : 
     198                 :            : static GEN
     199                 :     104304 : Zq_Z_div_safe(GEN a, GEN b, GEN T, GEN q, GEN p, long e)
     200                 :            : {
     201                 :            :   long v;
     202         [ +  + ]:     104304 :   if (e==1) return Fq_div(a, b, T, q);
     203                 :       1008 :   v = Z_pvalrem(b, p, &b);
     204         [ +  + ]:       1008 :   if (v>0)
     205                 :            :   {
     206                 :         35 :     long w = Z_pval(Q_content(a), p);
     207         [ -  + ]:         35 :     if (v>w) pari_err_INV("Zq_div",b);
     208                 :         35 :     a = Zq_divexact(a, powiu(p,v));
     209                 :            :   }
     210                 :     104304 :   return Fq_Fp_mul(a, Fp_inv(b, q), T, q);
     211                 :            : }
     212                 :            : 
     213                 :            : /*Gives the first precS terms of the Weierstrass series related to */
     214                 :            : /*E: y^2 = x^3 + a4x + a6.  Assumes (precS-2)*(2precS+3) < ULONG_MAX, i.e.
     215                 :            :  * precS < 46342 in 32-bit machines */
     216                 :            : static GEN
     217                 :      11450 : find_coeff(GEN a4, GEN a6, GEN T, GEN p, long precS, GEN pp, long e)
     218                 :            : {
     219                 :            :   GEN res, den;
     220                 :            :   long k, h;
     221         [ +  + ]:      11450 :   if (e > 1) { p = sqri(p); e *= 2; }
     222                 :      11450 :   res = cgetg(precS+1, t_VEC);
     223                 :      11450 :   den = cgetg(precS+1, t_VECSMALL);
     224         [ -  + ]:      11450 :   if (precS == 0) return res;
     225                 :      11450 :   gel(res, 1) = Fq_div(a4, stoi(-5), T, p);
     226                 :      11450 :   den[1] = 0;
     227         [ -  + ]:      11450 :   if (precS == 1) return res;
     228                 :      11450 :   gel(res, 2) = Fq_div(a6, stoi(-7), T, p);
     229                 :      11450 :   den[2] = 0;
     230         [ +  + ]:     115754 :   for (k = 3; k <= precS; ++k)
     231                 :            :   {
     232                 :     104304 :     pari_sp btop = avma;
     233                 :     104304 :     GEN a = gen_0, d;
     234                 :     104304 :     long v=0;
     235         [ +  + ]:     104304 :     if (e > 1)
     236         [ +  + ]:       9422 :       for (h = 1; h <= k-2; h++)
     237                 :       8414 :         v = maxss(v, den[h]+den[k-1-h]);
     238         [ +  + ]:     935176 :     for (h = 1; h <= k-2; h++)
     239                 :            :     {
     240                 :     830872 :       GEN b = Fq_mul(gel(res, h), gel(res, k-1-h), T, p);
     241         [ +  + ]:     830872 :       if (v)
     242                 :       1638 :         b = Fq_Fp_mul(b, powiu(pp, v-(den[h]+den[k-1-h])), T, p);
     243                 :     830872 :       a = Fq_add(a, b, T, p);
     244                 :            :     }
     245                 :     104304 :     v += Z_pvalrem(utoi((k-2) * (2*k + 3)), pp, &d);
     246                 :     104304 :     a = Zq_div(gmulgs(a, 3), d, T, p, pp, e);
     247                 :     104304 :     gel(res, k) = gerepileupto(btop, a);
     248                 :     104304 :     den[k] = v;
     249                 :            :   }
     250                 :      11450 :   return mkvec2(res, den);
     251                 :            : }
     252                 :            : 
     253                 :            : /****************************************************************************/
     254                 :            : /*               SIMPLE ELLIPTIC CURVE OVER Fq                              */
     255                 :            : /****************************************************************************/
     256                 :            : 
     257                 :            : static GEN
     258                 :       1582 : Fq_ellj(GEN a4, GEN a6, GEN T, GEN p)
     259                 :            : {
     260                 :       1582 :   pari_sp ltop=avma;
     261                 :       1582 :   GEN a43 = Fq_mulu(Fq_powu(a4, 3, T, p), 4, T, p);
     262                 :       1582 :   GEN j   = Fq_div(Fq_mulu(a43, 1728, T, p),
     263                 :            :                    Fq_add(a43, Fq_mulu(Fq_sqr(a6, T, p), 27, T, p), T, p), T, p);
     264                 :       1582 :   return gerepileupto(ltop, j);
     265                 :            : }
     266                 :            : 
     267                 :            : static GEN
     268                 :         35 : Zq_ellj(GEN a4, GEN a6, GEN T, GEN p, GEN pp, long e)
     269                 :            : {
     270                 :         35 :   pari_sp ltop=avma;
     271                 :         35 :   GEN a43 = Fq_mulu(Fq_powu(a4, 3, T, p), 4, T, p);
     272                 :         35 :   GEN j   = Zq_div(Fq_mulu(a43, 1728, T, p),
     273                 :            :                    Fq_add(a43, Fq_mulu(Fq_sqr(a6, T, p), 27, T, p), T, p), T, p, pp, e);
     274                 :         35 :   return gerepileupto(ltop, j);
     275                 :            : }
     276                 :            : /****************************************************************************/
     277                 :            : /*                              EIGENVALUE                                  */
     278                 :            : /****************************************************************************/
     279                 :            : 
     280                 :            : struct eigen_ellinit
     281                 :            : {
     282                 :            :   GEN a4, h, T, p;
     283                 :            :   GEN RHS, DRHS, X12, Gr, nGr,O;
     284                 :            :   ulong pp;
     285                 :            : };
     286                 :            : 
     287                 :            : static void
     288                 :       3972 : init_eigen(struct eigen_ellinit *Edat, GEN a4, GEN a6, GEN h, GEN T, GEN p)
     289                 :            : {
     290                 :       3972 :   pari_sp ltop = avma;
     291                 :       3972 :   GEN RHS  = FqX_rem(mkpoln(4, gen_1, gen_0, a4, a6), h, T, p);
     292                 :       3972 :   GEN DRHS = FqX_rem(mkpoln(3, utoi(3), gen_0, a4), h, T, p);
     293                 :       3972 :   GEN lambda = FqXQ_div(DRHS, FqX_mulu(RHS, 4, T, p), h, T, p);
     294                 :       3972 :   GEN C = FqX_sub(FqXQ_mul(lambda, DRHS, h, T, p), monomial(gen_2,1,0), T, p);
     295                 :       3972 :   GEN D = FqXQ_mul(FqX_mulu(lambda, 2, T, p),FqX_sub(pol_x(0), C, T, p), h, T, p);
     296                 :       3972 :   GEN X12 = mkvec2(C, FqX_Fq_add(D, gen_m1, T, p));
     297         [ +  + ]:       7871 :   GEN Gr = T ? FpXQXQ_halfFrobenius(RHS, h, T, p):
     298                 :       3899 :                FpXQ_pow(RHS, shifti(p, -1), h, p);
     299                 :       3972 :   GEN nGr = FqX_neg(Gr, T, p);
     300                 :       3972 :   gerepileall(ltop, 5, &RHS, &DRHS, &X12, &Gr, &nGr);
     301                 :       3972 :   Edat->a4    = gcopy(a4);
     302                 :       3972 :   Edat->h     = gcopy(h);
     303                 :       3972 :   Edat->T     = T;
     304                 :       3972 :   Edat->p     = p;
     305                 :       3972 :   Edat->pp    = 0;
     306                 :       3972 :   Edat->RHS   = RHS;
     307                 :       3972 :   Edat->DRHS  = DRHS;
     308                 :       3972 :   Edat->X12   = X12;
     309                 :       3972 :   Edat->Gr    = Gr;
     310                 :       3972 :   Edat->nGr   = nGr;
     311                 :       3972 :   Edat->O     = mkvec2(pol_x(0), pol_1(0));
     312                 :       3972 : }
     313                 :            : 
     314                 :            : static void
     315                 :        228 : init_eigenu(struct eigen_ellinit *Edat, GEN a4, GEN a6, GEN h, GEN T, ulong p)
     316                 :            : {
     317                 :        228 :   pari_sp ltop = avma;
     318                 :        228 :   long vT = get_Flx_var(T);
     319                 :        228 :   GEN g1 = pol1_Flx(vT), g0 = pol0_Flx(vT);
     320                 :        228 :   GEN RHS  = FlxqX_rem(mkpoln(4, g1, g0, a4, a6), h, T, p);
     321                 :        228 :   GEN DRHS = FlxqX_rem(mkpoln(3, Fl_to_Flx(3, T[1]), g0, a4), h, T, p);
     322                 :        228 :   GEN lambda = FlxqXQ_div(DRHS, FlxX_Fl_mul(RHS, 4, p), h, T, p);
     323                 :        228 :   GEN C = FlxX_sub(FlxqXQ_mul(lambda, DRHS, h, T, p), monomial(Fl_to_Flx(2,vT),1,0), p);
     324                 :        228 :   GEN D = FlxqXQ_mul(FlxX_double(lambda, p),FlxX_sub(pol_x(0), C, p), h, T, p);
     325                 :        228 :   GEN X12 = mkvec2(C, FlxX_Flx_add(D, Fl_to_Flx(p-1,vT), p));
     326                 :        228 :   GEN Gr = FlxqXQ_halfFrobenius(RHS,h,T,p);
     327                 :        228 :   GEN nGr = FlxX_neg(Gr, p);
     328                 :        228 :   GEN O = mkvec2(monomial(g1,1,0), monomial(g1,0,0));
     329                 :        228 :   gerepileall(ltop, 6, &RHS, &DRHS, &X12, &Gr, &nGr, &O);
     330                 :        228 :   Edat->a4    = gcopy(a4);
     331                 :        228 :   Edat->h     = gcopy(h);
     332                 :        228 :   Edat->T     = T;
     333                 :        228 :   Edat->p     = NULL;
     334                 :        228 :   Edat->pp    = p;
     335                 :        228 :   Edat->RHS   = RHS;
     336                 :        228 :   Edat->DRHS  = DRHS;
     337                 :        228 :   Edat->X12   = X12;
     338                 :        228 :   Edat->Gr    = Gr;
     339                 :        228 :   Edat->nGr   = nGr;
     340                 :        228 :   Edat->O     = O;
     341                 :        228 : }
     342                 :            : static GEN
     343                 :       5232 : eigen_elldbl(void *E, GEN P)
     344                 :            : {
     345                 :       5232 :   pari_sp ltop = avma;
     346                 :       5232 :   struct eigen_ellinit *Edat=(struct eigen_ellinit *)E;
     347                 :       5232 :   GEN T = Edat->T, p = Edat->p, h = Edat->h, x, y;
     348         [ -  + ]:       5232 :   if (ell_is_inf(P)) return gcopy(P);
     349                 :       5232 :   x = gel(P,1), y = gel(P,2);
     350 [ +  + ][ +  - ]:       5232 :   if (FqX_equal(x, pol_x(0)) && FqX_equal(y, pol_1(0)))
     351                 :       3713 :     return Edat->X12;
     352                 :            :   else
     353                 :            :   {
     354                 :       1519 :     GEN t1 = FqX_Fq_add(FqX_mulu(FqXQ_sqr(x,h,T,p),3,T, p), Edat->a4, T, p);
     355                 :       1519 :     GEN t2 = FqXQ_mul(FqX_mulu(y, 2, T, p), Edat->RHS, h, T, p);
     356                 :       1519 :     GEN lambda = FqXQ_div(t1, t2, h, T, p);
     357                 :       1519 :     GEN C = FqX_sub(FqXQ_mul(FqXQ_sqr(lambda, h, T, p), Edat->RHS, h, T, p),
     358                 :            :                     FqX_mulu(x, 2, T, p), T, p);
     359                 :       1519 :     GEN D = FqX_sub(FqXQ_mul(lambda, FqX_sub(x, C, T, p), h, T, p), y, T, p);
     360                 :       5232 :     return gerepilecopy(ltop, mkvec2(C,D));
     361                 :            :   }
     362                 :            : }
     363                 :            : 
     364                 :            : /* Returns the addition of [P[1], P[2]*Y] and of [Q[1], Q[2]*Y]
     365                 :            :  * Computations are done modulo Y^2 - (X^3 + a4X + a6)
     366                 :            :  * An inversion is equivalent to 4M, so that this function requires about 7M
     367                 :            :  * which is the same as with the method using ell-division polynomials
     368                 :            :  * Working in mixed projective coordinates would require 11M */
     369                 :            : static GEN
     370                 :      13284 : eigen_elladd(void *E, GEN P, GEN Q)
     371                 :            : {
     372                 :      13284 :   pari_sp ltop = avma;
     373                 :      13284 :   struct eigen_ellinit *Edat=(struct eigen_ellinit *)E;
     374                 :            :   GEN Px, Py, Qx, Qy;
     375                 :      13284 :   GEN T = Edat->T, p = Edat->p, h = Edat->h, lambda, C, D;
     376         [ -  + ]:      13284 :   if (ell_is_inf(P)) return gcopy(Q);
     377         [ -  + ]:      13284 :   if (ell_is_inf(Q)) return gcopy(P);
     378                 :      13284 :   Px = gel(P,1); Py = gel(P,2);
     379                 :      13284 :   Qx = gel(Q,1); Qy = gel(Q,2);
     380         [ +  + ]:      13284 :   if (FqX_equal(Px, Qx))
     381                 :            :   {
     382         [ +  - ]:       2105 :     if (FqX_equal(Py, Qy))
     383                 :       2105 :       return eigen_elldbl(E, P);
     384                 :            :     else
     385                 :          0 :       return ellinf();
     386                 :            :   }
     387                 :      11179 :   lambda = FqXQ_div(FqX_sub(Py, Qy, T, p), FqX_sub(Px, Qx, T, p), h, T, p);
     388                 :      11179 :   C = FqX_sub(FqX_sub(FqXQ_mul(FqXQ_sqr(lambda, h, T, p), Edat->RHS, h, T, p), Px, T, p), Qx, T, p);
     389                 :      11179 :   D = FqX_sub(FqXQ_mul(lambda, FqX_sub(Px, C, T, p), h, T, p), Py, T, p);
     390                 :      13284 :   return gerepilecopy(ltop, mkvec2(C,D));
     391                 :            : }
     392                 :            : 
     393                 :            : static GEN
     394                 :        291 : eigenu_elldbl(void *E, GEN P)
     395                 :            : {
     396                 :        291 :   pari_sp ltop = avma;
     397                 :        291 :   struct eigen_ellinit *Edat=(struct eigen_ellinit *)E;
     398                 :        291 :   GEN T = Edat->T, h = Edat->h, x, y;
     399                 :        291 :   long vT = get_Flx_var(T);
     400                 :        291 :   ulong p = Edat->pp;
     401         [ -  + ]:        291 :   if (ell_is_inf(P)) return gcopy(P);
     402                 :        291 :   x = gel(P,1), y = gel(P,2);
     403 [ +  + ][ +  - ]:        291 :   if (FlxX_equal(x, monomial(pol1_Flx(vT),1,0)) && FlxX_equal(y, monomial(pol1_Flx(vT),0,0)))
     404                 :        193 :     return Edat->X12;
     405                 :            :   else
     406                 :            :   {
     407                 :         98 :     GEN t1 = FlxX_Flx_add(FlxX_triple(FlxqXQ_sqr(x,h,T,p),p), Edat->a4, p);
     408                 :         98 :     GEN t2 = FlxqXQ_mul(FlxX_double(y, p), Edat->RHS, h, T, p);
     409                 :         98 :     GEN lambda = FlxqXQ_div(t1, t2, h, T, p);
     410                 :         98 :     GEN C = FlxX_sub(FlxqXQ_mul(FlxqXQ_sqr(lambda, h, T, p), Edat->RHS, h, T, p),
     411                 :            :                      FlxX_double(x, p), p);
     412                 :         98 :     GEN D = FlxX_sub(FlxqXQ_mul(lambda, FlxX_sub(x, C, p), h, T, p), y, p);
     413                 :        291 :     return gerepilecopy(ltop, mkvec2(C,D));
     414                 :            :   }
     415                 :            : }
     416                 :            : 
     417                 :            : /* Returns the addition of [P[1], P[2]*Y] and of [Q[1], Q[2]*Y]
     418                 :            :  * Computations are done modulo Y^2 - (X^3 + a4X + a6)
     419                 :            :  * An inversion is equivalent to 4M, so that this function requires about 7M
     420                 :            :  * which is the same as with the method using ell-division polynomials
     421                 :            :  * Working in mixed projective coordinates would require 11M */
     422                 :            : static GEN
     423                 :        996 : eigenu_elladd(void *E, GEN P, GEN Q)
     424                 :            : {
     425                 :        996 :   pari_sp ltop = avma;
     426                 :        996 :   struct eigen_ellinit *Edat=(struct eigen_ellinit *)E;
     427                 :            :   GEN Px, Py, Qx, Qy;
     428                 :        996 :   GEN T = Edat->T, h = Edat->h, lambda, C, D;
     429                 :        996 :   ulong p = Edat->pp;
     430         [ -  + ]:        996 :   if (ell_is_inf(P)) return gcopy(Q);
     431         [ -  + ]:        996 :   if (ell_is_inf(Q)) return gcopy(P);
     432                 :        996 :   Px = gel(P,1); Py = gel(P,2);
     433                 :        996 :   Qx = gel(Q,1); Qy = gel(Q,2);
     434         [ +  + ]:        996 :   if (FlxX_equal(Px, Qx))
     435                 :            :   {
     436         [ +  - ]:        121 :     if (FlxX_equal(Py, Qy))
     437                 :        121 :       return eigenu_elldbl(E, P);
     438                 :            :     else
     439                 :          0 :       return ellinf();
     440                 :            :   }
     441                 :        875 :   lambda = FlxqXQ_div(FlxX_sub(Py, Qy, p), FlxX_sub(Px, Qx, p), h, T, p);
     442                 :        875 :   C = FlxX_sub(FlxX_sub(FlxqXQ_mul(FlxqXQ_sqr(lambda, h, T, p), Edat->RHS, h, T, p), Px, p), Qx, p);
     443                 :        875 :   D = FlxX_sub(FlxqXQ_mul(lambda, FlxX_sub(Px, C, p), h, T, p), Py, p);
     444                 :        996 :   return gerepilecopy(ltop, mkvec2(C,D));
     445                 :            : }
     446                 :            : 
     447                 :            : static GEN
     448                 :       2289 : eigen_ellmulu(struct eigen_ellinit *E, GEN z, ulong n)
     449                 :            : {
     450                 :       2289 :   pari_sp av = avma;
     451 [ +  - ][ -  + ]:       2289 :   if (!n || ell_is_inf(z)) return mkvec(gen_0);
     452         [ +  + ]:       2289 :   if (n == 1) return gcopy(z);
     453         [ +  + ]:       1680 :   if (E->pp)
     454                 :         72 :     return gerepileupto(av, gen_powu(z, n, E, &eigenu_elldbl, &eigenu_elladd));
     455                 :            :   else
     456                 :       2289 :     return gerepileupto(av, gen_powu(z, n, E, &eigen_elldbl, &eigen_elladd));
     457                 :            : }
     458                 :            : 
     459                 :            : static GEN
     460                 :        498 : Fq_to_Flx(GEN a4, GEN T, ulong p)
     461                 :            : {
     462         [ -  + ]:        498 :   return typ(a4)==t_INT ? Z_to_Flx(a4, p, get_FpX_var(T)): ZX_to_Flx(a4, p);
     463                 :            : }
     464                 :            : 
     465                 :            : /*Finds the eigenvalue of the Frobenius given E, ell odd prime, h factor of the
     466                 :            :  *ell-division polynomial, p and tr the possible values for the trace
     467                 :            :  *(useful for primes with one root)*/
     468                 :            : static ulong
     469                 :       3227 : find_eigen_value(GEN a4, GEN a6, ulong ell, GEN h, GEN T, GEN p, GEN tr)
     470                 :            : {
     471                 :       3227 :   pari_sp ltop = avma;
     472                 :            :   GEN BP, Dr;
     473                 :            :   ulong t;
     474                 :            :   struct eigen_ellinit Edat;
     475         [ +  + ]:       3227 :   ulong pp = T ?itou_or_0(p): 0;
     476         [ +  + ]:       3227 :   if (pp)
     477                 :        215 :     init_eigenu(&Edat, Fq_to_Flx(a4, T, pp), Fq_to_Flx(a6, T, pp),
     478                 :            :                        ZXX_to_FlxX(h,pp, get_FpX_var(T)), ZXT_to_FlxT(T,pp), pp);
     479                 :            :   else
     480                 :       3012 :     init_eigen(&Edat, a4, a6, h, T, p);
     481                 :       3227 :   Dr = BP = Edat.O;
     482                 :            :   /*[0,Gr], BP, Dr are not points on the curve. */
     483                 :            :   /*To obtain the corresponding points, multiply the y-coordinates by Y */
     484         [ +  + ]:       3227 :   if (!tr)
     485                 :            :   {
     486                 :       2884 :     pari_sp btop = avma;
     487         [ +  - ]:      14364 :     for (t = 1; t <= (ell>>1); t++)
     488                 :            :     {
     489         [ +  + ]:      14364 :       if (gequal(gel(Dr,2), Edat.Gr))  { avma = ltop; return t; }
     490         [ +  + ]:      12376 :       if (gequal(gel(Dr,2), Edat.nGr)) { avma = ltop; return ell-t; }
     491         [ +  + ]:      11480 :       Dr = pp ? eigenu_elladd(&Edat, Dr, BP): eigen_elladd(&Edat, Dr, BP);
     492                 :      11480 :       Dr = gerepileupto(btop, Dr);
     493                 :            :     }
     494                 :            :   }
     495                 :            :   else
     496                 :            :   {
     497                 :        343 :     t = Fl_div(tr[1], 2, ell);
     498         [ +  + ]:        343 :     if (t < (ell>>1)) t = ell - t;
     499                 :        343 :     Dr = eigen_ellmulu(&Edat, BP, t);
     500         [ +  + ]:        343 :     if (gequal(gel(Dr,2), Edat.Gr)) { avma = ltop; return t; }
     501         [ +  - ]:        238 :     if (gequal(gel(Dr,2), Edat.nGr)) { avma = ltop; return ell - t; }
     502                 :            :   }
     503                 :       3227 :   pari_err_BUG("find_eigen_value"); return 0; /* NOT REACHED */
     504                 :            : }
     505                 :            : 
     506                 :            : /*Finds the eigenvalue of the Frobenius modulo ell^k given E, ell, k, h factor
     507                 :            :  *of the ell-division polynomial, lambda the previous eigen value and p */
     508                 :            : static ulong
     509                 :        973 : find_eigen_value_power(GEN a4, GEN a6, ulong ell, long k, GEN h, ulong lambda, GEN T, GEN p)
     510                 :            : {
     511                 :        973 :   pari_sp ltop = avma;
     512                 :            :   pari_sp btop;
     513                 :            :   struct eigen_ellinit Edat;
     514                 :            :   GEN BP, Dr, Gr, nGr;
     515                 :            :   /*[0,Gr], BP, Dr are not points on the curve. */
     516                 :            :   /*To obtain the corresponding points, multiply the y-coordinates by Y */
     517                 :        973 :   ulong t, ellk1 = upowuu(ell, k-1), ellk = ell*ellk1;
     518         [ +  + ]:        973 :   ulong pp = T ?itou_or_0(p): 0;
     519         [ +  + ]:        973 :   if (pp)
     520                 :         13 :     init_eigenu(&Edat, Fq_to_Flx(a4, T, pp), Fq_to_Flx(a6, T, pp),
     521                 :            :         ZXX_to_FlxX(h, pp, get_FpX_var(T)), ZXT_to_FlxT(T,pp), pp);
     522                 :            :   else
     523                 :        960 :     init_eigen(&Edat, a4, a6, h, T, p);
     524                 :        973 :   BP = eigen_ellmulu(&Edat, Edat.O, ellk1);
     525                 :        973 :   Dr = eigen_ellmulu(&Edat, Edat.O, lambda);
     526                 :        973 :   Gr = Edat.Gr; nGr = Edat.nGr;
     527                 :            : 
     528                 :        973 :   btop = avma;
     529         [ +  - ]:       1904 :   for (t = 0; t < ellk; t += ellk1)
     530                 :            :   {
     531         [ +  + ]:       1904 :     if (gequal(gel(Dr,2), Gr))  { avma = ltop; return t+lambda; }
     532         [ -  + ]:        931 :     if (gequal(gel(Dr,2), nGr)) { avma = ltop; return ellk-(t+lambda); }
     533         [ +  + ]:        931 :     Dr = pp ? eigenu_elladd(&Edat, Dr, BP): eigen_elladd(&Edat, Dr, BP);
     534         [ -  + ]:        931 :     if (gc_needed(btop, 1))
     535                 :          0 :       Dr = gerepileupto(btop, Dr);
     536                 :            :   }
     537                 :          0 :   pari_err_BUG("find_eigen_value_power");
     538                 :        973 :   return 0; /* NOT REACHED */
     539                 :            : }
     540                 :            : 
     541                 :            : /*Finds the kernel polynomial h, dividing the ell-division polynomial from the
     542                 :            :   isogenous curve Eb and trace term pp1. Uses CCR algorithm and returns h.
     543                 :            :   Return NULL if E and Eb are *not* isogenous. */
     544                 :            : static GEN
     545                 :       5725 : find_kernel(GEN a4, GEN a6, ulong ell, GEN a4t, GEN a6t, GEN pp1, GEN T, GEN p, GEN pp, long e)
     546                 :            : {
     547                 :       5725 :   const long ext = 2;
     548                 :       5725 :   pari_sp ltop = avma, btop;
     549                 :            :   GEN P, v, tlist, h;
     550                 :            :   long i, j, k;
     551                 :       5725 :   long deg = (ell - 1)/2, dim = 2 + deg + ext;
     552                 :       5725 :   GEN psi2 = Fq_elldivpol2(a4, a6, T, p);
     553                 :       5725 :   GEN Dpsi2 = Fq_elldivpol2d(a4, T, p);
     554                 :       5725 :   GEN C  = find_coeff(a4, a6, T, p, dim, pp, e);
     555                 :       5725 :   GEN Ct = find_coeff(a4t, a6t, T, p, dim, pp, e);
     556                 :       5725 :   GEN V = cgetg(dim+1, t_VEC);
     557         [ +  + ]:      69327 :   for (k = 1; k <= dim; k++)
     558                 :            :   {
     559                 :      63602 :     long v = mael(C,2,k);
     560                 :      63602 :     GEN z = gmul(gsub(gmael(Ct,1,k), gmael(C,1,k)), shifti(mpfact(2*k), -1));
     561 [ +  - ][ -  + ]:      63602 :     if (signe(z) && Zq_pval(z, pp) < v) return NULL;
     562                 :      63602 :     gel(V, k) = Zq_divexact(z, powiu(pp, v));
     563                 :            :   }
     564                 :       5725 :   btop = avma;
     565                 :       5725 :   v = zerovec(dim);
     566                 :       5725 :   gel(v, 1) = utoi(deg);
     567                 :       5725 :   gel(v, 2) = pp1;
     568                 :       5725 :   P = pol_x(0);
     569         [ +  + ]:      57877 :   for (k = 3; k <= dim; k++)
     570                 :            :   {
     571                 :      52152 :     GEN s, r = FqX_Fq_mul(Dpsi2, gel(P, 3), T, p);
     572         [ +  + ]:     415436 :     for (j = 4; j < lg(P); j++)
     573                 :            :     {
     574                 :     363284 :       long o = j - 2;
     575                 :     363284 :       GEN D = FqX_add(RgX_shift_shallow(Dpsi2, 1), FqX_mulu(psi2, o-1, T, p), T, p);
     576                 :     363284 :       GEN E = FqX_Fq_mul(D, Fq_mulu(gel(P, j), o, T, p), T, p);
     577                 :     363284 :       r = FqX_add(r, RgX_shift_shallow(E, o-2), T, p);
     578                 :            :     }
     579                 :      52152 :     P = r;
     580                 :      52152 :     s = Fq_mul(gel(P, 2), gel(v, 1), T, p);
     581         [ +  + ]:     467588 :     for (j = 3; j < lg(P)-1; j++)
     582                 :     415436 :       s = Fq_add(s, Fq_mul(gel(P, j), gel(v, j-1), T, p), T, p);
     583                 :      52152 :     gel(v, k) = Zq_Z_div_safe(Fq_sub(gel(V, k-2), s, T, p), gel(P, j), T, p, pp, e);
     584         [ -  + ]:      52152 :     if (gc_needed(btop, 1))
     585                 :            :     {
     586         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"find_kernel");
     587                 :          0 :       gerepileall(btop, 2, &v, &P);
     588                 :            :     }
     589                 :            :   }
     590                 :       5725 :   tlist = cgetg(dim, t_VEC);
     591                 :       5725 :   gel(tlist, dim-1) = gen_1;
     592         [ +  + ]:      57877 :   for (k = 1; k <= dim-2; k++)
     593                 :            :   {
     594                 :      52152 :     pari_sp btop = avma;
     595                 :      52152 :     GEN s = gel(v, k+1);
     596         [ +  + ]:     415436 :     for (i = 1; i < k; i++)
     597                 :     363284 :       s = Fq_add(s, Fq_mul(gel(tlist, dim-i-1), gel(v, k-i+1), T, p), T, p);
     598                 :      52152 :     gel(tlist, dim-k-1) = gerepileupto(btop, Zq_Z_div_safe(s, stoi(-k), T, p, pp, e));
     599                 :            :   }
     600         [ +  + ]:      15119 :   for (i = 1; i <= ext; i++)
     601         [ +  + ]:      10422 :     if (signe(Fq_red(gel(tlist, i),T, pp))) { avma = ltop; return NULL; }
     602                 :       4697 :   h = FqX_red(RgV_to_RgX(vecslice(tlist, ext+1, dim-1), 0),T,p);
     603         [ +  - ]:       5725 :   return signe(Fq_elldivpolmod(a4, a6, ell, h, T, pp)) ? NULL: h;
     604                 :            : }
     605                 :            : 
     606                 :            : static GEN
     607                 :       4318 : 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)
     608                 :            : {
     609                 :       4318 :   pari_sp ltop = avma;
     610                 :       4318 :   GEN dxxgj = FqX_eval(Dxxg, j, T, p);
     611                 :       4318 :   GEN dxJgj = FqX_eval(DxJg, j, T, p);
     612                 :       4318 :   GEN dJJgj = FqX_eval(DJJg, j, T, p);
     613                 :       4318 :   GEN E42 = Fq_sqr(E4, T, p), E6ovE4 = Zq_div(E6, E4, T, p, pp, e);
     614                 :       4318 :   GEN a = Fq_mul(gprime, dxxgj, T, p);
     615                 :       4318 :   GEN b = Fq_mul(Fq_mul(Fq_mulu(j,2*q, T, p), dxJgj, T, p), E6ovE4, T, p);
     616                 :       4318 :   GEN c = Fq_mul(Zq_div(Fq_sqr(E6ovE4, T, p), gprime, T, p, pp, e), j, T, p);
     617                 :       4318 :   GEN d = Fq_mul(Fq_mul(c,sqru(q), T, p), Fq_add(pJ, Fq_mul(j, dJJgj, T, p), T, p), T, p);
     618                 :       4318 :   GEN f = Fq_sub(Fq_div(E6ovE4,utoi(3), T, p),
     619                 :            :                  Zq_div(E42, Fq_mulu(E6,2,T, p), T, p, pp, e), T, p);
     620                 :       4318 :   GEN g = Fq_sub(Fq_sub(b,a,T,p), d, T, p);
     621                 :       4318 :   return gerepileupto(ltop, Fq_add(Zq_div(g,px,T,p,pp,e), Fq_mulu(f,q,T,p), T, p));
     622                 :            : }
     623                 :            : 
     624                 :            : /* Finds the isogenous EC, and the sum of the x-coordinates of the points in
     625                 :            :  * the kernel of the isogeny E -> Eb
     626                 :            :  * E: elliptic curve, ell: a prime, meqn: Atkin modular equation
     627                 :            :  * g: root of meqn defining isogenous curve Eb. */
     628                 :            : static GEN
     629                 :       1680 : find_isogenous_from_Atkin(GEN a4, GEN a6, ulong ell, struct meqn *MEQN, GEN g, GEN T, GEN pp, long e)
     630                 :            : {
     631                 :       1680 :   pari_sp ltop = avma, btop;
     632                 :       1680 :   GEN meqn = MEQN->eq, meqnx, Roots, gprime, u1;
     633                 :       1680 :   long k, vJ = MEQN->vy;
     634         [ +  + ]:       1680 :   GEN p = e==1 ? pp: powiu(pp, e);
     635                 :       1680 :   GEN E4 = Fq_div(a4, stoi(-3), T, p);
     636                 :       1680 :   GEN E6 = Fq_neg(Fq_halve(a6, T, p), T, p);
     637                 :       1680 :   GEN E42 = Fq_sqr(E4, T, p);
     638                 :       1680 :   GEN E43 = Fq_mul(E4, E42, T, p);
     639                 :       1680 :   GEN E62 = Fq_sqr(E6, T, p);
     640                 :       1680 :   GEN delta = Fq_div(Fq_sub(E43, E62, T, p), utoi(1728), T, p);
     641                 :       1680 :   GEN j = Zq_div(E43, delta, T, p, pp, e);
     642                 :       1680 :   GEN Dx = RgX_deriv(meqn);
     643                 :       1680 :   GEN DJ = deriv(meqn, vJ);
     644                 :       1680 :   GEN Dxg = FpXY_Fq_evaly(Dx, g, T, p, vJ);
     645                 :       1680 :   GEN px = FqX_eval(Dxg, j, T, p), dx = Fq_mul(px, g, T, p);
     646                 :       1680 :   GEN DJg = FpXY_Fq_evaly(DJ, g, T, p, vJ);
     647                 :       1680 :   GEN pJ = FqX_eval(DJg, j, T, p), dJ = Fq_mul(pJ, j, T, p);
     648                 :       1680 :   GEN Dxx = RgX_deriv(Dx);
     649                 :       1680 :   GEN DxJg = FqX_deriv(Dxg, T, p);
     650                 :            : 
     651                 :       1680 :   GEN Dxxg = FpXY_Fq_evaly(Dxx, g, T, p, vJ);
     652                 :       1680 :   GEN DJJg = FqX_deriv(DJg, T, p);
     653                 :            :   GEN a, b;
     654 [ +  + ][ +  + ]:       1680 :   if (!signe(dJ) || !signe(dx))
     655                 :            :   {
     656 [ -  + ][ #  # ]:         28 :     if (DEBUGLEVEL>0) err_printf("[A: d%c=0]",signe(dJ)? 'x': 'J');
     657                 :         28 :     avma = ltop; return NULL;
     658                 :            :   }
     659                 :       1652 :   a = Fq_mul(dJ, Fq_mul(g, E6, T, p), T, p);
     660                 :       1652 :   b = Fq_mul(E4, dx, T, p);
     661                 :       1652 :   gprime = Zq_div(a, b, T, p, pp, e);
     662                 :            : 
     663                 :       1652 :   u1 = compute_u(gprime, Dxxg, DxJg, DJJg, j, pJ, px, 1, E4, E6, T, p, pp, e);
     664                 :       1652 :   meqnx = FpXY_Fq_evaly(meqn, g, T, p, vJ);
     665                 :       1652 :   Roots = FqX_roots(meqnx, T, pp);
     666                 :            : 
     667                 :       1652 :   btop = avma;
     668         [ +  - ]:       2680 :   for (k = lg(Roots)-1; k >= 1; k--, avma = btop)
     669                 :            :   {
     670         [ +  + ]:       2680 :     GEN jt = e==1 ? gel(Roots, k): ZpXQX_liftroot(meqnx, gel(Roots, k), T, pp, e);
     671 [ +  + ][ +  + ]:       2680 :     if (signe(jt) == 0 || signe(Fq_sub(jt, utoi(1728), T, p)) == 0)
     672                 :            :     {
     673 [ -  + ][ #  # ]:         14 :       if (DEBUGLEVEL>0) err_printf("[A: jt=%ld]",signe(jt)? 1728: 0);
     674                 :         14 :       avma = ltop; return NULL;
     675                 :            :     }
     676                 :            :     else
     677                 :            :     {
     678                 :       2666 :       GEN pxstar = FqX_eval(Dxg, jt, T, p);
     679                 :       2666 :       GEN dxstar = Fq_mul(pxstar, g, T, p);
     680                 :       2666 :       GEN pJstar = FqX_eval(DJg, jt, T, p);
     681                 :       2666 :       GEN dJstar = Fq_mul(Fq_mulu(jt, ell, T, p), pJstar, T, p);
     682                 :       2666 :       GEN u = Fq_mul(Fq_mul(dxstar, dJ, T, p), E6, T, p);
     683                 :       2666 :       GEN v = Fq_mul(Fq_mul(dJstar, dx, T, p), E4, T, p);
     684                 :       2666 :       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);
     685                 :       2666 :       GEN E6t = Zq_div(Fq_mul(u, E4t, T, p), v, T, p, pp, e);
     686                 :       2666 :       GEN u2 = compute_u(gprime, Dxxg, DxJg, DJJg, jt, pJstar, pxstar, ell, E4t, E6t, T, p, pp, e);
     687                 :       2666 :       GEN pp1 = Fq_mulu(Fq_sub(u1, u2, T, p), 3*ell, T, p);
     688                 :       2666 :       GEN a4t = Fq_mul(mulsi(-3, powuu(ell,4)), E4t, T, p);
     689                 :       2666 :       GEN a6t = Fq_mul(mulsi(-2, powuu(ell,6)), E6t, T, p);
     690                 :       2666 :       GEN h = find_kernel(a4, a6, ell, a4t, a6t, pp1, T, p, pp, e);
     691         [ +  + ]:       2666 :       if (h) return gerepilecopy(ltop, mkvec3(a4t, a6t, h));
     692                 :            :     }
     693                 :            :   }
     694                 :          0 :   pari_err_BUG("find_isogenous_from_Atkin, kernel not found");
     695                 :       1680 :   return NULL;
     696                 :            : }
     697                 :            : 
     698                 :            : /* Finds E' ell-isogenous to E and the trace term p1 from canonical modular
     699                 :            :  *   equation meqn
     700                 :            :  * E: elliptic curve, ell: a prime, meqn: canonical modular equation
     701                 :            :  * g: root of meqn defining isogenous curve Eb. */
     702                 :            : static GEN
     703                 :       3073 : find_isogenous_from_canonical(GEN a4, GEN a6, ulong ell, struct meqn *MEQN, GEN g, GEN T, GEN pp, long e)
     704                 :            : {
     705                 :       3073 :   pari_sp ltop = avma;
     706                 :       3073 :   GEN meqn = MEQN->eq;
     707                 :       3073 :   long vJ = MEQN->vy;
     708         [ -  + ]:       3073 :   GEN p = e==1 ? pp: powiu(pp, e);
     709                 :            :   GEN h;
     710                 :       3073 :   GEN E4 = Fq_div(a4, stoi(-3), T, p);
     711                 :       3073 :   GEN E6 = Fq_neg(Fq_halve(a6, T, p), T, p);
     712                 :       3073 :   GEN E42 = Fq_sqr(E4, T, p);
     713                 :       3073 :   GEN E43 = Fq_mul(E4, E42, T, p);
     714                 :       3073 :   GEN E62 = Fq_sqr(E6, T, p);
     715                 :       3073 :   GEN delta = Fq_div(Fq_sub(E43, E62, T, p), utoi(1728), T, p);
     716                 :       3073 :   GEN j = Zq_div(E43, delta, T, p, pp, e);
     717                 :       3073 :   GEN Dx = RgX_deriv(meqn);
     718                 :       3073 :   GEN DJ = deriv(meqn, vJ);
     719                 :       3073 :   GEN Dxg = FpXY_Fq_evaly(Dx, g, T, p, vJ);
     720                 :       3073 :   GEN px  = FqX_eval(Dxg, j, T, p), dx  = Fq_mul(px, g, T, p);
     721                 :       3073 :   GEN DJg = FpXY_Fq_evaly(DJ, g, T, p, vJ);
     722                 :       3073 :   GEN pJ = FqX_eval(DJg, j, T, p), dJ = Fq_mul(j, pJ, T, p);
     723                 :       3073 :   GEN Dxx = RgX_deriv(Dx);
     724                 :       3073 :   GEN DxJg = FqX_deriv(Dxg, T, p);
     725                 :            : 
     726                 :       3073 :   GEN ExJ = FqX_eval(DxJg, j, T, p);
     727                 :       3073 :   ulong tis = ugcd(12, ell-1), is = 12 / tis;
     728                 :       3073 :   GEN itis = Fq_inv(stoi(-tis), T, p);
     729                 :       3073 :   GEN deltal = Fq_div(Fq_mul(delta, Fq_powu(g, tis, T, p), T, p), powuu(ell, 12), T, p);
     730                 :            :   GEN E4l, E6l, a4tilde, a6tilde, p_1;
     731         [ +  + ]:       3073 :   if (signe(dx)==0)
     732                 :            :   {
     733         [ -  + ]:          7 :     if (DEBUGLEVEL>0) err_printf("[C: dx=0]");
     734                 :          7 :     avma = ltop; return NULL;
     735                 :            :   }
     736         [ -  + ]:       3066 :   if (signe(dJ)==0)
     737                 :            :   {
     738                 :            :     GEN jl;
     739         [ #  # ]:          0 :     if (DEBUGLEVEL>0) err_printf("[C: dJ=0]");
     740                 :          0 :     E4l = Fq_div(E4, sqru(ell), T, p);
     741                 :          0 :     jl  = Zq_div(Fq_powu(E4l, 3, T, p), deltal, T, p, pp, e);
     742                 :          0 :     E6l = Zq_sqrt(Fq_mul(Fq_sub(jl, utoi(1728), T, p), deltal, T, p), T, p, pp, e);
     743                 :          0 :     p_1 = gen_0;
     744                 :            :   }
     745                 :            :   else
     746                 :            :   {
     747                 :            :     GEN jl, f, fd, Dgs, Djs, jld;
     748                 :       3066 :     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);
     749                 :       3066 :     GEN gd = Fq_mul(Fq_mul(E2s, itis, T, p), g, T, p);
     750                 :       3066 :     GEN jd = Zq_div(Fq_mul(Fq_neg(E42, T, p), E6, T, p), delta, T, p, pp, e);
     751                 :       3066 :     GEN E0b = Zq_div(E6, Fq_mul(E4, E2s, T, p), T, p, pp, e);
     752                 :       3066 :     GEN Dxxgj = FqXY_eval(Dxx, g, j, T, p);
     753                 :       3066 :     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);
     754                 :       3066 :     GEN DJgJj = FqX_eval(FqX_deriv(DJg, T, p), j, T, p);
     755                 :       3066 :     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);
     756                 :       3066 :     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);
     757                 :       3066 :     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);
     758                 :       3066 :     jl = Zq_div(Fq_powu(E4l, 3, T, p), deltal, T, p, pp, e);
     759         [ +  + ]:       3066 :     if (signe(jl)==0)
     760                 :            :     {
     761         [ -  + ]:          7 :       if (DEBUGLEVEL>0) err_printf("[C: jl=0]");
     762                 :          7 :       avma = ltop; return NULL;
     763                 :            :     }
     764                 :       3059 :     f =  Zq_div(powuu(ell, is), g, T, p, pp, e);
     765                 :       3059 :     fd = Fq_neg(Fq_mul(Fq_mul(E2s, f, T, p), itis, T, p), T, p);
     766                 :       3059 :     Dgs = FqXY_eval(Dx, f, jl, T, p);
     767                 :       3059 :     Djs = FqXY_eval(DJ, f, jl, T, p);
     768                 :       3059 :     jld = Zq_div(Fq_mul(Fq_neg(fd, T, p), Dgs, T, p), Fq_mulu(Djs, ell, T, p), T, p, pp, e);
     769                 :       3059 :     E6l = Zq_div(Fq_mul(Fq_neg(E4l, T, p), jld, T, p), jl, T, p, pp, e);
     770                 :       3059 :     p_1 = Fq_neg(Fq_halve(Fq_mulu(E2s, ell, T, p), T, p),T,p);
     771                 :            :   }
     772                 :       3059 :   a4tilde = Fq_mul(Fq_mul(stoi(-3), powuu(ell,4), T, p), E4l, T, p);
     773                 :       3059 :   a6tilde = Fq_mul(Fq_mul(stoi(-2), powuu(ell,6), T, p), E6l, T, p);
     774                 :       3059 :   h = find_kernel(a4, a6, ell, a4tilde, a6tilde, p_1, T, p, pp, e);
     775         [ -  + ]:       3059 :   if (!h) return NULL;
     776                 :       3073 :   return gerepilecopy(ltop, mkvec3(a4tilde, a6tilde, h));
     777                 :            : }
     778                 :            : 
     779                 :            : static GEN
     780                 :       4753 : find_isogenous(GEN a4,GEN a6, ulong ell, struct meqn *MEQN, GEN g, GEN T,GEN p)
     781                 :            : {
     782                 :       4753 :   ulong pp = itou_or_0(p);
     783 [ +  + ][ +  + ]:       4753 :   long e = (pp && pp <= 2*ell+3) ? 2+factorial_lval(ell, pp): 1;
     784         [ +  + ]:       4753 :   if (e > 1)
     785                 :            :   {
     786                 :         35 :     GEN pe = powiu(p, e);
     787                 :         35 :     GEN meqnj = FqXY_evalx(MEQN->eq, Zq_ellj(a4, a6, T, pe, p, e), T, pe);
     788                 :         35 :     g = ZpXQX_liftroot(meqnj, g, T, p, e);
     789                 :            :   }
     790                 :       4753 :   return (MEQN->type == 'C')
     791                 :            :     ? find_isogenous_from_canonical(a4, a6, ell, MEQN, g, T, p, e)
     792         [ +  + ]:       4753 :     : find_isogenous_from_Atkin(a4, a6, ell, MEQN, g, T, p, e);
     793                 :            : }
     794                 :            : 
     795                 :            : static GEN
     796                 :       4389 : FqX_homogenous_eval(GEN P, GEN A, GEN B, GEN T, GEN p)
     797                 :            : {
     798                 :       4389 :   long d = degpol(P), i, v = varn(A);
     799                 :       4389 :   GEN s =  scalar_ZX_shallow(gel(P, d+2), v), Bn = pol_1(v);
     800         [ +  + ]:      17115 :   for (i = d-1; i >= 0; i--)
     801                 :            :   {
     802                 :      12726 :     Bn = FqX_mul(Bn, B, T, p);
     803                 :      12726 :     s = FqX_add(FqX_mul(s, A, T, p), FqX_Fq_mul(Bn, gel(P,i+2), T, p), T, p);
     804                 :            :   }
     805                 :       4389 :   return s;
     806                 :            : }
     807                 :            : 
     808                 :            : static GEN
     809                 :        973 : FqX_homogenous_div(GEN P, GEN Q, GEN A, GEN B, GEN T, GEN p)
     810                 :            : {
     811                 :        973 :   GEN z = cgetg(3, t_RFRAC);
     812                 :        973 :   long d = degpol(Q)-degpol(P);
     813                 :        973 :   gel(z, 1) = FqX_homogenous_eval(P, A, B, T, p);
     814                 :        973 :   gel(z, 2) = FqX_homogenous_eval(Q, A, B, T, p);
     815         [ -  + ]:        973 :   if (d > 0)
     816                 :          0 :     gel(z, 1) = FqX_mul(gel(z, 1), FqX_powu(B, d, T, p), T, p);
     817         [ +  - ]:        973 :   else if (d < 0)
     818                 :        973 :     gel(z, 2) = FqX_mul(gel(z, 2), FqX_powu(B, -d, T, p), T, p);
     819                 :        973 :   return z;
     820                 :            : }
     821                 :            : 
     822                 :            : static GEN
     823                 :        980 : find_kernel_power(GEN Eba4, GEN Eba6, GEN Eca4, GEN Eca6, ulong ell, struct meqn *MEQN, GEN kpoly, GEN Ib, GEN T, GEN p)
     824                 :            : {
     825                 :        980 :   pari_sp ltop = avma, btop;
     826                 :            :   GEN a4t, a6t, gtmp;
     827                 :        980 :   GEN num_iso = FqX_numer_isog_abscissa(kpoly, Eba4, Eba6, T, p, 0);
     828                 :        980 :   GEN mpoly = FqXY_evalx(MEQN->eq, Fq_ellj(Eca4, Eca6, T, p), T, p);
     829                 :        980 :   GEN mroots = FqX_roots(mpoly, T, p);
     830                 :        980 :   GEN kpoly2 = FqX_sqr(kpoly, T, p);
     831                 :        980 :   long i, l1 = lg(mroots);
     832                 :        980 :   btop = avma;
     833         [ +  - ]:       1477 :   for (i = 1; i < l1; i++)
     834                 :            :   {
     835                 :            :     GEN h;
     836                 :       1477 :     GEN tmp = find_isogenous(Eca4, Eca6, ell, MEQN, gel(mroots, i), T, p);
     837         [ +  + ]:       1477 :     if (!tmp) { avma = ltop; return NULL; }
     838                 :       1470 :     a4t =  gel(tmp, 1);
     839                 :       1470 :     a6t =  gel(tmp, 2);
     840                 :       1470 :     gtmp = gel(tmp, 3);
     841                 :            : 
     842                 :            :     /*check that the kernel kpoly is the good one */
     843                 :       1470 :     h = FqX_homogenous_eval(gtmp, num_iso, kpoly2, T, p);
     844         [ +  + ]:       1470 :     if (signe(Fq_elldivpolmod(Eba4, Eba6, ell, h, T, p)))
     845                 :            :     {
     846                 :        973 :       GEN Ic = FqX_homogenous_div(num_iso, kpoly2, numer(Ib), denom(Ib), T, p);
     847                 :        973 :       GEN kpoly_new = FqX_homogenous_eval(gtmp, numer(Ic), denom(Ic), T, p);
     848                 :        973 :       return gerepilecopy(ltop, mkvecn(5, a4t, a6t, kpoly_new, gtmp, Ic));
     849                 :            :     }
     850                 :        497 :     avma = btop;
     851                 :            :   }
     852                 :          0 :   pari_err_BUG("failed to find kernel polynomial");
     853                 :        980 :   return NULL; /*NOT REACHED*/
     854                 :            : }
     855                 :            : 
     856                 :            : /****************************************************************************/
     857                 :            : /*                                  TRACE                                   */
     858                 :            : /****************************************************************************/
     859                 :            : enum mod_type {MTpathological, MTAtkin, MTElkies, MTone_root, MTroots};
     860                 :            : 
     861                 :            : static GEN
     862                 :        306 : Flxq_study_eqn(long ell, GEN mpoly, GEN T, ulong p, long *pt_dG, long *pt_r)
     863                 :            : {
     864                 :        306 :   GEN Xq = FlxqX_Frobenius(mpoly, T, p);
     865                 :        306 :   GEN G  = FlxqX_gcd(FlxX_sub(Xq, pol_x(0), p), mpoly, T, p);
     866                 :        306 :   *pt_dG = degpol(G);
     867         [ +  + ]:        306 :   if (!*pt_dG)
     868                 :            :   {
     869                 :         84 :     GEN L = FlxqXQ_matrix_pow(Xq, ell+1, ell+1, mpoly, T, p);
     870                 :         84 :     long vT = get_Flx_var(T);
     871                 :         84 :     long s = ell + 1 - FlxqM_rank(FlxM_Flx_add_shallow(L, Fl_to_Flx(p-1, vT), p), T, p);
     872                 :         84 :     *pt_r = (ell + 1)/s;
     873                 :         84 :     return NULL;
     874                 :            :   }
     875                 :        306 :   return G;
     876                 :            : }
     877                 :            : 
     878                 :            : static GEN
     879                 :       5768 : Fp_study_eqn(long ell, GEN mpoly, GEN p, long *pt_dG, long *pt_r)
     880                 :            : {
     881                 :       5768 :   GEN XP = FpX_Frobenius(mpoly, p);
     882                 :       5768 :   GEN G  = FpX_gcd(FpX_sub(XP, pol_x(0), p), mpoly, p);
     883                 :       5768 :   *pt_dG = degpol(G);
     884         [ +  + ]:       5768 :   if (!*pt_dG)
     885                 :            :   {
     886                 :       2716 :     GEN L = FpXQ_matrix_pow(XP, ell+1, ell+1, mpoly, p);
     887                 :       2716 :     long s = ell + 1 - FpM_rank(RgM_Rg_add_shallow(L, gen_m1), p);
     888                 :       2716 :     *pt_r = (ell + 1)/s;
     889                 :       2716 :     return NULL;
     890                 :            :   }
     891                 :       5768 :   return FpX_oneroot(G, p);
     892                 :            : }
     893                 :            : 
     894                 :            : static GEN
     895                 :        427 : FpXQ_study_eqn(long ell, GEN mpoly, GEN T, GEN p, long *pt_dG, long *pt_r)
     896                 :            : {
     897                 :            :   GEN G;
     898         [ +  + ]:        427 :   if (lgefint(p)==3)
     899                 :            :   {
     900                 :        306 :     ulong pp = p[2];
     901                 :        306 :     GEN Tp = ZXT_to_FlxT(T,pp);
     902                 :        306 :     GEN mpolyp = ZXX_to_FlxX(mpoly,pp,get_FpX_var(T));
     903                 :        306 :     G = Flxq_study_eqn(ell, mpolyp, Tp, pp, pt_dG, pt_r);
     904         [ +  + ]:        306 :     if (!G) return NULL;
     905                 :        222 :     G = FlxX_to_ZXX(G);
     906                 :            :   }
     907                 :            :   else
     908                 :            :   {
     909                 :        121 :     GEN Xq = FpXQX_Frobenius(mpoly, T, p);
     910                 :        121 :     G  = FpXQX_gcd(FpXX_sub(Xq, pol_x(0), p), mpoly, T, p);
     911                 :        121 :     *pt_dG = degpol(G);
     912         [ +  + ]:        121 :     if (!*pt_dG)
     913                 :            :     {
     914                 :         49 :       GEN L = FpXQXQ_matrix_pow(Xq, ell+1, ell+1, mpoly, T, p);
     915                 :         49 :       long s = ell + 1 - FqM_rank(RgM_Rg_add(L, gen_m1), T, p);
     916                 :         49 :       *pt_r = (ell + 1)/s;
     917                 :         49 :       return NULL;
     918                 :            :     }
     919                 :            :   }
     920                 :        427 :   return gel(FqX_roots(G, T, p), 1);
     921                 :            : }
     922                 :            : 
     923                 :            : /* Berlekamp variant */
     924                 :            : static GEN
     925                 :       6202 : study_modular_eqn(long ell, GEN mpoly, GEN T, GEN p, enum mod_type *mt, long *ptr_r)
     926                 :            : {
     927                 :       6202 :   pari_sp ltop = avma;
     928                 :       6202 :   GEN g = gen_0;
     929                 :       6202 :   *ptr_r = 0; /*gcc -Wall*/
     930         [ +  + ]:       6202 :   if (degpol(FqX_gcd(mpoly, FqX_deriv(mpoly, T, p), T, p)) > 0)
     931                 :          7 :     *mt = MTpathological;
     932                 :            :   else
     933                 :            :   {
     934                 :            :     long dG;
     935                 :       6195 :     g = T ? FpXQ_study_eqn(ell, mpoly, T, p, &dG, ptr_r)
     936         [ +  + ]:       6195 :             : Fp_study_eqn(ell, mpoly, p, &dG, ptr_r);
     937   [ +  +  +  + ]:       6195 :     switch(dG)
     938                 :            :     {
     939                 :       2849 :       case 0:  *mt = MTAtkin; break;
     940                 :        378 :       case 1:  *mt = MTone_root; break;
     941                 :       2898 :       case 2:  *mt = MTElkies;   break;
     942         [ +  - ]:       6195 :       default: *mt = (dG == ell + 1)? MTroots: MTpathological;
     943                 :            :     }
     944                 :            :   }
     945 [ -  + ][ #  #  :       6202 :   if (DEBUGLEVEL) switch(*mt)
             #  #  #  # ]
     946                 :            :   {
     947                 :          0 :     case MTone_root: err_printf("One root\t"); break;
     948                 :          0 :     case MTElkies: err_printf("Elkies\t"); break;
     949                 :          0 :     case MTroots: err_printf("l+1 roots\t"); break;
     950                 :          0 :     case MTAtkin: err_printf("Atkin\t"); break;
     951                 :          0 :     case MTpathological: err_printf("Pathological\n"); break;
     952                 :            :   }
     953         [ +  + ]:       6202 :   return g ? gerepilecopy(ltop, g): NULL;
     954                 :            : }
     955                 :            : 
     956                 :            : /*Returns the trace modulo ell^k when ell is an Elkies prime */
     957                 :            : static GEN
     958                 :       3276 : 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)
     959                 :            : {
     960                 :       3276 :   pari_sp ltop = avma, btop;
     961                 :            :   GEN tmp, Eba4, Eba6, Eca4, Eca6, Ib, kpoly;
     962                 :       3276 :   ulong lambda, ellk = upowuu(ell, k), pellk = umodiu(q, ellk);
     963                 :            :   long cnt;
     964                 :            : 
     965         [ -  + ]:       3276 :   if (DEBUGLEVEL) { err_printf("Trace mod %ld", ell); }
     966                 :       3276 :   Eba4 = a4;
     967                 :       3276 :   Eba6 = a6;
     968                 :       3276 :   tmp = find_isogenous(a4,a6, ell, MEQN, g, T, p);
     969         [ +  + ]:       3276 :   if (!tmp) { avma = ltop; return NULL; }
     970                 :       3227 :   Eca4 =  gel(tmp, 1);
     971                 :       3227 :   Eca6 =  gel(tmp, 2);
     972                 :       3227 :   kpoly = gel(tmp, 3);
     973                 :       3227 :   Ib = pol_x(0);
     974                 :       3227 :   lambda = find_eigen_value(a4, a6, ell, kpoly, T, p, tr);
     975         [ -  + ]:       3227 :   if (DEBUGLEVEL>1) err_printf(" [%ld ms]", timer_delay(ti));
     976 [ -  + ][ #  # ]:       3227 :   if (smallfact && ell>smallfact)
     977                 :            :   {
     978                 :          0 :     ulong pell = pellk%ell;
     979                 :          0 :     ulong ap = Fl_add(lambda, Fl_div(pell, lambda, ell), ell);
     980         [ #  # ]:          0 :     if (Fl_sub(pell, ap, ell)==ell-1) { avma = ltop; return mkvecsmall(ap); }
     981                 :            :   }
     982                 :       3227 :   btop = avma;
     983         [ +  + ]:       4200 :   for (cnt = 2; cnt <= k; cnt++)
     984                 :            :   {
     985                 :            :     GEN tmp;
     986         [ -  + ]:        980 :     if (DEBUGLEVEL) err_printf(", %Ps", powuu(ell, cnt));
     987                 :        980 :     tmp = find_kernel_power(Eba4, Eba6, Eca4, Eca6, ell, MEQN, kpoly, Ib, T, p);
     988         [ +  + ]:        980 :     if (!tmp) { avma = ltop; return NULL; }
     989                 :        973 :     lambda = find_eigen_value_power(a4, a6, ell, cnt, gel(tmp,3), lambda, T, p);
     990                 :        973 :     Eba4 = Eca4;
     991                 :        973 :     Eba6 = Eca6;
     992                 :        973 :     Eca4 = gel(tmp,1);
     993                 :        973 :     Eca6 = gel(tmp,2);
     994                 :        973 :     kpoly = gel(tmp,4);
     995                 :        973 :     Ib = gel(tmp, 5);
     996         [ -  + ]:        973 :     if (gc_needed(btop, 1))
     997                 :            :     {
     998         [ #  # ]:          0 :       if(DEBUGMEM>1) pari_warn(warnmem,"find_trace_Elkies_power");
     999                 :          0 :       gerepileall(btop, 6, &Eba4, &Eba6, &Eca4, &Eca6, &kpoly, &Ib);
    1000                 :            :     }
    1001         [ -  + ]:        973 :     if (DEBUGLEVEL>1) err_printf(" [%ld ms]", timer_delay(ti));
    1002                 :            :   }
    1003                 :       3220 :   avma = ltop;
    1004                 :       3276 :   return mkvecsmall(Fl_add(lambda, Fl_div(pellk, lambda, ellk), ellk));
    1005                 :            : }
    1006                 :            : 
    1007                 :            : /*Returns the possible values of the trace when ell is an Atkin prime, */
    1008                 :            : /*given r the splitting degree of the modular equation at J = E.j */
    1009                 :            : static GEN
    1010                 :       2849 : find_trace_Atkin(ulong ell, long r, GEN q)
    1011                 :            : {
    1012                 :       2849 :   pari_sp ltop = avma;
    1013                 :       2849 :   long nval = 0;
    1014                 :       2849 :   ulong teta, pell = umodiu(q, ell), invp = Fl_inv(pell, ell);
    1015                 :       2849 :   GEN val_pos = cgetg(1+ell, t_VECSMALL), P = gel(factoru(r), 1);
    1016                 :       2849 :   GEN S = mkvecsmall4(0, pell, 0, 1);
    1017                 :       2849 :   GEN U = mkvecsmall3(0, ell-1, 0);
    1018                 :       2849 :   pari_sp btop = avma;
    1019 [ +  + ][ +  - ]:       2849 :   if (r==2 && krouu(ell-pell, ell) < 0)
    1020                 :        532 :     val_pos[++nval] = 0;
    1021         [ +  + ]:      58597 :   for (teta = 1; teta < ell; teta++, avma = btop)
    1022                 :            :   {
    1023                 :      55748 :     ulong disc = Fl_sub(Fl_sqr(teta,ell), Fl_mul(4UL,pell,ell), ell);
    1024                 :            :     GEN a;
    1025         [ +  + ]:      55748 :     if (krouu(disc, ell) >= 0) continue;
    1026                 :      27398 :     S[3] = Fl_neg(teta, ell);
    1027                 :      27398 :     U[3] = Fl_mul(invp, teta, ell);
    1028                 :      27398 :     a = Flxq_powu(U, r/P[1], S, ell);
    1029 [ +  + ][ +  + ]:      27398 :     if (!Flx_equal1(a) && Flx_equal1(Flxq_powu(a, P[1], S, ell)))
    1030                 :            :     {
    1031                 :      17458 :       pari_sp av = avma;
    1032                 :      17458 :       long i, l=lg(P);
    1033         [ +  + ]:      29400 :       for (i = 2; i < l; i++, avma = av)
    1034         [ +  + ]:      15232 :         if (Flx_equal1(Flxq_powu(U, r/P[i], S, ell))) break;
    1035         [ +  + ]:      17458 :       if (i==l) val_pos[++nval] = teta;
    1036                 :            :     }
    1037                 :            :   }
    1038                 :       2849 :   return gerepileupto(ltop, vecsmall_shorten(val_pos, nval));
    1039                 :            : }
    1040                 :            : 
    1041                 :            : /*Returns the possible traces when there is only one root */
    1042                 :            : static GEN
    1043                 :        378 : find_trace_one_root(ulong ell, GEN q)
    1044                 :            : {
    1045                 :        378 :   ulong a = Fl_double(Fl_sqrt(umodiu(q,ell), ell), ell);
    1046                 :        378 :   return mkvecsmall2(a, ell - a);
    1047                 :            : }
    1048                 :            : 
    1049                 :            : static GEN
    1050                 :         70 : find_trace_lp1_roots(long ell, GEN q)
    1051                 :            : {
    1052                 :         70 :   ulong ell2 = ell * ell, pell = umodiu(q, ell2);
    1053                 :         70 :   ulong a  = Fl_sqrt(pell%ell, ell);
    1054                 :         70 :   ulong pa = Fl_add(Fl_div(pell, a, ell2), a, ell2);
    1055                 :         70 :   return mkvecsmall2(pa, ell2 - pa);
    1056                 :            : }
    1057                 :            : 
    1058                 :            : /*trace modulo ell^k: [], [t] or [t1,...,td] */
    1059                 :            : static GEN
    1060                 :       6202 : find_trace(GEN a4, GEN a6, GEN j, ulong ell, GEN q, GEN T, GEN p, long *ptr_kt,
    1061                 :            :   ulong smallfact, long vx, long vy)
    1062                 :            : {
    1063                 :       6202 :   pari_sp ltop = avma;
    1064                 :            :   GEN g, meqnj, tr, tr2;
    1065                 :       6202 :   long k = 1, kt, r;
    1066                 :            :   enum mod_type mt;
    1067                 :            :   struct meqn MEQN;
    1068                 :            :   pari_timer ti;
    1069                 :            : 
    1070         [ +  + ]:       6202 :   if (ell <= 13)
    1071                 :            :   {
    1072                 :       3010 :     long lp = expi(q);
    1073   [ +  +  +  + ]:       3010 :     switch(ell)
    1074                 :            :     {
    1075         [ +  + ]:        602 :       case 3: k = 3 + (lp > 160) + (lp > 350); break;
    1076         [ -  + ]:        602 :       case 5: k = 2 + (lp > 260); break;
    1077         [ -  + ]:        602 :       case 7: k = 2 + (lp > 390); break;
    1078         [ -  + ]:       1204 :       default:k = 1 + (lp > 260);
    1079                 :            :     }
    1080                 :            :   }
    1081                 :       6202 :   kt = k;
    1082         [ -  + ]:       6202 :   if (!get_modular_eqn(&MEQN, ell, vx, vy)) err_modular_eqn(ell);
    1083         [ -  + ]:       6202 :   if (DEBUGLEVEL)
    1084                 :          0 :   { err_printf("Process prime %5ld. ", ell); timer_start(&ti); }
    1085                 :       6202 :   meqnj = FqXY_evalx(MEQN.eq, j, T, p);
    1086                 :       6202 :   g = study_modular_eqn(ell, meqnj, T, p, &mt, &r);
    1087                 :            :   /* If l is an Elkies prime, search for a factor of the l-division polynomial.
    1088                 :            :   * Then deduce the trace by looking for eigenvalues of the Frobenius by
    1089                 :            :   * computing modulo this factor */
    1090   [ +  +  +  +  :       6202 :   switch (mt)
                      + ]
    1091                 :            :   {
    1092                 :            :   case MTone_root:
    1093                 :        378 :     tr2 = find_trace_one_root(ell, q);
    1094                 :        378 :     kt = k = 1;
    1095                 :            :     /* Must take k = 1 because we can't apply Hensel: no guarantee that a
    1096                 :            :      * root mod ell^2 exists */
    1097                 :        378 :     tr = find_trace_Elkies_power(a4,a6,ell, k, &MEQN, g, tr2, q, T, p, smallfact, &ti);
    1098         [ +  + ]:        378 :     if (!tr) tr = tr2;
    1099                 :        378 :     break;
    1100                 :            :   case MTElkies:
    1101                 :            :     /* Contrary to MTone_root, may look mod higher powers of ell */
    1102         [ +  + ]:       2898 :     if (cmpiu(p, 2*ell+3) <= 0)
    1103                 :         28 :       kt = k = 1; /* Not implemented in this case */
    1104                 :       2898 :     tr = find_trace_Elkies_power(a4,a6,ell, k, &MEQN, g, NULL, q, T, p, smallfact, &ti);
    1105         [ +  + ]:       2898 :     if (!tr)
    1106                 :            :     {
    1107         [ -  + ]:         21 :       if (DEBUGLEVEL) err_printf("[fail]\n");
    1108                 :         21 :       avma = ltop; return NULL;
    1109                 :            :     }
    1110                 :       2877 :     break;
    1111                 :            :   case MTroots:
    1112                 :         70 :     tr = find_trace_lp1_roots(ell, q);
    1113                 :         70 :     kt = 2;
    1114                 :         70 :     break;
    1115                 :            :   case MTAtkin:
    1116                 :       2849 :     tr = find_trace_Atkin(ell, r, q);
    1117         [ -  + ]:       2849 :     if (lg(tr)==1) pari_err_PRIME("ellsea",p);
    1118                 :       2849 :     kt = 1;
    1119                 :       2849 :     break;
    1120                 :            :   default: /* case MTpathological: */
    1121                 :          7 :     avma = ltop; return NULL;
    1122                 :            :   }
    1123         [ -  + ]:       6174 :   if (DEBUGLEVEL) {
    1124                 :          0 :     long n = lg(tr)-1;
    1125 [ #  # ][ #  # ]:          0 :     if (n > 1 || mt == MTAtkin)
    1126                 :            :     {
    1127                 :          0 :       err_printf("%3ld trace(s)",n);
    1128         [ #  # ]:          0 :       if (DEBUGLEVEL>1) err_printf(" [%ld ms]", timer_delay(&ti));
    1129                 :            :     }
    1130                 :          0 :     err_printf("\n");
    1131                 :            :   }
    1132                 :       6174 :   *ptr_kt = kt;
    1133                 :       6202 :   return gerepileupto(ltop, tr);
    1134                 :            : }
    1135                 :            : 
    1136                 :            : /* A partition of compile_atkin in baby and giant is represented as the binary
    1137                 :            :    developpement of an integer; if the i-th bit is 1, the i-th prime in
    1138                 :            :    compile-atkin is a baby. The optimum is obtained when the ratio between
    1139                 :            :    the number of possibilities for traces modulo giants (p_g) and babies (p_b)
    1140                 :            :    is near 3/4. */
    1141                 :            : static long
    1142                 :        567 : separation(GEN cnt)
    1143                 :            : {
    1144                 :            :   pari_sp btop;
    1145                 :        567 :   long k = lg(cnt)-1, l = (1L<<k)-1, best_i, i, j;
    1146                 :            :   GEN best_r, P, P3, r;
    1147                 :            : 
    1148                 :        567 :   P = gen_1;
    1149         [ +  + ]:       2779 :   for (j = 1; j <= k; ++j) P = mulis(P, cnt[j]);
    1150                 :            :   /* p_b * p_g = P is constant */
    1151                 :        567 :   P3 = mulsi(3, P);
    1152                 :        567 :   btop = avma;
    1153                 :        567 :   best_i = 0;
    1154                 :        567 :   best_r = P3;
    1155         [ +  + ]:      25557 :   for (i = 1; i < l; i++)
    1156                 :            :   {
    1157                 :            :     /* scan all possibilities */
    1158                 :      25032 :     GEN p_b = gen_1;
    1159         [ +  + ]:     224798 :     for (j = 0; j < k; j++)
    1160         [ +  + ]:     199766 :       if (i & (1L<<j)) p_b = mulis(p_b, cnt[1+j]);
    1161                 :      25032 :     r = subii(shifti(sqri(p_b), 2), P3); /* (p_b/p_g - 3/4)*4*P */
    1162         [ +  + ]:      25032 :     if (!signe(r)) { best_i = i; break; }
    1163         [ +  + ]:      24990 :     if (absi_cmp(r, best_r) < 0) { best_i = i; best_r = r; }
    1164         [ -  + ]:      24990 :     if (gc_needed(btop, 1))
    1165                 :          0 :       best_r = gerepileuptoint(btop, best_r);
    1166                 :            :   }
    1167                 :        567 :   return best_i;
    1168                 :            : }
    1169                 :            : 
    1170                 :            : /* x VEC defined modulo P (= *P), y VECSMALL modulo q, (q,P) = 1. */
    1171                 :            : /* Update in place:
    1172                 :            :  *   x to vector mod q P congruent to x mod P (resp. y mod q). */
    1173                 :            : /*   P ( <-- qP ) */
    1174                 :            : static void
    1175                 :       1078 : multiple_crt(GEN x, GEN y, GEN q, GEN P)
    1176                 :            : {
    1177                 :       1078 :   pari_sp ltop = avma, av;
    1178                 :       1078 :   long i, j, k, lx = lg(x)-1, ly = lg(y)-1;
    1179                 :            :   GEN  a1, a2, u, v, A2X;
    1180                 :       1078 :   (void)bezout(P,q,&u,&v);
    1181                 :       1078 :   a1 = mulii(P,u);
    1182                 :       1078 :   a2 = mulii(q,v); A2X = ZC_Z_mul(x, a2);
    1183                 :       1078 :   av = avma; affii(mulii(P,q), P);
    1184         [ +  + ]:      56938 :   for (i = 1, k = 1; i <= lx; i++, avma = av)
    1185                 :            :   {
    1186                 :      55860 :     GEN a2x = gel(A2X,i);
    1187         [ +  + ]:     983164 :     for (j = 1; j <= ly; ++j)
    1188                 :            :     {
    1189                 :     927304 :       GEN t = Fp_add(Fp_mulu(a1, y[j], P), a2x, P);
    1190                 :     927304 :       affii(t, gel(x, k++));
    1191                 :            :     }
    1192                 :            :   }
    1193                 :       1078 :   setlg(x, k); avma = ltop;
    1194                 :       1078 : }
    1195                 :            : 
    1196                 :            : /****************************************************************************/
    1197                 :            : /*                              MATCH AND SORT                              */
    1198                 :            : /****************************************************************************/
    1199                 :            : 
    1200                 :            : static GEN
    1201                 :       1134 : possible_traces(GEN compile, GEN mask, GEN *P, int larger)
    1202                 :            : {
    1203                 :       1134 :   GEN V, Pfinal = gen_1, C = shallowextract(compile, mask);
    1204                 :       1134 :   long i, lfinal = 1, lC = lg(C), lP;
    1205                 :       1134 :   pari_sp av = avma;
    1206                 :            : 
    1207         [ +  + ]:       3346 :   for (i = 1; i < lC; i++)
    1208                 :            :   {
    1209                 :       2212 :     GEN c = gel(C,i), t;
    1210                 :       2212 :     Pfinal = mulii(Pfinal, gel(c,1));
    1211                 :       2212 :     t = muluu(lfinal, lg(gel(c,2))-1);
    1212                 :       2212 :     lfinal = itou(t);
    1213                 :            :   }
    1214                 :       1134 :   Pfinal = gerepileuptoint(av, Pfinal);
    1215         [ +  + ]:       1134 :   if (larger)
    1216                 :        567 :     lP = lgefint(shifti(Pfinal,1));
    1217                 :            :   else
    1218                 :        567 :     lP = lgefint(Pfinal);
    1219                 :       1134 :   lfinal++;
    1220                 :            :   /* allocate room for final result */
    1221                 :       1134 :   V = cgetg(lfinal, t_VEC);
    1222         [ +  + ]:     877100 :   for (i = 1; i < lfinal; i++) gel(V,i) = cgeti(lP);
    1223                 :            : 
    1224                 :            :   {
    1225                 :       1134 :     GEN c = gel(C,1), v = gel(c,2);
    1226                 :       1134 :     long l = lg(v);
    1227         [ +  + ]:       5656 :     for (i = 1; i < l; i++) affsi(v[i], gel(V,i));
    1228                 :       1134 :     setlg(V, l); affii(gel(c,1), Pfinal); /* reset Pfinal */
    1229                 :            :   }
    1230         [ +  + ]:       2212 :   for (i = 2; i < lC; i++)
    1231                 :            :   {
    1232                 :       1078 :     GEN c = gel(C,i);
    1233                 :       1078 :     multiple_crt(V, gel(c,2), gel(c,1), Pfinal); /* Pfinal updated! */
    1234                 :            :   }
    1235                 :       1134 :   *P = Pfinal; return V;
    1236                 :            : }
    1237                 :            : 
    1238                 :            : static GEN
    1239                 :     163996 : cost(long mask, GEN cost_vec)
    1240                 :            : {
    1241                 :     163996 :   pari_sp ltop = avma;
    1242                 :            :   long i;
    1243                 :     163996 :   GEN c = gen_1;
    1244         [ +  + ]:    1853089 :   for (i = 1; i < lg(cost_vec); i++)
    1245         [ +  + ]:    1689093 :     if (mask&(1L<<(i-1)))
    1246                 :     732886 :       c = mulis(c, cost_vec[i]);
    1247                 :     163996 :   return gerepileuptoint(ltop, c);
    1248                 :            : }
    1249                 :            : 
    1250                 :            : static GEN
    1251                 :     132258 : value(long mask, GEN atkin, long k)
    1252                 :            : {
    1253                 :     132258 :   pari_sp ltop = avma;
    1254                 :            :   long i;
    1255                 :     132258 :   GEN c = gen_1;
    1256         [ +  + ]:    1493793 :   for (i = 1; i <= k; i++)
    1257         [ +  + ]:    1361535 :     if (mask&(1L<<(i-1)))
    1258                 :     591066 :       c = mulii(c, gmael(atkin, i, 1));
    1259                 :     132258 :   return gerepileuptoint(ltop, c);
    1260                 :            : }
    1261                 :            : 
    1262                 :            : static void
    1263                 :      64932 : set_cost(GEN B, long b, GEN cost_vec, long *pi)
    1264                 :            : {
    1265                 :      64932 :   pari_sp av = avma;
    1266                 :      64932 :   GEN costb = cost(b, cost_vec);
    1267                 :      64932 :   long i = *pi;
    1268         [ +  + ]:      86352 :   while (cmpii(costb, cost(B[i], cost_vec)) < 0) --i;
    1269                 :      64932 :   B[++i] = b;
    1270                 :      64932 :   *pi = i; avma = av;
    1271                 :      64932 : }
    1272                 :            : 
    1273                 :            : static GEN
    1274                 :       1169 : get_lgatkin(GEN compile_atkin, long k)
    1275                 :            : {
    1276                 :       1169 :   GEN v = cgetg(k+1, t_VECSMALL);
    1277                 :            :   long j;
    1278         [ +  + ]:       6125 :   for (j = 1; j <= k; ++j) v[j] = lg(gmael(compile_atkin, j, 2))-1;
    1279                 :       1169 :   return v;
    1280                 :            : }
    1281                 :            : 
    1282                 :            : static GEN
    1283                 :        602 : champion(GEN atkin, long k, GEN bound_champ)
    1284                 :            : {
    1285                 :        602 :   const long two_k = 1L<<k;
    1286                 :        602 :   pari_sp ltop = avma;
    1287                 :            :   long i, j, n, i1, i2;
    1288                 :        602 :   GEN B, Bp, cost_vec, res = NULL;
    1289                 :            : 
    1290                 :        602 :   cost_vec = get_lgatkin(atkin, k);
    1291         [ +  + ]:        602 :   if (k == 1) return mkvec2(gen_1, utoipos(cost_vec[1]));
    1292                 :            : 
    1293                 :        595 :   B  = zero_zv(two_k);
    1294                 :        595 :   Bp = zero_zv(two_k);
    1295                 :        595 :   Bp[2] = 1;
    1296         [ +  + ]:       2737 :   for (n = 2, j = 2; j <= k; j++)
    1297                 :            :   {
    1298                 :            :     long b;
    1299                 :       2142 :     i = 1;
    1300         [ +  + ]:      60704 :     for (i1 = 2, i2 = 1; i1 <= n; )
    1301                 :            :     {
    1302                 :      58562 :       pari_sp av = avma;
    1303                 :      58562 :       long b1 = Bp[i1], b2 = Bp[i2]|(1L<<(j-1));
    1304         [ +  + ]:      58562 :       if (cmpii(value(b1, atkin, k), value(b2, atkin, k)) < 0)
    1305                 :      58562 :         { b = b1; i1++; } else { b = b2; i2++; }
    1306                 :      58562 :       avma = av;
    1307                 :      58562 :       set_cost(B, b, cost_vec, &i);
    1308                 :            :     }
    1309         [ +  + ]:       8512 :     for ( ; i2 <= n; i2++)
    1310                 :            :     {
    1311                 :       6370 :       b = Bp[i2]|(1L<<(j-1));
    1312                 :       6370 :       set_cost(B, b, cost_vec, &i);
    1313                 :            :     }
    1314                 :       2142 :     n = i;
    1315         [ +  + ]:      47796 :     for (i = 1; i <= n; i++)
    1316                 :      45654 :       Bp[i] = B[i];
    1317                 :            :   }
    1318         [ +  + ]:     224679 :   for (i = 1; i <= two_k; i++)
    1319         [ +  + ]:     224084 :     if (B[i])
    1320                 :            :     {
    1321                 :      12712 :       GEN b = cost (B[i], cost_vec);
    1322                 :      12712 :       GEN v = value(B[i], atkin, k);
    1323         [ +  + ]:      12712 :       if (cmpii(v, bound_champ) <=0) continue;
    1324 [ +  + ][ +  - ]:       1484 :       if (res && gcmp(b, gel(res, 2)) >=0) continue;
    1325                 :        595 :       res = mkvec2(utoi(B[i]), b);
    1326                 :            :     }
    1327                 :        602 :   return gerepilecopy(ltop, res);
    1328                 :            : }
    1329                 :            : 
    1330                 :            : static GEN
    1331                 :       1134 : compute_diff(GEN v)
    1332                 :            : {
    1333                 :       1134 :   pari_sp av = avma;
    1334                 :       1134 :   long i, l = lg(v) - 1;
    1335                 :       1134 :   GEN diff = cgetg(l, t_VEC);
    1336         [ +  + ]:     875966 :   for (i = 1; i < l; i++) gel(diff, i) = subii(gel(v, i+1), gel(v, i));
    1337                 :       1134 :   return gerepileupto(av, ZV_sort_uniq(diff));
    1338                 :            : }
    1339                 :            : 
    1340                 :            : static int
    1341                 :      10192 : cmp_atkin(void*E, GEN a, GEN b)
    1342                 :            : {
    1343                 :      10192 :   long ta=typ(a)==t_INT, tb=typ(b)==t_INT, c;
    1344                 :            :   (void) E;
    1345 [ +  - ][ +  + ]:      10192 :   if (ta || tb) return ta-tb;
    1346                 :       3395 :   c = lg(gel(a,2)) - lg(gel(b,2));
    1347         [ +  + ]:       3395 :   if (c) return c;
    1348                 :      10192 :   return cmpii(gel(b,1), gel(a,1));
    1349                 :            : }
    1350                 :            : 
    1351                 :            : static void
    1352                 :       2422 : add_atkin(GEN atkin, GEN trace, long *nb)
    1353                 :            : {
    1354                 :       2422 :   long l = lg(atkin)-1;
    1355                 :       2422 :   long i, k = gen_search(atkin, trace, 1, NULL, cmp_atkin);
    1356 [ +  - ][ -  + ]:       4844 :   if (k==0 || k > l) return;
    1357         [ +  + ]:      47054 :   for (i = l; i > k; i--)
    1358                 :      44632 :     gel(atkin,i) = gel(atkin,i-1);
    1359         [ +  - ]:       2422 :   if (typ(gel(atkin,l))==t_INT) (*nb)++;
    1360                 :       2422 :   gel(atkin,k) = trace;
    1361                 :            : }
    1362                 :            : 
    1363                 :            : /* V = baby / giant, P = Pb / Pg */
    1364                 :            : static GEN
    1365                 :       1134 : BSGS_pre(GEN *pdiff, GEN V, GEN P, void *E, const struct bb_group *grp)
    1366                 :            : {
    1367                 :       1134 :   GEN diff = compute_diff(V);
    1368                 :       1134 :   GEN pre = cgetg(lg(diff), t_VEC);
    1369                 :       1134 :   long i, l = lg(diff);
    1370                 :       1134 :   gel(pre, 1) = grp->pow(E, P, gel(diff, 1));
    1371                 :            :   /* what we'd _really_ want here is a hashtable diff[i] -> pre[i]  */
    1372         [ +  + ]:      28609 :   for (i = 2; i < l; i++)
    1373                 :            :   {
    1374                 :      27475 :     pari_sp av = avma;
    1375                 :      27475 :     GEN d = subii(gel(diff, i), gel(diff, i-1));
    1376                 :      27475 :     GEN Q = grp->mul(E, gel(pre, i-1), grp->pow(E, P, d));
    1377                 :      27475 :     gel(pre, i) = gerepilecopy(av, Q);
    1378                 :            :   }
    1379                 :       1134 :   *pdiff = diff; return pre;
    1380                 :            : }
    1381                 :            : 
    1382                 :            : /* u = trace_elkies, Mu = prod_elkies. Let caller collect garbage */
    1383                 :            : /* Match & sort: variant from Lercier's thesis, section 11.2.3 */
    1384                 :            : /* baby/giant/table updated in place: this routines uses
    1385                 :            :  *   size(baby)+size(giant)+size(table)+size(table_ind) + O(log p)
    1386                 :            :  * bits of stack */
    1387                 :            : static GEN
    1388                 :        574 : match_and_sort(GEN compile_atkin, GEN Mu, GEN u, GEN q, void *E, const struct bb_group *grp)
    1389                 :            : {
    1390                 :            :   pari_sp av1, av2;
    1391                 :        574 :   GEN baby, giant, SgMb, Mb, Mg, den, Sg, dec_inf, div, pp1 = addis(q,1);
    1392                 :            :   GEN P, Pb, Pg, point, diff, pre, table, table_ind;
    1393                 :        574 :   long best_i, i, lbaby, lgiant, k = lg(compile_atkin)-1;
    1394                 :        574 :   GEN bound = sqrti(shifti(q, 2)), card;
    1395                 :        574 :   const long lcard = 100;
    1396                 :        574 :   long lq = lgefint(q), nbcard;
    1397                 :            :   pari_timer ti;
    1398                 :            : 
    1399         [ +  + ]:        574 :   if (k == 1)
    1400                 :            :   { /*only one Atkin prime, check the cardinality with random points */
    1401                 :          7 :     GEN r = gel(compile_atkin, 1), r1 = gel(r,1), r2 = gel(r,2);
    1402                 :          7 :     long l = lg(r2), j;
    1403                 :          7 :     GEN card = cgetg(l, t_VEC), Cs2, C, U;
    1404                 :          7 :     Z_chinese_pre(Mu, r1, &C,&U, NULL);
    1405                 :          7 :     Cs2 = shifti(C, -1);
    1406         [ +  + ]:         21 :     for (j = 1, i = 1; i < l; i++)
    1407                 :            :     {
    1408                 :         14 :       GEN t = Z_chinese_post(u, stoi(r2[i]), C, U, NULL);
    1409                 :         14 :       t = Fp_center(t, C, Cs2);
    1410         [ +  - ]:         14 :       if (absi_cmp(t, bound) <= 0)
    1411                 :         14 :         gel(card, j++) = subii(pp1, t);
    1412                 :            :     }
    1413                 :          7 :     setlg(card, j);
    1414                 :          7 :     return gen_select_order(card, E, grp);
    1415                 :            :   }
    1416         [ -  + ]:        567 :   if (DEBUGLEVEL>=2) timer_start(&ti);
    1417                 :        567 :   av1 = avma;
    1418                 :        567 :   best_i = separation( get_lgatkin(compile_atkin, k) );
    1419                 :        567 :   avma = av1;
    1420                 :            : 
    1421                 :        567 :   baby  = possible_traces(compile_atkin, stoi(best_i), &Mb, 1);
    1422                 :        567 :   giant = possible_traces(compile_atkin, subis(int2n(k), best_i+1), &Mg, 0);
    1423                 :        567 :   lbaby = lg(baby);
    1424                 :        567 :   lgiant = lg(giant);
    1425                 :        567 :   den = Fp_inv(Fp_mul(Mu, Mb, Mg), Mg);
    1426                 :        567 :   av2 = avma;
    1427         [ +  + ]:     511245 :   for (i = 1; i < lgiant; i++, avma = av2)
    1428                 :     510678 :     affii(Fp_mul(gel(giant,i), den, Mg), gel(giant,i));
    1429                 :        567 :   gen_sort_inplace(giant, (void*)&cmpii, &cmp_nodata, NULL);
    1430                 :        567 :   Sg = Fp_mul(negi(u), den, Mg);
    1431                 :        567 :   den = Fp_inv(Fp_mul(Mu, Mg, Mb), Mb);
    1432                 :        567 :   dec_inf = divii(mulii(Mb,addii(Mg,shifti(Sg,1))), shifti(Mg,1));
    1433                 :        567 :   togglesign(dec_inf); /* now, dec_inf = ceil(- (Mb/2 + Sg Mb/Mg) ) */
    1434                 :        567 :   div = mulii(truedivii(dec_inf, Mb), Mb);
    1435                 :        567 :   av2 = avma;
    1436         [ +  + ]:     365855 :   for (i = 1; i < lbaby; i++, avma = av2)
    1437                 :            :   {
    1438                 :     365288 :     GEN b = addii(Fp_mul(Fp_sub(gel(baby,i), u, Mb), den, Mb), div);
    1439         [ +  + ]:     365288 :     if (cmpii(b, dec_inf) < 0) b = addii(b, Mb);
    1440                 :     365288 :     affii(b, gel(baby,i));
    1441                 :            :   }
    1442                 :        567 :   gen_sort_inplace(baby, (void*)&cmpii, &cmp_nodata, NULL);
    1443                 :            : 
    1444                 :        567 :   SgMb = mulii(Sg, Mb);
    1445                 :        567 :   card = cgetg(lcard+1,t_VEC);
    1446         [ +  + ]:      57267 :   for (i = 1; i <= lcard; i++) gel(card,i) = cgetipos(lq+1);
    1447                 :            : 
    1448                 :        567 :   av2 = avma;
    1449                 :            : MATCH_RESTART:
    1450                 :        567 :   avma = av2;
    1451                 :        567 :   nbcard = 0;
    1452                 :        567 :   P = grp->rand(E);
    1453                 :        567 :   point = grp->pow(E,P, Mu);
    1454                 :        567 :   Pb = grp->pow(E,point, Mg);
    1455                 :        567 :   Pg = grp->pow(E,point, Mb);
    1456                 :            :   /* Precomputation for babies */
    1457                 :        567 :   pre = BSGS_pre(&diff, baby, Pb, E, grp);
    1458                 :            : 
    1459                 :            :   /*Now we compute the table of babies, this table contains only the */
    1460                 :            :   /*lifted x-coordinate of the points in order to use less memory */
    1461                 :        567 :   table = cgetg(lbaby, t_VECSMALL);
    1462                 :        567 :   av1 = avma;
    1463                 :            :   /* (p+1 - u - Mu*Mb*Sg) P - (baby[1]) Pb */
    1464                 :        567 :   point = grp->pow(E,P, subii(subii(pp1, u), mulii(Mu, addii(SgMb, mulii(Mg, gel(baby,1))))));
    1465                 :        567 :   table[1] = grp->hash(gel(point,1));
    1466         [ +  + ]:     365288 :   for (i = 2; i < lbaby; i++)
    1467                 :            :   {
    1468                 :     364721 :     GEN d = subii(gel(baby, i), gel(baby, i-1));
    1469                 :     364721 :     point =  grp->mul(E, point, grp->pow(E, gel(pre, ZV_search(diff, d)), gen_m1));
    1470                 :     364721 :     table[i] = grp->hash(gel(point,1));
    1471         [ +  + ]:     364721 :     if (gc_needed(av1,3))
    1472                 :            :     {
    1473         [ -  + ]:         22 :       if(DEBUGMEM>1) pari_warn(warnmem,"match_and_sort, baby = %ld", i);
    1474                 :         22 :       point = gerepileupto(av1, point);
    1475                 :            :     }
    1476                 :            :   }
    1477                 :        567 :   avma = av1;
    1478                 :            :   /* Precomputations for giants */
    1479                 :        567 :   pre = BSGS_pre(&diff, giant, Pg, E, grp);
    1480                 :            : 
    1481                 :            :   /* Look for a collision among the x-coordinates */
    1482                 :        567 :   table_ind = vecsmall_indexsort(table);
    1483                 :        567 :   table = perm_mul(table,table_ind);
    1484                 :            : 
    1485                 :        567 :   av1 = avma;
    1486                 :        567 :   point = grp->pow(E, Pg, gel(giant, 1));
    1487                 :        567 :   for (i = 1; ; i++)
    1488                 :            :   {
    1489                 :            :     GEN d;
    1490                 :     510678 :     long h = grp->hash(gel(point, 1));
    1491                 :     510678 :     long s = zv_search(table, h);
    1492         [ +  + ]:     510678 :     if (s) {
    1493 [ +  + ][ +  - ]:       1134 :       while (table[s] == h && s) s--;
    1494 [ +  + ][ +  + ]:       1134 :       for (s++; s < lbaby && table[s] == h; s++)
    1495                 :            :       {
    1496                 :        567 :         GEN B = gel(baby,table_ind[s]), G = gel(giant,i);
    1497                 :        567 :         GEN GMb = mulii(G, Mb), BMg = mulii(B, Mg);
    1498                 :        567 :         GEN Be = subii(subii(pp1, u), mulii(Mu, addii(SgMb, BMg)));
    1499                 :        567 :         GEN Bp = grp->pow(E,P, Be);
    1500                 :            :         /* p+1 - u - Mu (Sg Mb + GIANT Mb + BABY Mg) */
    1501         [ +  - ]:        567 :         if (gequal(gel(Bp,1),gel(point,1)))
    1502                 :            :         {
    1503                 :        567 :           GEN card1 = subii(Be, mulii(Mu, GMb));
    1504                 :        567 :           GEN card2 = addii(card1, mulii(mulsi(2,Mu), GMb));
    1505         [ -  + ]:        567 :           if (DEBUGLEVEL>=2) timer_printf(&ti,"match_and_sort");
    1506         [ +  + ]:        567 :           if (absi_cmp(subii(pp1, card1), bound) <= 0)
    1507                 :        490 :             affii(card1, gel(card, ++nbcard));
    1508         [ -  + ]:        567 :           if (nbcard >= lcard) goto MATCH_RESTART;
    1509         [ +  + ]:        567 :           if (absi_cmp(subii(pp1, card2), bound) <= 0)
    1510                 :        154 :             affii(card2, gel(card, ++nbcard));
    1511         [ -  + ]:        567 :           if (nbcard >= lcard) goto MATCH_RESTART;
    1512                 :            :         }
    1513                 :            :       }
    1514                 :            :     }
    1515         [ +  + ]:     510678 :     if (i==lgiant-1) break;
    1516                 :     510111 :     d = subii(gel(giant, i+1), gel(giant, i));
    1517                 :     510111 :     point = grp->mul(E,point, gel(pre, ZV_search(diff, d)));
    1518         [ +  + ]:     510111 :     if (gc_needed(av1,3))
    1519                 :            :     {
    1520         [ -  + ]:         22 :       if(DEBUGMEM>1) pari_warn(warnmem,"match_and_sort, giant = %ld", i);
    1521                 :         22 :       point = gerepileupto(av1, point);
    1522                 :            :     }
    1523                 :     510111 :   }
    1524                 :        567 :   setlg(card, nbcard+1);
    1525                 :        574 :   return gen_select_order(card, E, grp);
    1526                 :            : }
    1527                 :            : 
    1528                 :            : static GEN
    1529                 :        602 : get_bound_bsgs(long lp)
    1530                 :            : {
    1531                 :            :   GEN B;
    1532         [ +  + ]:        602 :   if (lp <= 160)
    1533                 :        574 :     B = divru(powru(dbltor(1.048), lp), 9);
    1534         [ +  + ]:         28 :   else if (lp <= 192)
    1535                 :         21 :     B = divrr(powru(dbltor(1.052), lp), dbltor(16.65));
    1536                 :            :   else
    1537                 :          7 :     B = mulrr(powru(dbltor(1.035), minss(lp,307)), dbltor(1.35));
    1538                 :        602 :   return mulru(B, 1000000);
    1539                 :            : }
    1540                 :            : 
    1541                 :            : /*FIXME: the name of the function does not quite match what it does*/
    1542                 :            : static const struct bb_group *
    1543                 :        574 : get_FqE_group(void ** pt_E, GEN a4, GEN a6, GEN T, GEN p)
    1544                 :            : {
    1545         [ +  + ]:        574 :   if (!T) return get_FpE_group(pt_E,a4,a6,p);
    1546         [ +  + ]:         28 :   else if (lgefint(p)==3)
    1547                 :            :   {
    1548                 :         21 :     ulong pp = uel(p,2);
    1549                 :         21 :     return get_FlxqE_group(pt_E, Fq_to_Flx(a4, T, pp), Fq_to_Flx(a6, T, pp),
    1550                 :            :                            ZXT_to_FlxT(T,pp),pp);
    1551                 :            :   }
    1552                 :        574 :   return get_FpXQE_group(pt_E,a4,a6,T,p);
    1553                 :            : }
    1554                 :            : 
    1555                 :            : /* E is an elliptic curve defined over Z or over Fp in ellinit format, defined
    1556                 :            :  * by the equation E: y^2 + a1*x*y + a2*y = x^3 + a2*x^2 + a4*x + a6
    1557                 :            :  * p is a prime number
    1558                 :            :  * set smallfact to stop whenever a small factor > smallfact of the order is
    1559                 :            :  * detected. Useful when searching for a good curve for cryptographic
    1560                 :            :  * applications */
    1561                 :            : GEN
    1562                 :        602 : Fq_ellcard_SEA(GEN a4, GEN a6, GEN q, GEN T, GEN p, long smallfact)
    1563                 :            : {
    1564                 :        602 :   const long MAX_ATKIN = 21;
    1565                 :        602 :   pari_sp ltop = avma, btop;
    1566                 :            :   long ell, i, nb_atkin, vx,vy;
    1567                 :            :   GEN TR, TR_mod, compile_atkin, bound, bound_bsgs, champ;
    1568                 :        602 :   GEN prod_atkin = gen_1, max_traces = gen_0;
    1569                 :            :   GEN j;
    1570                 :        602 :   double bound_gr = 1.;
    1571                 :        602 :   const double growth_factor = 1.26;
    1572                 :            :   forprime_t TT;
    1573                 :            :   void *E;
    1574                 :            : 
    1575 [ +  + ][ -  + ]:        602 :   if (!modular_eqn && !get_seadata(0)) return NULL;
    1576                 :        602 :   j = Fq_ellj(a4, a6, T, p);
    1577 [ +  + ][ +  - ]:        602 :   if (typ(j)==t_INT && (signe(j)==0 || equaliu(j, umodui(1728, p))))
                 [ -  + ]
    1578                 :          0 :     return T ? FpXQ_ellcard(Fq_to_FpXQ(a4, T, p), Fq_to_FpXQ(a6, T, p), T, p)
    1579         [ #  # ]:          0 :              : Fp_ellcard(a4, a6, p);
    1580                 :            :   /*First compute the trace modulo 2 */
    1581      [ +  +  + ]:        602 :   switch(FqX_nbroots(mkpoln(4, gen_1, gen_0, a4, a6), T, p))
    1582                 :            :   {
    1583                 :            :   case 3: /* bonus time: 4 | #E(Fq) = q+1 - t */
    1584         [ +  + ]:         56 :     i = mod4(q)+1; if (i > 2) i -= 4;
    1585                 :         56 :     TR_mod = utoipos(4);
    1586                 :         56 :     TR = stoi(i); break;
    1587                 :            :   case 1:
    1588                 :        315 :     TR_mod = gen_2;
    1589                 :        315 :     TR = gen_0; break;
    1590                 :            :   default : /* 0 */
    1591                 :        231 :     TR_mod = gen_2;
    1592                 :        231 :     TR = gen_1; break;
    1593                 :            :   }
    1594 [ -  + ][ #  # ]:        602 :   if (smallfact == 1 && !mpodd(TR))
    1595                 :            :   {
    1596         [ #  # ]:          0 :     if (DEBUGLEVEL) err_printf("Aborting: #E(Fq) divisible by 2\n");
    1597                 :          0 :     avma = ltop; return gen_0;
    1598                 :            :   }
    1599                 :        602 :   vy = fetch_var();
    1600                 :        602 :   vx = fetch_var_higher();
    1601                 :            : 
    1602                 :            :   /* compile_atkin is a vector containing informations about Atkin primes,
    1603                 :            :    * informations about Elkies primes lie in Mod(TR, TR_mod). */
    1604                 :        602 :   u_forprime_init(&TT, 3, 1000); /* way beyond what seadata provides */
    1605                 :        602 :   bound = sqrti(shifti(q, 4));
    1606                 :        602 :   bound_bsgs = get_bound_bsgs(expi(q));
    1607                 :        602 :   compile_atkin = zerovec(MAX_ATKIN); nb_atkin = 0;
    1608                 :        602 :   btop = avma;
    1609         [ +  - ]:       6209 :   while ( (ell = u_forprime_next(&TT)) )
    1610                 :            :   {
    1611                 :       6209 :     long ellkt, kt = 1, nbtrace;
    1612                 :            :     GEN trace_mod;
    1613         [ +  + ]:       6209 :     if (equalui(ell, p)) continue;
    1614                 :       6202 :     trace_mod = find_trace(a4, a6, j, ell, q, T, p, &kt, smallfact, vx,vy);
    1615         [ +  + ]:       6202 :     if (!trace_mod) continue;
    1616                 :            : 
    1617                 :       6174 :     nbtrace = lg(trace_mod) - 1;
    1618                 :       6174 :     ellkt = (long)upowuu(ell, kt);
    1619         [ +  + ]:       6174 :     if (nbtrace == 1)
    1620                 :            :     {
    1621                 :       3752 :       long t_mod_ellkt = trace_mod[1];
    1622 [ -  + ][ #  # ]:       3752 :       if (smallfact && ell > smallfact)
    1623                 :            :       { /* does ell divide q + 1 - t ? */
    1624                 :          0 :         long card_mod_ell = (umodiu(q,ell) + 1 - t_mod_ellkt) % ell ;
    1625         [ #  # ]:          0 :         if (!card_mod_ell)
    1626                 :            :         {
    1627         [ #  # ]:          0 :           if (DEBUGLEVEL)
    1628                 :          0 :             err_printf("\nAborting: #E(Fq) divisible by %ld\n",ell);
    1629                 :          0 :           avma = ltop; return gen_0;
    1630                 :            :         }
    1631                 :            :       }
    1632                 :       3752 :       (void)Z_incremental_CRT(&TR, t_mod_ellkt, &TR_mod, ellkt);
    1633                 :            :     }
    1634                 :            :     else
    1635                 :            :     {
    1636                 :       2422 :       add_atkin(compile_atkin, mkvec2(utoipos(ellkt), trace_mod), &nb_atkin);
    1637                 :       2422 :       prod_atkin = value(-1, compile_atkin, nb_atkin);
    1638                 :            :     }
    1639         [ +  + ]:       6174 :     if (cmpii(mulii(TR_mod, prod_atkin), bound) > 0)
    1640                 :            :     {
    1641                 :            :       GEN bound_tr;
    1642         [ +  + ]:        644 :       if (!nb_atkin) return gerepileuptoint(ltop, subii(addis(q, 1), TR));
    1643                 :        616 :       bound_tr = mulrr(bound_bsgs, dbltor(bound_gr));
    1644                 :        616 :       bound_gr *= growth_factor;
    1645         [ +  + ]:        616 :       if (signe(max_traces))
    1646                 :            :       {
    1647                 :         42 :         max_traces = divis(muliu(max_traces,nbtrace), ellkt);
    1648         [ -  + ]:         42 :         if (DEBUGLEVEL>=3)
    1649                 :          0 :           err_printf("At least %Ps remaining possibilities.\n",max_traces);
    1650                 :            :       }
    1651         [ +  + ]:        616 :       if (cmpir(max_traces, bound_tr) < 0)
    1652                 :            :       {
    1653                 :        602 :         GEN bound_atkin = truedivii(bound, TR_mod);
    1654                 :        602 :         champ = champion(compile_atkin, nb_atkin, bound_atkin);
    1655                 :        602 :         max_traces = gel(champ,2);
    1656         [ -  + ]:        602 :         if (DEBUGLEVEL>=2)
    1657                 :          0 :           err_printf("%Ps remaining possibilities.\n", max_traces);
    1658         [ +  + ]:        602 :         if (cmpir(max_traces, bound_tr) < 0)
    1659                 :            :         {
    1660                 :        574 :           GEN res, cat = shallowextract(compile_atkin, gel(champ,1));
    1661                 :            :           const struct bb_group *grp;
    1662         [ -  + ]:        574 :           if (DEBUGLEVEL)
    1663                 :          0 :             err_printf("Match and sort for %Ps possibilities.\n", max_traces);
    1664                 :        574 :           grp = get_FqE_group(&E,a4,a6,T,p);
    1665                 :        574 :           res = match_and_sort(cat, TR_mod, TR, q, E, grp);
    1666                 :        574 :           return gerepileuptoint(ltop, res);
    1667                 :            :         }
    1668                 :            :       }
    1669                 :            :     }
    1670         [ -  + ]:       5572 :     if (gc_needed(btop, 1))
    1671                 :       6209 :       gerepileall(btop,5, &TR,&TR_mod, &compile_atkin, &max_traces, &prod_atkin);
    1672                 :            :   }
    1673                 :        602 :   return NULL;/*not reached*/
    1674                 :            : }
    1675                 :            : 
    1676                 :            : GEN
    1677                 :        567 : Fp_ellcard_SEA(GEN a4, GEN a6, GEN p, long smallfact)
    1678                 :            : {
    1679                 :        567 :   return Fq_ellcard_SEA(a4, a6, p, NULL, p, smallfact);
    1680                 :            : }
    1681                 :            : 
    1682                 :            : GEN
    1683                 :          0 : ellsea(GEN E, GEN p, long smallfact)
    1684                 :            : {
    1685                 :          0 :   pari_sp av = avma;
    1686                 :          0 :   GEN a4 = modii(mulis(Rg_to_Fp(gel(E,10), p), -27), p);
    1687                 :          0 :   GEN a6 = modii(mulis(Rg_to_Fp(gel(E,11), p), -54), p);
    1688                 :          0 :   GEN card = Fp_ellcard_SEA(a4, a6, p, smallfact);
    1689         [ #  # ]:          0 :   if (!card) pari_err_PACKAGE("seadata");
    1690                 :          0 :   return gerepileuptoint(av, subii(addis(p,1),card));
    1691                 :            : }

Generated by: LCOV version 1.9