Benyamin Gholami on Sun, 26 Mar 2017 17:11:26 +0200


[Date Prev] [Date Next] [Thread Prev] [Thread Next] [Date Index] [Thread Index]

details


Sure Sir !
here you are:
note that i created a rank 1 surface with torsion Z/2Z*Z/4Z  and i'm searching on this family for curves with rank 6 ( or fine if greater) with codes that you gave me before. my problem is in most of cases it just give "1" as lower bound and the selmer bounds are not indeed accurate. i upload denis simon script and then use sieving. any helps would be appreciated .thanks to you and Prof.Cremona (mr bellabas helped me to write S(E,N) , so thank to him so)


























####################################################################################################################################################################
                                               Denis simon Gp script
##########################################################################################################################################################################################################################################################








\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\       Copyright (C) 2017 Denis Simon
\\
\\ Distributed under the terms of the GNU General Public License (GPL)
\\
\\    This code is distributed in the hope that it will be useful,
\\    but WITHOUT ANY WARRANTY; without even the implied warranty of
\\    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
\\    General Public License for more details.
\\
\\ The full text of the GPL is available at:
\\
\\                 http://www.gnu.org/licenses/
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

/*
  Auteur :
  Denis SIMON -> denis.simon@unicaen.fr
  adresse du fichier :
  www.math.unicaen.fr/~simon/ellQ.gp

  *********************************************
  *          VERSION 10/01/2017               *
  *********************************************


  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
  \\             English                      \\
  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

  This package provides functions to compute the rank of elliptic
  curves over Q using 2-descent.
  This package requires a version of gp >= 2.9
  This package requires the other package ellcommon.gp downloadable at
  www.math.unicaen.fr/~simon/ellcommon.gp

  They can be run under GP by the commands 
  gp > \r ellcommon.gp
  gp > \r ellQ.gp

  The main function is ellrank(), which takes as an argument
  any elliptic curve in the form [a1,a2,a3,a4,a6]
  the result is a vector [r,s,v], where
    r is a lower bound for the rank,
    s is the rank of the 2-Selmer group
    v is a set of independant points in E(Q)/2E(Q).

  Example:

  gp > ell = [1,2,3,4,5];
  gp > ellrank(ell)
  %1 = [1, 1, [[1,2]]
  In this example, the rank is exactly 1, and [1,2] has infinite order.

  more details on the computations may be obtained by setting
  DEBUGLEVEL_ell = 1 (the higher value, the more details)

  Other functions: 

  ell2descent_complete, ell2descent_twist, ell2descent_viaisog,
  ellhalf, ellredgen, ellsort, elltors2, elltorseven,
  locallysoluble, ratpoint, redquartic,
  bnfpSelmer, reducemodsquares

  Quick information is obtained by typing
  gp > ?NameOfTheFunction

  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
  \\    Comment utiliser ce programme ?       \\
  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

  Ce module contient des fonctions pour calculer le rang des courbes
  elliptiques sur Q, en utilisant la 2-descente.
  langage : GP

  Ce module utilise le module complementaire suivant :
    ellcommon.gp
  qui est disponibles a l'adresse
    www.math.unicaen.fr/~simon/ellcommon.gp

  Pour l'utiliser, lancer gp, puis taper
  \r ellcommon.gp
  \r ellQ.gp

  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
  \\  Description des principales fonctions   \\
  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

  Explications succintes :
  La fonction ellrank() accepte toutes les courbes sous la forme
  [a1,a2,a3,a4,a6]
  Les coefficients peuvent etre entiers ou non.
  L'algorithme utilise est celui de la 2-descente.
  La 2-torsion peut etre quelconque.
  Il suffit de taper : 

  gp > ell = [a1,a2,a3,a4,a6];
  gp > ellrank(ell)

  Retourne un vecteur [r,s,v] ou
    r est le rang probable (c'est toujours une minoration du rang),
    s est le 2-rang du groupe de Selmer,
    v est une liste de points independants dans E(Q)/2E(Q).

  Exemple :

  gp > ell = [1,2,3,4,5];
  gp > ellrank(ell)
  %1 = [1, 1, [[1,2]]
  Ici, le rang est exactement 1, et le point [1,2] est d'ordre infini.

  Courbes de la forme : k*y^2 = x^3+A*x^2+B*x+C
  sans 2-torsion, A,B,C entiers.
  gp > bnf = bnfinit(x^3+A*x^2+B*x+C);
  gp > ell = ellinit([0,A,0,B,C],1);
  gp > rank = ell2descent_twist(ell,bnf,k);

  Courbes avec #E[2](Q) >= 2 :
  ell doit etre sous la forme 
  y^2 = x^3 + A*x^2 + B*x
  avec A et B entiers.
  gp > ell = [0,A,0,B,0]
  gp > ell2descent_viaisog(ell)
  = algorithme de la 2-descente par isogenies
  Attention A et B doivent etre entiers

  Courbes avec #E[2](Q) = 4 : y^2 = (x-e1)*(x-e2)*(x-e3)
  gp > ell2descent_complete(e1,e2,e3)
  = algorithme de la 2-descente complete
  Attention : les ei doivent etre entiers.

  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
  \\     Autres fonctions utiles              \\
  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

  elltors2(E)      : determine le groupe E[2](Q)
  elltorseven(E)   : determine le groupe E[2^*](Q)
  ellhalf(E,P)     : liste les points Q tels que 2Q = P
  ellredgen(E,v)   : reduction des points de v sur E

  locallysoluble(pol): teste si y^2=pol(x) est ELS
  ratpoint(pol,lim): cherche un point sur y^2=pol(x)
  redquartic(pol): reduction de la quartique pol


  Aide en ligne :
  ?NomDeLaFonction


  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
  \\       Affichage des calculs              \\
  \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

  On peut avoir plus ou moins de details de calculs avec
  DEBUGLEVEL_ell = 0;
  DEBUGLEVEL_ell = 1; 2; 3;...

*/

