Ticket #6955: trac_6955-simon-update-extcode.patch

File trac_6955-simon-update-extcode.patch, 59.6 KB (added by John Cremona, 13 years ago)

Apply in SAGE_ROOT/data/extcode before other patch

  • pari/simon/ell.gp

    # HG changeset patch
    # User John Cremona <john.cremona@gmail.com>
    # Date 1270319608 -3600
    # Node ID 499528b85ed0854fbf340df609d751282aff8aba
    # Parent  4efecfcb95fa3d95349589215e9204a89665a1fc
    #6955: update pari/simon scripts
    
    diff -r 4efecfcb95fa -r 499528b85ed0 pari/simon/ell.gp
    a b  
    2020\\ www.math.unicaen.fr/~simon/ell.gp
    2121\\
    2222\\  *********************************************
    23 \\  *          VERSION 13/11/2007               *
     23\\  *          VERSION 25/03/2009               *
    2424\\  *********************************************
    2525\\
    2626\\ Programme de calcul du rang des courbes elliptiques
     
    127127    return(simplify(subst(lift(polsu),variable(lift(polsu)),subsx)) )
    128128  , return(simplify(lift(polsu))));
    129129}
     130if( DEBUGLEVEL_ell >= 4, print("nfsign"));
     131{
     132nfsign(nf,a,i) =
     133\\ return the sign of the algebraic number a in the i-th real embedding.
     134local(nf_roots,ay,def);
     135
     136  if( a == 0, return(0));
     137
     138  a = lift(a);
     139  if( type(a) != "t_POL",
     140    return(sign(a)));
     141
     142  nf_roots = nf.roots;
     143  def = default(realprecision);
     144
     145  ay = 0;
     146  while( ay == 0 || precision(ay) < 10,
     147
     148    ay = subst(a,variable(a),nf_roots[i]);
     149
     150    if( ay == 0 || precision(ay) < 10,
     151if( DEBUGLEVEL_ell >= 3,
     152  print(" **** Warning: doubling the real precision in nfsign **** ",
     153        2*default(realprecision)));
     154      default(realprecision,2*default(realprecision));
     155      nf_roots = real(polroots(nf.pol))
     156    )
     157  );
     158  default(realprecision,def);
     159
     160  return(sign(ay));
     161}
    130162if( DEBUGLEVEL_ell >= 4, print("degre"));
    131163{
    132164degre(idegre) =
    133 local(ideg,jdeg);
     165local(ideg = idegre, jdeg = 0);
    134166
    135   ideg = idegre; jdeg = 0;
    136167  while( ideg >>= 1, jdeg++);
    137168  return(jdeg);
    138169}
     
    144175{
    145176nfsqrt( nf, a) =
    146177\\ si a est un carre, renvoie [sqrt(a)], sinon [].
    147 local(alift,ta,res,pfact,r1,py);
     178local(alift,ta,res,pfact);
    148179
    149180if( DEBUGLEVEL_ell >= 5, print("entree dans nfsqrt ",a));
    150181  if( a==0 || a==1,
     
    166197  if( ta == "t_POL", a = Mod(a,nf.pol));
    167198
    168199\\ tous les plgements reels doivent etre >0
    169 \\   
    170   r1 = nf.sign[1];
    171   for( i = 1, r1,
    172     py = mysubst(alift,nf.roots[i]);
    173     if( sign(py) < 0,
     200
     201  for( i = 1, nf.r1,
     202    if( nfsign(nf,a,i) < 0,
    174203if( DEBUGLEVEL_ell >= 5, print("fin de nfsqrt"));
    175204      return([])));
     205
    176206\\ factorisation sur K du polynome X^2-a :
     207
    177208  if( variable(nf.pol) == x,
    178209    py = subst(nf.pol,x,y);
    179210    pfact = lift(factornf(x^2-mysubst(alift,Mod(y,py)),py)[1,1])
     
    191222  sqrtint(numerator(a))/sqrtint(denominator(a));
    192223}
    193224
     225
    194226\\
    195227\\ Fonctions propres a ell.gp
    196228\\
     
    213245\\ ideal doit etre sous la forme primedec
    214246  if( #nf.zk == 1,
    215247if( DEBUGLEVEL_ell >= 5, print("fin de nfmodid2"));
    216     return(a*Mod(1,ideal[1])));
     248    return(a*Mod(1,ideal.p)));
    217249  a = mynfeltmod(nf,a,nfbasistoalg(nf,ideal[2]));
    218   if( gcd(denominator(content(lift(a))),ideal[1]) == 1,
     250  if( gcd(denominator(content(lift(a))),ideal.p) == 1,
    219251if( DEBUGLEVEL_ell >= 5, print("fin de nfmodid2"));
    220     return(a*Mod(1,ideal[1])));
     252    return(a*Mod(1,ideal.p)));
    221253if( DEBUGLEVEL_ell >= 5, print("fin de nfmodid2"));
    222254  return(a);
    223255}
     
    241273local(alpha,beta,sig,aux,aux2,rep);
    242274
    243275if( DEBUGLEVEL_ell >= 5, print("entree dans mynfhilbertp ",p));
    244   if( a==0 || b==0, print("0 argument in mynfhilbertp"));
    245   if( p[1] == 2,
     276  if( a == 0 || b == 0, print("0 argument in mynfhilbertp"));
     277  if( p.p == 2,
    246278if( DEBUGLEVEL_ell >= 5, print("fin de mynfhilbertp"));
    247279    return(nfhilb2(nf,a,b,p)));
    248280  if( type(a) != "t_POLMOD", a = Mod(a,nf.pol));
     
    275307    for( j = 1, #S1, k = #Slist;
    276308      for( k = 1, #Slist,
    277309        if( Slist[k] == S1[j], test = 0; break));
    278       if( test, Slist=concat(Slist,[S1[j]]), test = 1);
     310      if( test, Slist = concat(Slist,[S1[j]]), test = 1);
    279311     ));
    280312if( DEBUGLEVEL_ell >= 5, print("fin de ideallistfactor"));
    281313  return(Slist);
     
    287319\\ =1 si l'equation X^2-aY^2-bZ^2=0 a une solution non triviale,
    288320\\ =-1 sinon,
    289321\\ a et b doivent etre non nuls.
    290 local(al,bl,r1,i,S);
     322local(al,bl,S);
    291323
    292324if( DEBUGLEVEL_ell >= 4, print("entree dans mynfhilbert ",[a,b]));
    293325  if( a == 0 || b == 0, error("mynfhilbert : argument = 0"));
    294326  al = lift(a); bl = lift(b);
    295327
    296 \\ solutions locales aux places infinies reelles
     328\\ solutions locales aux places reelles
    297329
    298   r1 = nf.sign[1];
    299   for( i = 1, r1,
    300     if( sign(mysubst(al,nf.roots[i])) < 0,
    301       if( sign(mysubst(bl,nf.roots[i])) < 0,
     330  for( i = 1, nf.r1,
     331    if( nfsign(nf,al,i) < 0 && nfsign(nf,bl,i) < 0,
    302332if( DEBUGLEVEL_ell >= 3, print("mynfhilbert non soluble a l'infini"));
    303333if( DEBUGLEVEL_ell >= 4, print("fin de mynfhilbert"));
    304         return(-1))));
     334      return(-1))
     335  );
    305336
    306337  if( type(a) != "t_POLMOD", a = Mod(a,nf.pol));
    307338  if( type(b) != "t_POLMOD", b = Mod(b,nf.pol));
    308339
    309340\\  solutions locales aux places finies (celles qui divisent 2ab)
    310341
    311   S = ideallistfactor (nf,[2,a,b]);
     342  S = ideallistfactor(nf,[2,a,b]);
    312343  forstep ( i = #S, 2, -1,
    313344\\ d'apres la formule du produit on peut eviter un premier
    314345    if( mynfhilbertp(nf,a,b, S[i]) == -1,
     
    336367  pp=[ p, nfbasistoalg(nf,p[2]), idval, 0, repres(nf,p) ];
    337368  if( idval,
    338369    pp[4] = idealstar(nf,idealpow(nf,p,1+2*idval)),
    339     pp[4] = p[1]^p[4]\2 );
     370    pp[4] = p.p^p.f\2 );
    340371if( DEBUGLEVEL_ell >= 5, print("fin de initp"));
    341372  return(pp);
    342373}
     
    412443  for( i = 1, #mat,
    413444    if( mat[i,i] != 1, fond = concat(fond,nf.zk[i])));
    414445  f = #fond;
    415   pp = p[1];
     446  pp = p.p;
    416447  rep = vector(pp^f,i,0);
    417448  rep[1] = 0;
    418449  ppi = 1;
     
    449480    if( f > q, aaa = nfbasistoalg(nf,p[2])^((q+1)>>1), aaa = 0);
    450481if( DEBUGLEVEL_ell >= 4, print("fin de nfissquarep"));
    451482    return(aaa));
    452   if( f, aaa = a*nfbasistoalg(nf,p[5]/p[1])^f, aaa = a);
     483  if( f, aaa = a*nfbasistoalg(nf,p[5]/p.p)^f, aaa = a);
    453484  if( pherm[1,1] != 2,
    454485\\ cas ou p ne divise pas 2
    455486\\ algorithme de Shanks
    456487    n = nfrandintmodid(nf,pherm);
    457488    while( nfpsquareodd(nf,n,p), n = nfrandintmodid(nf,pherm));
    458     pp = Mod(1,p[1]);
     489    pp = Mod(1,p.p);
    459490    n *= pp;
    460491    qq = idealnorm(nf,pherm)\2;
    461492    e = 1; while( !(qq%2), e++; qq \= 2);
     
    499530    for( i = 1, #zlog,
    500531      expo = zlog[i];
    501532      if( expo,
    502         if( !expo%2, expo = expo>>1, aux = zinit[2][i]; expo = expo*((aux+1)>>1)%aux);
    503         if( expo == 1,
    504           xx *= nfbasistoalg(nf,(zinit[2][3][i]))
    505         , xx *= nfbasistoalg(nf,(zinit[2][3][i]))^expo);
     533        if( !expo%2,
     534          expo = expo>>1
     535        , aux = zinit[2][i];
     536          expo = expo*((aux+1)>>1)%aux
     537        );
     538        xx *= nfbasistoalg(nf,zinit[2][3][i])^expo
    506539      )
    507540    );
    508     if( f, xx *= nfbasistoalg(nf,p[2])^(f>>1); id = idealpow(nf,p,q));
     541    if( f,
     542      xx *= nfbasistoalg(nf,p[2])^(f>>1);
     543      id = idealpow(nf,p,q));
    509544    xx = mynfeltreduce(nf,xx,id);
    510545  );
    511546if( DEBUGLEVEL_ell >= 4, print("fin de nfissquarep ",xx));
     
    526561  if( v%2,
    527562if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquareodd"));
    528563    return(0));
    529   ap = a*(1/nfbasistoalg(nf,p[2])^v);
     564  ap = a/nfbasistoalg(nf,p[2])^v;
    530565
    531566  norme = idealnorm(nf,p)\2;
    532   den = denominator(content(lift(ap)))%p[1];
    533   if(sign(den), ap*=Mod(1,p[1]));
     567  den = denominator(content(lift(ap)))%p.p;
     568  if(sign(den), ap*=Mod(1,p.p));
    534569  ap = ap^norme-1;
    535570  if( ap == 0,
    536571if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquareodd"));
     
    554589if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquare"));
    555590    return(1));
    556591
    557   if( (p[1]%2),
     592  if( p.p != 2,
    558593if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquare"));
    559594    return(nfpsquareodd(nf,a,p)));
    560595
     
    563598if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquare"));
    564599    return(0));
    565600  if( valap,
    566     zlog = ideallog(nf,a*(nfbasistoalg(nf,p[5])/p[1])^valap,zinit)
     601    zlog = ideallog(nf,a*(nfbasistoalg(nf,p[5])/p.p)^valap,zinit)
    567602  ,
    568603    zlog = ideallog(nf,a,zinit));
    569604  for( i = 1, #zinit[2][2],
     
    644679if( DEBUGLEVEL_ell >= 5, print("fin de nflemma7"));
    645680      return(-1));
    646681    q = mu+nu-lambda;
    647     if( q>2*v,
     682    if( q > 2*v,
    648683if( DEBUGLEVEL_ell >= 5, print("fin de nflemma7"));
    649684      return(-1));
    650685    if( nfpsquareq(nf,gx*nfbasistoalg(nf,p[5]/2)^lambda,p,q),
    651686if( DEBUGLEVEL_ell >= 5, print("fin de nflemma7"));
    652687      return(1))
    653688  ,
    654     if( lambda>= 2*nu,
     689    if( lambda >= 2*nu,
    655690if( DEBUGLEVEL_ell >= 5, print("fin de nflemma7"));
    656691      return(0));
    657692    if( lambda%2,
     
    739774nfqpsolublebig( nf, pol, p,ap=0,b=1) =
    740775local(deg,i,xx,z,Px,j,cont,pi,pol2,Roots);
    741776
    742 if( DEBUGLEVEL_ell >= 4, print("entree dans nfqpsolublebig avec ",p[1]));
     777if( DEBUGLEVEL_ell >= 4, print("entree dans nfqpsolublebig avec ",p.p));
    743778  deg = poldegree(pol);
    744779
    745780  if( nfpsquareodd(nf,polcoeff(pol,0),p),
     
    753788  cont = idealval(nf,polcoeff(pol,0),p);
    754789  for( i = 1, deg,
    755790    if( cont, cont = min(cont,idealval(nf,polcoeff(pol,i),p))));
    756   if( cont, pi = nfbasistoalg(nf,p[5]/p[1]));
     791  if( cont, pi = nfbasistoalg(nf,p[5]/p.p));
    757792  if( cont > 1, pol *= pi^(2*(cont\2)));
    758793   
    759794\\ On essaye des valeurs de x au hasard
     
    790825\\ p est un ideal premier de nf, sous la forme idealprimedec
    791826local(factlist,sol);
    792827
    793   factlist = nffactormod(nf,pol,p);
    794 \\ CETTE LIGNE NE DOIT PAS RESTER TRES LONGTEMPS
    795   if( type(factlist) == "t_VEC", factlist = factlist[1], factlist = factlist[,1]);
     828  factlist = nffactormod(nf,pol,p)[,1];
    796829  sol = [];
    797830  for( i = 1, #factlist,
    798831    if( poldegree(factlist[i]) == 1,
     
    823856if( DEBUGLEVEL_ell >= 4, print("nflocallysoluble"));
    824857{
    825858nflocallysoluble( nf, pol, r=0,a=1,b=1) =
    826 local(pol0,plist,add,ff,p,r1,Delta,vecpol,ar,er,Deltar,vecpolr,Sturmr);
     859local(pol0,plist,add,ff,p,Delta,vecpol,vecpolr,Sturmr);
    827860
    828861if( DEBUGLEVEL_ell >= 4, print("entree dans nflocallysoluble ",[pol,r,a,b]));
    829862  pol0 = pol;
     
    844877      for( i = 1, #plist,
    845878        p =  plist[i]; 
    846879if( DEBUGLEVEL_ell >= 3, print("p = ",p));
    847         if( p[1] < LIMBIGPRIME,
     880        if( p.p < LIMBIGPRIME,
    848881          if( !nfqpsoluble(nf,pol,initp(nf,p)),
    849882if( DEBUGLEVEL_ell >= 2, print(" non ELS en ",p));
    850883if( DEBUGLEVEL_ell >= 4, print("fin de nflocallysoluble"));
    851884            return(0)),
    852885          if( !nfqpsolublebig(nf,pol,p,r/a,b),
    853 if( DEBUGLEVEL_ell >= 2, print(" non ELS en ",p[1]," ( = grand premier  )"));
     886if( DEBUGLEVEL_ell >= 2, print(" non ELS en ",p.p," ( = grand premier  )"));
    854887if( DEBUGLEVEL_ell >= 4, print("fin de nflocallysoluble"));
    855888            return(0))));
    856889);
    857890\\ places reelles
    858   r1 = nf.sign[1];
    859   if( r1,
     891  if( nf.r1,
    860892    Delta = poldisc(pol); vecpol = Vec(pol);
    861     for( i = 1, r1,
    862       ar = mysubst(pollead(pol),nf.roots[i]);
    863       if( ar > 0, next);
    864       er = mysubst(polcoeff(pol,0),nf.roots[i]);
    865       if( er > 0, next);
    866       Deltar = mysubst(Delta,nf.roots[i]);
    867       if( Deltar < 0, next);
    868       vecpolr = vector(poldegree(pol)+1,j,mysubst(vecpol[j],nf.roots[i]));
     893    for( i = 1, nf.r1,
     894      if( nfsign(nf,pollead(pol),i) > 0, next);
     895      if( nfsign(nf,polcoeff(pol,0),i) > 0, next);
     896      if( nfsign(nf,Delta,i) < 0, next);
     897      vecpolr = vector(#vecpol,j,mysubst(vecpol[j],nf.roots[i]));
    869898      Sturmr = polsturm(Pol(vecpolr));
    870899      if( Sturmr == 0,
    871900if( DEBUGLEVEL_ell >= 2, print(" non ELS a l'infini"));
     
    10051034  oddclass = 0;
    10061035  while( !oddclass,
    10071036    KS2gen = bnfsunit(bnf,idealfactor(bnf,KS2prod)[,1]~);
    1008     oddclass = (KS2gen[5][1]%2);
     1037    oddclass = KS2gen[5][1]%2;
    10091038    if( !oddclass,
    10101039      KS2prod = idealmul(bnf,KS2prod,(KS2gen[5][3][1])));
    10111040 );
    10121041  KS2gen = KS2gen[1];
     1042\\  A CHANGER : KS2gen = matbasistoalg(bnf,KS2gen);
    10131043  for( i = 1, #KS2gen,
    10141044    KS2gen[i] = nfbasistoalg(bnf, KS2gen[i]));
    10151045  KS2gen = concat(Mod(lift(bnf.tufu),bnf.pol),KS2gen);
     
    10421072    if( !oddclass,
    10431073      KS2prod = idealmul(bnf,KS2prod,(KS2gen[5][3][1]))));
    10441074  KS2gen = KS2gen[1];
     1075\\ A CHANGER KS2gen = matbasistoalg(bnf,KS2gen);
    10451076  for( i = 1, #KS2gen,
    10461077    KS2gen[i] = nfbasistoalg(bnf, KS2gen[i]));
    10471078  KS2gen = concat(Mod(lift(bnf.tufu),bnf.pol),KS2gen);
     
    11701201
    11711202if( DEBUGLEVEL_ell >= 4, print("entree dans nfchinremain"));
    11721203  l = #fact[,1];
    1173   fact2 = vector(l,i,0);
    1174   for( i = 1, l,
    1175     fact2[i] = idealdiv(nf,b,idealpow(nf,fact[i,1],fact[i,2])));
     1204  fact2 = vector(l,i,idealdiv(nf,b,idealpow(nf,fact[i,1],fact[i,2])));
     1205\\  for( i = 1, l,
     1206\\    fact2[i] = idealdiv(nf,b,idealpow(nf,fact[i,1],fact[i,2])));
    11761207  fact2 = idealaddtoone(nf,fact2);
     1208\\ A CHANGER : fact2 = matbasistoalg(nf,fact2);
    11771209  for( i = 1, l,
    11781210    fact2[i] = nfbasistoalg(nf,fact2[i]));
    11791211if( DEBUGLEVEL_ell >= 4, print("fin de nfchinremain"));
     
    12121244if( DEBUGLEVEL_ell >= 4, print("bbbnf.clgp = ",bbbnf.clgp));
    12131245  for( i = 1, #bbbnf.clgp[2],
    12141246    if( bbbnf.clgp[2][i]%2 == 0,
    1215       SL0 = idealmul(bbbnf,SL0,bbbnf.clgp[3][i])));
     1247      SL0 = idealmul(bbbnf,SL0,bbbnf.clgp[3][i][1,1])));
    12161248  SL1 = idealmul(bbbnf,SL0,rnfeltup(rrrnf,bleg));
    12171249  SL = idealfactor(bbbnf,SL1)[,1]~;
    12181250  sunL = bnfsunit(bbbnf,SL);
    1219   fondsunL = concat(bbbnf.futu,nfbasistoalg(bbbnf,sunL[1]));
     1251\\ A CHANGER : fondsunL = concat(bbbnf.futu,matbasistoalg(bbbnf,sunL[1]));
     1252  fondsunL = concat(bbbnf.futu,vector(#sunL[1],i,nfbasistoalg(bbbnf,sunL[1][i])));
    12201253  normfondsunL = norm(rnfeltabstorel( rrrnf,fondsunL));
    12211254  SK = idealfactor(bnf,idealnorm(bbbnf,SL1))[,1]~;
    12221255  sunK = bnfsunit(bnf,SK);
    1223   fondsunK = concat(bnf.futu,nfbasistoalg(bnf,sunK[1]));
     1256\\ A CHANGER :  fondsunK = concat(bnf.futu,matbasistoalg(bnf,sunK[1]));
     1257  fondsunK = concat(bnf.futu,vector(#sunK[1],i,nfbasistoalg(bnf,sunK[1][i])));
    12241258  vecbleg = bnfissunit(bnf,sunK,bleg);
    12251259  matnorm = matrix(#fondsunK,#normfondsunL,i,j,0);
    12261260  for( i = 1, #normfondsunL,
     
    15271561if( DEBUGLEVEL_ell >= 3, print("KS2gen = ",KS2gen[1]));
    15281562
    15291563  LS2genunit = lift(Lrnf.futu);
    1530   LS2genunit = concat(LS2genunit,lift(nfbasistoalg(Lrnf,LS2gen[1])));
     1564\\ A CHANGER : LS2genunit = concat(LS2genunit,lift(matbasistoalg(Lrnf,LS2gen[1])));
     1565  LS2genunit = concat(LS2genunit,vector(#LS2gen[1],i,lift(nfbasistoalg(Lrnf,LS2gen[1][i]))));
     1566
    15311567
    15321568  LS2genunit = subst(LS2genunit,x,ttheta);
    15331569  LS2genunit = LS2genunit*Mod(1,polrel);
     
    17221758        found = ispointtriv = 1
    17231759      )
    17241760    );
    1725    
     1761
    17261762\\ \\\\\\\\\\\
    17271763\\ Construction de la quartique
    17281764\\ \\\\\\\\\\\
     
    18471883      rang = m1+1)
    18481884  );
    18491885if( DEBUGLEVEL_ell >= 1, print("listpointsmwr = ",listpointsmwr));
    1850   for( i = 1, #listpointsmwr,
     1886  for( i = 1, #listpointsmwr,
     1887    if( #listpointsmwr[i] == 3,
     1888      listpointsmwr[i] = vecextract(listpointsmwr[i],3));
    18511889    if( !ellisoncurve(ellnf,listpointsmwr[i]),
    18521890      error("bnfell2descent : MAUVAIS POINT ")));
    18531891if( DEBUGLEVEL_ell >= 4, print("fin de bnfell2descent_gen"));
     
    18621900\\ attention bnf a un polynome en y.
    18631901\\ si bigflag !=0, on reduit les quartiques
    18641902\\ si flag3 != 0, on utilise bnfqfsolve2
    1865 local(urst,urst1,den,eqtheta,rnfeq,bbnf,ext,rang);
     1903local(urst,urst1,den,factden,eqtheta,rnfeq,bbnf,ext,rang);
    18661904
    18671905if( DEBUGLEVEL_ell >= 3, print("entree dans bnfellrank"));
    18681906  if( #ell <= 5, ell = ellinit(ell,1));
     
    18771915
    18781916\\ removes denominators
    18791917  while( (den = idealinv(bnf,idealadd(bnf,idealadd(bnf,1,ell.a2),idealadd(bnf,ell.a4,ell.a6))))[1,1] > 1,
    1880     den = idealfactor(bnf,den); den[,2] = vectorv(#den[,2],i,1);
    1881     den = factorback(den,bnf)[1,1];
     1918    factden = idealfactor(bnf,den)[,1];
     1919    den = 1;     
     1920    for( i = 1, #factden,
     1921      den = idealmul(bnf,den,factden[i]));
     1922    den = den[1,1];
    18821923    urst1 = [1/den,0,0,0];
    18831924    ell = ellchangecurve(ell,urst1);
    18841925    urst = ellcomposeurst(urst,urst1);
     
    19521993      KS2prod = idealmul(bnf,KS2prod,(KS2gen[5][3][1])));
    19531994  );
    19541995  KS2gen = KS2gen[1];
     1996\\ A CHANGER : KS2gen = matbasistoalg(bnf,KS2gen);
    19551997  for( i = 1, #KS2gen,
    19561998    KS2gen[i] = nfbasistoalg(bnf, KS2gen[i]));
    19571999  KS2gen = concat(Mod(lift(bnf.tufu),bnf.pol),KS2gen);
    1958 if( DEBUGLEVEL_ell >= 2, 
     2000if( DEBUGLEVEL_ell >= 2,
    19592001  print("#K(S,2)gen          = ",#KS2gen);
    19602002  print("K(S,2)gen = ",KS2gen));
    19612003
  • pari/simon/ellQ.gp

    diff -r 4efecfcb95fa -r 499528b85ed0 pari/simon/ellQ.gp
    a b  
    11\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
    2 \\       Copyright (C) 2007 Denis Simon
     2\\       Copyright (C) 2008 Denis Simon
    33\\
    44\\ Distributed under the terms of the GNU General Public License (GPL)
    55\\
     
    2020\\ www.math.unicaen.fr/~simon/ellQ.gp
    2121\\
    2222\\  *********************************************
    23 \\  *          VERSION 13/11/2007               *
     23\\  *          VERSION 29/04/2008               *
    2424\\  *********************************************
    2525\\
    2626\\ Programme de calcul du rang des courbes elliptiques sur Q.
     
    8383\\ Variables globales usuelles
    8484\\
    8585
    86   DEBUGLEVEL_ell = 2; \\ pour avoir plus ou moins de details
     86  DEBUGLEVEL_ell = 1; \\ pour avoir plus ou moins de details
    8787  LIM1 = 5;       \\ limite des points triviaux sur les quartiques
    8888  LIM3 = 50;      \\ limite des points sur les quartiques ELS
    89   LIMTRIV = 10;   \\ limite des points triviaux sur la courbe elliptique
     89  LIMTRIV = 50;   \\ limite des points triviaux sur la courbe elliptique
    9090
    9191\\
    9292\\  Variables globales techniques
     
    224224\\ Fonctions propres a ellQ.gp
    225225\\
    226226
     227if( DEBUGLEVEL_ell >= 4, print("ellhalf"));
     228{
     229ellhalf(ell,P)=
     230\\ renvoie tous les points Q sur ell tels que 2Q = P.
     231local(pol2,ratroots,half,x2,y2,P2);
     232
     233  if(#ell < 13, ell=ellinit(ell,1));
     234
     235  pol2 = Pol([4,ell.b2,2*ell.b4,ell.b6]); \\ polynome de 2-division
     236
     237  if( P == [0],
     238    ratroots = polratroots(pol2);
     239    half = vector(#ratroots,i,[ratroots[i],ellordinate(ell,ratroots[i])[1]]);
     240    half = concat( [[0]], half);
     241    return(half)
     242  );
     243
     244  x2=Pol([1,0,-ell.b4,-2*ell.b6,-ell.b8]); \\ x(2P) = x2/pol2
     245
     246  half = [];
     247  ratroots = polratroots(x2-P[1]*pol2); 
     248  if( #ratroots == 0, return(half));
     249  for( i = 1, #ratroots,
     250    y2 = ellordinate(ell,ratroots[i]);
     251    for( j = 1, #y2,
     252      P2 = [ratroots[i],y2[j]];
     253      if( ellpow(ell,P2,2) == P, half = concat(half,[P2]))
     254    )
     255  );
     256  return(half);
     257}
     258if( DEBUGLEVEL_ell >= 4, print("elltors2"));
     259{
     260elltors2(ell)=
     261\\ Calcule le sous-groupe de 2-torsion de la courbe elliptique ell.
     262local(pol2,ratroots,tors2);
     263
     264if( DEBUGLEVEL_ell >= 4, print("calcul de la 2-torsion"));
     265  if(#ell < 13, ell=ellinit(ell,1));
     266  tors2 = ellhalf(ell,[0]);
     267  if( #tors2 == 1,
     268    tors2 = [1, [], []],
     269  if( #tors2 == 2,
     270    tors2 = [2, [2], [tors2[2]]]
     271  , tors2 = [4, [2,2], [tors2[2],tors2[3]]]
     272  ));
     273if( DEBUGLEVEL_ell >= 4, print("E[2] = ",tors2));
     274  return(tors2);
     275}
     276if( DEBUGLEVEL_ell >= 4, print("elltorseven"));
     277{
     278elltorseven(ell)=
     279\\ Calcule le 2-Sylow sous-groupe de torsion de la courbe elliptique ell.
     280local(torseven,P2);
     281
     282if( DEBUGLEVEL_ell >= 4, print("calcul de la 2^n-torsion"));
     283  if(#ell < 13, ell=ellinit(ell,1));
     284  torseven = elltors2(ell);
     285
     286  while( torseven[1] != 1,
     287    P2 = ellhalf(ell,torseven[3][1]);
     288    if( #P2 > 0,
     289       torseven[1] *= 2;
     290       torseven[2][1] *= 2;
     291       torseven[3][1] = P2[1];
     292       next
     293    );
     294    if( #torseven[3] == 1, break());
     295
     296    P2 = ellhalf(ell,torseven[3][2]);
     297    if( #P2 > 0,
     298       torseven[1] *= 2;
     299       torseven[2][2] *= 2;
     300       torseven[3][2] = P2[1];
     301       next
     302    );
     303    P2 = ellhalf(ell,elladd(ell,torseven[3][1],torseven[3][2]));
     304    if( #P2 > 0,
     305       torseven[1] *= 2;
     306       torseven[2][1] *= 2;
     307       torseven[3][1] = P2[1];
     308       next
     309    );
     310    break()
     311  );
     312 
     313if( DEBUGLEVEL_ell >= 4, print("E[2^n] = ",torseven));
     314  return(torseven);
     315}
    227316if( DEBUGLEVEL_ell >= 4, print("polratroots"));
    228317{
    229318polratroots(pol) =
     
    237326}
    238327if( DEBUGLEVEL_ell >= 4, print("ratpoint"));
    239328{
    240 ratpoint(pol,lim,singlepoint=1,tryhard=0) =
     329ratpoint(pol,lim=1,singlepoint=1,tryhard=0) =
    241330\\ Recherche de points sur y^2=pol(x).
    242331\\ Les coeff de pol sont entiers.
    243332\\ Si singlepoint >= 1, cherche un seul point, sinon plusieurs.
     
    286375        tab5[xx+1,zz+1] = !issquare([xx^4,xx^3*zz,xx^2*zz^2,xx*zz^3,zz^4]*pol5)))
    287376  );
    288377 
     378  lead = pollead(pol);
     379  pol0 = polcoeff(pol,0);
     380
    289381  if( odd,
    290382    vecz = vector(lim,i,i^2);
    291     vecx = vector(lim,i,i)
    292383  ,
    293 \\ si le degre de pol est pair, il faut que le coeff constant soit
    294 \\ un carre mod xx. Idem pour le dominant mod zz.
     384\\ si le degre de pol est pair, il faut que le coeff dominant soit
     385\\ un carre mod zz.
    295386    vecz = vector(lim);
    296     vecx = vector(lim);
    297     lead = pollead(pol);
    298     pol0 = polcoeff(pol,0);
    299     zz = 0; xx = 0;
     387    zz = 0;
    300388    for( i = 1, lim,
    301       xx++; while( !issquare(Mod(pol0,xx)),xx++); vecx[i] = xx;
    302389      zz++; while( !issquare(Mod(lead,zz)),zz++); vecz[i] = zz
    303390  ));
     391\\ le coeff constant doit etre un carre mod xx.
     392  vecx = vector(lim);
     393  xx = 0;
     394  for( i = 1, lim,
     395    xx++; while( !issquare(Mod(pol0,xx)),xx++); vecx[i] = xx);
    304396
    305 if( DEBUGLEVEL_ell >= 4, print("xmax = ",xx));
    306 if( DEBUGLEVEL_ell >= 4, print("zmax = ",zz));
     397if( DEBUGLEVEL_ell >= 4, print("xmax = ",vecx[lim]));
     398if( DEBUGLEVEL_ell >= 4, print("zmax = ",vecz[lim]));
    307399
    308400if( DEBUGLEVEL_ell >= 5, print("vecx = ",vecx));
    309401if( DEBUGLEVEL_ell >= 5, print("vecz = ",vecz));
     
    313405    for( ix = max(1,somme-lim), min(lim,somme-1),
    314406      xx = vecx[ix]; iz = somme-ix; zz = vecz[iz];
    315407      if( gcd(zz,xx) > 1, next);
     408      if( odd && !issquare(lead*Mod(xx,zz)), next);
    316409      for( eps = 1, 2, if( eps == 2, zz = -zz);
    317410      if( deg4 &&
    318411        (tab16[xx%16+1,zz%16+1] || tab9[xx%9+1,zz%9+1] || tab5[xx%5+1,zz%5+1])
     
    331424  );
    332425
    333426\\
    334 \\ Try another strategy when pol has a nontrivial content
     427\\ Essaye une autre strategie quand pol a un content non trivial
    335428\\
    336429
    337430  if( !odd && tryhard,
    338 if( DEBUGLEVEL_ell >= 4, print(" Try hard *************"));
     431if( DEBUGLEVEL_ell >= 4, print(" Autre strategie dans ratpoint **********"));
    339432    K = content(pol);
    340433    if( K != 1,
    341434      pol /= K;
     
    347440        e = factK[,2]%2;
    348441        ind = #e; while( !e[ind], ind--);
    349442        p = factK[ind,1];
    350         if( valuation( pollead(pol), p) >= 1
    351          && valuation( pollead(pol), p) <= 2,
     443        if( valuation( pollead(pol), p) == 1 ||
     444            ( valuation( pollead(pol), p) >= 2 && valuation( polcoeff(pol,poldegree(pol)-1), p) == 0),
    352445if( DEBUGLEVEL_ell >= 4, print(" utilise une racine de pol mod p = ",p));
    353446          sol = ratpoint(K/p^2*subst(polrecip(pol),variable(pol),p*variable(pol)),lim,singlepoint,1);
    354447          if( #sol > 0,
     
    365458          r = -centerlift(polcoeff(factpol[i],0));
    366459          if( valuation(subst(pol,variable(pol),r),p) > 2, next);
    367460          M = [p,r;0,1];
    368           U = redquartique(subst(pol,variable(pol),p*variable(pol)+r));
     461          U = redquartique(subst(K*pol,variable(pol),p*variable(pol)+r));
    369462          if( content(U[1]) != p, next);
    370463          sol = ratpoint(K/p^2*U[1],lim,singlepoint,1);
    371464          if( #sol > 0,
     
    459552    return(0));
    460553  pnup = pnu*p;
    461554  nu++;
    462   if( p< LIMBIGPRIME || !LIMBIGPRIME, 
     555  if( p< LIMBIGPRIME || !LIMBIGPRIME,
    463556    for( i = 0, p-1,
    464557      if( zpsoluble(pol,p,nu,pnup,x0+pnu*i),
    465558if( DEBUGLEVEL_ell >= 5, print("fin de zpsoluble"));
     
    572665    if( !test, default(realprecision, prec *= 2))
    573666  );
    574667
     668\\ On n'utilise plus
    575669\\  q = Vec(sum( i = 1, d, norm(x-r[i])));
    576   q = Vec(sum( i = 1, d, norm(x-r[i]) / normderiv[i]^(1/(d-2)))); \\ uses Cremona-Stoll normalization
     670\\ mais la normalisation de Cremona-Stoll
     671  q = Vec(sum( i = 1, d, norm(x-r[i]) / normderiv[i]^(1/(d-2))));
    577672  M = QfbReduce([q[1],q[2]/2;q[2]/2,q[3]]);
    578673  pol = subst(pol,variable(pol),Pol(M[1,])/Pol(M[2,]))*Pol(M[2,])^poldegree(pol);
    579674
     
    601696
    602697  return(delta*subst(Pol(reduc),x,xx)^2);
    603698}
     699if( DEBUGLEVEL_ell >= 4, print("ellsort"));
     700{
     701ellsort(listpts) =
     702\\ tri des points listpts sur une courbe elliptique
     703\\ suivant la hauteur naive.
     704local(n,v,aux,ord);
     705
     706  v = vector(n = #listpts);
     707  for( i = 1, n,
     708    if( listpts[i] == [0], v[i] = [0,0,0]; next);
     709    aux = denominator(listpts[i][2])/denominator(listpts[i][1]);
     710    v[i] = vecsort(abs([listpts[i][1]*aux^2, listpts[i][2]*aux^3,aux]),,4);
     711  );
     712  ord = vecsort(v,,3);
     713  return(vector(n,i,listpts[ord[i]]));
     714}
    604715if( DEBUGLEVEL_ell >= 4, print("ellredgen"));
    605716{
    606717ellredgen(ell,listgen,K=1) =
    607 \\ reduction des generateurs de listgen en utilisant hauteur naive
     718\\ reduction des generateurs de listgen
    608719\\ sur la courbe ell = [a1,a2,a3,a4,a6]
    609720\\ ou K*y^2 = x^3 + a2*x^2 + a4*x + a6 (lorsque a1 = a3 = 0);
    610721local(d,sqrtK,urst,M,U,limgoodrelations,listgen2);
     
    683794    listgen = listgen2
    684795  );
    685796
     797  listgen = ellsort(listgen);
     798if( DEBUGLEVEL_ell >= 4, print("generateurs tries = ",listgen));
     799
    686800  listgen = ellchangepointinverse(listgen,urst);
    687801  if( K != 1,
    688802    for( i = 1, d,
    689       for( j = 1, d,
     803      for( j = 1, 2,
    690804        listgen[i][j] /= K^j)));
    691805
    692806\\ on ne garde que les points (x,y) avec y >= 0
     
    779893if( DEBUGLEVEL_ell >= 3, print("KS2gen = ",KS2gen));
    780894
    781895  LS2genunit = ext.tufu;
    782   LS2genunit = concat(lift(nfbasistoalg(ext,LS2gen[1])),LS2genunit);
     896  LS2genunit = concat(LS2gen[1],LS2genunit);
    783897
    784898  LS2genunit = subst(LS2genunit,x,ttheta);
    785899  LS2genunit = LS2genunit*Mod(1,polrel);
     
    10231137ellrank(ell,help=[]) =
    10241138\\ Algorithme de la 2-descente sur la courbe elliptique ell.
    10251139\\ help est une liste de points connus sur ell.
    1026 local(urst,urst1,den,tors,eqtheta,bnf,rang,time1);
     1140local(urst,urst1,den,eqell,tors2,bnf,rang,time1);
    10271141
    10281142if( DEBUGLEVEL_ell >= 3, print("entree dans ellrank"));
    10291143  if( #ell < 13, ell = ellinit(ell,1));
    10301144
    1031 \\ removes the coefficients a1 and a3
     1145\\ supprime les coefficients a1 et a3
    10321146  urst = [1,0,0,0];
    10331147  if( ell.a1 != 0 || ell.a3 != 0,
    10341148    urst1 = [1,0,-ell.a1/2,-ell.a3/2];
     
    10361150    urst = ellcomposeurst(urst,urst1)
    10371151  );
    10381152
    1039 \\ removes denominators
     1153\\ supprime les denominateurs
    10401154  while( (den = denominator([ell.a2,ell.a4,ell.a6])) > 1,
    10411155    den = factor(den); den[,2] = vectorv(#den[,2],i,1);
    10421156    den = factorback(den);
     
    10461160  );
    10471161
    10481162  help = ellchangepoint(help,urst);
     1163  eqell = Pol([1,ell.a2,ell.a4,ell.a6]);
     1164if( DEBUGLEVEL_ell >= 1, print("courbe elliptique : Y^2 = ",eqell));
    10491165
    10501166\\ choix de l'algorithme suivant la 2-torsion
    1051   eqtheta = Pol([1,ell.a2,ell.a4,ell.a6]);
    1052 if( DEBUGLEVEL_ell >= 1, print("courbe elliptique : Y^2 = ",eqtheta));
    1053   f = polratroots(eqtheta);
    10541167
    1055   if( #f == 0,                                  \\ cas 1: 2-torsion triviale
     1168  tors2 = ellhalf(ell,[0]);
     1169if( DEBUGLEVEL_ell >= 1, print("E[2] = ",tors2)); 
     1170
     1171  if( #tors2 == 1,                              \\ cas 1: 2-torsion triviale
    10561172if( DEBUGLEVEL_ell >= 3, print1("bnfinit "));
    1057 if( DEBUGLEVEL_ell >= 2, gettime());
    1058     bnf = bnfinit(eqtheta,1);
    1059 if( DEBUGLEVEL_ell >= 2, time1 = gettime());
     1173if( DEBUGLEVEL_ell >= 4, gettime());
     1174    bnf = bnfinit(eqell,1);
     1175if( DEBUGLEVEL_ell >= 4, time1 = gettime());
    10601176if( DEBUGLEVEL_ell >= 3, print("ok"));
    10611177    rang = ell2descent_gen(ell,bnf,1,help);
    1062 if( DEBUGLEVEL_ell >= 2, print("temps dans bnfinit = ",time1));
    1063 if( DEBUGLEVEL_ell >= 2, print("temps pour le reste = ",gettime()));
     1178if( DEBUGLEVEL_ell >= 4, print("temps dans bnfinit = ",time1));
     1179if( DEBUGLEVEL_ell >= 4, print("temps pour le reste = ",gettime()));
    10641180  ,
    1065   if( #f >= 1,                                  \\ cas 2: 2-torsion >= Z/2Z
    1066     if( f[1] != 0,
    1067       urst1 = [1,f[1],0,0];
     1181  if( #tors2 >= 2,                              \\ cas 2: 2-torsion >= Z/2Z
     1182    if( ell.a6 != 0,
     1183      urst1 = [1,tors2[2][1],0,0];
    10681184      ell = ellchangecurve(ell,urst1);
    10691185      urst = ellcomposeurst(urst,urst1)
    10701186    );
    1071     rang = ell2descent_viaisog(ell)
     1187    eqell = Pol([1,ell.a2,ell.a4,ell.a6]);
     1188if( DEBUGLEVEL_ell >= 1, print("courbe elliptique : Y^2 = ",eqell));
     1189
     1190    rang = ell2descent_viaisog(ell,help)
    10721191  ,                                             \\ cas 3: 2-torsion = Z/2Z*Z/2Z
    1073 \\    rang = ell2descent_complete(f[1],f[2],f[3])
     1192\\    rang = ell2descent_complete(tors2[2][1],tors2[3][2],tors2[4][3])
    10741193  ));
    10751194
    10761195  rang[3] = ellchangepointinverse(rang[3],urst);
     
    12121331, print("rang        = ",rang));
    12131332  if( rang, print("points = ",listepoints));
    12141333);
    1215   ell = ellinit([0,-(e1+e2+e3),0,e1*e2+e2*e3+e3*e1,e1*e2*e2],1);
    1216   listepoints = ellredgen(ell,listepoints);
    1217   listepoints = concat(elltors(ell,#ell<19)[3],listepoints);
     1334  ell = ellinit([0,-(e1+e2+e3),0,e1*e2+e2*e3+e3*e1,-e1*e2*e3]);
     1335  if( ELLREDGENFLAG, listepoints = ellredgen(ell,listepoints));
     1336  listepoints = concat(ellsort(elltorseven(ell)[3]),listepoints);
    12181337
    12191338return([rang,selmer,listepoints]);
    12201339}
     
    12931412}
    12941413if( DEBUGLEVEL_ell >= 4, print("ell2descent_viaisog"));
    12951414{
    1296 ell2descent_viaisog(ell) =
     1415ell2descent_viaisog(ell,help=[]) =
    12971416\\ Calcul du rang des courbes elliptiques avec 2-torsion
    12981417\\ par la methode des 2-isogenies.
    12991418\\
     
    13181437  P = Pol([1,ell.a2,ell.a4]);
    13191438  Pfact = factor(P)[,1];
    13201439  tors = #Pfact;
    1321   if( #Pfact > 1,
    1322     listpointstriv=[[0,0],[-polcoeff(Pfact[1],0),0],[-polcoeff(Pfact[2],0),0]]
    1323   , listpointstriv = [[0,0]]);
    1324 
     1440  listpointstriv = concat(help,elltorseven(ell)[3]);
     1441 
    13251442  apinit = -2*ell.a2; bpinit = ell.a2^2-4*ell.a4;
    13261443
    13271444  plist = factor(6*ell.disc)[,1];
     
    13661483  print("K(a^2-4b,2)gen     = ",KS2gen));
    13671484
    13681485  P = Pol([1,apinit,bpinit]);
    1369   Pfact = factor(P)[,1];
    1370   if( #Pfact > 1,
    1371     listpointstriv = [[0,0],[-polcoeff(Pfact[1],0),0],[-polcoeff(Pfact[2],0),0]]
    1372   , listpointstriv = [[0,0]]);
    1373 
     1486  listpointstriv = elltorseven([0,apinit,0,bpinit,0])[3];
     1487   
    13741488if( DEBUGLEVEL_ell >= 3, print(" Recherche de points triviaux sur la courbe"));
    13751489  P *= x;
    13761490if( DEBUGLEVEL_ell >= 3, print("Y^2 = ",P));
     
    14081522  if( certain && certainp, print1(" "), print1("<"));
    14091523  print("= ",1<<(n2+np2-n1-np1));
    14101524
    1411   print("#E(K)[2]            = ",1<<tors);
     1525  print("#E(Q)[2]            = ",1<<tors);
    14121526);
    14131527  rang = n1+np1-2;
    14141528if( DEBUGLEVEL_ell >= 1,
     
    14741588
    14751589\\ fin de strange
    14761590 
    1477   pointgen = ellredgen(ell,pointgen);
    1478   pointgen = concat(elltors(ell,#ell<19)[3],pointgen);
     1591  if( ELLREDGENFLAG, pointgen = ellredgen(ell,pointgen));
     1592  pointgen = concat(ellsort(elltorseven(ell)[3]),pointgen);
    14791593if( DEBUGLEVEL_ell >= 1, print("points = ",pointgen));
    14801594if( DEBUGLEVEL_ell >= 3, print("fin de ell2descent_viaisog"));
    14811595  return([rang,n2+np2-2+tors,pointgen]);
  • pari/simon/qfsolve.gp

    diff -r 4efecfcb95fa -r 499528b85ed0 pari/simon/qfsolve.gp
    a b  
    11\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
    2 \\       Copyright (C) 2005 Denis Simon
     2\\       Copyright (C) 2007 Denis Simon
    33\\
    44\\ Distributed under the terms of the GNU General Public License (GPL)
    55\\
     
    2020\\ www.math.unicaen.fr/~simon/qfsolve.gp
    2121\\
    2222\\  *********************************************         
    23 \\  *          VERSION 25/10/2005               *
     23\\  *          VERSION 02/10/2009               *
    2424\\  *********************************************         
    2525\\
    2626\\ Programme de resolution des equations quadratiques
     
    5454\\
    5555
    5656\\
    57 global(DEBUGLEVEL);
    58 DEBUGLEVEL = 0;
     57DEBUGLEVEL_qfsolve = 0;
    5958
    6059\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
    6160\\
    6261\\  DEBUT DES SOURCES                \\
    6362\\
    6463\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
    65 {QfbReduce(M)=
    66 \\ M=[a,b;b;c] est a coefficients entiers.
     64{QfbReduce(M) =
     65\\ M = [a,b;b;c] est a coefficients entiers.
    6766\\ Reduction de la forme quadratique binaire
    68 \\   qf=(a,b,c)=a*X^2+2*b*X*Y+c*Y^2
     67\\   qf = (a,b,c)=a*X^2+2*b*X*Y+c*Y^2
    6968\\ Renvoit la matrice de reduction de det = +1.
    7069
    7170local(a,b,c,H,test,di,q,r,nexta,nextb,nextc,aux);
    7271
    73 if( DEBUGLEVEL >= 5, print("entree dans QfbReduce avec ",M));
     72if( DEBUGLEVEL_qfsolve >= 5, print("entree dans QfbReduce avec ",M));
    7473
    7574  a = M[1,1]; b = M[1,2]; c = M[2,2];
    7675 
    7776  H = matid(2); test = 1;
    78   while(test && a,
     77  while( test && a,
    7978    di = divrem(b,a); q = di[1]; r = di[2];
    8079    if( 2*r > abs(a), r -= abs(a); q += sign(a));
    8180    H[,2] -= q*H[,1];
    8281    nextc = a; nextb = -r; nexta= (nextb-b)*q+c;
    8382 
    84     if( test = abs(nexta)<abs(a),
     83    if( test = abs(nexta) < abs(a),
    8584      c = nextc; b = nextb; a = nexta;
    86       aux = H[,1]; H[,1] = -H[,2]; H[,2] = aux;
     85      aux = H[,1]; H[,1] = -H[,2]; H[,2] = aux
    8786    )
    8887  );
    8988
    90 if( DEBUGLEVEL >= 5, print("sortie de QfbReduce avec ",H));
    91 return(H)
     89if( DEBUGLEVEL_qfsolve >= 5, print("sortie de QfbReduce avec ",H));
     90return(H);
    9291}
    93 {IndefiniteLLL(G,c=1,base=0)=
     92{IndefiniteLLL(G,c=1,base=0) =
     93\\ Performs first a LLL reduction on a positive definite
     94\\ quadratic form QD bounding the indefinite G.
     95\\ Then finishes the reduction with IndefiniteLLL2.
     96
     97local(n,M,QD,M1,S,red);
     98
     99  n = length(G);
     100  M = matid(n);
     101  QD = G;
     102  for( i = 1, n-1,
     103    if( !QD[i,i],
     104return(IndefiniteLLL2(G,c,base))
     105    );
     106    M1 = matid(n);
     107    M1[i,] = -QD[i,]/QD[i,i];
     108    M1[i,i] = 1;
     109    M = M*M1;
     110    QD = M1~*QD*M1
     111  );
     112  M = M^(-1);
     113  QD = M~*abs(QD)*M;
     114  S = qflllgram(QD/content(QD));
     115  red = IndefiniteLLL2(S~*G*S,c,base);
     116  if( type(red) == "t_COL",
     117return(S*red));
     118  if( length(red) == 3,
     119return([red[1],S*red[2],S*red[3]]));
     120return([red[1],S*red[2]]);
     121}
     122{IndefiniteLLL2(G,c=1,base=0) =
    94123\\ following Cohen's book p. 86
    95124\\ but without b and bstar: works on G
    96 \\ returns [H~*G*H,H] where det(H)=1 and H~*G*H is reduced.
     125\\ returns [H~*G*H,H] where det(H) = 1 and H~*G*H is reduced.
    97126\\ Exit with a norm 0 vector if one such is found.
    98 \\ If base==1 and norm 0 is obtained, returns [H~*G*H,H,sol] where
     127\\ If base == 1 and norm 0 is obtained, returns [H~*G*H,H,sol] where
    99128\\   sol is a norm 0 vector and is the 1st column of H.
    100129
    101 local(n,H,M,A,aux,sol,k,nextk,swap,q,di,HM,au1,aux2,Mkk1,bk1new,Mkk1new,newG);
     130local(n,H,M,A,aux,sol,k,nextk,swap,q,di,HM,aux1,aux2,Mkk1,bk1new,Mkk1new,newG);
    102131
    103132  n = length(G);
    104 if( DEBUGLEVEL >= 3, print("LLL dim ",n," avec G=",log(vecmax(abs(G)))/log(10)));
    105 if( DEBUGLEVEL >= 4, print("LLL avec ");printp(G));
     133if( DEBUGLEVEL_qfsolve >= 3, print("LLL dim ",n," avec G=",log(vecmax(abs(G)))/log(10)));
     134if( DEBUGLEVEL_qfsolve >= 4, print("LLL avec ");printp(G));
    106135
    107136 if( n <= 1, return([G,matid(n)]));
    108137
     
    111140\\ compute Gram-Schmidt 
    112141
    113142  for( i = 1, n,
    114     if( !(A[i,i] = G[i,i]) ,
     143    if( !(A[i,i] = G[i,i]),
    115144      if( base,
    116145        aux = H[,1]; H[,1] = H[,i]; H[,i] = -aux;
    117146        return([H~*G*H,H,H[,1]])
     
    137166
    138167    swap = 1;
    139168    while( swap,
    140       swap=0;
     169      swap = 0;
    141170
    142171\\ red(k,k-1);
    143172      if( q = round(M[k,k-1]),
     
    161190        sol /= content(sol);
    162191        if( base,
    163192          H = H*completebasis(sol,1);
    164           aux = H[,1]; H[,1] = H[,n]; H[,n] = -aux;     
     193          aux = H[,1]; H[,1] = H[,n]; H[,n] = -aux;
    165194          return([H~*G*H,H,H[,1]])
    166195        , return(H*sol)
    167196        )
     
    195224
    196225        M[k,k-1] = Mkk1new;
    197226
    198 if( DEBUGLEVEL >=4, newG=H~*G*H;print(vector(n,i,matdet(vecextract(newG,1<<i-1,1<<i-1)))));
     227if( DEBUGLEVEL_qfsolve >=4, newG=H~*G*H;print(vector(n,i,matdet(vecextract(newG,1<<i-1,1<<i-1)))));
    199228
    200         if( k !=2, k--)
     229        if( k != 2, k--)
    201230      )
    202231    );
    203232
     
    214243    );
    215244    k++
    216245  );
    217 return([H~*G*H,H])
     246return([H~*G*H,H]);
    218247}
    219 {kermodp(M,p)=
     248{kermodp(M,p) =
    220249\\ Compute the kernel of M mod p.
    221250\\ returns [d,U], where
    222251\\ d = dim (ker M mod p)
     
    229258  d = length(U);
    230259  U = completebasis(U);
    231260  U = matrix(n,n,i,j,U[i,n+1-j]);
    232 return([d,U])
     261return([d,U]);
    233262}
    234 {Qfparam(G,sol,fl=3)=
     263{Qfparam(G,sol,fl=3) =
    235264\\ G est une matrice symetrique 3*3, et sol une solution de sol~*G*sol=0.
    236265\\ Renvoit une parametrisation des solutions avec de bons invariants,
    237266\\ sous la forme d'une matrice 3*3, dont chaque ligne contient
     
    241270
    242271local(U,G1,G2);
    243272
    244 if( DEBUGLEVEL >= 5, print("entree dans Qfparam"));
    245   sol /= content(sol); 
    246 \\ construction de U telle que U[,3] est sol, et det(U)=+-1
     273if( DEBUGLEVEL_qfsolve >= 5, print("entree dans Qfparam"));
     274  sol /= content(sol);
     275\\ construction de U telle que U[,3] est sol, et det(U) = +-1
    247276  U = completebasis(sol,1);
    248277  G1 = U~*G*U; \\ G1 a un 0 en bas a droite.
    249278  G2 = [-2*G1[1,3],-2*G1[2,3],0;
     
    258287         U[2,1]^2,2*U[2,1]*U[2,2],U[2,2]^2];
    259288    sol = sol*U
    260289  );
    261 if( DEBUGLEVEL >= 5, print("sortie de Qfparam"));
    262 return(sol)
     290if( DEBUGLEVEL_qfsolve >= 5, print("sortie de Qfparam"));
     291return(sol);
    263292}
    264 {LLLgoon3(G,c=1)=
     293{LLLgoon3(G,c=1) =
    265294\\ reduction LLL de la forme quadratique G (matrice de Gram)
    266295\\ en dim 3 seulement avec detG = -1 et sign(G) = [1,2];
    267296
     
    276305  U2 = [bez[1],G2[3,2]/bez[3],0;bez[2],-G2[3,1]/bez[3],0;0,0,-1];
    277306  G3 = U2~*G2*U2;
    278307\\ la matrice G3 a des 0 sous l'anti-diagonale.
    279   cc = G3[1,1]%2; 
    280   U3 = [1,0,0;  cc,1,0; 
     308  cc = G3[1,1]%2;
     309  U3 = [1,0,0;  cc,1,0;
    281310        round(-(G3[1,1]+cc*(2*G3[1,2]+G3[2,2]*cc))/2/G3[1,3]),
    282311        round(-(G3[1,2]+cc*G3[2,2])/G3[1,3]),1];
    283 return([U3~*G3*U3,red[2]*U1*U2*U3])
     312return([U3~*G3*U3,red[2]*U1*U2*U3]);
    284313}
    285 {completebasis(v,redflag=0)=
     314{completebasis(v,redflag=0) =
    286315\\ Donne une matrice unimodulaire dont la derniere colonne est v.
    287316\\ Si redflag <> 0, alors en plus on reduit le reste.
    288317
     
    292321  n = length(v~);
    293322  if( n==1 || !redflag, return(U));
    294323  re = qflll(vecextract(U,1<<n-1,1<<(n-1)-1));
    295 return( U*matdiagonalblock([re,Mat(1)]))
     324return( U*matdiagonalblock([re,Mat(1)]));
    296325}
    297 {LLLgoon(G,c=1)=
     326{LLLgoon(G,c=1) =
    298327\\ reduction LLL de la forme quadratique G (matrice de Gram)
    299328\\ ou l'on continue meme si on trouve un vecteur isotrope
    300329
     
    305334  if( length(red) == 2, return(red));
    306335\\ sinon :
    307336  U1 = red[2];
    308   G2 = red[1]; \\ On a G2[1,1]=0
     337  G2 = red[1]; \\ On a G2[1,1] = 0
    309338  U2 = mathnf(Mat(G2[1,]),4)[2];
    310339  G3 = U2~*G2*U2;
    311340\\ la matrice de G3 a des 0 sur toute la 1ere ligne,
     
    331360  red = LLLgoon(matrix(n-2,n-2,i,j,G5[i+1,j+1]),c);
    332361  U5 = matdiagonalblock([Mat(1),red[2],Mat(1)]);
    333362  G6 = U5~*G5*U5;
    334 return([G6,U1*U2*U3*U4*U5])
     363return([G6,U1*U2*U3*U4*U5]);
    335364}
    336 {QfWittinvariant(G,p)=
     365{QfWittinvariant(G,p) =
    337366\\ calcule l'invariant c (=invariant de Witt) d'une forme quadratique,
    338367\\ p-adique (reel si p = -1)
    339368
     
    348377  c = prod( i = 1, n,
    349378        prod( j = i+1, n,
    350379          hilbert( diag[i], diag[j], p)));
    351 return(c)
     380return(c);
    352381}
    353 {Qflisteinvariants(G,fa=[])=
     382{Qflisteinvariants(G,fa=[]) =
    354383\\ G est une forme quadratique, ou une matrice symetrique,
    355384\\ ou un vecteur contenant des formes quadratiques de meme discriminant.
    356385\\ fa = factor(-abs(2*matdet(G)))[,1] est facultatif.
    357386
    358387local(l,sol,n,det);
    359388
    360 if( DEBUGLEVEL >= 4, print("entree dans Qflisteinvariants",G));
     389if( DEBUGLEVEL_qfsolve >= 4, print("entree dans Qflisteinvariants",G));
    361390  if( type(G) != "t_VEC", G = [G]);
    362391  l = length(G);
    363392  for( j = 1, l,   
     
    371400\\ En dimension 2, chaque invariant est un unique symbole de Hilbert.
    372401    det = -matdet(G[1]);
    373402    sol = matrix(length(fa),l,i,j,hilbert(G[j][1,1],det,fa[i])<0);
    374 if( DEBUGLEVEL >= 4, print("sortie de Qflisteinvariants"));
     403if( DEBUGLEVEL_qfsolve >= 4, print("sortie de Qflisteinvariants"));
    375404    return([fa,sol])
    376405  );
    377406
     
    380409    n = length(G[j]);
    381410\\ En dimension n, on calcule un produit de n symboles de Hilbert.
    382411    det = vector(n+1, i, matdet(matrix(i-1,i-1,k,m,G[j][k,m])));
    383     for(i = 1, length(fa),
     412    for( i = 1, length(fa),
    384413      sol[i,j] = prod( k = 1, n-1, hilbert(-det[k],det[k+1],fa[i]))*hilbert(det[n],det[n+1],fa[i]) < 0;
    385414    )
    386415  );
    387 if( DEBUGLEVEL >= 4, print("sortie de Qflisteinvariants"));
    388 return([fa,sol])
     416if( DEBUGLEVEL_qfsolve >= 4, print("sortie de Qflisteinvariants"));
     417return([fa,sol]);
    389418}
    390 {Qfsolvemodp(G,p)=
     419{Qfsolvemodp(G,p) =
    391420\\ p a prime number.
    392421\\ finds a solution mod p for the quadatic form G
    393422\\ such that det(G) !=0 mod p and dim G = n>=3;
     
    413442    if( issquare( N2 = -vdet[3]/vdet[1]),         s = sqrt(N2); sol = s*x2+x3; break);
    414443    if( issquare( N3 = -vdet[2]*vdet[3]/vdet[1]), s = sqrt(N3); sol = s*x1+x3; break);
    415444    r = 1;
    416     while ( !issquare( s = (1-N1*r^2)/N3), r = random(p));
     445    while( !issquare( s = (1-N1*r^2)/N3), r = random(p));
    417446    s = sqrt(s); sol = x1+r*x2+s*x3; break
    418447  );
    419448  sol = vectorv(n, j, if( j <= 3, sol[j]));
    420 return(sol)
     449return(sol);
    421450}
    422 {Qfminim(G,factdetG=0)=
     451{Qfminim(G,factdetG=0) =
    423452\\ Minimisation de la forme quadratique entiere non degeneree G,
    424453\\   de dimension n >=2. On suppose que G est symetrique a coefficients entiers.
    425454\\ Renvoit [G',U,factd] avec U \in GLn(Q) telle que G'=U~*G*U*constante est entiere
     
    439468  i = 1;
    440469  while(i <= length(factdetG[,1]),
    441470    p = factdetG[i,1];
    442     if (p == -1, i++; next);
    443 if( DEBUGLEVEL >= 4, print("p=",p,"^",factdetG[i,2]));
     471    if( p == -1, i++; next);
     472if( DEBUGLEVEL_qfsolve >= 4, print("p=",p,"^",factdetG[i,2]));
    444473    vp = factdetG[i,2];
    445474    if( vp == 0, i++; next);
    446475\\ Le cas vp=1 n'est minimisable que si n est impair.
     
    450479    );
    451480    Ker = kermodp(G,p); dimKer = Ker[1]; Ker = Ker[2];
    452481\\ Rem: on a toujours dimKer <= vp
    453 if( DEBUGLEVEL >= 4, print("dimKer = ",dimKer));
     482if( DEBUGLEVEL_qfsolve >= 4, print("dimKer = ",dimKer));
    454483\\ cas trivial: G est divisible par p.
    455484    if( dimKer == n,
    456485      G /= p;
     
    462491\\ 1er cas: la dimension du noyau est plus petite que la valuation
    463492\\ alors le noyau mod p contient un noyau mod p^2
    464493    if( dimKer < vp,
    465 if( DEBUGLEVEL >= 4, print("cas 1"));
     494if( DEBUGLEVEL_qfsolve >= 4, print("cas 1"));
    466495      Ker2 = kermodp(matrix(dimKer,dimKer,j,k,G[j,k]/p),p);
    467496      dimKer2 = Ker2[1]; Ker2 = Ker2[2];
    468497      for( j = 1, dimKer2, Ker2[,j] /= p);
     
    470499      G = Ker2~*G*Ker2;
    471500      U = U*Ker2;
    472501      factdetG[i,2] -= 2*dimKer2;
    473 if( DEBUGLEVEL >= 4, print("fin cas 1"));
    474       next;
     502if( DEBUGLEVEL_qfsolve >= 4, print("fin cas 1"));
     503      next
    475504    );
    476505
    477506\\ Maintenant, vp = dimKer
     
    480509       (dimKer == 2 && issquare(di=Mod((G[1,2]^2-G[1,1]*G[2,2])/p^2,p))),
    481510\\ on cherche dans le noyau un elt de norme p^2...
    482511      if( dimKer > 2,
    483 if( DEBUGLEVEL >= 4, print("cas 2.1"));
     512if( DEBUGLEVEL_qfsolve >= 4, print("cas 2.1"));
    484513        dimKer = 3;
    485514        sol = Qfsolvemodp(matrix(3,3,j,k,G[j,k]/p),p)
    486515      , 
    487 if( DEBUGLEVEL >= 4, print("cas 2.2"));
     516if( DEBUGLEVEL_qfsolve >= 4, print("cas 2.2"));
    488517        if( G[1,1]%p^2 == 0,
    489           sol = [1,0]~;
    490         , sol = [-G[1,2]/p+sqrt(di),Mod(G[1,1]/p,p)]~;
     518          sol = [1,0]~
     519        , sol = [-G[1,2]/p+sqrt(di),Mod(G[1,1]/p,p)]~
    491520        )
    492521      );
    493522      sol = centerlift(sol); sol /= content(sol);
    494 if( DEBUGLEVEL >= 4, print("sol=",sol));
     523if( DEBUGLEVEL_qfsolve >= 4, print("sol=",sol));
    495524      Ker = vectorv(n, j, if( j<= dimKer, sol[j], 0)); \\ on complete avec des 0
    496525      Ker = completebasis(Ker,1);
    497526      Ker[,n] /= p;
    498527      G = Ker~*G*Ker;
    499528      U = U*Ker;
    500529      factdetG[i,2] -= 2;
    501 if( DEBUGLEVEL >= 4, print("fin cas 2"));
    502       next;
     530if( DEBUGLEVEL_qfsolve >= 4, print("fin cas 2"));
     531      next
    503532    );
    504533
    505534\\ Maintenant, vp = dimKer <= 2
     
    508537\\ Dans certains cas, on peut echanger le noyau et l'image
    509538\\   pour continuer a reduire.
    510539
    511     m = (n-1)\2-1; 
     540    m = (n-1)\2-1;
    512541    if( ( vp == 1 && issquare(Mod(-(-1)^m*matdet(G)/G[1,1],p)))
    513542     || ( vp == 2 && n%2 == 1 && n >= 5)
    514543     || ( vp == 2 && n%2 == 0 && !issquare(Mod((-1)^m*matdet(G)/p^2,p)))
    515544    ,
    516 if( DEBUGLEVEL >= 4, print("cas 3"));
     545if( DEBUGLEVEL_qfsolve >= 4, print("cas 3"));
    517546      Ker = matid(n);
    518547      for( j = dimKer+1, n, Ker[j,j] = p);
    519548      G = Ker~*G*Ker/p;
    520549      U = U*Ker;
    521550      factdetG[i,2] -= 2*dimKer-n;
    522 if( DEBUGLEVEL >= 4, print("fin cas 3"));
    523       next;
     551if( DEBUGLEVEL_qfsolve >= 4, print("fin cas 3"));
     552      next
    524553    );
    525554
    526555\\ On n'a pas pu minimiser.
    527556\\ Si n == 3 ou n == 4 cela demontre la non-solubilite locale en p.
    528557    if( n == 3 || n == 4,
    529 if( DEBUGLEVEL >= 1, print("pas de solution locale en ",p));
     558if( DEBUGLEVEL_qfsolve >= 1, print("pas de solution locale en ",p));
    530559      return(p));
    531560
    532 if( DEBUGLEVEL >= 4, print("plus de minimisation possible"));
     561if( DEBUGLEVEL_qfsolve >= 4, print("plus de minimisation possible"));
    533562    factd = concat(factd~,Mat([p,vp])~)~;
    534     i++;
     563    i++
    535564  );
    536565\\print("Un=",log(vecmax(abs(U))));
    537566  aux = qflll(U);
    538567\\print("Ur=",log(vecmax(abs(U*aux))));
    539 return([aux~*G*aux,U*aux,factd])
     568return([aux~*G*aux,U*aux,factd]);
    540569}
    541 {mymat(qfb)=qfb=Vec(qfb);[qfb[1],qfb[2]/2;qfb[2]/2,qfb[3]]
     570{mymat(qfb) = qfb = Vec(qfb);[qfb[1],qfb[2]/2;qfb[2]/2,qfb[3]];
    542571}
    543 {Qfbsqrtgauss(G,factdetG)=
     572{Qfbsqrtgauss(G,factdetG) =
    544573\\ pour l'instant, ca ne marche que pour detG sans facteur carre
    545574\\ sauf en 2, ou la valuation est 2 ou 3.
    546575\\ factdetG contient la factorisation de 2*abs(disc G)
    547576local(a,b,c,d,m,n,p,aux,Q1,M);
    548 if( DEBUGLEVEL >=3, print("entree dans Qfbsqrtgauss avec",G,factdetG));
     577if( DEBUGLEVEL_qfsolve >=3, print("entree dans Qfbsqrtgauss avec",G,factdetG));
    549578  G = Vec(G);
    550579  a = G[1]; b = G[2]/2; c = G[3]; d = a*c-b^2;
    551580
    552581\\ 1ere etape: on resout m^2 = a (d), m*n = -b (d), n^2 = c (d)
    553582  m = n = Mod(1,1);
    554 factdetG[1,2]-=3;
     583  factdetG[1,2] -= 3;
    555584  for( i = 1, length(factdetG[,1]),
    556     if(!factdetG[i,2],next);
     585    if( !factdetG[i,2], next);
    557586    p = factdetG[i,1];
    558     if(gcd(a,p)==1,
     587    if( gcd(a,p) == 1,
    559588      aux = sqrt(Mod(a,p));
    560589      m = chinese(m,aux);
    561590      n = chinese(n,-b/aux)
     
    566595    )
    567596  );
    568597  m = centerlift(m);  n = centerlift(n);
    569 if( DEBUGLEVEL >=4, print("m=",m);print("n=",n)); 
     598if( DEBUGLEVEL_qfsolve >=4, print("m=",m); print("n=",n));
    570599 
    571 \\ 2eme etape: on construit Q1 de det -1 tq Q1(x,y,0)=G(x,y)
     600\\ 2eme etape: on construit Q1 de det -1 tq Q1(x,y,0) = G(x,y)
    572601  Q1 = [(n^2-c)/d, (m*n+b)/d, n ;
    573602        (m*n+b)/d, (m^2-a)/d, m ;
    574603        n,         m,         d ];
     
    576605
    577606\\ 3eme etape: on reduit Q1 jusqu'a [0,0,-1;0,1,0;-1,0,0]
    578607  M = LLLgoon3(Q1)[2][3,];
    579   if(M[1]<0,M=-M);
    580 if( DEBUGLEVEL >=3, print("fin de Qfbsqrtgauss "));
     608  if( M[1] < 0, M = -M);
     609if( DEBUGLEVEL_qfsolve >=3, print("fin de Qfbsqrtgauss "));
    581610  if( M[1]%2,
    582611    return(Qfb(M[1],2*M[2],2*M[3]))
    583   , return(Qfb(M[3],-2*M[2],2*M[1])))
     612  , return(Qfb(M[3],-2*M[2],2*M[1])));
    584613}
    585 {class2(D,factdetG,Winvariants,U2)=
     614{class2(D,factdetG,Winvariants,U2) =
    586615\\ On suit l'algorithme de Shanks/Bosma-Stevenhagen
    587616\\ pour construire tout le 2-groupe de classes
    588617\\ seulement pour D discriminant fondamental.
    589 \\ lorsque D=1(4), on travaille avec 4D.
     618\\ lorsque D = 1(4), on travaille avec 4D.
    590619\\ Si on connait la factorisation de abs(2*D),
    591620\\ on peut la donner dans factdetG, et dans ce cas le reste
    592621\\ de l'algorithme est polynomial.
     
    595624
    596625local(factD,n,rang,m,listgen,vD,p,vp,aux,invgen,im,Ker,Kerim,listgen2,G2,struct,E,compteur,red);
    597626
    598 if( DEBUGLEVEL >= 1, print("Construction du 2-groupe de classe de discriminant ",D));
     627if( DEBUGLEVEL_qfsolve >= 1, print("Construction du 2-groupe de classe de discriminant ",D));
    599628  if( D%4 == 2 || D%4 == 3, print("Discriminant not congruent to 0,1 mod 4");return(0));
    600629
    601   if(D==-4, return([[1],[Qfb(1,0,1)]]));
     630  if( D==-4, return([[1],[Qfb(1,0,1)]]));
    602631
    603   if(!factdetG, factdetG = factor(2*abs(D)));
     632  if( !factdetG, factdetG = factor(2*abs(D)));
    604633  factD = concat([-1],factdetG[,1]~);
    605   if( D%4 == 1, D*=4; factdetG[1,2]+=2);
     634  if( D%4 == 1, D *= 4; factdetG[1,2] += 2);
    606635 
    607   n = length(factD);  rang = n-3;
    608   if(D>0, m = rang+1, m = rang); 
    609 if( DEBUGLEVEL >= 3, print("factD = ",factD));
    610   listgen = vector(m);
     636  n = length(factD); rang = n-3;
     637  if(D>0, m = rang+1, m = rang);
     638if( DEBUGLEVEL_qfsolve >= 3, print("factD = ",factD));
     639  listgen = vector(max(0,m));
    611640 
    612641  if( vD = valuation(D,2),
    613     E = Qfb(1,0,-D/4);
     642    E = Qfb(1,0,-D/4)
    614643  , E = Qfb(1,1,(1-D)/4)
    615644  );
    616 if( DEBUGLEVEL >= 3, print("E = ",E));
     645if( DEBUGLEVEL_qfsolve >= 3, print("E = ",E));
    617646
    618   if( type(Winvariants)=="t_COL" && (Winvariants == 0 || length(matinverseimage(U2*Mod(1,2),Winvariants))>0), return([[1],[E]]));
     647  if( type(Winvariants) == "t_COL" && (Winvariants == 0 || length(matinverseimage(U2*Mod(1,2),Winvariants))>0), return([[1],[E]]));
    619648 
    620649  for( i = 1, m, \\ on  ne regarde pas factD[1]=-1, ni factD[2]=2
    621650    p = factD[i+2];
    622651    vp = valuation(D,p);
    623652    aux = p^vp;
    624     if (vD,
     653    if( vD,
    625654      listgen[i] = Qfb(aux,0,-D/4/aux)
    626     , listgen[i] = Qfb(aux,aux,(aux-D/aux)/4));
     655    , listgen[i] = Qfb(aux,aux,(aux-D/aux)/4))
    627656  );
    628657  if( vD == 2 && D%16 != 4,
    629     m++;rang++; listgen = concat(listgen,[Qfb(2,2,(4-D)/8)]));
     658    m++; rang++; listgen = concat(listgen,[Qfb(2,2,(4-D)/8)]));
    630659  if( vD == 3,
    631     m++;rang++; listgen = concat(listgen,[Qfb(2^(vD-2),0,-D/2^vD)]));
     660    m++; rang++; listgen = concat(listgen,[Qfb(2^(vD-2),0,-D/2^vD)]));
    632661 
    633 if( DEBUGLEVEL >= 3, print("listgen = ",listgen));
    634 if( DEBUGLEVEL >= 2, print("rang = ",rang));
     662if( DEBUGLEVEL_qfsolve >= 3, print("listgen = ",listgen));
     663if( DEBUGLEVEL_qfsolve >= 2, print("rang = ",rang));
    635664
    636665  if( !rang, return([[1],[E]]));
    637666
    638667 invgen = Qflisteinvariants(listgen,factD)[2]*Mod(1,2);
    639 if( DEBUGLEVEL >= 3, printp("invgen = ",lift(invgen)));
     668if( DEBUGLEVEL_qfsolve >= 3, printp("invgen = ",lift(invgen)));
    640669
    641670  struct = vector(m,i,2);
    642671  im = lift(matinverseimage(invgen,matimage(invgen)));
    643672  while( (length(im) < rang)
    644   || (type(Winvariants) == "t_COL" && length(matinverseimage(concat(invgen,U2),Winvariants)==0)),
     673  || (type(Winvariants) == "t_COL" && length(matinverseimage(concat(invgen,U2),Winvariants) == 0)),
    645674    Ker = lift(matker(invgen));
    646675    Kerim = concat(Ker,im);
    647676    listgen2 = vector(m);
     
    658687    listgen = listgen2;
    659688    invgen = invgen*Kerim;
    660689
    661 if( DEBUGLEVEL >= 4, print("listgen = ",listgen));
    662 if( DEBUGLEVEL >= 4, printp("invgen = ",lift(invgen)));
     690if( DEBUGLEVEL_qfsolve >= 4, print("listgen = ",listgen));
     691if( DEBUGLEVEL_qfsolve >= 4, printp("invgen = ",lift(invgen)));
    663692
    664693    for( i = 1, length(Ker),
    665694      G2 = Qfbsqrtgauss(listgen[i],factdetG);
    666       struct[i]<<=1;
     695      struct[i] <<= 1;
    667696      listgen[i] = G2;
    668697      invgen[,i] = Qflisteinvariants(G2,factD)[2][,1]*Mod(1,2)
    669698    );
    670699
    671 if( DEBUGLEVEL >= 3, print("listgen = ",listgen));
    672 if( DEBUGLEVEL >= 3, printp("invgen = ",lift(invgen)));
    673 if( DEBUGLEVEL >= 3, print("struct = ",struct));
     700if( DEBUGLEVEL_qfsolve >= 3, print("listgen = ",listgen));
     701if( DEBUGLEVEL_qfsolve >= 3, printp("invgen = ",lift(invgen)));
     702if( DEBUGLEVEL_qfsolve >= 3, print("struct = ",struct));
    674703
    675704    im = lift(matinverseimage(invgen,matimage(invgen)))
    676   ); 
     705  );
    677706 
    678707  listgen2 = vector(rang);
    679708  for( i = 1, rang,
     
    690719\\ listgen = vector(rang,i,listgen[m-rang+i]);
    691720  struct = vector(rang,i,struct[m-rang+i]);
    692721
    693 if( DEBUGLEVEL >= 2, print("listgen = ",listgen));
    694 if( DEBUGLEVEL >= 2, print("struct = ",struct));
     722if( DEBUGLEVEL_qfsolve >= 2, print("listgen = ",listgen));
     723if( DEBUGLEVEL_qfsolve >= 2, print("struct = ",struct));
    695724
    696 return([struct,listgen])
     725return([struct,listgen]);
    697726}
    698 {Qfsolve(G,factD)=
     727{Qfsolve(G,factD) =
    699728\\ Resolution de la forme quadratique X^tGX=0 de dimension n >= 3.
    700729\\ On suppose que G est entiere et primitive.
    701730\\ La solution peut etre un vectorv ou une matrice.
     
    705734\\ Si on connait la factorisation de -abs(2*matdet(G)),
    706735\\ on peut la passer par le parametre factD pour gagner du temps.
    707736
    708 local(n,M,signG,d,Min,U,codim,aux,G1,detG1,M1,subspace1,G2,subspace2,M2,solG2,Winvariants,dQ,factd,U2,clgp2,V,detG2,dimseti,solG1,sol);
     737local(n,M,signG,d,Min,U,codim,aux,G1,detG1,M1,subspace1,G2,subspace2,M2,solG2,Winvariants,dQ,factd,U2,clgp2,V,detG2,dimseti,solG1,sol,Q);
    709738
    710 if( DEBUGLEVEL >= 1, print("entree dans Qfsolve"));
     739if( DEBUGLEVEL_qfsolve >= 1, print("entree dans Qfsolve"));
    711740\\
    712741\\ 1ere reduction des coefficients de G
    713742\\
    714743
    715744  n = length(G);
    716745  M = IndefiniteLLL(G);
    717   if(type(M) == "t_COL",
    718 if( DEBUGLEVEL >= 1, print("solution triviale",M));
     746  if( type(M) == "t_COL",
     747if( DEBUGLEVEL_qfsolve >= 1, print("solution ",M));
    719748    return(M));
    720749  G = M[1]; M = M[2];
    721750
    722751\\ Solubilite reelle
    723752  signG = qfsign(G);
    724753  if( signG[1] == 0 || signG[2] == 0,
    725 if( DEBUGLEVEL >= 1, print("pas de solution reelle"));
     754if( DEBUGLEVEL_qfsolve >= 1, print("pas de solution reelle"));
    726755    return(-1));
    727756  if( signG[1] < signG[2], G = -G; signG = signG*[0,1;1,0]);
    728757
    729758\\ Factorisation du determinant
    730759  d = matdet(G);
    731760  if( !factD,
    732 if( DEBUGLEVEL >= 1, print("factorisation du determinant"));
     761if( DEBUGLEVEL_qfsolve >= 1, print("factorisation du determinant"));
    733762    factD = factor(-abs(2*d));
    734 if( DEBUGLEVEL >= 1, print(factD));
     763if( DEBUGLEVEL_qfsolve >= 1, print(factD))
    735764  );
    736765  factD[1,2] = 0;
    737766  factD[2,2] --;
     
    740769\\ Minimisation et solubilite locale.
    741770\\
    742771
    743 if( DEBUGLEVEL >= 1, print("minimisation du determinant"));
     772if( DEBUGLEVEL_qfsolve >= 1, print("minimisation du determinant"));
    744773  Min = Qfminim(G,factD);
    745   if (type(Min) == "t_INT",
    746 if( DEBUGLEVEL >= 1, print("pas de solution locale en ",Min));
     774  if( type(Min) == "t_INT",
     775if( DEBUGLEVEL_qfsolve >= 1, print("pas de solution locale en ",Min));
    747776    return(Min));
    748777
    749778  M = M*Min[2];
    750779  G = Min[1];
    751780\\  Min[3] contient la factorisation de abs(matdet(G));
    752781
    753 if( DEBUGLEVEL >= 4, print("G minime = ",G));
    754 if( DEBUGLEVEL >= 4, print("d=",d));
     782if( DEBUGLEVEL_qfsolve >= 4, print("G minime = ",G));
     783if( DEBUGLEVEL_qfsolve >= 4, print("d=",d));
    755784
    756785\\ Maintenant, on sait qu'il y a des solutions locales
    757786\\ (sauf peut-etre en 2 si n==4),
     
    762791\\ Reduction de G et recherche de solution triviales
    763792\\ par exemple quand det G=+-1, il y en a toujours.
    764793
    765 if( DEBUGLEVEL >= 1, print("reduction"));
     794if( DEBUGLEVEL_qfsolve >= 1, print("reduction"));
    766795  U = IndefiniteLLL(G);
    767796  if(type(U) == "t_COL",
    768 if( DEBUGLEVEL >= 1, print("solution ",M*U));
     797if( DEBUGLEVEL_qfsolve >= 1, print("solution ",M*U));
    769798    return(M*U));
    770799  G = U[1]; M = M*U[2];
    771800
     
    775804\\
    776805
    777806  if( n >= 6 && n%2 == 0 && matsize(Min[3])[1] != 0,
    778 if( DEBUGLEVEL >= 1, print("On passe en dimension ",n+1));
     807if( DEBUGLEVEL_qfsolve >= 1, print("On passe en dimension ",n+1));
    779808    codim = 1; n++;
    780809\\ On calcule le plus grand diviseur carre de d.   
    781810    aux = prod( i = 1, matsize(Min[3])[1], if( Min[3][i,1] == 2, Min[3][i,1], 1));
     
    793822    Min = Qfminim(G1,factD);
    794823    G1 = Min[1];
    795824    M1 = Min[2];
    796     subspace1 = matrix(n,n-1,i,j,i==j)
     825    subspace1 = matrix(n,n-1,i,j, i == j)
    797826  , codim = 0;
    798827    G1 = G;
    799828    subspace1 = M1 = matid(n)
     
    805834\\
    806835
    807836  if( matsize(Min[3])[1] == 0, \\ if( abs(d) == 1,
    808 if( DEBUGLEVEL >= 2, print(" detG2 = 1"));
     837if( DEBUGLEVEL_qfsolve >= 2, print(" detG2 = 1"));
    809838     G2 = G1;
    810839     subspace2 = M2 = matid(n);
    811840     solG2 = LLLgoon(G2,1)
    812841  ,
    813 if( DEBUGLEVEL >= 1, print("On passe en dimension ",n+2));
     842if( DEBUGLEVEL_qfsolve >= 1, print("On passe en dimension ",n+2));
    814843    codim += 2;
    815844    subspace2 = matrix( n+2, n, i, j, i == j);
    816845    d = prod( i = 1, matsize(Min[3])[1],Min[3][i,1]);    \\ d = abs(matdet(G1));
    817846    if( signG[2]%2 == 1, d = -d);                        \\ d = matdet(G1);
    818847    if( Min[3][1,1] == 2, factD = [-1], factD = [-1,2]); \\ si d est pair ...
    819848    factD = concat(factD, Min[3][,1]~);
    820 if( DEBUGLEVEL >= 4, print("factD=",factD));
     849if( DEBUGLEVEL_qfsolve >= 4, print("factD=",factD));
    821850
    822851\\ Solubilite locale en 2 (c'est le seul cas qui restait a traiter !!)
    823852    if( n == 4 && d%8 == 1,
    824853      if( QfWittinvariant(G,2) == 1,
    825 if( DEBUGLEVEL >= 1, print("pas de solution locale en 2"));
     854if( DEBUGLEVEL_qfsolve >= 1, print("pas de solution locale en 2"));
    826855        return(2)));
    827856 
    828857\\
     
    839868    if( n >= 5, dQ *= 8);
    840869   
    841870\\ invariants p-adiques
    842 \\ pour p=2, on ne choisit pas.
     871\\ pour p = 2, on ne choisit pas.
    843872    if( n == 4,
    844 if( DEBUGLEVEL >= 1, print("calcul des invariants locaux de G1"));
     873if( DEBUGLEVEL_qfsolve >= 1, print("calcul des invariants locaux de G1"));
    845874      aux = Qflisteinvariants(-G1,factD)[2][,1];
    846875      for( i = 3, length(factD), Winvariants[i] = aux[i])
    847876    ,
     
    850879    );
    851880    Winvariants[2] = sum( i = 1, length(factD), Winvariants[i])%2;
    852881
    853 if( DEBUGLEVEL >= 1,
     882if( DEBUGLEVEL_qfsolve >= 1,
    854883  print("Recherche d'un forme binaire de discriminant = ",dQ);
    855884  print("et d'invariants de Witt = ",Winvariants));
    856885
     
    866895    factd[1,2]++;
    867896    U2 = matrix(length(factD), n == 4, i,j, hilbert(2,dQ,factD[i])<0);
    868897    clgp2 = class2(dQ,factd,Winvariants,U2);
    869 if( DEBUGLEVEL >= 4, print("clgp2=",clgp2));
     898if( DEBUGLEVEL_qfsolve >= 4, print("clgp2=",clgp2));
    870899
    871900    clgp2 = clgp2[2];
    872901    U = Qflisteinvariants(clgp2,factD)[2];
    873902    if( n == 4, U = concat(U,U2));
    874 if( DEBUGLEVEL >= 4, printp("U=",U));
     903if( DEBUGLEVEL_qfsolve >= 4, printp("U=",U));
    875904
    876905    V = lift(matinverseimage(U*Mod(1,2),Winvariants*Mod(1,2)));
    877906    if( !length(V), next);
    878 if( DEBUGLEVEL >= 4, print("V=",V));
     907if( DEBUGLEVEL_qfsolve >= 4, print("V=",V));
    879908
    880909    if( dQ%2 == 1, Q = qfbprimeform(4*dQ,1), Q = qfbprimeform(dQ,1));
    881910    for( i = 1, length(clgp2),
     
    883912    Q = mymat(Q);
    884913    if( norml2(V) > 1, aux = QfbReduce(Q); Q = aux~*Q*aux);
    885914    if( n == 4 && V[length(V)], Q*=  2);
    886 if( DEBUGLEVEL >= 2, print("Q=",Q));
    887 if( DEBUGLEVEL >= 3, print("invariants de Witt de Q=",Qflisteinvariants(Q,factD)));
     915if( DEBUGLEVEL_qfsolve >= 2, print("Q=",Q));
     916if( DEBUGLEVEL_qfsolve >= 3, print("invariants de Witt de Q=",Qflisteinvariants(Q,factD)));
    888917
    889918\\
    890919\\ Construction d'une forme de dim=n+2 potentiellement unimodulaire
    891920\\
    892921
    893922    G2 = matdiagonalblock([G1,-Q]);
    894 if( DEBUGLEVEL >= 4, print("G2=",G2));
     923if( DEBUGLEVEL_qfsolve >= 4, print("G2=",G2));
    895924
    896 if( DEBUGLEVEL >= 2, print("minimisation de la forme quadratique de dimension ",length(G2)));
     925if( DEBUGLEVEL_qfsolve >= 2, print("minimisation de la forme quadratique de dimension ",length(G2)));
    897926\\ Minimisation de G2
    898927    detG2 = matdet(G2);
    899928    factd = matrix(length(factD)-1,2);
    900929    for( i = 1, length(factD)-1,
    901930      factd[i,2] = valuation(detG2, factd[i,1] = factD[i+1]));
    902 if( DEBUGLEVEL >= 3, print("det(G2)=",factd));
     931if( DEBUGLEVEL_qfsolve >= 3, print("det(G2) = ",factd));
    903932    Min = Qfminim(G2,factd);
    904933    M2 = Min[2]; G2 = Min[1];
    905934if( abs(matdet(G2)) > 2, print("******* ERREUR dans Qfsolve: det(G2) <> +-1 *******",matdet(G2));return(0));
    906 if( DEBUGLEVEL >= 4, print("G2=",G2));
     935if( DEBUGLEVEL_qfsolve >= 4, print("G2=",G2));
    907936
    908937\\ Maintenant det(G2) = +-1
    909938
    910939\\ On construit un seti pour G2 (Sous-Espace Totalement Isotrope)
    911 if( DEBUGLEVEL >= 2, print("recherche d'un espace de solutions pour G2"));
     940if( DEBUGLEVEL_qfsolve >= 2, print("recherche d'un espace de solutions pour G2"));
    912941    solG2 = LLLgoon(G2,1);
    913942    if( matrix(codim+1,codim+1,i,j,solG2[1][i,j]) != 0,
    914 if( DEBUGLEVEL >= 2, print(" pas assez de solutions dans G2"));return(0));
     943if( DEBUGLEVEL_qfsolve >= 2, print(" pas assez de solutions dans G2"));return(0))
    915944  );
    916945
    917946\\ G2 doit avoir un espace de solutions de dimension > codim
     
    921950print("dimseti = ",dimseti," <= codim = ",codim);
    922951print("************* ERREUR : pas assez de solutions pour G2"); return(0));   
    923952  solG2 = matrix(length(G2),dimseti,i,j,solG2[2][i,j]);
    924 if( DEBUGLEVEL >= 3, print("solG2=",solG2));
     953if( DEBUGLEVEL_qfsolve >= 3, print("solG2=",solG2));
    925954
    926955\\ La solution de G1 se trouve a la fois dans solG2 et dans subspace2
    927 if( DEBUGLEVEL >= 1, print("Reconstruction de la solution de G1"));
     956if( DEBUGLEVEL_qfsolve >= 1, print("Reconstruction de la solution de G1"));
    928957  solG1 = matintersect(subspace2,M2*solG2);
    929958  solG1 = subspace2~*solG1;
    930 if( DEBUGLEVEL >= 3, print("solG1 = ",solG1));
     959if( DEBUGLEVEL_qfsolve >= 3, print("solG1 = ",solG1));
    931960
    932961\\ La solution de G se trouve a la fois dans solG et dans subspace1
    933 if( DEBUGLEVEL >= 1, print("Reconstruction de la solution de G"));
     962if( DEBUGLEVEL_qfsolve >= 1, print("Reconstruction de la solution de G"));
    934963  sol = matintersect(subspace1,M1*solG1);
    935964  sol = subspace1~*sol;
    936965  sol = M*sol;
    937966  sol /= content(sol);
    938967  if( length(sol) == 1, sol = sol[,1]);
    939 if( DEBUGLEVEL >= 3, print("sol = ",sol));
    940 if( DEBUGLEVEL >= 1, print("fin de Qfsolve"));
     968if( DEBUGLEVEL_qfsolve >= 3, print("sol = ",sol));
     969if( DEBUGLEVEL_qfsolve >= 1, print("fin de Qfsolve"));
    941970  return(sol);
    942971}
    943 {matdiagonalblock(v)=
     972{matdiagonalblock(v) =
    944973local(lv,lt,M);
    945974  lv = length(v);
    946975  lt = sum( i = 1, lv, length(v[i]));
     
    950979    for( j = 1, length(v[i]),
    951980      for( k = 1, length(v[i]),
    952981        M[lt+j,lt+k] = v[i][j,k]));
    953     lt += length(v[i]);
     982    lt += length(v[i])
    954983  );
    955 return(M)
     984return(M);
    956985}
    957986
    958987