{
\\
\\ Usual global variables
\\

global(DEBUGLEVEL_ell, LIM1, LIM3, LIMTRIV, ELLREDGENFLAG, COMPLETE):small;

  DEBUGLEVEL_ell = 0; \\ From 0 to 5: choose a higher value to have
                      \\ more details printed.
  LIM1 = 5;           \\ Limit for the search of trivial points on quartics
  LIM3 = 50;          \\ Limit for the search of points on ELS quartics
  LIMTRIV = 3;        \\ Limit for the search of trivial points on the elliptic curve
  ELLREDGENFLAG = 1;  \\ to reduce the generators at the end
  COMPLETE = 0;       \\ Use Complete 2-descent when full 2-torsion,
                      \\ otherwise 2-descent via isogenies.

\\
\\  Technical global variables
\\

global(MAXPROB, LIMBIGPRIME):small;

  MAXPROB = 20;
  LIMBIGPRIME = 30; \\ for primes larger than this limit: use a probabilistic test
                    \\ LIMBIGPRIME = 0 means: only deterministic tests
}

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\          SCRIPT                             \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\    MANIPULATION OF GLOBAL VARIABLES         \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{default_ellQ(
  DEBUGLEVEL_ell_val:small = 0,
  LIM1_val:small = 5,
  LIM3_val:small = 50,
  LIMTRIV_val:small = 3,
  ELLREDGENFLAG_val:small = 1,
  COMPLETE_val:small = 0,
  MAXPROB_val:small = 20,
  LIMBIGPRIME_val:small = 30
  ) = 

  DEBUGLEVEL_ell = DEBUGLEVEL_ell_val;
  print("  DEBUGLEVEL_ell = ",DEBUGLEVEL_ell);

  LIM1 = LIM1_val;
  print("  LIM1 = ",LIM1);

  LIM3 = LIM3_val;
  print("  LIM3 = ",LIM3);

  LIMTRIV = LIMTRIV_val;
  print("  LIMTRIV = ",LIMTRIV);

  ELLREDGENFLAG = ELLREDGENFLAG_val;
  print("  ELLREDGENFLAG = ",ELLREDGENFLAG);

  COMPLETE = COMPLETE_val;
  print("  COMPLETE = ",COMPLETE);

  MAXPROB = MAXPROB_val;
  print("  MAXPROB = ",MAXPROB);

  LIMBIGPRIME = LIMBIGPRIME_val;
  print("  LIMBIGPRIME = ",LIMBIGPRIME);
}

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\    FUNCTIONS FOR POLYNOMIALS                \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{ratpoint( pol, lim:small = 1, singlepoint = 1, infinity = 1) =
\\ Search for points on y^2=pol(x).
\\ The coeff of pol must be integers.
\\ If singlepoint >= 1, stop after a first point is found.

my(listpoints,point1,odd,deg4,pol16,tab16,pol9,tab9,pol5,tab5,pol0,vecz,vecx,lead,zz,xx,evpol,iz,factpol,deg,vz,epsmax);

if( DEBUGLEVEL_ell >= 4,
  print("    Starting ratpoint with pol = ",pol); 
  print("    lim = ",lim););

  if( !singlepoint, listpoints = []);
  point1 = [];

\\
\\          trivial solutions
\\

\\ the leading coeff is a square
  if( infinity && issquare(pollead(pol)),
    point1 = [ 1, sqrtrat(pollead(pol)), 0];
if( DEBUGLEVEL_ell >= 3, print("   trivial solution: lead(pol) is a square"));
    if( singlepoint,
if( DEBUGLEVEL_ell >= 4, print("    end of ratpoint"));
      return(point1));
    listpoints = concat(listpoints,[point1]));

\\ the constant coeff is a square
  if( issquare(polcoeff(pol,0)),
    point1 = [ 0, sqrtrat(polcoeff(pol,0)) ];
if( DEBUGLEVEL_ell >= 3, print("   trivial solution: pol(0) is a square"));
    if( singlepoint,
if( DEBUGLEVEL_ell >= 4, print("    end of ratpoint"));
      return(point1));
    listpoints = concat(listpoints,[point1]));

\\ roots of pol ?
  factpol = nfroots(,pol);
  if( #factpol, 
if( DEBUGLEVEL_ell >= 3, print("   trivial solution: roots of pol",factpol));
    if( singlepoint,
if( DEBUGLEVEL_ell >= 4, print("    end of ratpoint"));
      return([factpol[1],0]));
    listpoints = concat(listpoints,vector(#factpol,i,[factpol[i],0]))
  );

\\
\\       Sieve
\\

\\ initialisation of the sieve modulo 16, 9 and 5
\\ used only with even degree when lim is large

  deg = poldegree(pol);
  odd = deg%2;
  deg4 = ( !odd && lim > 20);
  if( deg4,

    pol16 = (Vec(pol)*Mod(1,16))~;
    tab16 = matrix(16,16);
    for(xx = 0, 16-1, 
      for(zz = 0, 16-1,
        tab16[xx+1,zz+1] = !issquare(vector(deg+1,i,xx^(deg+1-i)*zz^(i-1))*pol16)));
    pol9 = (Vec(pol)~)*Mod(1,9);
    tab9 = matrix(9,9);
    for(xx = 0, 9-1, 
      for(zz = 0, 9-1,
        tab9[xx+1,zz+1] = !issquare(vector(deg+1,i,xx^(deg+1-i)*zz^(i-1))*pol9)));
    pol5 = (Vec(pol)~)*Mod(1,5);
    tab5 = matrix(5,5);
    for(xx = 0, 5-1, 
      for(zz = 0, 5-1,
        tab5[xx+1,zz+1] = !issquare(vector(deg+1,i,xx^(deg+1-i)*zz^(i-1))*pol5)));
  );

\\ if the degree is odd, search only for square denominators
  if( odd, 
    vecz = vector(lim,i,i^2);
  ,
\\ if the degree is even, the leading coeff must be
\\ a square modulo zz.
    lead = pollead(pol);
    vecz = vector(lim);
    zz = 0;
    for( i = 1, lim,
      zz++; while( !issquare(Mod(lead,zz)),zz++); vecz[i] = zz
  ));

\\ the constant coeff must be a square modulo xx.
  pol0 = polcoeff(pol,0);
  vecx = vector(lim);
  xx = 0;
  for( i = 1, lim,
    xx++; while( !issquare(Mod(pol0,xx)),xx++); vecx[i] = xx);

if( DEBUGLEVEL_ell >= 4, print("    xmax = ",vecx[lim]));
if( DEBUGLEVEL_ell >= 4, print("    zmax = ",vecz[lim]));

if( DEBUGLEVEL_ell >= 5, print("    vecx = ",vecx));
if( DEBUGLEVEL_ell >= 5, print("    vecz = ",vecz));

  if( deg4,
    vz = vector(lim,i,Pol(
      vector(deg+1,j,polcoeff(pol,deg+1-j)*vecz[i]^(j-1))));
  );

\\ if deg is even and pol has no odd terms,
\\ it is enough to consider xx > 0
  if( !odd && subst(pol,x,-x)==pol, epsmax = 1, epsmax = 2);

\\ loop over x = xx/zz
\\ the loop on [xx,zz] is done diagonally
\\ to start with the smallest values of both xx and zz.
  for( somme = 2, 2*lim,
    for( ix = max(1,somme-lim), min(lim,somme-1),
      xx = vecx[ix]; iz = somme-ix; zz = vecz[iz];
      if( gcd(zz,xx) != 1, next);
      for( eps = 1, epsmax,
\\ when eps = 1, xx > 0; when eps = 2, xx < 0.
        if( deg4,
          if( tab16[xx%16+1,zz%16+1], xx=-xx;next);
          if( tab9[xx%9+1,zz%9+1],    xx=-xx;next);
          if( tab5[xx%5+1,zz%5+1],    xx=-xx;next);
          evpol = subst(vz[iz],'x,xx)
        ,
          evpol = subst(pol,variable(pol),xx/zz)
        );
        if( issquare(evpol),
          point1 = [xx/zz,sqrtrat(subst(pol,variable(pol),xx/zz))];
if( DEBUGLEVEL_ell >= 4, print("    point found by ratpoint = ",point1));
          if( singlepoint, break(3));
          listpoints = concat(listpoints,[point1])
        );
        xx = -xx
  )));

  if( point1 != [], 
if( DEBUGLEVEL_ell >= 3, print("   point found by ratpoint = ",point1));
if( DEBUGLEVEL_ell >= 4, print("    end of ratpoint "));
    if( singlepoint, return(point1), return(listpoints))
  );

return([]);
}
{redquartic(pol) =
\\ reduction of the quartic polynomial.
\\ (also works with other degrees)
my(mylocalprec,prec0,d,disc2,test,r,normderiv,disc2v,q,M);

if( DEBUGLEVEL_ell >= 4, print("    starting redquartic"));
if( DEBUGLEVEL_ell >= 3, print("   reduction of the quartic ",pol));

\\ choice of the real precision used in the computation
  mylocalprec = prec0 = default(realprecision);
  d = poldegree(pol);
  disc2 = poldisc(pol)^2;
  test = 0;
  while( test == 0,
if( DEBUGLEVEL_ell >= 4, print("    precision = ",mylocalprec));
    r = polroots(pol);
    normderiv = vector( d, i, norm(subst(pol',variable(pol),r[i])));
    disc2v = prod( i = 1, d, normderiv[i]) * pollead(pol)^(2*d-4);
    test = abs(disc2v-disc2) < 10^(-mylocalprec\2);
    if( !test, default(realprecision, mylocalprec *= 2))
  );

\\ former choice of the quadratic form
\\  q = Vec(sum( i = 1, d, norm(x-r[i])));
\\ Now, uses the quadratic form normalized as in Cremona-Stoll
  q = Vec(sum( i = 1, d, norm('x-r[i]) / normderiv[i]^(1/(d-2))));
  M = qfbreduce([q[1],q[2]/2;q[2]/2,q[3]]);
  pol = subst(pol,variable(pol),Pol(M[1,])/Pol(M[2,]))*Pol(M[2,])^poldegree(pol);

  if( mylocalprec != prec0, default(realprecision,prec0));

if( DEBUGLEVEL_ell >= 3, print("   reduced quartic = ",pol));
if( DEBUGLEVEL_ell >= 4, print("    end of redquartic"));

  return([pol,M]);
}
{listratpoint( pol, redflag = 0) =
my(list,i,K,ff,C,p,M,U,newpol,factpol,ll,listf,rr);

if( DEBUGLEVEL_ell >= 5, print("     Starting listratpoint with ",pol));
  list = [[pol,matid(2),1,1]];
  i = 1; 
  while( i <= #list,

    pol = list[i][1];

    K = abs(content(pol));
    if( K != 1,
      pol = (list[i][1] /= K);
      list[i][3] *= K
    );

    K = list[i][3];
    if( K == 1, i++; next);

    ff = factor(K);
    if( vecmax(ff[,2]) > 1,
      ff[,2] \= 2;
      C = factorback(ff);
      list[i][4] *= C;
      K = ( list[i][3] /= C^2);
      if( K == 1, i++; next);
      ff = factor(K);
    );

    p = ff[1,1];
    M = list[i][2];
    C = list[i][4];

    if( pollead(pol)%p == 0,
      U = M*[1,0;0,p];
      if( content(U) == 1,
        newpol = subst(pol,'x,'x/p)*p^(poldegree(pol)-1);
        list = concat(list, [[newpol,U,K/p,C*p]])
      )
    );

    factpol = centerlift(polrootsmod(pol,p));
    for( j = 1, #factpol,
      U = M*[p,factpol[j];0,1];
      if( content(U) == 1,
        newpol = subst(pol,'x,p*'x+factpol[j])/p;
        list = concat(list, [[newpol,U,K/p,C*p]])
    ));

    i++;
  );

  ll = sum( i = 1, #list, list[i][3] == 1);
  listf = vector(ll);
  i = 1;
  for( j = 1, #list,
    if( list[j][3] == 1,
      listf[i] = list[j]; i++));

  if( redflag,
    for( i = 1, #listf,
      rr = redquartic(listf[i][1]);
      listf[i][1] = rr[1];
      listf[i][2] = listf[i][2]*rr[2]
    )
  );

if( DEBUGLEVEL_ell >= 5, print("     Output of listratpoint = ",listf));

return(listf);
}
{ratpoint2( pol, lim:small = 1, singlepoint = 1, redflag = 0) =
my(listpoints,list,rr,y2,aux);

  listpoints = [];
  list = listratpoint(pol,redflag);
  for( i = 1, #list,
    rr = ratpoint(list[i][1],lim,singlepoint);
    if( singlepoint && #rr, rr=[rr]);
    for( j = 1, #rr,
      y2 = rr[j][2]*list[i][4];
      if( #rr[j] == 2,
        aux = [rr[j][1],1]~
      , aux = [rr[j][1],rr[j][3]]~
      );
      aux = list[i][2] * aux;
      if( aux[2] == 0,
        rr[j] = [aux[1],y2,0]
      , rr[j] = [aux[1]/aux[2],y2/aux[2]^(poldegree(pol)\2)]
      );
    );
    if( singlepoint && #rr, return(rr[1]));
    listpoints = concat(listpoints,rr);
  );
  listpoints = vecsort(listpoints,,2);
return(listpoints);
}
{polrealrootsisolate(pol) =
\\ pol is a squarefree polynomial in Z[x].
\\ Returns a list of vectors [a,b] with a and b rationals
\\ such that the intervals ]a,b] are disjoints and contain
\\ all the real roots of pol, and excatly one in each interval.
my(st,a,res,ind,b,c,stab,stac);

if( DEBUGLEVEL_ell >= 5, print("     starting polrealrootsisolate with pol = ",pol));
  st = polsturm(pol);
  if( !st, return([]));
  a = 1;
  while( polsturm(pol,[-a,a]) < st, a <<= 1);
  res = [[-a,a,st]];
  ind = 1;
  while( #res < st,
    while( res[ind][3] == 1, ind++);
    a = res[ind][1]; b = res[ind][2]; stab = res[ind][3];
    c = (a+b)/2;
    stac = polsturm(pol,[a,c]);
    if( stac == 0, res[ind][1] = c; next);
    if( stac == stab, res[ind][2] = c; next);
    res[ind] = [a,c,stac];
    res = concat(res,[[c,b,stab-stac]]);
  );
  res = vector(st,i,[res[i][1],res[i][2]]);
  res = vecsort(res,1);
if( DEBUGLEVEL_ell >= 5, print("     end of polrealrootsisolate with res = ",res));
  return(res);
}
{polrealrootsimprove(pol,v) =
\\ pol is a polynomial and v is a vector v=[a,b]
\\ such that pol contains exactly one root in the interval ]a,b].
\\ Returns another interval with the same property, but with half length.
\\ (dichotomy)
my(c,v2,vc);

  c = (v[1]+v[2])/2;
  v2 = subst(pol,variable(pol),v[2]);
  if( v2 == 0, return([c,v[2]]));
  vc = subst(pol,variable(pol),c);
  if( sign(v2)*sign(vc) >= 0, v[2] = c, v[1] = c);
  return(v);
}
{polrootsmodpn(pol,p) =
\\ Compute a list v. Each element of v is of the form
\\ [a,b], with maximal b <= valuation(poldisc(pol),p)
\\ and a is a root of pol modulo p^b.
my(vd,rac,i,pol2,r,newrac);

if( DEBUGLEVEL_ell >= 5, print("     starting polrootsmodpn ",p,": ",pol));

  vd = valuation(poldisc(pol),p);
  rac = [[0,0]];
  i = 1;
  while (i <= #rac,
\\    if( rac[i][2] > vd, i++; next);
    if( rac[i][2] >= vd, i++; next);
    pol2 = subst(pol,'x,rac[i][1]+'x*p^rac[i][2]);
    pol2 /= content(pol2);
    r = lift(polrootsmod(pol2,p));
    if( #r == 0, i++; next);
    newrac = vector(#r,j,[rac[i][1] + p^rac[i][2]*r[j],rac[i][2]+1]);
    rac = concat(rac, vector(#r-1,j,newrac[j+1]));
    rac[i] = newrac[1];
  );
if( DEBUGLEVEL_ell >= 5, print("     end of polrootsmodpn ",rac));
  return(rac);
}

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\    FUNCTIONS FOR LOCAL COMPUTATIONS         \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{ppinit( nf, p) =
\\ a little more structure than idealprimedec()
my(pdec,pp);

  pdec = idealprimedec(nf,p);
  pp = vector(#pdec,i,
    [ pdec[i]
    , nfbasistoalg(nf,pdec[i][2])
    , if( p == 2, idealstar(nf,idealpow(nf,pdec[i],1+2*pdec[i].e)),0)
    , nfmodprinit(nf,pdec[i])
    ]);
  return(pp);
}
{nfpsquareodd( nf, a, pr) =
\\ pr is a prime ideal of nf as output by nfmodprinit
\\ a is an element of nf.
\\ Returns 1 if a is a square in the p-adics, 0 otherwise
\\ works only for p odd.
my(p,v,ap,den,norme);

if( DEBUGLEVEL_ell >= 5, print("     starting nfpsquareodd(",a,",",pr));
  if( a == 0,
if( DEBUGLEVEL_ell >= 5, print("     end of nfpsquareodd"));
    return(1));

  p = pr[3];
  v = idealval(nf,lift(a),p);
  if( v%2,
if( DEBUGLEVEL_ell >= 5, print("     end of nfpsquareodd"));
    return(0));
  ap = nfalgtobasis(nf,a/nfbasistoalg(nf,p[2])^v);
  den = valuation(denominator(content(ap)),p.p);
  if( den,
    den += den%2;
    ap = p.p^den*nfeltmul(nf,ap,nfeltpow(nf,p[2],-den*p.e))
  );

  norme = (p.p^p.f-1)/2;
  ap = nfeltpowmodpr(nf,ap,norme,pr);
  ap[1] -= 1;
  if( ap == 0,
if( DEBUGLEVEL_ell >= 5, print("     end of nfpsquareodd"));
    return(1));
  if( idealval(nf,ap,p) > 0,
if( DEBUGLEVEL_ell >= 5, print("     end of nfpsquareodd"));
    return(1));
if( DEBUGLEVEL_ell >= 5, print("     end of nfpsquareodd"));
  return(0);
}
{psquare( a, p) =
\\ a is an integer.
\\ p is a prime integer.
\\ Returns 1 if a is a square in the p-adics, 0 otherwise.
my(v,ap);

if( DEBUGLEVEL_ell >= 5, print("     starting psquare ",[a,p]));

  if( a == 0, 
if( DEBUGLEVEL_ell >= 5, print("     end of psquare 1"));
    return(1));

  v = valuation(a,p);
  if( v%2,
if( DEBUGLEVEL_ell >= 5, print("     end of psquare 0"));
    return(0));
  if( p == 2, 
    ap = (a>>v)%8-1,
    ap = kronecker(a/p^v,p)-1
  );
if( DEBUGLEVEL_ell >= 5, print("     end of psquare ", !ap));
  return(!ap);
}
{lemma6(pol, p, nu, xx) =
\\ technical lemma for local solubility of quartics
\\ Only for p <> 2.
my(gx,gpx,lambda,mu);

  gx = subst( pol, variable(pol), xx);
  if( psquare(gx,p), return(1));
  gpx = subst( pol', variable(pol), xx);
  lambda = valuation(gx,p); mu = valuation(gpx,p);

  if( (mu<+oo) && (lambda > 2*mu), return(1));
\\  if( (lambda >= mu+nu) && (nu > mu), return(1));
  if( (lambda >= 2*nu) && (mu >= nu), return(0));
  return(-1);
}
{lemma7( pol, nu, xx) =
\\ technical lemma for local solubility of quartics
\\ at p = 2.
my(gx,gpx,lambda,mu,q);

  gx = subst( pol, variable(pol), xx);
  if( psquare(gx,2), return(1));
  gpx = subst( pol', variable(pol), xx);
  lambda = valuation(gx,2); mu = valuation(gpx,2);
  if( (mu<+oo) && (lambda > 2*mu), return(1));
  if( nu > mu,
    if( lambda%2, return(-1));
    q = mu+nu-lambda;
    if( q == 1, return(1));
    if( q == 2 && (gx>>lambda)%4 == 1, return(1));
    return(-1));
  q = lambda-2*nu;
  if( q >= 0, return(0));
  if( q == -2 && (gx>>lambda)%4 == 1, return(0));
  return(-1);
}
{zp_soluble(pol, p, nu, pnu, x0, pnup) =
my(result,pol2,fact,x1);

if( DEBUGLEVEL_ell >= 5, print("     starting zp_soluble ",[pol,p,x0,nu]));
  if( p == 2,
    result = lemma7(pol,nu,x0),
    result = lemma6(pol,p,nu,x0));
  if( result == +1,
if( DEBUGLEVEL_ell >= 5, print("     end of zp_soluble 1 lemma"));
    return(1));
  if( result == -1,
if( DEBUGLEVEL_ell >= 5, print("     end of zp_soluble 0 lemma"));
    return(0));
  pnup = pnu*p;
  nu++;
  if( p < LIMBIGPRIME || !LIMBIGPRIME,
    for( i = 0, p-1,
      if( zp_soluble(pol,p,nu,pnup,x0+pnu*i),
if( DEBUGLEVEL_ell >= 5, print("     end of zp_soluble"));
        return(1)))
  ,
    pol2 = subst(pol,variable(pol),x0+pnu*variable(pol));
    pol2 /= content(pol2);
    pol2 = pol2*Mod(1,p);
    if( !poldegree(pol2), return(0));
    fact = factormod(pol2,p)[,1];
    for( i = 1, #fact,
      x1 = -centerlift(polcoeff(fact[i],0));
      if( zp_soluble(pol,p,nu,pnup,x0+pnu*x1),
if( DEBUGLEVEL_ell >= 5, print("     end of zp_soluble"));
        return(1)));
    for( i = 1, MAXPROB,
      x1 = random(p);
      if( zp_soluble(pol,p,nu,pnup,x0+pnu*x1),
if( DEBUGLEVEL_ell >= 5, print("     end of zp_soluble"));
        return(1)))
  );
if( DEBUGLEVEL_ell >= 2,
  if( p >= LIMBIGPRIME, 
    print("  ******* probabilistic test at p = ",p,"*******")));
if( DEBUGLEVEL_ell >= 5, print("     end of zp_soluble"));
  return(0);
}
{qp_soluble(pol, p) = 
if( DEBUGLEVEL_ell >= 5, 
  print("     starting qp_soluble ",p);
  print("     pol = ",pol));
  if( psquare(pollead(pol),p),
if( DEBUGLEVEL_ell >= 5, print("     end of qp_soluble 1"));
    return(1));
  if( psquare(polcoeff(pol,0),p),
if( DEBUGLEVEL_ell >= 5, print("     end of qp_soluble 1"));
    return(1));
  if( zp_soluble(pol,p,0,1,0),
if( DEBUGLEVEL_ell >= 5, print("     end of qp_soluble 1"));
    return(1));
  if( zp_soluble(polrecip(pol),p,1,p,0),
if( DEBUGLEVEL_ell >= 5, print("     end of qp_soluble 1"));
    return(1));
if( DEBUGLEVEL_ell >= 5, print("     end of qp_soluble 0"));
  return(0);
}
{locallysoluble(pol) =
\\ Determines if y^2 = pol(x,z) is everywhere locally soluble
my(c,disc0,plist,p,vc);

if( DEBUGLEVEL_ell >= 4, print("    starting locallysoluble: ",pol));

\\ real place
  if( !(poldegree(pol)%2) && sign(pollead(pol)) < 0 
         && sign(polcoeff(pol,0)) < 0 && polsturm(pol) == 0,
if( DEBUGLEVEL_ell >= 3, print("   not ELS at infinity"));
if( DEBUGLEVEL_ell >= 4, print("    end of locallysoluble"));
    return(0));

\\
\\ finite places
\\
  pol *= denominator(content(pol))^2;
  c = content(pol);

  disc0 = poldisc(pol);
  plist = factor (abs(2*disc0));
if( DEBUGLEVEL_ell >= 4, print("    list of bad primes = ",plist));
  for( i = 1, #plist[,1],
    p = plist[i,1];
if( DEBUGLEVEL_ell >= 4, print("    p = ",p));
    vc = valuation(c,p);
    if( vc >= 2, 
      pol /= p^(2*(vc\2));
      plist[i,2] -= 2*(vc\2)*(2*poldegree(pol)-2)
    );
    if( poldegree(pol) == 4 && p != 2 && plist[i,2] < 2, next);
    if( !qp_soluble(pol,p),
if( DEBUGLEVEL_ell >= 3, print("   not ELS at ",p));
if( DEBUGLEVEL_ell >= 4, print("    end of locallysoluble"));
      return(0)));

if( DEBUGLEVEL_ell >= 2, print("  quartic ELS: Y^2 = ",pol));
if( DEBUGLEVEL_ell >= 4, print("    end of locallysoluble"));
  return(1);
}
{LS2localimage(nf,gen,pp) =
my(p,LS2image,ph,ival,delta);

if( DEBUGLEVEL_ell >= 4, print("    starting LS2localimage"));

  p = pp[1][1].p;
  LS2image = matrix( if( p == 2, sum(i=1,#pp,1+#pp[i][3].cyc), 2*#pp), #gen);

  for( j = 1, #gen,
    ph = [];
    for( i = 1, #pp,
      ival = idealval(nf,gen[j],pp[i][1]);
      ph = concat(ph,[ival]); 
      delta = gen[j]/pp[i][2]^ival;
      if( p == 2,
        ph = concat(ph,ideallog(nf,delta,pp[i][3])~);
      , ph = concat(ph,[1-nfpsquareodd(nf,delta,pp[i][4])]);
      )
    );
    LS2image[,j] = ph~
  );
  LS2image *= Mod(1,2);

if( DEBUGLEVEL_ell >= 4, print("    LS2image = ",lift(LS2image)));
if( DEBUGLEVEL_ell >= 4, print("    end of LS2localimage"));
  return(LS2image);
}

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\    GENERIC FUNCTIONS FOR ELLIPTIC CURVES    \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{ellhalf(ell,P)=
\\ returns all the points Q on ell such that 2Q = P.
my(pol2,ratroots,half,x2,y2,P2);

  if(#ell < 13, ell=ellinit(ell,1));

  pol2 = Pol([4,ell.b2,2*ell.b4,ell.b6]); \\ 2-division polynomial

  if( P == [0], 
    ratroots = nfroots(,pol2);
    half = vector(#ratroots,i,[ratroots[i],-(ell.a1*ratroots[i]+ell.a3)/2]);
    half = concat( [[0]], half);
    return(half)
  );

  x2 = Pol([1,0,-ell.b4,-2*ell.b6,-ell.b8]); \\ x(2P) = x2/pol2 

  half = [];
  ratroots = nfroots(,x2-P[1]*pol2);  
  if( #ratroots == 0, return(half));
  for( i = 1, #ratroots,
    y2 = ellordinate(ell,ratroots[i]);
    for( j = 1, #y2,
      P2 = [ratroots[i],y2[j]];
      if( ellmul(ell,P2,2) == P, half = concat(half,[P2]))
    )
  );

  return(half);
}
{elltors2(ell)=
\\ Compute the 2-torsion subgroup of the elliptic curve ell.
my(tors2);

if( DEBUGLEVEL_ell >= 3, print("   computing the 2-torsion"));

  tors2 = ellhalf(ell,[0]);
  if( #tors2 == 1, 
    tors2 = [1, [], []],
  if( #tors2 == 2, 
    tors2 = [2, [2], [tors2[2]]]
  , tors2 = [4, [2,2], [tors2[2],tors2[3]]]
  ));
if( DEBUGLEVEL_ell >= 3, print("   E[2] = ",tors2));
  return(tors2);
}
{elltorseven(ell)=
\\ Compute the 2-Sylow subgroup of the torsion of the elliptic curve ell.
my(torseven,P2);

if( DEBUGLEVEL_ell >= 4, print("    computing the 2^n-torsion"));
  if(#ell < 13, ell=ellinit(ell,1));
  torseven = elltors2(ell);

  while( torseven[1] != 1,
    P2 = ellhalf(ell,torseven[3][1]);
    if( #P2 > 0,
       torseven[1] *= 2;
       torseven[2][1] *= 2;
       torseven[3][1] = P2[1];
       next
    );
    if( #torseven[3] == 1, break());

    P2 = ellhalf(ell,torseven[3][2]);
    if( #P2 > 0,
       torseven[1] *= 2;
       torseven[2][2] *= 2;
       torseven[3][2] = P2[1];
       next
    );
    P2 = ellhalf(ell,elladd(ell,torseven[3][1],torseven[3][2]));
    if( #P2 > 0,
       torseven[1] *= 2;
       torseven[2][1] *= 2;
       torseven[3][1] = P2[1];
       next
    );
    break()
  );
  
if( DEBUGLEVEL_ell >= 4, print("    E[2^n] = ",torseven)); 
  return(torseven);
}
{ellsort(listpts) =
\\ Sorting the points listpts on an elliptic curve
\\ using the naive height.
my(n,v,aux,ord);

  v = vector(n = #listpts);
  for( i = 1, n, 
    if( listpts[i] == [0], v[i] = [0,0,0]; next);
    aux = denominator(listpts[i][2])/denominator(listpts[i][1]);
    v[i] = vecsort(abs([listpts[i][1]*aux^2, listpts[i][2]*aux^3,aux]),,4);
  );
  ord = Vec(vecsort(v,,3)); \\ ord = vecsort(v,,3);
  return(vector(n,i,listpts[ord[i]]));
}
{ellremovetorsion(ell,listgen) =
\\ Extracting the points of infinite order from listgen
my(d,extra);

if( DEBUGLEVEL_ell >= 5, print("     removing torsion from ",listgen));
  d = #listgen;
  extra = 0;
  for( i = 1, d,
\\ points of order 1 or 2
    if( listgen[i] == [0]
     || listgen[i] == ellmul(ell,listgen[i],-1)
    , extra += 1<<(i-1);
      next
    );
\\ detection of infinite order points by looking at
\\ 8*9*5*7*P modulo the prime 1048583
    if( ell.disc%1048583 != 0
     && denominator(listgen[i])%1048583 != 0
     && ellmul(ell,listgen[i]*Mod(1,1048583),2520) != [0]
    , next
    );
\\ detection of torsion points by ellorder()
    if( ellorder(ell,listgen[i]),
      extra += 1<<(i-1)
    )
  );
  if( extra,
    listgen = vecextract(listgen,1<<#listgen-1-extra);
  );
if( DEBUGLEVEL_ell >= 5, print("     without torsion = ",listgen));
  return(listgen);
}
{ellredgen(ell0,listgen,K=1) =
\\ reduction of the generators of points in listgen
\\ on the elliptic curve ell = [a1,a2,a3,a4,a6]
\\ or K*y^2 = x^3 + a2*x^2 + a4*x + a6 (when a1 = a3 = 0);
\\ using the canonical height.
my(d,ell=ell0,sqrtK,urst,extra,M,U,listgen2,tors2,vt);

if( DEBUGLEVEL_ell >= 3, print("   Reduction of the generators ",listgen));
if( DEBUGLEVEL_ell >= 5, print("     ell=",ell));          
  d = #listgen;
  if( d == 0, return([]));

\\ removing torsion points from listgen
  listgen = ellremovetorsion(ell0,listgen);
  d = #listgen;
  if( d == 0, return([]));  

  if( #ell < 13, ell = ellinit(ell,1));

  if( K != 1,
    if( ell.a1 != 0 || ell.a3 != 0, error(" ellredgen: a1*a3 != 0"));
    ell[2] *= K; ell[4] *= K^2; ell[5] *= K^3;
    ell[6] *= K; ell[7] *= K^2; ell[8] *= K^3; ell[9] *= K^4;
    ell[10] *= K^2; ell[11] *= K^3; ell[12] *= K^6;
    sqrtK = sqrt(K);
    if( #ell == 19,
      ell[14] *= K;
      ell[15] /= sqrtK; ell[16] /= sqrtK;
      ell[17] *= sqrtK; ell[18] *= sqrtK;
      ell[19] /= K
    );

    for( i = 1, d,
      for( j = 1, #listgen[i],
        listgen[i][j] *= K^j))
  );

if( d == 1,
  urst = [1,0,0,0];
,
  if( #ell < 19, ell = ellinit(ell));
  ell = ellminimalmodel(ell,&urst);
  listgen = ellchangepoint(listgen,urst);
if( DEBUGLEVEL_ell >= 5, print("     ell = ",ell));
if( DEBUGLEVEL_ell >= 5, print("     listgen = ",listgen));

\\ Looking for relations between the points in listgen
\\ using LLL on the height matrix

  extra = 1;
  while( extra,
    M = ellheightmatrix(ell,listgen);
if( DEBUGLEVEL_ell >= 4, print("    height matrix = ",M));
    if( abs(matdet(M)) > 10^(-default(realprecision)+10),break);
    U = qflll(round(M*10^(default(realprecision)-10)),4);
    U = concat(U[1],U[2]);
if( DEBUGLEVEL_ell >= 4, print("    change of basis proposed by LLL = ",U));
\\ the columns of U that have very small coefficients
\\ are either exact relations or reductions (coeff <= 20)
\\ the other ones are irrelevant.
    extra = 0;
    for( i = 1, d,
      if( vecmax(abs(U[,i])) > 20, extra += 1<<(i-1))
    );
    U = vecextract(U,1<<d-1-extra);
    U = completebasis(U);
if( DEBUGLEVEL_ell >= 4, print("    change of basis 1 = ",U));

    listgen2 = vector(d);
    for( i = 1, d,
      listgen2[i] = [0];
      for( j = 1, d,
        listgen2[i] = elladd(ell,listgen2[i],ellmul(ell,listgen[j],U[j,i]))));
    listgen = listgen2;
  );

\\ Extracting the points of infinite order

\\ removing torsion points from listgen
    listgen = ellremovetorsion(ell,listgen);
    d = #listgen;
    if( d == 0, return([]));  
  );

if( DEBUGLEVEL_ell >= 3, print("   infinite order points = ",listgen));

\\ Now, the points should be of infinite order and independant
\\ Reducing the points of infinite order

  if( d > 1,
    M = ellheightmatrix(ell,listgen);
if( DEBUGLEVEL_ell >= 4, print("    height matrix = ",M));
    U = qflllgram(M);
if( DEBUGLEVEL_ell >= 4, print("    change of basis 2 = ",U));

    listgen2 = vector(d);
    for( i = 1, d,
      listgen2[i] = [0];
      for( j = 1, d,
        listgen2[i] = elladd(ell,listgen2[i],ellmul(ell,listgen[j],U[j,i]))));
    listgen = listgen2;
  );

if( DEBUGLEVEL_ell >= 3, print("   infinite order points = ",listgen));

  listgen = ellchangepoint(listgen,ellinverturst(urst));

\\ Reducing modulo the 2-torsion

  tors2 = elltorseven(ell0);
  if( tors2[1] > 1,
    vt = vector(tors2[2][1],j,ellmul(ell0,tors2[3][1],j-1));
    if( #tors2[2] == 2,
      vt = concat(vt,vector(#vt,j,elladd(ell0,vt[j],tors2[3][2])))
    );
    for( i = 1, d,
      listgen[i] = ellsort(vector(#vt,j,elladd(ell0,listgen[i],vt[j])))[1];
    );
  );

if( DEBUGLEVEL_ell >= 3, print("   infinite order points = ",listgen));
  
  if( K != 1,
    for( i = 1, d,
      for( j = 1, 2,
        listgen[i][j] /= K^j)));

\\ keep only the points (x,y) with y >= 0

  if( ell0.a1 == 0 && ell0.a3 == 0, 
    for( i = 1, d,
      if( #listgen[i] == 2,
        listgen[i][2] = abs(listgen[i][2]))));

if( DEBUGLEVEL_ell >= 2, print("  reduced generators = ",listgen));
  return(listgen);
}

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\    FUNCTIONS FOR NUMBER FIELDS              \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{reducemodsquares(delta,d) =
\\ Uses LLL to find z such that delta*z^2 has a small coefficient in x^d.
\\ delta must be a t_POLMOD
my(deg,xx,z,qd,Qd,reduc);

  deg = poldegree(component(delta,1));  \\ deg = poldegree(delta.mod);
  xx = Mod('x,component(delta,1)); \\ xx = Mod(x,delta.mod);
  z = subst(Pol(vector(deg,i,eval(Str("a"i)))),'x,xx);
  qd = polcoeff(lift(delta*z^2),d,'x);
  Qd = simplify(matrix(deg,deg,i,j,deriv(deriv(qd,eval(Str("a"i))),eval(Str("a"j)))/2));

  reduc = qflllgram_indef(Qd);
  if( #reduc == 2, reduc = reduc[2][,1]);

  return(delta*subst(Pol(reduc),'x,xx)^2);
}
{bnfpSelmer(bnf,S=1,p) =
\\ p is a prime integer and bnf a big number field.
\\ Compute the p-Selmer group of the number field bnf
\\ relative to the prime ideals dividing S.
\\ This group is denoted by LS2 in the sequel.
\\ Returns [gen,S'] where gen is a vector containing the generators
\\ of the p-Selmer group, represented as elements of bnf modulo p-powers,
\\ and S' is the support of gen.
my(S1,oddclass,multS,Slist,LS2gen,newprimes,newprimesval,kerval);

if( DEBUGLEVEL_ell >= 3, print("   Constructing the field Selmer group: L(S,",p,")"));
  S1 = idealhnf(bnf,S);

  oddclass = 0; multS = 1;
  while( !oddclass,
    if( multS != 1, S1 = idealmul(bnf,S1,multS));
    Slist = idealfactor(bnf,S1)[,1]~;
if( DEBUGLEVEL_ell >= 4, print("    constructing the S-units "));
if( DEBUGLEVEL_ell >= 4, print("    S1 = ",Slist));
    LS2gen = bnfsunit(bnf,Slist);

\\ If the class group is divisible by p,
\\ need to enlarge S1.
    oddclass = LS2gen[5].no % p;
    if( !oddclass,
if( DEBUGLEVEL_ell >= 4, print("    class group divisible by p = ",LS2gen[5].no));
      multS = idealmul(bnf,S,LS2gen[5].gen[1]);
    )
  );
  LS2gen = Mod(LS2gen[1],bnf.pol);

\\ The valuation of the generators must be divisible by p outside S.
  newprimes = [];
  for( i = 1, #Slist,
    if( idealadd(bnf,S,Slist[i]) == 1,
      newprimes = concat(newprimes,[Slist[i]])));
if( DEBUGLEVEL_ell >= 4, print("    newprimes = ",newprimes));
  newprimesval = matrix(#newprimes,#LS2gen,i,j,
    idealval(bnf,LS2gen[j],newprimes[i]));
if( DEBUGLEVEL_ell >= 4, print("    newprimesval = ",newprimesval));
  kerval = lift(matker(newprimesval*Mod(1,p)));
if( DEBUGLEVEL_ell >= 4, print("    kerval = ",kerval));
  LS2gen = vector(#kerval,i,
    prod( j = 1, #LS2gen,
      LS2gen[j]^kerval[j,i]));

\\ Add the units
  LS2gen = concat(Mod(bnf[8][5],bnf.pol),LS2gen); \\ LS2gen = concat(bnf.fu,LS2gen);
\\ Add also the torsion unit if its order is divisible by p.
  if( bnf[8][4][1]%p == 0, \\ if( bnf.tu[1]%p == 0,
    LS2gen = concat( [Mod(bnf[8][4][2],bnf.pol)], LS2gen)); \\ LS2gen = concat( [Mod(bnf.tu[2],bnf.pol)], LS2gen));
if( DEBUGLEVEL_ell >= 3, print("   #LS2gen = ",#LS2gen));
if( DEBUGLEVEL_ell >= 4, print("    LS2gen = ",LS2gen));
  return([LS2gen,Slist]);
}
{kersign(gen,rootapprox) =
\\ Determine the kernel of the sign map
\\ restricted to the subgroup generated by gen,
\\ and relative to the real embedding corresponding to
\\ the root of pol contained in the interval rootapprox.
my(signs,elt,elt2,d,st,kers,compt);

if( DEBUGLEVEL_ell >= 3, print("   Computing the kernel of the sign ",rootapprox));

\\ determination of the signs

  signs = vector(#gen);
  for( i = 1, #gen,
    elt = lift(gen[i]);
if( DEBUGLEVEL_ell >= 5, print("     Computing the sign of elt = ",elt));
    if( poldegree(elt) == 0, signs[i] = sign(simplify(elt)) < 0; next);
    d = poldisc(elt);
    if( poldegree(elt) == 2,
      if( d <= 0, signs[i] = sign(pollead(elt)) < 0; next));
    elt2 = if( d == 0, elt/gcd(elt,elt'), elt);
    st = 1;
    compt = 0;
    while( st,
      st = polsturm(elt2,[rootapprox[1],rootapprox[2]]);
      if( st,
        rootapprox = polrealrootsimprove(component(gen[i],1),rootapprox);
\\        rootapprox = polrealrootsimprove(gen[i].mod,rootapprox);
\\ if the sign of elt is too difficult to determine, 
\\ try a reduction modulo squares
        if( compt++ == 5, gen[i] = reducemodsquares(gen[i]); i--; next(2));
\\ if the sign of elt is still too difficult to determine, 
\\ try the sign of 1/elt.
        if( compt%5 == 0, gen[i] = 1/gen[i]; i--; next(2))
    ));
    signs[i] = sign(subst(elt,variable(elt),rootapprox[2])) < 0
  );
if( DEBUGLEVEL_ell >= 4, print("    signs = ",signs));

\\ construction of the kernel
  kers = matker(Mat(signs*Mod(1,2)))*Mod(1,2);
if( DEBUGLEVEL_ell >= 4, print("    kers = ",lift(kers)));
  return(kers);
}
{kernorm(gen,S,p) =
\\ gen is a generating set for a subgroup of K^* / K^p.
\\ Compute the kernel of the norm map.
\\ Uses the fact that all valuations are 0 mod p,
\\ except maybe at primes in S.
my(normgen,normmap,kern);
  
if( DEBUGLEVEL_ell >= 3, print("   Computing the kernel of the norm map"));

  if( p == 2, S = concat([-1],S));
  normgen = norm(gen);
if( DEBUGLEVEL_ell >= 4, print("    normgen = ",normgen));

\\ matrix of the norm map
  normmap = matrix(#S,#normgen,i,j,
    if( i == 1 && p == 2,
      sign(normgen[j]) < 0
    , valuation(normgen[j],S[i])));
if( DEBUGLEVEL_ell >= 4, print("    normmap = ",normmap));

\\ construction of the kernel
  kern = matker(normmap*Mod(1,p))*Mod(1,p);
if( DEBUGLEVEL_ell >= 4, print("    ker = ",lift(kern)));
  return(kern);
}

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\    FUNCTIONS FOR 2-DESCENT                  \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{elllocalimage( nf, pp, K = 1) =
\\ pol is the cubic polynomial defining the elliptic curve,
\\ nf is nfinit(pol),
\\ p is a prime integer, and pp = ppinit(p).
\\ Returns the image of the p-adic points 
\\ E(Qp)/2E(Qp) in Kp/Kp^2.
\\ The algorithm consists of choosing random p-adic points in E(Qp)
\\ until the number of images is equal to #E(Qp)[2] / |2|_p
my(X,p,prank,rac,pts,bound,essai,mrank,r,xx,delta,ph,delta2,mylocalprec,ival);

if( DEBUGLEVEL_ell >= 4, print("    starting elllocalimage"));

  X = Mod('x,nf.pol);
  p = pp[1][1][1];
  prank  = #pp - (p != 2);
if( DEBUGLEVEL_ell >= 4, print("    prank = ",prank));

  rac = polrootsmodpn(K*nf.pol,p);
if( DEBUGLEVEL_ell >= 5, print("     rac = ",rac));

  pts = matrix(0,0);
  bound = p+6;
  essai = 0;
  mrank = 0;
  while( mrank < prank,

    essai ++;
    if( essai%16 == 0,
      pts = matimage(pts);
      bound *= p;
    );

    r = random(#rac)+1; mylocalprec = random(rac[r][2]+3)-2;
    xx = rac[r][1]+p^mylocalprec*random(bound);
if( DEBUGLEVEL_ell >= 5, print("     xx = ",xx));
    delta = K*(xx-X);

\\ rem: K*pol(xx) = norm(delta) ( = y^2 for a point on the elliptic curve)
    if( !psquare(K*subst(nf.pol,'x,xx),p), next);
    ph = [];
    for( i = 1, #pp,
      ph = concat(ph,[ ival = idealval(nf,delta,pp[i][1])]);
      delta2 = delta/pp[i][2]^ival;
      if( p == 2,
        ph = concat(ph,ideallog(nf,delta2,pp[i][3])~);
      , ph = concat(ph,[1-nfpsquareodd(nf,delta2,pp[i][4])])
      )
    );
if( DEBUGLEVEL_ell >= 5, print("     ph = ",ph));

    pts = concat(pts,ph~*Mod(1,2));
    mrank = matrank(pts*Mod(1,2));
if( DEBUGLEVEL_ell >= 5, print("     pts = ",lift(pts)));
if( DEBUGLEVEL_ell >= 5, print("    matrank = ",mrank));
  );

  pts = matimage(pts);
if( DEBUGLEVEL_ell >= 5, print("     essai = ",essai));
if( DEBUGLEVEL_ell >= 4, print("    end of elllocalimage"));
  return(pts);
}
{ell2descent_twist(ell,bnf,K=1,help=[],redflag=0) =
\\ This algorithm performs 2-descent on the elliptic curve ell twisted by K
\\ K*y^2=x^3+A*x^2+B*x+C
\\ when ell has trivial 2-torsion.

\\ ell=[0,A,0,B,C], with K,A,B,C integers.
\\ bnf is bnfinit(x^3+A*x^2+B*x+C,1).

\\
\\ help is a list of known points (maybe empty) on the twisted curve.
\\ if redflag != 0, reduces the elements 
\\ of the field Selmer group modulo squares.

my(A,B,C,polrel,polprime,ttheta,badprimes,S,LS2,selmer,rootapprox,p,pp,locimage,LS2image,listpointstriv,listpoints,iwhile,expo,zc,liftzc,den,point,idealfactorzc,idealzc,baseidealzc,q2,sol,param,q1,pol,redq,q0,pointxx,point2,rang);

if( DEBUGLEVEL_ell >= 4, print("    starting ell2descent_twist"));

  if( #ell < 13, ell = ellinit(ell,1));

  if( ell.a1 != 0 || ell.a3 != 0,
    error(" ell2descent_twist: the curve is not of the form [0,a,0,b,c]"));
  if( denominator(ell.a2) > 1 || denominator(ell.a4) > 1 || denominator(ell.a6) >1,
    error(" ell2descent_twist: non integral coefficients"));

  A = ell.a2; if( DEBUGLEVEL_ell >= 2, print("  A = ",A));
  B = ell.a4; if( DEBUGLEVEL_ell >= 2, print("  B = ",B));
  C = ell.a6; if( DEBUGLEVEL_ell >= 2, print("  C = ",C));
  
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\      Construction of L(S,2)      \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

if( DEBUGLEVEL_ell >= 2, print(); print("  Computing L(S,2)"));

  polrel = Pol([1,A,B,C]);
  polprime = polrel';
  ttheta = Mod('x,polrel);

  if( !bnf, 
if( DEBUGLEVEL_ell >= 3, print("   bnfinit(",polrel,")"));
    bnf = bnfinit(polrel,1));

  badprimes = abs(K*idealadd(bnf,polprime,bnf.index));
if( DEBUGLEVEL_ell >= 5, print("     badprimes = ",badprimes[1,1]));
  S = bnfpSelmer(bnf,badprimes,2);
  LS2 = S[1]; S = S[2];

if( DEBUGLEVEL_ell >= 2, print("  L(S,2) = ",LS2));

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\   Construction of the Selmer group \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

if( DEBUGLEVEL_ell >= 2, print(); print("  Computing the Selmer group"));

\\ elements with square norm
  selmer = kernorm(LS2,vector(#S,i,S[i].p),2);
if( DEBUGLEVEL_ell >= 3, print("   selmer = ",lift(selmer)));

\\ the first real embedding must be > 0
\\ since the norm is a square, this is automatic
\\ if there is a single real embedding.
  if( bnf.r1 == 3,
    rootapprox = polrealrootsisolate(polrel)[1];
    selmer = matintersect(selmer,kersign(LS2,rootapprox))*Mod(1,2);
if( DEBUGLEVEL_ell >= 3, print("   selmer = ",lift(selmer)));
  );

\\ p-adic points
if( DEBUGLEVEL_ell >= 3, print("   p-adic points"));
  badprimes = factorint(badprimes[1,1]*2)[,1];
if( DEBUGLEVEL_ell >= 2, print("  badprimes = ",badprimes));
  for( i = 1, #badprimes,
    p = badprimes[i];
if( DEBUGLEVEL_ell >= 4, print("    p = ",p));
    pp = ppinit(bnf.nf,p);
    locimage = elllocalimage(bnf.nf,pp,K);
    LS2image = LS2localimage(bnf.nf,LS2,pp);
    locimage = matintersect(LS2image,locimage);
    selmer = matintersect(
      selmer,
      concat(
        matker(LS2image),
        matinverseimage(LS2image,locimage)*Mod(1,2)));
    selmer = matimage(selmer*Mod(1,2));
if( DEBUGLEVEL_ell >= 4, print("    selmer = ",lift(selmer)));
  if( !#selmer, break);
  );

if( DEBUGLEVEL_ell >= 2, print("  selmer = ",lift(selmer)));
if( DEBUGLEVEL_ell >= 2, print("  Selmer rank = ",#selmer));

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\   Search for trivial points      \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

  if( #selmer,
if( DEBUGLEVEL_ell >= 2, print(); print("  Search for trivial points on the curve"));
    listpointstriv = ratpoint(K^3*subst(polrel,'x,'x/K),LIMTRIV,0);
    for( i = 1, #listpointstriv, 
      if( #listpointstriv[i] == 3,
        listpointstriv[i] = [0]
      , for( j = 1, 2, listpointstriv[i][j] /= K^j))
     );
    listpointstriv = concat(help,listpointstriv);
if( DEBUGLEVEL_ell >= 2, print("  Trivial points on the curve = ",listpointstriv));
  );

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\   Run through the Selmer group   \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

if( DEBUGLEVEL_ell >= 2, print(); print("  Run through the Selmer group"));

  listpoints = [];
  selmer = lift(selmer);
  iwhile = 1;
  while( iwhile < 1<<#selmer,
if( DEBUGLEVEL_ell >= 2, print());
if( DEBUGLEVEL_ell >= 4, print("   iwhile = ",iwhile));

\\ the next element zc as an algebraic number modulo squares

    expo = selmer*vectorv(#selmer,i,bittest(iwhile,i-1));
    zc = prod( i = 1, #LS2, LS2[i]^expo[i]);
if( DEBUGLEVEL_ell >= 2, print("  zc = ",zc));
    liftzc = lift(zc);

\\ Reduction modulo squares

    if( redflag, 
      zc = reducemodsquares(zc,2);
      liftzc = lift(zc);
      den = denominator(content(liftzc))^2;
      zc *= den; liftzc *= den;
if( DEBUGLEVEL_ell >= 2, print("  zc reduced = ",zc))
    );

\\ Does it come from a trivial point ?

    for( i = 1, #listpointstriv,
      point = listpointstriv[i];
      if( #point == 2,
        if( nfissquare(bnf.nf,K*(point[1]-'x)*liftzc),
if( DEBUGLEVEL_ell >= 2, print("  comes from the trivial point ",point));
          listpoints = concat(listpoints,[point]);
          iwhile = 1 << (degre(iwhile)+1);
          next(2)
    )));
    
if( DEBUGLEVEL_ell >= 2, print("  does not come from a trivial point"));

\\ Construction of the quadratic form q2
\\ Change the basis using the square factors of zc

    idealfactorzc = idealfactor(bnf,zc);
    idealfactorzc[,2] *= -1;
    idealfactorzc[,2] \= 2;
\\    idealzc = idealfactorback(bnf,idealfactorzc);
    idealzc = matid(3);
    for( i = 1, #idealfactorzc[,1],
      idealzc = idealmul(bnf,idealzc,idealpow(bnf,idealfactorzc[i,1],idealfactorzc[i,2]));
    );
    baseidealzc = vector(3,i,nfbasistoalg(bnf,idealzc[,i]));
    q2 = matrix(3,3,i,j,trace(zc*baseidealzc[i]*baseidealzc[j]/polprime));
if( DEBUGLEVEL_ell >= 4, print("    q2 = ",q2));
if( DEBUGLEVEL_ell >= 4, print("    q2/content(q2) = ",q2/content(q2)));

\\ Solution of the quadratic equation q2=0

    sol = qfsolve(q2/content(q2));
if( DEBUGLEVEL_ell >= 4,print("    sol = ",sol));
    if( type(sol) == "t_INT", 
      error(" ell2descent_twist: WRONG ELEMENT IN THE SELMER GROUP, please report"));


\\ Parametrizing the solutions of q2=0

    param = qfparam(q2,sol)*['x^2,'x,1]~;
    param /= content(param);

\\ Construction of the quartic

    q1 = -matrix(3,3,i,j,trace(zc*baseidealzc[i]*baseidealzc[j]*(ttheta+A)/polprime));
    pol = param~*q1*param;
if( DEBUGLEVEL_ell >= 2, print("  quartic: ",K,"*Y^2 = ",pol));
    redq = redquartic(pol);
    pol = redq[1];
    den = denominator(content(K*pol));
    pol *= den^2;
if( DEBUGLEVEL_ell >= 2, print("  reduced: ",K,"*Y^2 = ",pol));

\\ Search for points on the quartic

    point = ratpoint(K*pol,LIM1,1);
    if( point == [], point = ratpoint2(K*pol,LIM3,1));
    if( point == [], iwhile ++; next );

    if( #point == 2, point = concat(point,[1]));
if( DEBUGLEVEL_ell >= 2, print("  point on the reduced quartic = ",point));
    point = concat(redq[2]*[point[1],point[3]]~,[point[2]/den]~);
if( DEBUGLEVEL_ell >= 2, print("  point on the quartic = ",point));

\\ Construction of the point on the elliptic curve from the point on the quartic

    param = subst(param,'x,'x/'y)*'y^2;
    param = subst(subst(param,'x,point[1]),'y,point[2]);
    param *= K/point[3];
if( DEBUGLEVEL_ell >= 3, print("   reconstruction of the point on the curve"));
    q0 = matrix(3,3,i,j,trace(zc*baseidealzc[i]*baseidealzc[j]*(ttheta^2+A*ttheta+B)/polprime));
    pointxx = param~*q0*param/K;
    point2 = [ pointxx, sqrtrat(subst(polrel,'x,pointxx)/K)];
if( DEBUGLEVEL_ell >= 1, print(" point on the curve = ",point2));
    listpoints = concat(listpoints,[point2]);
    iwhile = 1 << (degre(iwhile)+1)
  );

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\      Conclusion report           \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

  rang = #listpoints;

if( DEBUGLEVEL_ell >= 2,
  print();
  print("  rank of found points     = ",#listpoints);
  print("  rank of the Selmer group = ",#selmer));

if( DEBUGLEVEL_ell >= 1, afficheselmer(rang,#selmer));
  
  if(  (#selmer - rang)%2,
    rang ++;
if( DEBUGLEVEL_ell >= 1, 
      print(" III should be a square, hence ");
      afficheselmer(rang,#selmer));
  );

\\ Verification

if( DEBUGLEVEL_ell >= 1, print("listpoints = ",listpoints));
  for( i = 1, #listpoints, 
    if( subst(polrel,'x,listpoints[i][1])-K*listpoints[i][2]^2,
      error(" ell2descent_twist: WRONG POINT = ",listpoints[i]," please report")));

\\ Reduction of the points

  listpoints = vecsort(listpoints,,2);
  if( #listpoints >= 2 && ELLREDGENFLAG,
    listpoints = ellredgen(ell,listpoints,K));

if( DEBUGLEVEL_ell >= 4, print("    end of ell2descent_twist"));
  return([rang,#selmer,listpoints]);
}
{afficheselmer(m1,m2,tors2) =

    print("#E(Q)[2]      = ",1<<tors2);
    print("#S(E/Q)[2]    = ",1<<m2);
  if( m1+tors2 == m2,
    print("#E(Q)/2E(Q)   = ",1<<(m1+tors2));
    print("#III(E/Q)[2]  = 1");
    print("rank(E/Q)     = ",m1);
  ,
    print("#E(Q)/2E(Q)  >= ",1<<(m1+tors2));
    print("#III(E/Q)[2] <= ",1<<(m2-m1-tors2));
    print("rank(E/Q)    >= ",m1)
  );
}
{ellrank(ell,help=[]) =
\\ Algorithm of 2-descent on the elliptic curve ell.
\\ help is a list of known points on ell.
my(urst,urst1,den,eqell,tors2,bnf,rang,time1);

if( DEBUGLEVEL_ell >= 3, print("   starting ellrank"));
  if( #ell < 13, ell = ellinit(ell,1));

\\ kill the coefficients a1 and a3
  urst = [1,0,0,0];
  if( ell.a1 != 0 || ell.a3 != 0,
    urst1 = [1,0,-ell.a1/2,-ell.a3/2];
    ell = ellchangecurve(ell,urst1);
    urst = ellcomposeurst(urst,urst1)
  );

\\ kill denominators
  while( (den = denominator([ell.a2,ell.a4,ell.a6])) > 1,
    den = factor(den); den[,2] = vectorv(#den[,2],i,1);
    den = factorback(den);
    urst1 = [1/den,0,0,0];
    ell = ellchangecurve(ell,urst1);
    urst = ellcomposeurst(urst,urst1)
  );

  help = ellchangepoint(help,urst);
  eqell = Pol([1,ell.a2,ell.a4,ell.a6]);
if( DEBUGLEVEL_ell >= 1, print(" Elliptic curve: Y^2 = ",eqell));

\\ Choice of the algorithm depending on the 2-torsion structure

  tors2 = ellhalf(ell,[0]);
if( DEBUGLEVEL_ell >= 1, print(" E[2] = ",tors2));

  if( #tors2 == 1,                              \\ case 1: 2-torsion trivial
if( DEBUGLEVEL_ell >= 3, print1("   bnfinit(",eqell,")"));
if( DEBUGLEVEL_ell >= 4, gettime());
    bnf = bnfinit(eqell,1);
if( DEBUGLEVEL_ell >= 4, time1 = gettime());
if( DEBUGLEVEL_ell >= 3, print(" done"));
    rang = ell2descent_twist(ell,bnf,1,help);
if( DEBUGLEVEL_ell >= 4, print("    time for bnfinit  = ",time1));
if( DEBUGLEVEL_ell >= 4, print("    time for the rest = ",gettime()));
  ,
  if( #tors2 == 2 || !COMPLETE,                 \\ case 2: 2-torsion >= Z/2Z
    if( ell.a6 != 0,
      urst1 = [1,tors2[2][1],0,0];
      ell = ellchangecurve(ell,urst1);
      urst = ellcomposeurst(urst,urst1)
    );
    eqell = Pol([1,ell.a2,ell.a4,ell.a6]);
if( DEBUGLEVEL_ell >= 1, print(" Elliptic curve: Y^2 = ",eqell));

    rang = ell2descent_viaisog(ell,help)
  ,                                             \\ case 3: 2-torsion = Z/2Z*Z/2Z
    rang = ell2descent_complete(tors2[2][1],tors2[3][1],tors2[4][1])
  ));

  rang[3] = ellchangepoint(rang[3],ellinverturst(urst));
if( DEBUGLEVEL_ell >= 3, print("   end of ellrank"));

  return(rang);
}
{ell2descent_complete(e1,e2,e3,help) =
\\ Compute the rank of the elliptic curve
\\ E: Y^2 = (x-e1)*(x-e2)*(x-e3)
\\ using the complete 2-descent algorithm (see J.Silverman).
\\ Returns [r,s,v] with
\\   r is a lower bound for the rank of E(Q)
\\   s is the rank of Selmer[2]
\\   v is a system of independant points on E(Q)/2E(Q)

\\ e1, e2, e3 must be integers.
\\ help is a list of known points on E.

my(ee,d32,d31,d21,G1,G2,G3,vect1,vect2,vect3,selmer,rang,listpoints,b1,b2,q1,sol1,param1,param1x,quart,point,z1,solx,soly,strange,ell);

if( DEBUGLEVEL_ell >= 2, print("  Algorithm of complete 2-descent"));

\\ sort the integers e1, e2, e3 in increasing order

  ee = vecsort([e1,e2,e3]);
  e1 = ee[1]; e2 = ee[2]; e3 = ee[3];

\\ Computation of the groups G1 and G2

  d32 = e3-e2; d31 = e3-e1; d21 = e2-e1;
  G1 = factor(d31*d21)[,1];   \\ (G1 > 0)
  G2 = factor(-d32*d21)[,1];  \\ (G2 < 0)
  G3 = d31*d32;

if( DEBUGLEVEL_ell >= 3, print("   G1 = ",G1));
if( DEBUGLEVEL_ell >= 3, print("   G2 = ",G2));

\\ Run through G1*G2

  vect1 = vector(#G1,i,[0,1]);
  vect2 = vector(#G2,i,[0,1]);
  selmer = 0;
  rang = 0;
  listpoints = [];

  forvec( X = vect1,
    b1 = prod( i = 1, #G1, G1[i]^X[i]);

\\ b1*b2*b3 must be a square, where b3 is a divisor of d32*d31
    vect3 = vect2;
    for( i = 2, #G2, 
      if( G3%G2[i] !=0,
        vect3[i] = [1,1]*valuation(b1,G2[i])));

    forvec( Y = vect3,
      b2 = prod( i = 1, #G2, G2[i]^Y[i]);

if( DEBUGLEVEL_ell >= 3, print("   [b1,b2] = ",lift([b1,b2])));

\\ Trivial points coming from the 2-torsion

      if( b1==1 && b2==1,
if( DEBUGLEVEL_ell >= 4, print("    trivial point [0]"));
        selmer++; rang++; next);
      if( issquare(-d21*b2) && issquare(d31*d21*b1),
if( DEBUGLEVEL_ell >= 3, print("   trivial point [e1,0]"));
        selmer++; rang++; listpoints = concat(listpoints,[[e1,0]]); next);
      if( issquare(d21*b1) && issquare(-d32*d21*b2),
if( DEBUGLEVEL_ell >= 3, print("   trivial point [e2,0]"));
        selmer++; rang++; listpoints = concat(listpoints,[[e2,0]]); next);
      if( issquare(d31*b1) && issquare(d32*b2),
if( DEBUGLEVEL_ell >= 3, print("   trivial point [e3,0]"));
        selmer++; rang++; listpoints = concat(listpoints,[[e3,0]]); next);

\\ Trivial points coming from help

      for( i = 1, #help,
        if( #help[i] != 2 || help[i][2] == 0, next);
        if( issquare(b1*(help[i][1]-e1)) && issquare(b2*(help[i][1]-e2)),
if( DEBUGLEVEL_ell >= 3, print("   trivial point from help ",help[i]));
          selmer++; rang++; 
          listpoints = concat(listpoints,[help[i]]); next(2));
      );

\\ If one can solve 2 quadratic equations
\\ (1) q1: b1*z1^2-b2*z2^2 = e2-e1
\\ (2) q2: b1*z1^2-b1*b2*z3^2 = e3-e1
\\ then (x,y) = (b1*z1^2+e1,b1*b2*z1*z2*z3) is a point on E
\\ we also have
\\ (3) q3 = q1-q2: b1*b2*z3^2-b2*z2^2=e2-e3

\\ Solution of the q1

      q1 = matdiagonal([b1,-b2,-d21]);
if( DEBUGLEVEL_ell >= 3, print("   q1 = ",q1));
      sol1 = qfsolve(q1);
      if( type(sol1) == "t_INT",
if( DEBUGLEVEL_ell >= 3, print("   q1 not ELS at ",sol1));
        next);
if( DEBUGLEVEL_ell >= 3, print("   solution of q1 = ",sol1));
      param1 = qfparam(q1,sol1,1);
      param1 /= content(param1);
if( DEBUGLEVEL_ell >= 3, print("   parametrization of q1 = ",param1));
      param1x = param1*['x^2,'x,1]~;

\\ Solution of the q2
\\ only useful to detect local non solubility

\\ my(q2,sol2);
\\      q2 = matdiagonal([b1,-b1*b2,-d31]);
\\if( DEBUGLEVEL_ell >= 3, print("   q2 = ",q2));
\\      sol2 = qfsolve(q2);
\\      if( type(sol2) == "t_INT",
\\if( DEBUGLEVEL_ell >= 3, print("   q2 not ELS at ",sol2));
\\        next);

\\ Construction of the quartic

      quart = b1*b2/gcd(b1,b2)^2*(b1*param1x[1]^2-d31*param1x[3]^2);
if( DEBUGLEVEL_ell >= 3, print("   quart = ",quart));

\\ Search for trivial points on the quartic

      point = [];
\\      point = ratpoint(quart,LIM1,1);

\\ Local solubility of the quartic
   
      if( point == [] && !locallysoluble(quart),
if( DEBUGLEVEL_ell >= 3, print("   quartic not ELS "));
        next);
if( DEBUGLEVEL_ell >= 2, print("  y^2 = ",quart));
      selmer++;

\\ Search for points on the quartic

      if( point == [], point = ratpoint2(quart,LIM3,1));
      if( point != [], 
if( DEBUGLEVEL_ell >= 2, print("  point found on the quartic !!"));
if( DEBUGLEVEL_ell >= 3, print("   ",point));
        if( #point == 2,
          z1 = subst(param1x[1],'x,point[1])/subst(param1x[3],'x,point[1])
        , z1 = param1[1,1]/param1[3,1]);
        solx = b1*z1^2+e1;
        soly = sqrtrat((solx-e1)*(solx-e2)*(solx-e3));
        listpoints = concat(listpoints,[[solx,soly]]);
if( DEBUGLEVEL_ell >= 1, print(" point on the elliptic curve = ",[solx,soly]));
        rang++
      ,
if( DEBUGLEVEL_ell >= 2, print("  no point found on the quartic"))
      )
    )
  );

\\ end

if( DEBUGLEVEL_ell >= 1, 
  print("#S^(2)      = ",selmer));
  if( rang > selmer/2, rang = selmer);
if( DEBUGLEVEL_ell >= 1,
  strange = rang != selmer;
  if( strange,
  print("#E[K]/2E[K]>= ",rang)
, print("#E[K]/2E[K] = ",rang));
  print("#E[2]       = 4");
);
  rang = ceil(log(rang)/log(2))-2;
  selmer = valuation(selmer,2);
if( DEBUGLEVEL_ell >= 1,
  if( strange,
  print(selmer-2," >= rank  >= ",rang)
, print("rank        = ",rang));
  if( rang, print("points = ",listpoints));
);
  ell = ellinit([0,-(e1+e2+e3),0,e1*e2+e2*e3+e3*e1,-e1*e2*e3],1);
  listpoints = vecsort(listpoints,,2);
  if( ELLREDGENFLAG, listpoints = ellredgen(ell,listpoints));
  listpoints = concat(ellsort(elltorseven(ell)[3]),listpoints);

  return([rang,selmer,listpoints]);
}
{ellcount( c, d, KS2gen, listpointstriv=[]) =
my(found,listgen,listpointscount,m1,m2,lastloc,mask,i,d1,iaux,j,triv,pol,point,qf,solqf,para,point1,v);

if( DEBUGLEVEL_ell >= 4, print("    starting ellcount ",[c,d]));
if( DEBUGLEVEL_ell >= 4, print("    KS2gen = ",KS2gen));
if( DEBUGLEVEL_ell >= 4, print("    listpointstriv = ",listpointstriv));

  found = 0;
  listgen = KS2gen;
  listpointscount = [];

  m1 = m2 = 0; lastloc = -1;

  mask = 1 << #KS2gen;
  i = 1;
  while( i < mask,
    d1 = 1; iaux = i; j = 1;
    while( iaux, 
      if( iaux%2, d1 *= listgen[j]);
      iaux >>= 1; j++);
if( DEBUGLEVEL_ell >= 3, print("   d1 = ",d1));
    triv = 0;
    for( j = 1, #listpointstriv,
      if( listpointstriv[j][1] && issquare(d1*listpointstriv[j][1]),
        listpointscount = concat(listpointscount,[listpointstriv[j]]);
if( DEBUGLEVEL_ell >= 2, print("  trivial point"));
        triv = 1; m1++;
        if( degre(i) > lastloc, m2++);
        found = 1; lastloc = -1; break));
    if( !triv,
    pol = Pol([d1,0,c,0,d/d1]);
if( DEBUGLEVEL_ell >= 3, print("   quartic = y^2 = ",pol));
    point = ratpoint(pol,LIM1,1);
    if( point != [],
if( DEBUGLEVEL_ell >= 2, print("  point on the quartic"));
if( DEBUGLEVEL_ell >= 3, print(   point));
      m1++;
      listpointscount = concat(listpointscount,[d1*point[1]*point]);
      if( degre(i) > lastloc, m2++);
      found = 1; lastloc = -1
    ,
      if( locallysoluble(pol),
        if( degre(i) > lastloc, m2++; lastloc = degre(i));
\\        point = ratpoint2(pol,LIM3,1);
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\ Instead of solving directly y^2 = d1*x^4+c*x^2+d/d1,
\\ we solve first y^2 = d1*X^2+c*X+d/d1, then solve the quartic X = x^2
\\ which gives a new quartic
        qf = 2*[d1,c/2,0;c/2,d/d1,0;0,0,-1];
        solqf = qfsolve(qf);
        para = qfparam(qf,solqf,2)*['x^2,'x,1]~;
if( DEBUGLEVEL_ell >= 3, 
  print("   the conic y^2 = ",Pol([d1,c,d/d1]));
  print("   is parametrized by [x,y] = "subst([para[1]/para[2],para[3]/para[2]],'x,'t)));
        point1 = ratpoint2(para[1]*para[2],LIM3,1);
        if( point1 != [],
          if(#point1 == 2,
            para = subst(para,'x,point1[1])
          , point1 = [1,point1[2]/point1[1]^2,0];
            para = vector(3,ii,polcoeff(para[ii],2))
          );
          point = [point1[2]/para[2],para[3]/para[2]];
        , point = [];
        );
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
        if( point != [],
if( DEBUGLEVEL_ell >= 2, print("  point on the quartic"));
if( DEBUGLEVEL_ell >= 3, print(   point));
          m1++;
          listpointscount = concat(listpointscount,[d1*point[1]*point]);
          if( degre(i) > lastloc, m2++);
          found = 1; lastloc = -1
        ,
if( DEBUGLEVEL_ell >= 2, print("  no point found on the quartic"))
          ))));
    if( found,
      found = 0;
      v = 0; iaux = (i>>1);
      while( iaux, iaux >>= 1; v++);
      mask >>= 1;
      listgen = vecextract(listgen,(1<<#listgen)-(1<<v)-1);
      i = (1<<v)
    , i++)
  );

  for( i = 1, #listpointscount,
   if( #listpointscount[i] > 1,
    if( subst('x^3+c*'x^2+d*'x,'x,listpointscount[i][1])-listpointscount[i][2]^2 != 0,
      error(" ellcount: WRONG POINT, please report ",i))));
if( DEBUGLEVEL_ell >= 4, print("    end of ellcount"));

  return([listpointscount,[m1,m2]]);
}
{ell2descent_viaisog(ell,help=[]) =
\\ Computation of the rank of the elliptic curve ell
\\ having rational 2-torsion, using the algorithm via 2-isogenies.
\\
\\ ell must be on the form
\\ y^2=x^3+ax^2+bx -> ell = [0,a,0,b,0]
\\ with a and b integers.

my(P,Pfact,tors,listpointstriv,KS2prod,KS2gen,listpoints,pointgen,n1,n2,certain,apinit,bpinit,np1,np2,listpoints2,aux1,aux2,certainp,rang,strange);

if( DEBUGLEVEL_ell >= 2, print("  Algorithm of 2-descent via isogenies"));
  if( #ell < 13, ell = ellinit(ell,1));

  if( ell.disc == 0,
    error(" ell2descent_viaisog: singular curve !!"));
  if( ell.a1 != 0 || ell.a3 != 0 || ell.a6 != 0, 
    error(" ell2descent_viaisog: the curve is not on the form [0,a,0,b,0]"));
  if( denominator(ell.a2) > 1 || denominator(ell.a4) > 1,
    error(" ell2descent_viaisog: non-integral coefficients"));

\\
\\ Working with the initial curve
\\

\\ Construction of trivial points: torsion

  P = Pol([1,ell.a2,ell.a4]);
  Pfact = factor(P)[,1];
  tors = #Pfact;
  listpointstriv = concat(help,elltorseven(ell)[3]);

\\ Construction of trivial points: small naive height

if( DEBUGLEVEL_ell >= 3, print("   Search for trivial points on the curve"));
  P *= 'x;
if( DEBUGLEVEL_ell >= 3, print("   Y^2 = ",P));
  listpointstriv = concat( listpointstriv, ratpoint(P,LIMTRIV,0));
if( DEBUGLEVEL_ell >= 1, print(" trivial points on E(Q) = ",listpointstriv); print());

  KS2prod = -abs(ell.a4);
  if( ell.a2^2-4*ell.a4 < 0, KS2prod *= -1);
  KS2gen = factor(KS2prod)[,1];

if( DEBUGLEVEL_ell >= 2,
  print("  #K(b,2)gen          = ",#KS2gen);
  print("  K(b,2)gen = ",KS2gen));

  listpoints = ellcount(ell.a2,ell.a4,KS2gen,listpointstriv);
  pointgen = listpoints[1];
if( DEBUGLEVEL_ell >= 1, print(" points on E(Q) = ",pointgen); print());
  n1 = listpoints[2][1]; n2 = listpoints[2][2];

  certain = (n1 == n2);
if( DEBUGLEVEL_ell >= 1,
  if( certain, 
    print("[E(Q):phi'(E'(Q))]  = ",1<<n1);
    print("#S^(phi')(E'/Q)     = ",1<<n2);
    print("#III(E'/Q)[phi']    = 1"); print()
  ,
    print("[E(Q):phi'(E'(Q))] >= ",1<<n1);
    print("#S^(phi')(E'/Q)     = ",1<<n2);
    print("#III(E'/Q)[phi']   <= ",1<<(n2-n1)); print())
);

\\
\\ Working with the isogeneous curve
\\

  apinit = -2*ell.a2; bpinit = ell.a2^2-4*ell.a4;
  KS2prod = -abs(bpinit);
  if( ell.a4 < 0, KS2prod *= -1);
  KS2gen = factor(KS2prod)[,1];

if( DEBUGLEVEL_ell >= 2,
  print("  #K(a^2-4b,2)gen     = ",#KS2gen);
  print("  K(a^2-4b,2)gen     = ",KS2gen));

\\ Construction of trivial points: torsion

  P = Pol([1,apinit,bpinit]);
  listpointstriv = elltorseven([0,apinit,0,bpinit,0])[3];
   
\\ Construction of trivial points: small naive height

if( DEBUGLEVEL_ell >= 3, print("   Search for trivial points on the curve"));
  P *= 'x;
if( DEBUGLEVEL_ell >= 3, print(" Y^2 = ",P));
  listpointstriv = concat( listpointstriv, ratpoint(P,LIMTRIV,0));
if( DEBUGLEVEL_ell >= 1, print(" trivial points on E'(Q) = ",listpointstriv); print());

  listpoints = ellcount(apinit,bpinit,KS2gen,listpointstriv);

if( DEBUGLEVEL_ell >= 1, print(" points on E'(Q) = ",listpoints[1]));
  np1 = listpoints[2][1]; np2 = listpoints[2][2];
  listpoints2 = vector(#listpoints[1],i,0);
  for( i = 1, #listpoints[1],
    listpoints2[i] = [0,0];
    aux1 = listpoints[1][i][1]^2;
    if( aux1 != 0,
      aux2 = listpoints[1][i][2];
      listpoints2[i][1] = aux2^2/aux1/4;
      listpoints2[i][2] = aux2*(bpinit-aux1)/aux1/8
    , listpoints2[i] = listpoints[1][i]));
if( DEBUGLEVEL_ell >= 1, print(" points on E(Q) = ",listpoints2); print());
  pointgen = concat(pointgen,listpoints2);

  certainp = (np1 == np2);
if( DEBUGLEVEL_ell >= 1,
  if( certainp,
    print("[E'(Q):phi(E(Q))]   = ",1<<np1);
    print("#S^(phi)(E/Q)       = ",1<<np2);
    print("#III(E/Q)[phi]      = 1"); print()
  ,
    print("[E'(Q):phi(E(Q))]  >= ",1<<np1);
    print("#S^(phi)(E/Q)       = ",1<<np2);
    print("#III(E/Q)[phi]     <= ",1<<(np2-np1)); print());

  if( !certain && (np2 > np1), print1(1<<(np2-np1)," <= "));
  print1("#III(E/Q)[2]       ");
  if( certain && certainp, print1(" "), print1("<"));
  print("= ",1<<(n2+np2-n1-np1));

  print("#E(Q)[2]            = ",1<<tors);
);
  rang = n1+np1-2;
if( DEBUGLEVEL_ell >= 1,
  if( certain && certainp,
    print("#E(Q)/2E(Q)         = ",(1<<(rang+tors)));
    print("rank                = ",rang); print()
  ,
    print("#E(Q)/2E(Q)        >= ",(1<<(rang+tors))); print();
    print("",rang," <= rank          <= ",n2+np2-2); print()
  ));

  strange = (n2+np2-n1-np1)%2;
  if( strange, 
if( DEBUGLEVEL_ell >= 1,
      print(" !!! III should be a square !!!"); print("hence"));
    if( certain, 
      np1++;
      certainp = (np1 == np2);
if( DEBUGLEVEL_ell >= 1,
        if( certainp,
          print("[E'(Q):phi(E(Q))]   = ",1<<np1);
          print("#S^(phi)(E/Q)       = ",1<<np2);
          print("#III(E/Q)[phi]      = 1"); print()
        ,
          print("[E'(Q):phi(E(Q))]  >= ",1<<np1);
          print("#S^(phi)(E/Q)       = ",1<<np2);
          print("#III(E/Q)[phi]     <= ",1<<(np2-np1)); print())
      )
    ,
    if( certainp,
      n1++;
      certain = (n1 == n2);
if( DEBUGLEVEL_ell >= 1, 
        if( certain,
          print("[E(Q):phi'(E'(Q))]   = ",1<<n1);
          print("#S^(phi')(E'/Q)       = ",1<<n2);
          print("#III(E'/Q)[phi']      = 1"); print()
        ,
          print("[E(Q):phi'(E'(Q))]  >= ",1<<n1);
          print("#S^(phi')(E'/Q)      = ",1<<n2);
          print("#III(E'/Q)[phi']    <= ",1<<(n2-n1)); print())
      )
    , n1++)
  );

if( DEBUGLEVEL_ell >= 1,
  print("#S^(2)(E/Q)           = ",1<<(n2+np2-1));
  if( !certain && (np2 > np1), print1(" ",1<<(np2-np1)," <= "));
  print1("#III(E/Q)[2]       ");
  if( certain && certainp, print1(" "), print1("<"));
  print("= ",1<<(n2+np2-n1-np1));
  print("#E(Q)[2]            = ",1<<tors);
);
  rang = n1+np1-2;
if( DEBUGLEVEL_ell >= 1,
  if( certain && certainp, 
    print("#E(Q)/2E(Q)         = ",(1<<(rang+tors))); print();
    print("rank                = ",rang); print()
  ,
    print("#E(Q)/2E(Q)        >= ",(1<<(rang+tors))); print();
    print(rang," <= rank          <= ",n2+np2-2); print())
  ));

\\ end of strange 
  
  pointgen = vecsort(pointgen,,2);
  if( ELLREDGENFLAG, pointgen = ellredgen(ell,pointgen));
  pointgen = concat(ellsort(elltorseven(ell)[3]),pointgen);
if( DEBUGLEVEL_ell >= 1, print("points = ",pointgen));
if( DEBUGLEVEL_ell >= 3, print("   end of ell2descent_viaisog"));

  return([rang,n2+np2-2+tors,pointgen]);
}

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\            HELP MESSAGES                \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{
\\                  functions for elliptic curves
  addhelp(ell2descent_complete,
    "ell2descent_complete(e1,e2,e3): Performs a complete 2-descent on the elliptic curve y^2 = (x-e1)*(x-e2)*(x-e3). See ?ellrank for the format of the output.");
  addhelp(ell2descent_twist,
    "ell2descent_twist((E,bnf,k=1,help=[]): E is a vector of the form [0,A,0,B,C], (or the result of ellinit of such a vector) A,B,C integers such that x^3+A*x^2+B*x+C; bnf is the corresponding bnfinit(,1); Performs 2-descent on the elliptic curve Ek: k*y^2=x^3+A*x^2+B*x+C. See ?ellrank for the format of the output.");
  addhelp(ell2descent_viaisog,
    "ell2descent_viaisog(E,help=[]): E is an elliptic curve of the form [0,a,0,b,0], with a, b integers. Performs a 2-descent via isogeny on E. See ?ellrank for the format of the output.");
  addhelp(ellrank,
    "ellrank(E,help=[]): E is any elliptic curve defined over Q. Returns a vector [r,s,v], where r is a lower bound for the rank of E, s is the rank of its 2-Selmer group and v is a list of independant points in E(Q)/2E(Q). If help is a vector of nontrivial points on E, the result might be faster. This function might be used in conjunction with elltors2(E). See also ?default_ellQ");
  addhelp(ellhalf,
    "ellhalf(E,P): returns the vector of all points Q on the elliptic curve E such that 2Q = P");
  addhelp(ellredgen,
    "ellredgen(E,v): returns a vector of smallest possible points on the elliptic curve E generating the same subgroup as v, up to torsion.");
  addhelp(ellsort,
    "ellsort(v): v being a vector of points on some elliptic curve, returns the vector v sorted according to the naive height.");
  addhelp(elltors2,
    "elltors2(E): for an elliptic curve E, returns the group E(K)[2], where K is the field of definition of the coefficients of E (Q, R, Qp or Fp)."); 
  addhelp(elltorseven,
    "elltorseven(E): for an elliptic curve E, returns 2-Sylow subgroup of E(K)_tors, where K is the field of definition of the coefficients of E: (Q, R, Qp or Fp).");


\\                  functions for polynomials
  addhelp(locallysoluble,
    "locallysoluble(pol): returns 1 if y^2=pol(x) is everywhere locally soluble, 0 otherwise.");
  addhelp(ratpoint,
    "ratpoint(pol,lim=1,singlepoint=1): search for rational points on y^2=pol(x), for about within the bounds given by lim. The coefficients of pol must be integral. If singlepoint=1, returns at most one point, otherwise as many as possible.");
  addhelp(redquartic,
    "redquartic(pol): reduction of the quartic pol using Cremona-Stoll algorithm. Returns [p,M], where p is the reduced quartic and M is the GL2(Z) transformation. Also works with other degree polynomials.");


\\                  functions for number fields
  addhelp(bnfpSelmer,
    "bnfpSelmer(K,S,p): K being a number field given by bnfinit, S an ideal of K, and p a prime number, computes a set of generators of the group K(S,p) = { x in K^* / K^*^p, v_P(x) = 0 (mod p) for all P coprime to S}");
  addhelp(reducemodsquares,
    "reducemodsquares(delta,d): delta being a t_POLMOD, returns another delta'=delta*z^2, such that delta' has a small coefficient in x^d.");


\\                  others
  addhelp(default_ellQ,
    "default_ellQ(DEBUGLEVEL_ell, LIM1, LIM3, LIMTRIV, ELLREDGENFLAG, COMPLETE, MAXPROB, LIMBIGPRIME): set the value of the global variables used for ellrank() and other related functions. DEBUGLEVEL_ell: 0-5: choose the quantity of information printed during the computation (default=0: print nothing); LIM1 (resp LIM3): search limit for easy (resp hard) points on quartics; LIMTRIV: search limit for trivial points on elliptic curves; ELLREDGENFLAG: if != 0, try to reduce the generators at the end; COMPLETE: if != 0 and full 2-torsion, use complete 2-descent, otherwise via 2-isogeny; MAXPROB, LIMBIGPRIME: technical.");
/*  addhelp(DEBUGLEVEL_ell,
    "DEBUGLEVEL_ell: Choose a higher value of this global variable to have more details of the computations printed during the 2-descent algorithm. 0 = don't print anything; 1 = (default) just print the result; 2 = print more details including the Selmer group and the nontrivial quartics.");
*/
}

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\    FUNCTIONS FROM QFSOLVE NOT IN IN GP      \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

DEBUGLEVEL_qfsolve=0;

\\ Gives a unimodular matrix with the last column equal to v.
\\ redflag = 0 or 1. If redflag = 1, then the n-#v first columns are reduced.
{completebasis(v,redflag=0) =
my(Mv,U,re);
my(n);

  if( type(v) != "t_COL" && type(v) != "t_MAT",
    error("wrong type in completebasis"));

  if( type(v) == "t_COL", Mv = Mat(v), Mv = v);
  n = length(Mv[,1]);
  if( n == length(Mv), return(Mv));
  U = (mathnf(Mv~,1)[2]~)^-1;
  if( n==1 || !redflag, return(U));
\\ extract the n-#v columns and LLL-reduce them
  re = qflll(vecextract(U,1<<n-1,1<<(n-#Mv)-1));
\\ apply the reduction
  re = U*matdiagonalblock([re,matid(#Mv)]);

return(re);
}

\\ M = [a,b;b;c] has integral coefficients.
\\ Gauss reduction of the binary quadratic form
\\   qf = (a,b,c)=a*X^2+2*b*X*Y+c*Y^2
\\ Returns the reduction matrix with det = +1.
{qfbreduce(M) =
my(H,di,aux);
my(a,b,c,q,r,nexta,nextb,nextc);
my(test);

if( DEBUGLEVEL_qfsolve >= 5, print("     starting qfbreduce with ",M));

  a = M[1,1]; b = M[1,2]; c = M[2,2];

  H = matid(2); test = 1;
  while( test && a,
    di = divrem(b,a); q = di[1]; r = di[2];
    if( 2*r > abs(a), r -= abs(a); q += sign(a));
    H[,2] -= q*H[,1];
    nextc = a; nextb = -r; nexta= (nextb-b)*q+c;

    test = abs(nexta) < abs(a);
    if( test,
      c = nextc; b = nextb; a = nexta;
      aux = H[,1]; H[,1] = -H[,2]; H[,2] = aux
    )
  );

if( DEBUGLEVEL_qfsolve >= 5, print("     end of qfbreduce with ",H));
return(H);
}

\\ Performs first a LLL reduction on a positive definite
\\ quadratic form QD bounding the indefinite G.
\\ Then finishes the reduction with qfsolvetriv().
{qflllgram_indef(G,c=1,base=0) =
my(M,QD,M1,S,red);
my(n);

if( DEBUGLEVEL_qfsolve >= 4, print("    qflllgram_indef with G = ",G));
  n = length(G);
  M = matid(n);
  QD = G;
  for( i = 1, n-1,
    if( !QD[i,i],
      return(qfsolvetriv(G,base))
    );
    M1 = matid(n);
    M1[i,] = -QD[i,]/QD[i,i];
    M1[i,i] = 1;
    M = (M*M1);
    QD = (M1~*QD*M1)
  );
  M = (M^(-1));
  QD = (M~*abs(QD)*M);
  S = qflllgram(QD/content(QD));
  if( #S < n, S = completebasis(S));
  red = qfsolvetriv(S~*G*S,base);
  if( type(red) == "t_COL",
    red = (S*red);
    return(red));
  red[2] = S*red[2];
  if( length(red) == 3,
    red[3] = S*red[3]);
return(red);
}

\\ Try to solve G = 0 with small coefficients
\\ This is proved to work if
\\ -  det(G) = 1, dim <= 6 and G is LLL reduced
\\
\\ Returns [G,matid] if no solution is found.
\\ Exit with a norm 0 vector if one such is found.
\\ If base == 1 and norm 0 is obtained, returns [H~*G*H,H,sol] where
\\   sol is a norm 0 vector and is the 1st column of H.
{qfsolvetriv(G,base=0) =
my(n);
my(H,sol,GG);

  n = length(G);
  H = matid(n);

\\ case 1: A basis vector is isotropic
  for( i = 1, n,
    if( G[i,i] == 0,
      sol = H[,i];
      if( base == 0, return(sol));
      H[,i] = H[,1]; H[,1] = sol;
      return([H~*G*H,H,sol])
    )
  );

\\ case 2: G has a block +- [1,0;0,-1] on the diagonal
  for( i = 2, n,
    if( G[i-1,i] == 0 && G[i-1,i-1]*G[i,i] == -1,
      H[i-1,i] = -1; sol = H[,i];
      if( base == 0, return(sol));
      H[,i] = H[,1]; H[,1] = sol;
      return([H~*G*H,H,sol])
    )
  );

\\ case 3: a principal minor is 0
  for( i = 1, n,
    GG = vecextract(G,1<<i-1,1<<i-1);
    if( matdet(GG) != 0, next);
    sol = matker(GG)[,1];
    sol /= content(sol);
    sol = concat(sol,vectorv(n-i));
    if( base == 0, return(sol));
    H = completebasis(sol);
    H[,n] = -H[,1]; H[,1] = sol;
    return([H~*G*H,H,sol])
  );

return([G,H]);
}







\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\       Copyright (C) 2017 Denis Simon
\\
\\ Distributed under the terms of the GNU General Public License (GPL)
\\
\\    This code is distributed in the hope that it will be useful,
\\    but WITHOUT ANY WARRANTY; without even the implied warranty of
\\    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
\\    General Public License for more details.
\\
\\ The full text of the GPL is available at:
\\
\\                 http://www.gnu.org/licenses/
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

/*
  Auteur :
  Denis SIMON -> denis.simon@unicaen.fr
  adresse du fichier :
  www.math.unicaen.fr/~simon/ellcommon.gp

  *********************************************
  *          VERSION 10/01/2017               *
  *********************************************

*/

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\          SCRIPT                             \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
\\    COMMON FUNCTIONS TO ell.gp AND ellQ.gp   \\
\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\

{ellcomposeurst(urst1,urst2) =
my(u1 = urst1[1], r1 = urst1[2], s1 = urst1[3], t1 = urst1[4],
      u2 = urst2[1], r2 = urst2[2], s2 = urst2[3], t2 = urst2[4]);
  [u1*u2,u1^2*r2+r1,u1*s2+s1,u1^3*t2+s1*u1^2*r2+t1];
}
{ellinverturst(urst) =
my(u = urst[1], r = urst[2], s = urst[3], t = urst[4]);
  [1/u,-r/u^2,-s/u,(r*s-t)/u^3];
}
{mysubst(polsu,subsx) = 
  if( type(lift(polsu)) == "t_POL",
    return(simplify(subst(lift(polsu),variable(lift(polsu)),subsx)))
  , return(simplify(lift(polsu))));
}
{degre(idegre) =
my(ideg = idegre, jdeg = 0);

  while( ideg >>= 1, jdeg++);
  return(jdeg);
}
{nfrealsign(nf,a,i) =
\\ return the sign of the algebraic number a in the i-th real embedding.
my(nf_roots,ay,prec0);

  if( a == 0, return(0));

  a = lift(a);
  if( type(a) != "t_POL",
    return(sign(a)));

  nf_roots = nf.roots;
  prec0 = default(realprecision);

  ay = 0;
  while( ay == 0 || precision(ay) < 10,

    ay = subst(a,variable(a),nf_roots[i]);

    if( ay == 0 || precision(ay) < 10,
if( DEBUGLEVEL_ell >= 3, 
  print(" **** Warning: doubling the real precision in nfrealsign **** ",
        2*default(realprecision)));
      default(realprecision,2*default(realprecision));
      nf_roots = real(polroots(nf.pol))
    )
  );
  default(realprecision,prec0);

  return(sign(ay));
}
{nfsqrt(nf, a) =
\\ if a is a square in the number field nf returns [sqrt(a)], otherwise [].
my(alift,ta,minpola,py,pfact);

  if( a==0 || a==1, return([a]));

  alift = lift(a);
  ta = type(a);
  if( !poldegree(alift), alift = polcoeff(alift,0));

  if( type(alift) != "t_POL",
    if( issquare(alift), return([sqrtrat(alift)])));

  if( poldegree(nf.pol) <= 1, return([]));
  if( ta == "t_POL", a = Mod(a,nf.pol));

\\ the norm should be a square

  if( !issquare(norm(a)), return([]));

\\ the real embeddings must all be >0

  minpola = minpoly(a);
  if( polsturm(minpola,[-oo,0]), return([]));

\\ factorization over nf of the polynomial X^2-a

  if( variable(nf.pol) == 'x,
    py = subst(nf.pol,'x,'y);
    pfact = lift(factornf('x^2-mysubst(alift,Mod('y,py)),py)[1,1])
  ,
    pfact = lift(factornf('x^2-a,nf.pol)[1,1]));
  if( poldegree(pfact) == 2, return([]));
  return([subst(polcoeff(pfact,0),'y,Mod(variable(nf.pol),nf.pol))]);
}
{nfissquare(nf, a) = #nfsqrt(nf,a) > 0;
}
{sqrtrat(a) = 
  sqrtint(numerator(a))/sqrtint(denominator(a));
}








A(t) = -t^4+96*t^2-1024
B(t) = -2048*t^4+196608*t^2-4194304
C(t) =1024*t^8-196608*t^6+13631488*t^4-402653184*t^2+4294967296

S(E, N) =
  { my (s = 0.0);
    forprime(p = 2, N, my(a = ellap(E,p)); s += (2-a)/(p+1-a)*log(p));
    return (s);  }
{
  for(a=1, 10000,
    for(b=1, 10000, if(gcd(a,b)==1,
      my(t=a/b);
      my(Et=ellinit([0,A(t),0,B(t),C(t)]));
      if(Et!=[],
        my(E=ellminimalmodel(Et));
        my (v = S(E,50));
         my (x = S(E,100));
       my (b = S(E,200));
       my (d = S(E,400));
       my (e = S(E,600));
       my (f = S(E,1000));
       my (k= S(E,5000));
        if(v >3 && x>8 && b>10 && d>13 && e>16 && f>19 && k >27 , print(t,"   ",floor(e),"    ",floor(f) ,"    ",floor(k),"          ","rank >= ",floor(ellrank(E)[1]),"      2-Selmer rank =",floor(ellrank(E)[2]));
  )))))
}