Ticket #6955: ell.gp

File ell.gp, 69.4 KB (added by John Cremona, 13 years ago)
Line 
1\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
2\\       Copyright (C) 2007 Denis Simon
3\\
4\\ Distributed under the terms of the GNU General Public License (GPL)
5\\
6\\    This code is distributed in the hope that it will be useful,
7\\    but WITHOUT ANY WARRANTY; without even the implied warranty of
8\\    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
9\\    General Public License for more details.
10\\
11\\ The full text of the GPL is available at:
12\\
13\\                 http://www.gnu.org/licenses/
14\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
15
16\\
17\\ Auteur:
18\\ Denis SIMON -> simon@math.unicaen.fr
19\\ adresse du fichier:
20\\ www.math.unicaen.fr/~simon/ell.gp
21\\
22\\  *********************************************
23\\  *          VERSION 25/03/2009               *
24\\  *********************************************
25\\
26\\ Programme de calcul du rang des courbes elliptiques
27\\ dans les corps de nombres.
28\\ langage: GP
29\\ pour l'utiliser, lancer gp, puis taper
30\\ \r ell.gp
31\\
32\\
33\\ Explications succintes :
34\\ definition du corps :
35\\ bnf=bnfinit(y^2+1);
36\\ (il est indispensable que la variable soit y).
37\\ on peut ensuite poser :
38\\ X = Mod(y,bnf.pol);
39\\
40\\ La fonction bnfellrank() accepte toutes les courbes sous la forme
41\\ [a1,a2,a3,a4,a6]
42\\ Les coefficients peuvent etre entiers ou non.
43\\ L'algorithme utilise est celui de la 2-descente.
44\\ La 2-torsion peut etre quelconque.
45\\ Il suffit de taper :
46\\
47\\ gp > ell = [a1,a2,a3,a4,a6];
48\\ gp > bnfellrank(bnf,ell)
49\\
50\\ Retourne un vecteur [r,s,vec]
51\\ ou r est le rang probable (c'est toujours une minoration du rang),
52\\ s est le 2-rang du groupe de Selmer,
53\\ vec est une liste de points dans E(K)/2E(K).
54\\
55\\ Courbes avec #E[2](K) >= 2:
56\\ ell doit etre sous la forme
57\\ y^2 = x^3 + A*^2 + B*x
58\\ avec A et B entiers algebriques
59\\ gp > ell = [0,A,0,B,0]
60\\ gp > bnfell2descent_viaisog(ell)
61\\ = algorithme de la 2-descente par isogenies
62\\ Attention A et B doivent etre entiers
63\\
64\\ Courbes avec #E[2](K) = 4: y^2 = (x-e1)*(x-e2)*(x-e3)
65\\ -> bnfell2descent_complete(bnf,e1,e2,e3);
66\\ = algorithme de la 2-descente complete
67\\ Attention: les ei doivent etre entiers algebriques.
68\\
69\\
70\\ On peut avoir plus ou moins de details de calculs avec
71\\ DEBUGLEVEL_ell = 0;
72\\ DEBUGLEVEL_ell = 1; 2; 3;...
73\\
74
75{
76\\
77\\ Variables globales usuelles
78\\
79
80  DEBUGLEVEL_ell = 1;   \\ pour avoir plus ou moins de details
81  LIM1 = 2;    \\ limite des points triviaux sur les quartiques
82  LIM3 = 4;    \\ limite des points sur les quartiques ELS
83  LIMTRIV = 2; \\ limite des points triviaux sur la courbe elliptique
84
85\\
86\\  Variables globales techniques
87\\
88
89  BIGINT = 32000;   \\ l'infini
90  MAXPROB = 20;
91  LIMBIGPRIME = 30; \\ pour distinguer un petit nombre premier d'un grand
92                    \\ utilise un test probabiliste pour les grands
93                    \\ si LIMBIGPRIME = 0, n'utilise aucun test probabiliste
94  NBIDEAUX = 10;
95
96}
97
98\\
99\\  Programmes
100\\
101
102\\
103\\ Fonctions communes ell.gp et ellQ.gp
104\\
105
106{
107ellinverturst(urst) =
108local(u = urst[1], r = urst[2], s = urst[3], t = urst[4]);
109  [1/u,-r/u^2,-s/u,(r*s-t)/u^3];
110}
111{
112ellchangecurveinverse(ell,v) = ellchangecurve(ell,ellinverturst(v));
113}
114{
115ellchangepointinverse(pt,v) = ellchangepoint(pt,ellinverturst(v));
116}
117{
118ellcomposeurst(urst1,urst2) =
119local(u1 = urst1[1], r1 = urst1[2], s1 = urst1[3], t1 = urst1[4],
120      u2 = urst2[1], r2 = urst2[2], s2 = urst2[3], t2 = urst2[4]);
121  [u1*u2,u1^2*r2+r1,u1*s2+s1,u1^3*t2+s1*u1^2*r2+t1];
122}
123if( DEBUGLEVEL_ell >= 4, print("mysubst"));
124{
125mysubst(polsu,subsx) =
126  if( type(lift(polsu)) == "t_POL",
127    return(simplify(subst(lift(polsu),variable(lift(polsu)),subsx)) )
128  , return(simplify(lift(polsu))));
129}
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}
162if( DEBUGLEVEL_ell >= 4, print("degre"));
163{
164degre(idegre) =
165local(ideg = idegre, jdeg = 0);
166
167  while( ideg >>= 1, jdeg++);
168  return(jdeg);
169}
170if( DEBUGLEVEL_ell >= 4, print("nfissquare"));
171{
172nfissquare(nf, a) = #nfsqrt(nf,a) > 0;
173}
174if( DEBUGLEVEL_ell >= 4, print("nfsqrt"));
175{
176nfsqrt( nf, a) =
177\\ si a est un carre, renvoie [sqrt(a)], sinon [].
178local(alift,ta,res,pfact);
179
180if( DEBUGLEVEL_ell >= 5, print("entree dans nfsqrt ",a));
181  if( a==0 || a==1,
182if( DEBUGLEVEL_ell >= 5, print("fin de nfsqrt"));
183    return([a]));
184
185  alift = lift(a);
186  ta = type(a);
187  if( !poldegree(alift), alift = polcoeff(alift,0));
188
189  if( type(alift) != "t_POL",
190    if( issquare(alift),
191if( DEBUGLEVEL_ell >= 5, print("fin de nfsqrt"));
192      return([sqrtrat(alift)])));
193
194  if( poldegree(nf.pol) <= 1,
195if( DEBUGLEVEL_ell >= 5, print("fin de nfsqrt"));
196    return([]));
197  if( ta == "t_POL", a = Mod(a,nf.pol));
198
199\\ tous les plgements reels doivent etre >0
200
201  for( i = 1, nf.r1,
202    if( nfsign(nf,a,i) < 0,
203if( DEBUGLEVEL_ell >= 5, print("fin de nfsqrt"));
204      return([])));
205
206\\ factorisation sur K du polynome X^2-a :
207
208  if( variable(nf.pol) == x,
209    py = subst(nf.pol,x,y);
210    pfact = lift(factornf(x^2-mysubst(alift,Mod(y,py)),py)[1,1])
211  ,
212    pfact = lift(factornf(x^2-a,nf.pol)[1,1]));
213  if( poldegree(pfact) == 2,
214if( DEBUGLEVEL_ell >= 5, print("fin de nfsqrt"));
215    return([]));
216if( DEBUGLEVEL_ell >= 5, print("fin de nfsqrt"));
217  return([subst(polcoeff(pfact,0),y,Mod(variable(nf.pol),nf.pol))]);
218}
219if( DEBUGLEVEL_ell >= 4, print("sqrtrat"));
220{
221sqrtrat(a) =
222  sqrtint(numerator(a))/sqrtint(denominator(a));
223}
224
225
226\\
227\\ Fonctions propres a ell.gp
228\\
229
230if( DEBUGLEVEL_ell >= 4, print("nfpolratroots"));
231{
232nfpolratroots(nf,pol) =
233local(f,ans);
234  f = nffactor(nf,lift(pol))[,1];
235  ans = [];
236  for( j = 1, #f,
237    if( poldegree(f[j]) == 1,
238      ans = concat(ans,[-polcoeff(f[j],0)/polcoeff(f[j],1)])));
239  return(ans);
240}
241if( DEBUGLEVEL_ell >= 4, print("nfmodid2"));
242{
243nfmodid2(nf,a,ideal) =
244if( DEBUGLEVEL_ell >= 5, print("entree dans nfmodid2"));
245\\ ideal doit etre sous la forme primedec
246  if( #nf.zk == 1,
247if( DEBUGLEVEL_ell >= 5, print("fin de nfmodid2"));
248    return(a*Mod(1,ideal.p)));
249  a = mynfeltmod(nf,a,nfbasistoalg(nf,ideal[2]));
250  if( gcd(denominator(content(lift(a))),ideal.p) == 1,
251if( DEBUGLEVEL_ell >= 5, print("fin de nfmodid2"));
252    return(a*Mod(1,ideal.p)));
253if( DEBUGLEVEL_ell >= 5, print("fin de nfmodid2"));
254  return(a);
255}
256if( DEBUGLEVEL_ell >= 4, print("nfhilb2"));
257{
258nfhilb2(nf,a,b,p) =
259local(res);
260
261if( DEBUGLEVEL_ell >= 5, print("entree dans nfhilb2"));
262  if( nfqpsoluble(nf,a*x^2+b,initp(nf,p)), res = 1, res = -1);
263if( DEBUGLEVEL_ell >= 5, print("fin de nfhilb2"));
264  return(res);
265}
266if( DEBUGLEVEL_ell >= 4, print("mynfhilbertp"));
267{
268mynfhilbertp(nf,a,b,p) =
269\\ calcule le symbole de Hilbert quadratique local (a,b)_p
270\\ * en l'ideal premier p du corps nf,
271\\ * a et b sont des elements non nuls de nf, sous la forme
272\\ * de polymods ou de polynomes, et p renvoye par primedec.
273local(alpha,beta,sig,aux,aux2,rep);
274
275if( DEBUGLEVEL_ell >= 5, print("entree dans mynfhilbertp ",p));
276  if( a == 0 || b == 0, print("0 argument in mynfhilbertp"));
277  if( p.p == 2,
278if( DEBUGLEVEL_ell >= 5, print("fin de mynfhilbertp"));
279    return(nfhilb2(nf,a,b,p)));
280  if( type(a) != "t_POLMOD", a = Mod(a,nf.pol));
281  if( type(b) != "t_POLMOD", b = Mod(b,nf.pol));
282 
283  alpha = idealval(nf,a,p); beta = idealval(nf,b,p);
284if( DEBUGLEVEL_ell >= 5, print("[alpha,beta] = ",[alpha,beta]));
285  if( (alpha%2 == 0) && (beta%2 == 0),
286if( DEBUGLEVEL_ell >= 5, print("fin de mynfhilbertp"));
287    return(1));
288  aux2 = idealnorm(nf,p)\2;
289  if( alpha%2 && beta%2 && aux2%2, sig = 1, sig = -1);
290  if( beta, aux = nfmodid2(nf,a^beta/b^alpha,p), aux = nfmodid2(nf,b^alpha,p));
291  aux = aux^aux2 + sig;
292  aux = lift(lift(aux));
293  if( aux == 0, rep = 1, rep = (idealval(nf,aux,p) >=  1) );
294if( DEBUGLEVEL_ell >= 5, print("fin de mynfhilbertp"));
295  if( rep, return(1), return(-1));
296}
297if( DEBUGLEVEL_ell >= 4, print("ideallistfactor"));
298{
299ideallistfactor(nf,listfact) =
300local(Slist,S1,test,i,j,k);
301
302if( DEBUGLEVEL_ell >= 5, print("entree dans ideallistfactor"));
303  Slist = []; test = 1;
304  for( i = 1, #listfact,
305    if( listfact[i] == 0, next);
306    S1 = idealfactor(nf,listfact[i])[,1];
307    for( j = 1, #S1, k = #Slist;
308      for( k = 1, #Slist,
309        if( Slist[k] == S1[j], test = 0; break));
310      if( test, Slist = concat(Slist,[S1[j]]), test = 1);
311     ));
312if( DEBUGLEVEL_ell >= 5, print("fin de ideallistfactor"));
313  return(Slist);
314}
315if( DEBUGLEVEL_ell >= 4, print("mynfhilbert"));
316{
317mynfhilbert(nf,a,b) =
318\\ calcule le symbole de Hilbert quadratique global (a,b):
319\\ =1 si l'equation X^2-aY^2-bZ^2=0 a une solution non triviale,
320\\ =-1 sinon,
321\\ a et b doivent etre non nuls.
322local(al,bl,S);
323
324if( DEBUGLEVEL_ell >= 4, print("entree dans mynfhilbert ",[a,b]));
325  if( a == 0 || b == 0, error("mynfhilbert : argument = 0"));
326  al = lift(a); bl = lift(b);
327
328\\ solutions locales aux places reelles
329
330  for( i = 1, nf.r1,
331    if( nfsign(nf,al,i) < 0 && nfsign(nf,bl,i) < 0,
332if( DEBUGLEVEL_ell >= 3, print("mynfhilbert non soluble a l'infini"));
333if( DEBUGLEVEL_ell >= 4, print("fin de mynfhilbert"));
334      return(-1))
335  );
336
337  if( type(a) != "t_POLMOD", a = Mod(a,nf.pol));
338  if( type(b) != "t_POLMOD", b = Mod(b,nf.pol));
339
340\\  solutions locales aux places finies (celles qui divisent 2ab)
341
342  S = ideallistfactor(nf,[2,a,b]);
343  forstep ( i = #S, 2, -1,
344\\ d'apres la formule du produit on peut eviter un premier
345    if( mynfhilbertp(nf,a,b, S[i]) == -1,
346if( DEBUGLEVEL_ell >= 3, print("mynfhilbert non soluble en : ",S[i]));
347if( DEBUGLEVEL_ell >= 4, print("fin de mynfhilbert"));
348      return(-1)));
349if( DEBUGLEVEL_ell >= 4, print("fin de mynfhilbert"));
350  return(1);
351}
352if( DEBUGLEVEL_ell >= 4, print("initp"));
353{
354initp( nf, p) =
355\\   pp[1] est l'ideal sous forme reduite
356\\   pp[2] est un entier de Zk avec une valuation 1 en p
357\\   pp[3] est la valuation de 2 en p
358\\   pp[4] sert a detecter les carres dans Qp
359\\     si p|2 il faut la structure de Zk/p^(1+2v) d'apres Hensel
360\\     sinon il suffit de calculer x^(N(p)-1)/2
361\\   pp[5] est un systeme de representants de Zk/p
362\\     c'est donc un ensemble de cardinal p^f .
363local(idval,pp);
364
365if( DEBUGLEVEL_ell >= 5, print("entree dans initp"));
366  idval = idealval(nf,2,p);
367  pp=[ p, nfbasistoalg(nf,p[2]), idval, 0, repres(nf,p) ];
368  if( idval,
369    pp[4] = idealstar(nf,idealpow(nf,p,1+2*idval)),
370    pp[4] = p.p^p.f\2 );
371if( DEBUGLEVEL_ell >= 5, print("fin de initp"));
372  return(pp);
373}
374if( DEBUGLEVEL_ell >= 4, print("deno"));
375{
376deno(num) =
377\\ calcule un denominateur du polynome num
378
379  if( num == 0, return(1));
380  if( type(num) == "t_POL",
381     return(denominator(content(num))));
382  return(denominator(num));
383}
384if( DEBUGLEVEL_ell >= 4, print("nfratpoint"));
385{
386nfratpoint(nf,pol,lim,singlepoint=1) =
387\\ Si singlepoint == 1, cherche un seul point, sinon plusieurs.
388local(compt1,compt2,deg,n,AA,point,listpoints,denoz,vectx,xx,evpol,sq);
389
390if( DEBUGLEVEL_ell >= 4, print("entree dans nfratpoint avec pol = ",pol); print("lim = ",lim));
391  compt1 = 0; compt2 = 0;
392  deg = poldegree(pol); n = poldegree(nf.pol);
393  AA = lim<<1;
394  if( !singlepoint, listpoints = []);
395
396\\ cas triviaux
397  sq = nfsqrt(nf,polcoeff(pol,0));
398  if( sq!= [],
399    point = [ 0, sq[1], 1];
400    if( singlepoint,
401if( DEBUGLEVEL_ell >= 4, print("fin de nfratpoint"));
402      return(point));
403    listpoints = concat(listpoints,[point])
404  );
405  sq = nfsqrt(nf,pollead(pol));
406  if( sq != [],
407    point = [ 1, sq[1], 0];
408    if( singlepoint,
409if( DEBUGLEVEL_ell >= 4, print("fin de nfratpoint"));
410      return(point));
411    listpoints = concat(listpoints,[point])
412  );
413 
414\\ boucle generale
415  point = [];
416  vectx = vector(n,i,[-lim,lim]);
417  for( denoz = 1, lim,
418    forvec( xx = vectx,
419      if( denoz == 1 || gcd(content(xx),denoz) == 1,
420        xpol = nfbasistoalg(nf,xx~);
421        evpol = subst(pol,x,xpol/denoz);
422        sq = nfsqrt(nf,evpol);
423        if( sq != [],
424          point = [xpol/denoz, sq[1], 1];
425          if( singlepoint, break(2));
426          listpoints = concat(listpoints,[point])));
427  ));
428
429  if( singlepoint, listpoints = point);
430if( DEBUGLEVEL_ell >= 4, print("sortie de nfratpoint"));
431if( DEBUGLEVEL_ell >= 3, print("points trouves par nfratpoint = ",listpoints));
432  return(listpoints);
433}
434if( DEBUGLEVEL_ell >= 4, print("repres"));
435{
436repres(nf,p) =
437\\ calcule un systeme de representants Zk/p
438local(fond,mat,i,j,k,f,rep,pp,ppi,pp2,jppi,gjf);
439
440if( DEBUGLEVEL_ell >= 5, print("entree dans repres"));
441  fond = [];
442  mat = idealhnf(nf,p);
443  for( i = 1, #mat,
444    if( mat[i,i] != 1, fond = concat(fond,nf.zk[i])));
445  f = #fond;
446  pp = p.p;
447  rep = vector(pp^f,i,0);
448  rep[1] = 0;
449  ppi = 1;
450  pp2 = pp\2;
451  for( i = 1, f,
452    for( j = 1, pp-1,
453      if( j <= pp2, gjf = j*fond[i], gjf = (j-pp)*fond[i]);
454      jppi = j*ppi;
455      for( k = 0, ppi-1, rep[jppi+k+1] = rep[k+1]+gjf ));
456    ppi *= pp);
457if( DEBUGLEVEL_ell >= 5, print("fin de repres"));
458  return(Mod(rep,nf.pol));
459}
460if( DEBUGLEVEL_ell >= 4, print("val"));
461{
462val(nf,num,p) =
463  if( num == 0, BIGINT, idealval(nf,lift(num),p));
464}
465if( DEBUGLEVEL_ell >= 4, print("nfissquarep"));
466{
467nfissquarep(nf,a,p,q) =
468\\ suppose que a est un carre modulo p^q
469\\ et renvoie sqrt(a) mod p^q (ou plutot p^(q/2))
470local(pherm,f,aaa,n,pp,qq,e,z,xx,yy,r,aux,b,m,vp,inv2x,zinit,zlog,expo);
471
472if( DEBUGLEVEL_ell >= 5, print("entree dans nfissquarep ",a,p,q));
473  if( a == 0 || a == 1,
474if( DEBUGLEVEL_ell >= 4, print("fin de nfissquarep"));
475    return(a));
476  pherm = idealhnf(nf,p);
477if( DEBUGLEVEL_ell >= 5, print("pherm = ",pherm));
478  f = idealval(nf,a,p);
479  if( f >= q,
480    if( f > q, aaa = nfbasistoalg(nf,p[2])^((q+1)>>1), aaa = 0);
481if( DEBUGLEVEL_ell >= 4, print("fin de nfissquarep"));
482    return(aaa));
483  if( f, aaa = a*nfbasistoalg(nf,p[5]/p.p)^f, aaa = a);
484  if( pherm[1,1] != 2,
485\\ cas ou p ne divise pas 2
486\\ algorithme de Shanks
487    n = nfrandintmodid(nf,pherm);
488    while( nfpsquareodd(nf,n,p), n = nfrandintmodid(nf,pherm));
489    pp = Mod(1,p.p);
490    n *= pp;
491    qq = idealnorm(nf,pherm)\2;
492    e = 1; while( !(qq%2), e++; qq \= 2);
493    z = mynfeltreduce(nf,lift(lift(n^qq)),pherm);
494    yy = z;r = e;
495    xx = mynfeltreduce(nf,lift(lift((aaa*pp)^(qq\2))),pherm);
496    aux = mynfeltreduce(nf,aaa*xx,pherm);
497    b = mynfeltreduce(nf,aux*xx,pherm);
498    xx = aux;
499    aux = b;m = 0;
500    while( !val(nf,aux-1,p), m++; aux = mynfeltreduce(nf,aux^2,pherm));
501    while( m,
502      if( m == r, error("nfissquarep : m = r"));
503      yy *= pp;
504      aux = mynfeltreduce(nf,lift(lift(yy^(1<<(r-m-1)))),pherm);
505      yy = mynfeltreduce(nf,aux^2,pherm);
506      r = m;
507      xx = mynfeltreduce(nf,xx*aux,pherm);
508      b = mynfeltreduce(nf,b*yy,pherm);
509      aux = b;m = 0;
510      while( !val(nf,aux-1,p), m++; aux = mynfeltreduce(nf,aux^2,pherm));
511    );
512\\ lift de Hensel
513\\
514    if( q > 1,
515      vp = idealval(nf,xx^2-aaa,p);
516      if( vp < q-f,
517        yy = 2*xx;
518        inv2x = nfbasistoalg(nf,idealaddtoone(nf,yy,p)[1])/yy;
519        while( vp < q, vp++; xx -= (xx^2-aaa)*inv2x);
520      );
521      if( f, xx *= nfbasistoalg(nf,p[2])^(f>>1));
522    );
523    xx = mynfeltreduce(nf,xx,idealpow(nf,p,q))
524  ,
525\\ cas ou p divise 2 */
526    if( q-f > 1, id = idealpow(nf,p,q-f), id = pherm);
527    zinit = idealstar(nf,id,2);
528    zlog = ideallog(nf,aaa,zinit);
529    xx = 1;
530    for( i = 1, #zlog,
531      expo = zlog[i];
532      if( 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
539      )
540    );
541    if( f,
542      xx *= nfbasistoalg(nf,p[2])^(f>>1);
543      id = idealpow(nf,p,q));
544    xx = mynfeltreduce(nf,xx,id);
545  );
546if( DEBUGLEVEL_ell >= 4, print("fin de nfissquarep ",xx));
547  return(xx);
548}
549if( DEBUGLEVEL_ell >= 4, print("nfpsquareodd"));
550{
551nfpsquareodd( nf, a, p) =
552\\ renvoie 1 si a est un carre dans ZK_p 0 sinon
553\\ seulement pour p premier avec 2
554local(v,ap,norme,den);
555
556if( DEBUGLEVEL_ell >= 5, print("entree dans nfpsquareodd"));
557  if( a == 0,
558if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquareodd"));
559    return(1));
560  v = idealval(nf,lift(a),p);
561  if( v%2,
562if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquareodd"));
563    return(0));
564  ap = a/nfbasistoalg(nf,p[2])^v;
565
566  norme = idealnorm(nf,p)\2;
567  den = denominator(content(lift(ap)))%p.p;
568  if(sign(den), ap*=Mod(1,p.p));
569  ap = ap^norme-1;
570  if( ap == 0,
571if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquareodd"));
572    return(1));
573  ap = lift(lift(ap));
574  if( idealval(nf,ap,p) > 0,
575if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquareodd"));
576    return(1));
577if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquareodd"));
578  return(0);
579}
580if( DEBUGLEVEL_ell >= 4, print("nfpsquare"));
581{
582nfpsquare( nf, a, p, zinit) =
583\\ a est un entier de K
584\\ renvoie 1 si a est un carre dans ZKp 0 sinon
585local(valap,zlog,i);
586
587if( DEBUGLEVEL_ell >= 5, print("entree dans nfpsquare ",[a,p,zinit]));
588  if( a == 0,
589if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquare"));
590    return(1));
591
592  if( p.p != 2,
593if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquare"));
594    return(nfpsquareodd(nf,a,p)));
595
596  valap = idealval(nf,a,p);
597  if( valap%2,
598if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquare"));
599    return(0));
600  if( valap,
601    zlog = ideallog(nf,a*(nfbasistoalg(nf,p[5])/p.p)^valap,zinit)
602  ,
603    zlog = ideallog(nf,a,zinit));
604  for( i = 1, #zinit[2][2],
605    if( !(zinit[2][2][i]%2) && (zlog[i]%2),
606if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquare"));
607      return(0)));
608if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquare"));
609  return(1);
610}
611if( DEBUGLEVEL_ell >= 4, print("nfpsquareq"));
612{
613nfpsquareq( nf, a, p, q) =
614\\ cette fonction renvoie 1 si a est un carre
615\\ ?inversible? modulo P^q et 0 sinon.
616\\ P divise 2, et ?(a,p)=1?.
617local(vala,zinit,zlog,i);
618
619if( DEBUGLEVEL_ell >= 5, print("entree dans nfpsquareq ",[a,p,q]));
620  if( a == 0,
621if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquareq"));
622    return(1));
623  vala = idealval(nf,a,p);
624  if( vala >= q,
625if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquareq"));
626    return(1));
627  if( vala%2,
628if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquareq"));
629    return(0));
630  zinit = idealstar(nf,idealpow(nf,p,q-vala),2);
631  zlog = ideallog(nf,a*nfbasistoalg(nf,p[5]/2)^vala,zinit);
632  for( i = 1, #zinit[2][2],
633    if( !(zinit[2][2][i]%2) && (zlog[i]%2),
634if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquareq"));
635      return(0)));
636if( DEBUGLEVEL_ell >= 5, print("fin de nfpsquareq"));
637  return(1);
638}
639if( DEBUGLEVEL_ell >= 4, print("nflemma6"));
640{
641nflemma6( nf, pol, p, nu, xx) =
642local(gx,gpx,lambda,mu);
643
644if( DEBUGLEVEL_ell >= 5, print("entree dans nflemma6"));
645  gx = subst( pol, x, xx);
646  if( nfpsquareodd(nf,gx,p),
647if( DEBUGLEVEL_ell >= 5, print("fin de nflemma6"));
648    return(1));
649  gpx = subst( pol', x, xx);
650  lambda = val(nf,gx,p);mu = val(nf,gpx,p);
651
652  if( lambda>2*mu,
653if( DEBUGLEVEL_ell >= 5, print("fin de nflemma6"));
654    return(1));
655  if( (lambda >= 2*nu)  && (mu >= nu),
656if( DEBUGLEVEL_ell >= 5, print("fin de nflemma6"));
657    return(0));
658if( DEBUGLEVEL_ell >= 5, print("fin de nflemma6"));
659  return(-1);
660}
661if( DEBUGLEVEL_ell >= 4, print("nflemma7"));
662{
663nflemma7( nf, pol, p, nu, xx, zinit) =
664local(gx,gpx,v,lambda,mu,q);
665
666if( DEBUGLEVEL_ell >= 5, print("entree dans nflemma7 ",[xx,nu]));
667  gx = subst( pol, x, xx);
668  if( nfpsquare(nf,gx,p,zinit),
669if( DEBUGLEVEL_ell >= 5, print("fin de nflemma7"));
670    return(1));
671  gpx = subst( pol', x, xx);
672  v = p[3];
673  lambda = val(nf,gx,p);mu = val(nf,gpx,p);
674  if( lambda>2*mu,
675if( DEBUGLEVEL_ell >= 5, print("fin de nflemma7"));
676    return(1));
677  if( nu > mu,
678    if( lambda%2,
679if( DEBUGLEVEL_ell >= 5, print("fin de nflemma7"));
680      return(-1));
681    q = mu+nu-lambda;
682    if( q > 2*v,
683if( DEBUGLEVEL_ell >= 5, print("fin de nflemma7"));
684      return(-1));
685    if( nfpsquareq(nf,gx*nfbasistoalg(nf,p[5]/2)^lambda,p,q),
686if( DEBUGLEVEL_ell >= 5, print("fin de nflemma7"));
687      return(1))
688  ,
689    if( lambda >= 2*nu,
690if( DEBUGLEVEL_ell >= 5, print("fin de nflemma7"));
691      return(0));
692    if( lambda%2,
693if( DEBUGLEVEL_ell >= 5, print("fin de nflemma7"));
694      return(-1));
695    q = 2*nu-lambda;
696    if( q > 2*v,
697if( DEBUGLEVEL_ell >= 5, print("fin de nflemma7"));
698      return(-1));
699    if( nfpsquareq(nf,gx*nfbasistoalg(nf,p[5]/2)^lambda,p,q),
700if( DEBUGLEVEL_ell >= 5, print("fin de nflemma7"));
701      return(0))
702  );
703if( DEBUGLEVEL_ell >= 5, print("fin de nflemma7"));
704  return(-1);
705}
706if( DEBUGLEVEL_ell >= 4, print("nfzpsoluble"));
707{
708nfzpsoluble( nf, pol, p, nu, pnu, x0) =
709local(result,pnup,lrep);
710
711if( DEBUGLEVEL_ell >= 5, print("entree dans nfzpsoluble ",[lift(x0),nu]));
712  if( p[3] == 0,
713    result = nflemma6(nf,pol,p[1],nu,x0),
714    result = nflemma7(nf,pol,p[1],nu,x0,p[4]));
715  if( result == +1,
716if( DEBUGLEVEL_ell >= 5, print("fin de nfzpsoluble"));
717    return(1));
718  if( result == -1,
719if( DEBUGLEVEL_ell >= 5, print("fin de nfzpsoluble"));
720    return(0));
721  pnup = pnu*p[2];
722  lrep = #p[5];
723  nu++;
724  for( i = 1, lrep,
725    if( nfzpsoluble(nf,pol,p,nu,pnup,x0+pnu*p[5][i]),
726if( DEBUGLEVEL_ell >= 5, print("fin de nfzpsoluble"));
727      return(1)));
728if( DEBUGLEVEL_ell >= 5, print("fin de nfzpsoluble"));
729  return(0);
730}
731if( DEBUGLEVEL_ell >= 4, print("mynfeltmod"));
732{
733mynfeltmod(nf,a,b) =
734local(qred);
735
736  qred = round(nfalgtobasis(nf,a/b));
737  qred = a-b*nfbasistoalg(nf,qred);
738  return(qred);
739}
740if( DEBUGLEVEL_ell >= 4, print("mynfeltreduce"));
741{
742mynfeltreduce(nf,a,id) =
743  nfbasistoalg(nf,nfeltreduce(nf,nfalgtobasis(nf,a),id));
744}
745if( DEBUGLEVEL_ell >= 4, print("nfrandintmodid"));
746{
747nfrandintmodid( nf, id) =
748local(res);
749
750if( DEBUGLEVEL_ell >= 5, print("entree dans nfrandintmodid"));
751  res = 0;
752  while( !res,
753    res = nfrandint(nf,0);
754    res = mynfeltreduce(nf,res,id));
755if( DEBUGLEVEL_ell >= 5, print("fin de nfrandintmodid"));
756  return(res);
757}
758if( DEBUGLEVEL_ell >= 4, print("nfrandint"));
759{
760nfrandint( nf, borne) =
761local(l,res,i);
762
763if( DEBUGLEVEL_ell >= 5, print("entree dans nfrandint"));
764  l = #nf.zk;
765  res = vectorv(l,i,0);
766  for( i = 1, l,
767      if( borne, res[i] = random(borne<<1)-borne, res[i] = random() ));
768  res = nfbasistoalg(nf,res);
769if( DEBUGLEVEL_ell >= 5, print("fin de nfrandint"));
770  return(res);
771}
772if( DEBUGLEVEL_ell >= 4, print("nfqpsolublebig"));
773{
774nfqpsolublebig( nf, pol, p,ap=0,b=1) =
775local(deg,i,xx,z,Px,j,cont,pi,pol2,Roots);
776
777if( DEBUGLEVEL_ell >= 4, print("entree dans nfqpsolublebig avec ",p.p));
778  deg = poldegree(pol);
779
780  if( nfpsquareodd(nf,polcoeff(pol,0),p),
781if( DEBUGLEVEL_ell >= 4, print("fin de nfqpsolublebig"));
782    return(1));
783  if( nfpsquareodd(nf,pollead(pol),p),
784if( DEBUGLEVEL_ell >= 4, print("fin de nfqpsolublebig"));
785    return(1));
786
787\\ on tient compte du contenu de pol
788  cont = idealval(nf,polcoeff(pol,0),p);
789  for( i = 1, deg,
790    if( cont, cont = min(cont,idealval(nf,polcoeff(pol,i),p))));
791  if( cont, pi = nfbasistoalg(nf,p[5]/p.p));
792  if( cont > 1, pol *= pi^(2*(cont\2)));
793   
794\\ On essaye des valeurs de x au hasard
795  if( cont%2,
796    pol2 = pol*pi
797  , pol2 = pol;
798    for( i = 1, MAXPROB,
799      xx = nfrandint(nf,0);
800      z = 0; while( !z, z = random());
801      xx = -ap*z+b*xx;
802      Px=polcoeff(pol,deg);
803      forstep (j=deg-1,0,-1,Px=Px*xx+polcoeff(pol,j));
804      Px *= z^(deg);
805      if( nfpsquareodd(nf,Px,p),
806if( DEBUGLEVEL_ell >= 4, print("fin de nfqpsolublebig"));
807        return(1));
808    )
809  );
810
811\\ On essaye les racines de pol
812  Roots = nfpolrootsmod(nf,pol2,p);
813  pi = nfbasistoalg(nf,p[2]);
814  for( i = 1, #Roots,
815    if( nfqpsolublebig(nf,subst(pol,x,pi*x+Roots[i]),p),
816      return(1)));
817
818if( DEBUGLEVEL_ell >= 4, print("fin de nfqpsolublebig"));
819  return(0);
820}
821if( DEBUGLEVEL_ell >= 4, print("nfpolrootsmod"));
822{
823nfpolrootsmod(nf,pol,p) =
824\\ calcule les racines modulo l'ideal p du polynome pol.
825\\ p est un ideal premier de nf, sous la forme idealprimedec
826local(factlist,sol);
827
828  factlist = nffactormod(nf,pol,p)[,1];
829  sol = [];
830  for( i = 1, #factlist,
831    if( poldegree(factlist[i]) == 1,
832      sol = concat(sol, [-polcoeff(factlist[i],0)/polcoeff(factlist[i],1)])));
833  return(sol);
834}
835if( DEBUGLEVEL_ell >= 4, print("nfqpsoluble"));
836{
837nfqpsoluble( nf, pol, p) =
838
839if( DEBUGLEVEL_ell >= 4, print("entree dans nfqpsoluble ",p));
840if( DEBUGLEVEL_ell >= 5, print("pol = ",pol));
841  if( nfpsquare(nf,pollead(pol),p[1],p[4]),
842if( DEBUGLEVEL_ell >= 5, print("fin de nfqpsoluble"));
843    return(1));
844  if( nfpsquare(nf,polcoeff(pol,0),p[1],p[4]),
845if( DEBUGLEVEL_ell >= 5, print("fin de nfqpsoluble"));
846    return(1));
847  if( nfzpsoluble(nf,pol,p,0,1,0),
848if( DEBUGLEVEL_ell >= 5, print("fin de nfqpsoluble"));
849    return(1));
850  if( nfzpsoluble(nf,polrecip(pol),p,1, p[2],0),
851if( DEBUGLEVEL_ell >= 5, print("fin de nfqpsoluble"));
852    return(1));
853if( DEBUGLEVEL_ell >= 5, print("fin de nfqpsoluble"));
854  return(0);
855}
856if( DEBUGLEVEL_ell >= 4, print("nflocallysoluble"));
857{
858nflocallysoluble( nf, pol, r=0,a=1,b=1) =
859local(pol0,plist,add,ff,p,Delta,vecpol,vecpolr,Sturmr);
860
861if( DEBUGLEVEL_ell >= 4, print("entree dans nflocallysoluble ",[pol,r,a,b]));
862  pol0 = pol;
863\\
864\\ places finies de plist */
865\\
866  pol *= deno(content(lift(pol)))^2;
867  for( ii = 1, 3,
868    if( ii == 1, plist = idealprimedec(nf,2));
869    if( ii == 2 && r, plist = idealfactor(nf,poldisc(pol0/pollead(pol0))/pollead(pol0)^6/2^12)[,1]);
870    if( ii == 2 && !r, plist = idealfactor(nf,poldisc(pol0))[,1]);
871    if( ii == 3,
872      add = idealadd(nf,a,b);
873      ff = factor(idealnorm(nf,add))[,1];
874      addprimes(ff);
875if( DEBUGLEVEL_ell >= 4, print("liste de premiers = ",ff));
876      plist = idealfactor(nf,add)[,1]);
877      for( i = 1, #plist,
878        p =  plist[i]; 
879if( DEBUGLEVEL_ell >= 3, print("p = ",p));
880        if( p.p < LIMBIGPRIME,
881          if( !nfqpsoluble(nf,pol,initp(nf,p)),
882if( DEBUGLEVEL_ell >= 2, print(" non ELS en ",p));
883if( DEBUGLEVEL_ell >= 4, print("fin de nflocallysoluble"));
884            return(0)),
885          if( !nfqpsolublebig(nf,pol,p,r/a,b),
886if( DEBUGLEVEL_ell >= 2, print(" non ELS en ",p.p," ( = grand premier  )"));
887if( DEBUGLEVEL_ell >= 4, print("fin de nflocallysoluble"));
888            return(0))));
889);
890\\ places reelles
891  if( nf.r1,
892    Delta = poldisc(pol); vecpol = Vec(pol);
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]));
898      Sturmr = polsturm(Pol(vecpolr));
899      if( Sturmr == 0,
900if( DEBUGLEVEL_ell >= 2, print(" non ELS a l'infini"));
901if( DEBUGLEVEL_ell >= 4, print("fin de nflocallysoluble"));
902        return(0));
903  ));
904if( DEBUGLEVEL_ell >= 2, print(" quartique ELS "));
905if( DEBUGLEVEL_ell >= 4, print("fin de nflocallysoluble"));
906  return(1);
907}
908if( DEBUGLEVEL_ell >= 4, print("nfellcount"));
909{
910nfellcount( nf, c, d, KS2gen, pointstriv) =
911local(found,listgen,listpointscount,m1,m2,lastloc,mask,i,d1,iaux,j,triv,pol,point,deuxpoints);
912
913if( DEBUGLEVEL_ell >= 4, print("entree dans nfellcount ",[c,d]));
914  found = 0;
915  listgen = KS2gen;
916  listpointscount = [];
917
918  m1 = m2 = 0; lastloc = -1;
919
920  mask = 1 << #KS2gen;
921  i = 1;
922  while( i < mask,
923    d1 = 1; iaux = i; j = 1;
924    while( iaux,
925      if( iaux%2, d1 *= listgen[j]);
926      iaux >>= 1; j++);
927if( DEBUGLEVEL_ell >= 2, print("d1 = ",d1));
928    triv = 0;
929    for( j = 1, #pointstriv,
930      if( pointstriv[j][3]*pointstriv[j][1]
931        && nfissquare(nf,d1*pointstriv[j][1]*pointstriv[j][3]),
932        listpointscount = concat(listpointscount,[pointstriv[j]]);
933if( DEBUGLEVEL_ell >= 2, print("point trivial"));
934        triv = 1; m1++;
935        if( degre(i) > lastloc, m2++);
936        found = 1; lastloc = -1; break));
937    if( !triv,
938    pol = Pol([d1,0,c,0,d/d1]);
939if( DEBUGLEVEL_ell >= 3, print("quartique = y^2 = ",pol));
940    point = nfratpoint(nf,pol,LIM1,1);
941    if( point != [],
942if( DEBUGLEVEL_ell >= 2, print("point sur la quartique"));
943if( DEBUGLEVEL_ell >= 3, print(point));
944      m1++;
945      if( point[3] != 0,
946        aux = d1*point[1]/point[3]^2;
947        deuxpoints = [ aux*point[1], aux*point[2]/point[3] ]
948      ,
949        deuxpoints = [0]);
950      listpointscount = concat(listpointscount,[deuxpoints]);   
951      if( degre(i) > lastloc, m2++);
952      found = 1; lastloc = -1
953    ,
954      if( nflocallysoluble(nf,pol),
955        if( degre(i) > lastloc, m2++; lastloc = degre(i));
956        point = nfratpoint(nf,pol,LIM3,1);
957        if( point != [],
958if( DEBUGLEVEL_ell >= 2, print("point sur la quartique"));
959if( DEBUGLEVEL_ell >= 3, print(point));
960          m1++;
961          aux = d1*point[1]/point[3]^2;
962          deuxpoints = [ aux*point[1], aux*point[2]/point[3] ];
963          listpointscount = concat(listpointscount,[deuxpoints]);
964          if( degre(i) > lastloc, m2++);
965          found = 1; lastloc = -1
966        ,
967if( DEBUGLEVEL_ell >= 2, print("pas de point trouve sur la quartique"));
968          ))));
969    if( found,
970      found = 0;
971      v = 0; iaux = (i>>1);
972      while( iaux, iaux >>= 1; v++);
973      mask >>= 1;
974      listgen = vecextract(listgen,(1<<#listgen)-(1<<v)-1);
975      i = (1<<v)
976    , i++)
977  );
978  for( i = 1, #listpointscount,
979   if( #listpointscount[i] > 1,
980    if( subst(x^3+c*x^2+d*x,x,listpointscount[i][1])-listpointscount[i][2]^2 != 0,
981      error("nfellcount : MAUVAIS POINT = ",listpointscount[i]))));
982if( DEBUGLEVEL_ell >= 4, print("fin de nfellcount"));
983  return([listpointscount,[m1,m2]]);
984}
985if( DEBUGLEVEL_ell >= 4, print("bnfell2descent_viaisog"));
986{
987bnfell2descent_viaisog( bnf, ell) =
988\\ Calcul du rang des courbes elliptiques avec 2-torsion
989\\ dans le corps de nombres bnf
990\\ par la methode des 2-isogenies.
991\\
992\\ ell = [a1,a2,a3,a4,a6]
993\\ y^2+a1xy+a3y=x^3+a2x^2+a4x+a6
994\\
995\\ ell doit etre sous la forme
996\\ y^2=x^3+ax^2+bx -> ell = [0,a,0,b,0]
997\\ avec a et b entiers.
998local(i,P,Pfact,tors,pointstriv,apinit,bpinit,plist,KS2prod,oddclass,KS2gen,listpoints,pointgen,n1,n2,certain,np1,np2,listpoints2,aux1,aux2,certainp,rang,strange);
999
1000if( DEBUGLEVEL_ell >= 2, print("Algorithme de la 2-descente par isogenies"));
1001if( DEBUGLEVEL_ell >= 3, print("entree dans bnfell2descent_viaisog"));
1002  if( variable(bnf.pol) != y,
1003    error(" bnfell2descent_viaisog : la variable du corps de nombres doit etre y "));
1004  ell = ellinit(Mod(lift(ell),bnf.pol),1);
1005
1006  if( ell.disc == 0,
1007    error(" bnfell2descent_viaisog : courbe singuliere !!"));
1008  if( ell.a1 != 0 || ell.a3 != 0 || ell.a6 != 0,
1009    error(" bnfell2descent_viaisog : la courbe n'est pas sous la forme [0,a,0,b,0]"));
1010  if( denominator(nfalgtobasis(bnf,ell.a2)) > 1 || denominator(nfalgtobasis(bnf,ell.a4)) > 1,
1011    error(" bnfell2descent_viaisog : coefficients non entiers"));
1012
1013  P = Pol([1,ell.a2,ell.a4])*Mod(1,bnf.pol);
1014  Pfact = factornf(P,bnf.pol)[,1];
1015  tors = #Pfact;
1016  if( #Pfact > 1,
1017    pointstriv = [[0,0,1],[-polcoeff(Pfact[1],0),0,1],[-polcoeff(Pfact[2],0),0,1]]
1018  , pointstriv = [[0,0,1]]);
1019
1020  apinit = -2*ell.a2; bpinit = ell.a2^2-4*ell.a4;
1021
1022\\ calcul des ideaux premiers de plist
1023\\ et de quelques renseignements associes
1024  plist = idealfactor(bnf,6*ell.disc)[,1];
1025
1026if( DEBUGLEVEL_ell >= 3, print(" Recherche de points triviaux sur la courbe"));
1027  P *= x;
1028if( DEBUGLEVEL_ell >= 3, print("Y^2 = ",P));
1029  pointstriv = concat( pointstriv, nfratpoint(bnf.nf,P,LIMTRIV,0));
1030if( DEBUGLEVEL_ell >= 1, print("points triviaux sur E(K) = ");
1031  print(lift(pointstriv)); print());
1032
1033  KS2prod = ell.a4;
1034  oddclass = 0;
1035  while( !oddclass,
1036    KS2gen = bnfsunit(bnf,idealfactor(bnf,KS2prod)[,1]~);
1037    oddclass = KS2gen[5][1]%2;
1038    if( !oddclass,
1039      KS2prod = idealmul(bnf,KS2prod,(KS2gen[5][3][1])));
1040 );
1041  KS2gen = KS2gen[1];
1042\\  A CHANGER : KS2gen = matbasistoalg(bnf,KS2gen);
1043  for( i = 1, #KS2gen,
1044    KS2gen[i] = nfbasistoalg(bnf, KS2gen[i]));
1045  KS2gen = concat(Mod(lift(bnf.tufu),bnf.pol),KS2gen);
1046if( DEBUGLEVEL_ell >= 2, 
1047  print("#K(b,2)gen          = ",#KS2gen);
1048  print("K(b,2)gen = ",KS2gen));
1049
1050  listpoints = nfellcount(bnf.nf,ell.a2,ell.a4,KS2gen,pointstriv);
1051  pointgen = listpoints[1];
1052if( DEBUGLEVEL_ell >= 1, print("points sur E(K) = ",lift(pointgen)); print());
1053  n1 = listpoints[2][1]; n2 = listpoints[2][2];
1054
1055  certain = (n1 == n2);
1056if( DEBUGLEVEL_ell >= 1,
1057  if( certain,
1058    print("[E(K):phi'(E'(K))]  = ",1<<n1);
1059    print("#S^(phi')(E'/K)     = ",1<<n2);
1060    print("#III(E'/K)[phi']    = 1"); print()
1061  ,
1062    print("[E(K):phi'(E'(K))] >= ",1<<n1);
1063    print("#S^(phi')(E'/K)     = ",1<<n2);
1064    print("#III(E'/K)[phi']   <= ",1<<(n2-n1)); print())
1065);
1066
1067  KS2prod = bpinit;
1068  oddclass = 0;
1069  while( !oddclass,
1070    KS2gen = bnfsunit(bnf,idealfactor(bnf,KS2prod)[,1]~);
1071    oddclass = (KS2gen[5][1]%2);
1072    if( !oddclass,
1073      KS2prod = idealmul(bnf,KS2prod,(KS2gen[5][3][1]))));
1074  KS2gen = KS2gen[1];
1075\\ A CHANGER KS2gen = matbasistoalg(bnf,KS2gen);
1076  for( i = 1, #KS2gen,
1077    KS2gen[i] = nfbasistoalg(bnf, KS2gen[i]));
1078  KS2gen = concat(Mod(lift(bnf.tufu),bnf.pol),KS2gen);
1079if( DEBUGLEVEL_ell >= 2, 
1080  print("#K(a^2-4b,2)gen     = ",#KS2gen);
1081  print("K(a^2-4b,2)gen     = ",KS2gen));
1082
1083  P = Pol([1,apinit,bpinit])*Mod(1,bnf.pol);
1084  Pfact= factornf(P,bnf.pol)[,1];
1085  if( #Pfact > 1,
1086    pointstriv=[[0,0,1],[-polcoeff(Pfact[1],0),0,1],[-polcoeff(Pfact[2],0),0,1]]
1087  , pointstriv = [[0,0,1]]);
1088
1089if( DEBUGLEVEL_ell >= 3, print(" Recherche de points triviaux sur la courbe"));
1090  P *= x;
1091if( DEBUGLEVEL_ell >= 3, print("Y^2 = ",P));
1092  pointstriv = concat( pointstriv, nfratpoint(bnf.nf,P,LIMTRIV,0));
1093if( DEBUGLEVEL_ell >= 1, print("points triviaux sur E'(K) = ");
1094  print(lift(pointstriv)); print());
1095
1096  listpoints = nfellcount(bnf.nf,apinit,bpinit,KS2gen,pointstriv);
1097if( DEBUGLEVEL_ell >= 1, print("points sur E'(K) = ",lift(listpoints[1])));
1098  np1 = listpoints[2][1]; np2 = listpoints[2][2];
1099  listpoints2 = vector(#listpoints[1],i,0);
1100  for( i = 1, #listpoints[1],
1101    listpoints2[i] = [0,0];
1102    aux1 = listpoints[1][i][1]^2;
1103    if( aux1 != 0,
1104      aux2 = listpoints[1][i][2];
1105      listpoints2[i][1] = aux2^2/aux1/4;
1106      listpoints2[i][2] = aux2*(bpinit-aux1)/aux1/8
1107    , listpoints2[i] = listpoints[1][i]));
1108if( DEBUGLEVEL_ell >= 1, print("points sur E(K) = ",lift(listpoints2)); print());
1109  pointgen = concat(pointgen,listpoints2);
1110
1111  certainp = (np1 == np2);
1112if( DEBUGLEVEL_ell >= 1,
1113  if( certainp,
1114    print("[E'(K):phi(E(K))]   = ",1<<np1);
1115    print("#S^(phi)(E/K)       = ",1<<np2);
1116    print("#III(E/K)[phi]      = 1"); print()
1117  ,
1118    print("[E'(K):phi(E(K))]  >= ",1<<np1);
1119    print("#S^(phi)(E/K)       = ",1<<np2);
1120    print("#III(E/K)[phi]     <= ",1<<(np2-np1)); print());
1121
1122  if( !certain && (np2>np1), print1(1<<(np2-np1)," <= "));
1123  print1("#III(E/K)[2]       ");
1124  if( certain && certainp, print1(" "), print1("<"));
1125  print("= ",1<<(n2+np2-n1-np1));
1126
1127  print("#E(K)[2]            = ",1<<tors);
1128);
1129  rang = n1+np1-2;
1130if( DEBUGLEVEL_ell >= 1,
1131  if( certain && certainp,
1132    print("#E(K)/2E(K)         = ",(1<<(rang+tors)));
1133    print("rang                = ",rang); print()
1134  ,
1135    print("#E(K)/2E(K)        >= ",(1<<(rang+tors))); print();
1136    print(rang," <= rang          <= ",n2+np2-2); print()
1137  ));
1138
1139  strange = (n2+np2-n1-np1)%2;
1140  if( strange,
1141if( DEBUGLEVEL_ell >= 1,
1142      print(" !!! III doit etre un carre !!!"); print("donc"));
1143    if( certain,
1144      np1++;
1145      certainp = (np1 == np2);
1146if( DEBUGLEVEL_ell >= 1,
1147        if( certainp,
1148          print("[E'(K):phi(E(K))]   = ",1<<np1);
1149          print("#S^(phi)(E/K)       = ",1<<np2);
1150          print("#III(E/K)[phi]      = 1"); print()
1151        ,
1152          print("[E'(K):phi(E(K))]  >= ",1<<np1);
1153          print("#S^(phi)(E/K)       = ",1<<np2);
1154          print("#III(E/K)[phi]     <= ",1<<(np2-np1)); print())
1155      )
1156    ,
1157    if( certainp,
1158      n1++;
1159      certain = ( n1 == n2);
1160if( DEBUGLEVEL_ell >= 1,
1161        if( certain,
1162          print("[E(K):phi'(E'(K))]   = ",1<<n1);
1163          print("#S^(phi')(E'/K)       = ",1<<n2);
1164          print("#III(E'/K)[phi']      = 1"); print()
1165        ,
1166          print("[E(K):phi'(E'(K))]  >= ",1<<n1);
1167          print("#S^(phi')(E'/K)      = ",1<<n2);
1168          print("#III(E'/K)[phi']    <= ",1<<(n2-n1)); print())
1169      )
1170    , n1++)       
1171  );
1172
1173if( DEBUGLEVEL_ell >= 1,
1174  if( !certain && (np2>np1), print1(1<<(np2-np1)," <= "));
1175  print1("#III(E/K)[2]       ");
1176  if( certain && certainp, print1(" "), print1("<"));
1177  print("= ",1<<(n2+np2-n1-np1));
1178  print("#E(K)[2]            = ",1<<tors);
1179);
1180  rang = n1+np1-2;
1181if( DEBUGLEVEL_ell >= 1,
1182  if( certain && certainp,
1183    print("#E(K)/2E(K)         = ",(1<<(rang+tors))); print();
1184    print("rang                = ",rang); print()
1185  ,
1186    print("#E(K)/2E(K)        >= ",(1<<(rang+tors))); print();
1187    print(rang," <= rang          <= ",n2+np2-2); print())
1188  ));
1189
1190\\ fin de strange
1191
1192if( DEBUGLEVEL_ell >= 1, print("points = ",pointgen));
1193if( DEBUGLEVEL_ell >= 3, print("fin de bnfell2descent_viaisog"));
1194  return([rang,n2+np2-2+tors,pointgen]);
1195}
1196if( DEBUGLEVEL_ell >= 4, print("nfchinremain"));
1197{
1198nfchinremain( nf, b, fact) =
1199\\ Chinese Remainder Theorem
1200local(l,fact2,i);
1201
1202if( DEBUGLEVEL_ell >= 4, print("entree dans nfchinremain"));
1203  l = #fact[,1];
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])));
1207  fact2 = idealaddtoone(nf,fact2);
1208\\ A CHANGER : fact2 = matbasistoalg(nf,fact2);
1209  for( i = 1, l,
1210    fact2[i] = nfbasistoalg(nf,fact2[i]));
1211if( DEBUGLEVEL_ell >= 4, print("fin de nfchinremain"));
1212  return(fact2);
1213}
1214if( DEBUGLEVEL_ell >= 4, print("bnfqfsolve2"));
1215{
1216bnfqfsolve2(bnf, aleg, bleg, auto=[y]) =
1217\\ Solves Legendre Equation x^2-aleg*Y^2=bleg*Z^2
1218\\ Using quadratic norm equations
1219\\ auto contient les automorphismes de bnf sous forme de polynomes
1220\\ en y, avec auto[1]=y .
1221local(aux,solvepolrel,auxsolve,solvepolabs,exprxy,rrrnf,bbbnf,SL0,i,SL1,SL,sunL,fondsunL,normfondsunL,SK,sunK,fondsunK,vecbleg,matnorm,matnormmod,expsolution,solution,reste,carre,verif);
1222
1223if( DEBUGLEVEL_ell >= 3, print("entree dans bnfqfsolve2"));
1224  solvepolrel = x^2-aleg;
1225if( DEBUGLEVEL_ell >= 4, print("aleg = ",aleg));
1226if( DEBUGLEVEL_ell >= 4, print("bleg = ",bleg));
1227
1228  if( #auto > 1,
1229if( DEBUGLEVEL_ell >= 4, print("factorisation du discriminant avec les automorhpismes de bnf"));
1230    for( i = 2, #auto,
1231      aux = abs(polresultant(lift(aleg)-subst(lift(aleg),y,auto[i]),bnf.pol));
1232      if( aux, addprimes(factor(aux)[,1]))));
1233
1234  auxsolve = rnfequation(bnf,solvepolrel,1);
1235  solvepolabs = auxsolve[1];
1236  exprxy = auxsolve[2];
1237  if( auxsolve[3],
1238if( DEBUGLEVEL_ell >=5, print(" CECI EST LE NOUVEAU CAS auxsolve[3] != 0")));
1239if( DEBUGLEVEL_ell >= 4, print(" bbbnfinit ",solvepolabs));
1240  rrrnf = rnfinit(bnf,solvepolrel);
1241  bbbnf = bnfinit(solvepolabs,1);
1242if( DEBUGLEVEL_ell >= 4, print(" done"));
1243  SL0 = 1;
1244if( DEBUGLEVEL_ell >= 4, print("bbbnf.clgp = ",bbbnf.clgp));
1245  for( i = 1, #bbbnf.clgp[2],
1246    if( bbbnf.clgp[2][i]%2 == 0,
1247      SL0 = idealmul(bbbnf,SL0,bbbnf.clgp[3][i][1,1])));
1248  SL1 = idealmul(bbbnf,SL0,rnfeltup(rrrnf,bleg));
1249  SL = idealfactor(bbbnf,SL1)[,1]~;
1250  sunL = bnfsunit(bbbnf,SL);
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])));
1253  normfondsunL = norm(rnfeltabstorel( rrrnf,fondsunL));
1254  SK = idealfactor(bnf,idealnorm(bbbnf,SL1))[,1]~;
1255  sunK = bnfsunit(bnf,SK);
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])));
1258  vecbleg = bnfissunit(bnf,sunK,bleg);
1259  matnorm = matrix(#fondsunK,#normfondsunL,i,j,0);
1260  for( i = 1, #normfondsunL,
1261    matnorm[,i] = lift(bnfissunit( bnf,sunK,normfondsunL[i] )));
1262  matnormmod = matnorm*Mod(1,2);
1263  expsolution = lift(matinverseimage( matnormmod, vecbleg*Mod(1,2)));
1264if( expsolution == []~, error(" bnfqfsolve2 : IL N'Y A PAS DE SOLUTION "));
1265  solution = prod( i = 1, #expsolution, fondsunL[i]^expsolution[i]);
1266  solution = rnfeltabstorel(rrrnf,solution);
1267  reste = (lift(vecbleg) - matnorm*expsolution)/2;
1268  carre = prod( i = 1, #vecbleg, fondsunK[i]^reste[i]);
1269  solution *= carre;
1270  x1=polcoeff(lift(solution),1,x);x0=polcoeff(lift(solution),0,x);
1271  verif = x0^2 - aleg*x1^2-bleg;
1272if( verif, error(" bnfqfsolve2 : MAUVAIS POINT"));
1273if( DEBUGLEVEL_ell >= 3, print("fin de bnfqfsolve2"));
1274  return([x0,x1,1]);
1275}
1276if( DEBUGLEVEL_ell >= 4, print("bnfqfsolve"));
1277{
1278bnfqfsolve(bnf, aleg, bleg, flag3, auto=[y]) =
1279\\ cette fonction resout l'equation X^2-aleg*Y^2=bleg*Z^2
1280\\ dans le corps de nombres nf.
1281\\ la solution est [X,Y,Z],
1282\\ [0,0,0] sinon.
1283local(nf,aa,bb,na,nb,maxna,maxnb,mat,resl,t,sq,pol,vecrat,alpha,xx,yy,borne,test,sun,fact,suni,k,f,l,aux,alpha2,maxnbiter,idbb,rem,nbiter,mask,oldnb,newnb,bor,testici,de,xxp,yyp,rap,verif);
1284
1285if( DEBUGLEVEL_ell >=5, print("entree dans bnfqfsolve"));
1286if( DEBUGLEVEL_ell >= 3, print("(a,b) = (",aleg,",",bleg,")"));
1287  nf = bnf.nf;
1288  aleg = Mod(lift(aleg),nf.pol); aa = aleg;
1289  bleg = Mod(lift(bleg),nf.pol); bb = bleg;
1290
1291  if( aa == 0,
1292if( DEBUGLEVEL_ell >= 5, print("fin de bnfqfsolve"));
1293    return([0,1,0]~));
1294  if( bb == 0,
1295if( DEBUGLEVEL_ell >= 5, print("fin de bnfqfsolve"));
1296    return([0,0,1]~));
1297
1298  na = abs(norm(aa)); nb = abs(norm(bb));
1299  if( na > nb, maxnb = na, maxnb = nb);
1300  maxnb <<= 20;
1301  mat = Mod(matid(3),nf.pol); borne = 1;
1302  test = 0; nbiter = 0;
1303
1304  while( 1,
1305    if( flag3 && bnf.clgp[1]>1, resl = bnfqfsolve2(bnf,aa,bb,auto)~; break);
1306if( DEBUGLEVEL_ell >= 4, print("(na,nb,a,b) = ",lift([na,nb,aa,bb,norm(aa),norm(bb)])));
1307if( DEBUGLEVEL_ell >= 5, print("***",nb,"*** "));
1308    if( nb >= maxnb,
1309      mat = Mod(matid(3),nf.pol);
1310      aa = aleg; bb = bleg; na = abs(norm(aleg)); nb = abs(norm(bleg)));
1311    if( aa == 1, resl = [1,1,0]~; break);
1312    if( bb == 1, resl = [1,0,1]~; break);
1313    if( aa+bb == 1, resl = [1,1,1]~; break);
1314    if( aa+bb == 0, resl = [0,1,1]~; break);
1315    if( aa == bb  && aa != 1,
1316      t = aa*mat[,1];
1317      mat[,1] = mat[,3]; mat[,3] = t;
1318      aa = -1; na = 1);
1319    if( issquare(na),
1320      sq = nfsqrt(nf,aa);
1321      if( sq != [], resl = [sq[1],1,0]~; break));
1322    if( issquare(nb),
1323      sq = nfsqrt(nf,bb);
1324      if( sq != [], resl = [sq[1],0,1]~; break));
1325    if( na > nb,
1326      t = aa; aa = bb; bb = t;
1327      t = na; na = nb; nb = t;
1328      t = mat[,3]; mat[,3] = mat[,2]; mat[,2] = t);
1329    if( nb == 1,
1330if( DEBUGLEVEL_ell >= 4, print("(a,b) = ",lift([aa,bb])));
1331if( DEBUGLEVEL_ell >= 4, print("(na,nb) = ",lift([na,nb])));
1332      if( aleg == aa && bleg == bb, mat = Mod(matid(3),nf.pol));
1333      if( flag3, resl = bnfqfsolve2(bnf,aa,bb,auto)~; break);
1334      pol = aa*x^2+bb; 
1335      vecrat = nfratpoint(nf,pol,borne++,1);
1336      if( vecrat != 0, resl=[vecrat[2],vecrat[1],vecrat[3]]~; break);
1337
1338      alpha = 0;
1339if( DEBUGLEVEL_ell >= 4, print("borne = ",borne));
1340      while( alpha==0,
1341        xx = nfrandint(nf,borne); yy = nfrandint(nf,borne);
1342        borne++;
1343        alpha = xx^2-aa*yy^2 );
1344      bb *= alpha; nb *= abs(norm(alpha));
1345      t = xx*mat[,1]+yy*mat[,2];
1346      mat[,2] = xx*mat[,2]+aa*yy*mat[,1];
1347      mat[,1] = t;
1348      mat[,3] *= alpha
1349    ,
1350      test = 1;
1351if( DEBUGLEVEL_ell >= 4, print("on factorise bb = ",bb));
1352      sun = bnfsunit(bnf,idealfactor(bnf,bb)[,1]~);
1353      fact = lift(bnfissunit(bnf,sun,bb));
1354if( DEBUGLEVEL_ell >= 4, print("fact = ",fact));
1355      suni = concat(bnf.futu,vector(#sun[1],i,nfbasistoalg(bnf,sun[1][i])));
1356      for( i = 1, #suni,
1357        if( (f = fact[i]>>1),
1358          test =0;
1359          for( k = 1, 3, mat[k,3] /= suni[i]^f);
1360          nb /= abs(norm(suni[i]))^(2*f);
1361          bb /= suni[i]^(2*f)));
1362if( DEBUGLEVEL_ell >= 4, print("on factorise bb = ",bb));
1363      fact = idealfactor(nf,bb);
1364if( DEBUGLEVEL_ell >= 4, print("fact = ",fact));
1365      l = #fact[,1];
1366
1367      if( test,
1368        aux = 1;
1369        for( i = 1, l,
1370          if( (f = fact[i,2]>>1) &&
1371               !(fact[i,1][1]%2) && !nfpsquareodd(nf,aa,fact[i,1]),
1372            aux=idealmul(nf,aux,idealpow(nf,fact[i,1],f))));
1373        if( aux != 1,
1374          test = 0;
1375          alpha = nfbasistoalg(nf,idealappr(nf,idealinv(nf,aux)));
1376          alpha2 = alpha^2;
1377          bb *= alpha2; nb *= abs(norm(alpha2));
1378          mat[,3] *= alpha));
1379      if( test,
1380        maxnbiter = 1<<l;
1381        sq = vector(l,i,nfissquarep(nf,aa,fact[i,1],fact[i,2]));
1382        l = #sq;
1383if( DEBUGLEVEL_ell >= 4, print("sq = ",sq); print("fact = ",fact); print("l = ",l));
1384        if( l > 1,
1385          idbb = idealhnf(nf,bb);
1386          rem = nfchinremain(nf,idbb,fact));
1387        test = 1; nbiter = 1;
1388        while( test && nbiter <= maxnbiter,
1389          if( l > 1,
1390            mask = nbiter; xx = 0;
1391            for( i = 1, l,
1392              if( mask%2, xx += rem[i]*sq[i], xx -= rem[i]*sq[i] ); mask >>= 1)
1393          ,
1394            test = 0; xx = sq[1]);
1395          xx = mynfeltmod(nf,xx,bb);
1396          alpha = xx^2-aa;
1397          if( alpha == 0, resl=[xx,1,0]~; break(2));
1398          t = alpha/bb;
1399if( DEBUGLEVEL_ell >= 4, print("[alpha,bb] = ",[alpha,bb]));
1400          oldnb = nb;
1401          newnb = abs(norm(t));
1402if( DEBUGLEVEL_ell >= 4, print("[oldnb,newnb,oldnb/newnb] = ",[oldnb,newnb,oldnb/newnb+0.]));
1403          while( nb > newnb,
1404            mat[,3] *= t;
1405            bb = t; nb = newnb;
1406            t = xx*mat[,1]+mat[,2];
1407            mat[,2] = aa*mat[,1] + xx*mat[,2];
1408            mat[,1] = t;
1409            xx = mynfeltmod(nf,-xx,bb);
1410            alpha = xx^2-aa;
1411            t = alpha/bb;
1412            newnb = abs(norm(t));
1413          );
1414        if( nb == oldnb, nbiter++, test = 0);
1415        );
1416        if( nb == oldnb,
1417          if( flag3, resl = bnfqfsolve2(bnf,aa,bb,auto)~; break);
1418          pol = aa*x^2+bb;
1419          vecrat =nfratpoint(nf,pol,borne++<<1,1);
1420          if( vecrat != 0, resl = [vecrat[2],vecrat[1],vecrat[3]]~; break);
1421
1422          bor = 1000; yy = 1; testici = 1;
1423          for( i = 1, 10000, de = nfbasistoalg(nf,vectorv(poldegree(nf.pol),j,random(bor)));
1424            if( idealadd(bnf,de,bb) != matid(poldegree(bnf.pol)),next);
1425            xxp = mynfeltmod(bnf,de*xx,bb); yyp = mynfeltmod(bnf,de*yy,bb);
1426            rap = (norm(xxp^2-aa*yyp^2)/nb^2+0.);
1427            if( abs(rap) < 1,
1428if( DEBUGLEVEL_ell >= 4, print("********** \n \n MIRACLE ",rap," \n \n ***"));
1429              t = (xxp^2-aa*yyp^2)/bb;
1430              mat[,3] *= t;
1431              bb = t; nb = abs(norm(bb));
1432if( DEBUGLEVEL_ell >= 4, print("newnb = ",nb));
1433              t = xxp*mat[,1]+yyp*mat[,2];
1434              mat[,2] = aa*yyp*mat[,1] + xxp*mat[,2];
1435              mat[,1] = t;
1436              xx = xxp; yy = -yyp; testici = 0;
1437              ));
1438
1439          if( testici, 
1440            alpha = 0;
1441            while( alpha == 0,
1442              xx = nfrandint(nf,4*borne); yy = nfrandint(nf,4*borne);
1443              borne++;
1444              alpha = xx^2-aa*yy^2);
1445            bb *= alpha; nb *= abs(norm(alpha));
1446            t = xx*mat[,1] + yy*mat[,2];
1447            mat[,2] = xx*mat[,2]+aa*yy*mat[,1];
1448            mat[,1] = t;
1449            mat[,3] *= alpha;)))));
1450  resl = lift(mat*resl);
1451if( DEBUGLEVEL_ell >= 5, print("resl1 = ",resl));
1452if( DEBUGLEVEL_ell >= 5, print("content = ",content(resl)));
1453  resl /= content(resl);
1454  resl = Mod(lift(resl),nf.pol);
1455if( DEBUGLEVEL_ell >=5, print("resl3 = ",resl));
1456  fact = idealadd(nf,idealadd(nf,resl[1],resl[2]),resl[3]);
1457  fact = bnfisprincipal(bnf,fact,3);
1458  resl *=1/nfbasistoalg(nf,fact[2]);
1459if( DEBUGLEVEL_ell >= 5, print("resl4 = ",resl));
1460if( DEBUGLEVEL_ell >= 3, print("resl = ",resl));
1461  verif = (resl[1]^2-aleg*resl[2]^2-bleg*resl[3]^2 == 0);
1462  if( !verif && DEBUGLEVEL_ell >= 0, error(" bnfqfsolve : MAUVAIS POINT"));
1463if( DEBUGLEVEL_ell >= 3, print("fin de bnfqfsolve"));
1464  return(resl);
1465}
1466if( DEBUGLEVEL_ell >= 4, print("bnfredquartique2"));
1467{
1468bnfredquartique2( bnf, pol, r,a,b) =
1469\\ reduction d'une quartique issue de la 2-descente
1470\\ en connaissant les valeurs de r, a et b.
1471local(gcc,princ,den,ap,rp,pol2);
1472
1473if( DEBUGLEVEL_ell >= 4, print("entree dans bnfredquartique2"));
1474if( DEBUGLEVEL_ell >= 4, print([r,a,b]));
1475if( DEBUGLEVEL_ell >= 3, print(" reduction de la quartique ",pol));
1476
1477  if( a == 0,
1478    rp = 0
1479  ,
1480    gcc = idealadd(bnf,b,a);
1481    if( gcc == 1,
1482      rp = nfbasistoalg(bnf,idealaddtoone(bnf.nf,a,b)[1])/a;
1483      rp = mynfeltmod(bnf,r*rp,b)
1484    ,
1485      princ = bnfisprincipal(bnf,gcc,3);
1486      if( princ[1] == 0, gcc = nfbasistoalg(bnf,princ[2])
1487      ,
1488if( DEBUGLEVEL_ell >= 3, print(" quartique non reduite"));
1489if( DEBUGLEVEL_ell >= 4, print("fin de bnfredquartique2"));
1490        return([pol,0,1]));
1491      rp = nfbasistoalg(bnf,idealaddtoone(bnf.nf,a/gcc,b/gcc)[1])/(a/gcc);
1492      rp = mynfeltmod(bnf,r*rp,b)/gcc;
1493      b /= gcc;
1494    )
1495  );
1496  pol2 = subst(pol/b,x,rp+b*x)/b^3;
1497if( DEBUGLEVEL_ell >= 3, print(" quartique reduite = ",pol2));
1498if( DEBUGLEVEL_ell >= 4, print("fin de bnfredquartique2"));
1499  return([pol2,rp,b]);
1500}
1501if( DEBUGLEVEL_ell >= 4, print("bnfell2descent_gen"));
1502{
1503bnfell2descent_gen( bnf, ell, ext, help=[], bigflag=1, flag3=1, auto=[y]) =
1504\\ bnf a un polynome en y.
1505\\ si ell= y^2=P(x), alors ext est
1506\\ ext[1] est une equation relative du corps (=P(x)),
1507\\ ext[2] est le resultat rnfequation(bnf,P,2);
1508\\ ext[3] est le buchinitfu (sur Q) de l'extension.
1509\\ dans la suite ext est note L = K(theta).
1510\\ help est une liste de points deja connus sur ell.
1511\\ si bigflag !=0 alors on applique bnfredquartique.
1512\\ si flag3 ==1 alors on utilise bnfqfsolve2 (equation aux normes) pour resoudre Legendre
1513\\ auto est une liste d'automorphismes connus de bnf
1514\\ (ca peut aider a factoriser certains discriminiants).
1515\\ ell est de la forme y^2=x^3+A*x^2+B*x+C
1516\\ ie ell=[0,A,0,B,C], avec A,B et C entiers.
1517\\
1518local(nf,unnf,ellnf,A,B,C,S,plist,Lrnf,SLprod,oddclass,SLlist,LS2gen,polrel,alpha,ttheta,KS2gen,LS2genunit,normLS2,normcoord,LS2coordtilda,LS2tilda,i,aux,j,listgen,listpoints,listpointstriv,listpointsmwr,list,m1,m2,loc,lastloc,maskwhile,iwhile,zc,iaux,liftzc,ispointtriv,point,c,b,a,sol,found,alphac,r,denc,dena,cp,alphacp,beta,mattr,vec,z1,ff,cont,d,e,polorig,pol,redq,transl,multip,UVW,pointxx,point2,v,rang);
1519
1520if( DEBUGLEVEL_ell >= 4, print("entree dans bnfell2descent_gen"));
1521
1522\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
1523\\      construction de L(S,2)         \\
1524\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
1525
1526  nf = bnf.nf;
1527  unnf = Mod(1,nf.pol);
1528  ellnf = ell*unnf;
1529  if( #ellnf <= 5, ellnf = ellinit(ellnf,1));
1530  A = ellnf.a2; if( DEBUGLEVEL_ell >= 2, print("A = ",A));
1531  B = ellnf.a4; if( DEBUGLEVEL_ell >= 2, print("B = ",B));
1532  C = ellnf.a6; if( DEBUGLEVEL_ell >= 2, print("C = ",C));
1533
1534  S = 6*lift(ellnf.disc);
1535  plist = idealfactor(nf,S)[,1];
1536  Lrnf = ext[3];
1537  SLprod = subst(lift(ext[1]'),y,lift(ext[2][2]));
1538if( DEBUGLEVEL_ell >= 3, print(ext[2]));
1539  oddclass = 0;
1540  while( !oddclass,
1541\\ Constructoin de S:
1542    SLlist = idealfactor(Lrnf,SLprod)[,1]~;
1543\\ Construction des S-unites
1544    LS2gen = bnfsunit(Lrnf,SLlist);
1545if( DEBUGLEVEL_ell >= 4, print("LS2gen = ",LS2gen));
1546\\ on ajoute la partie paire du groupe de classes.
1547    oddclass = LS2gen[5][1]%2;
1548    if( !oddclass,
1549if( DEBUGLEVEL_ell >= 3, print("2-class group ",LS2gen[5][3][1][1,1]));
1550      S *= LS2gen[5][3][1][1,1];
1551      SLprod = idealmul(Lrnf,SLprod,(LS2gen[5][3][1])));
1552  );
1553 
1554  polrel = ext[1];
1555  alpha = Mod(Mod(y,nf.pol),polrel); \\ alpha est l'element primitif de K
1556  ttheta = Mod(x,polrel);            \\ ttheta est la racine de P(x)
1557
1558  KS2gen = bnfsunit(bnf,idealfactor(nf,S)[,1]~);
1559 
1560if( DEBUGLEVEL_ell >= 3, print("#KS2gen = ",#KS2gen[1]));
1561if( DEBUGLEVEL_ell >= 3, print("KS2gen = ",KS2gen[1]));
1562
1563  LS2genunit = lift(Lrnf.futu);
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
1567
1568  LS2genunit = subst(LS2genunit,x,ttheta);
1569  LS2genunit = LS2genunit*Mod(1,polrel);
1570if( DEBUGLEVEL_ell >= 3, print("#LS2genunit = ",#LS2genunit));
1571if( DEBUGLEVEL_ell >= 3, print("LS2genunit = ",LS2genunit));
1572
1573\\ dans LS2gen, on ne garde que ceux dont la norme est un carre.
1574
1575  normLS2gen = norm(LS2genunit);
1576if( DEBUGLEVEL_ell >= 4, print("normLS2gen = ",normLS2gen));
1577
1578\\ matrice de l'application norme
1579
1580  normcoord = matrix(#KS2gen[1]+#bnf[8][5]+1,#normLS2gen,i,j,0);
1581  for( i = 1, #normLS2gen,
1582    normcoord[,i] = bnfissunit(bnf,KS2gen,normLS2gen[i]));
1583if( DEBUGLEVEL_ell >= 4, print("normcoord = ",normcoord));
1584
1585\\ construction du noyau de la norme
1586
1587  LS2coordtilda = lift(matker(normcoord*Mod(1,2)));
1588if( DEBUGLEVEL_ell >= 4, print("LS2coordtilda = ",LS2coordtilda));
1589  LS2tilda = vector(#LS2coordtilda[1,],i,0);
1590  for( i = 1, #LS2coordtilda[1,],
1591    aux = 1;
1592    for( j = 1, #LS2coordtilda[,i],
1593      if( sign(LS2coordtilda[j,i]),
1594        aux *= LS2genunit[j]));
1595    LS2tilda[i] = aux;
1596  );
1597
1598if( DEBUGLEVEL_ell >= 3, print("LS2tilda = ",LS2tilda));
1599if( DEBUGLEVEL_ell >= 3, print("norm(LS2tilda) = ",norm(LS2tilda)));
1600
1601\\ Fin de la construction de L(S,2)
1602
1603  listgen = LS2tilda;
1604if( DEBUGLEVEL_ell >= 2, print("LS2gen = ",listgen));
1605if( DEBUGLEVEL_ell >= 2, print("#LS2gen = ",#listgen));
1606  listpoints = [];
1607
1608if( DEBUGLEVEL_ell >= 3, print("(A,B,C) = ",[A,B,C]));
1609
1610\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
1611\\   Recherche de points triviaux   \\
1612\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
1613
1614if( DEBUGLEVEL_ell >= 2, print(" Recherche de points triviaux sur la courbe "));
1615  listpointstriv = nfratpoint(nf,x^3+A*x^2+B*x+C,LIMTRIV,0);
1616  listpointstriv = concat(help,listpointstriv);
1617if( DEBUGLEVEL_ell >= 1, print("points triviaux sur la courbe = ",listpointstriv));
1618
1619\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
1620\\          parcours de L(S,2)         \\
1621\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
1622
1623  listpointsmwr = [];
1624  list = [ 6, ellnf.disc, 0];
1625  m1 = 0; m2 = 0; lastloc = -1;
1626  maskwhile = 1<<#listgen;
1627  listELS = [0]; listnotELS = [];
1628  iwhile = 1;
1629  while( iwhile < maskwhile,
1630if( DEBUGLEVEL_ell >= 4, print("iwhile = ",iwhile); print("listgen = ",listgen));
1631
1632\\ utilise la structure de groupe pour detecter une eventuelle solubilite locale.
1633if( DEBUGLEVEL_ell >= 4, print("listELS = ",listELS);
1634                   print("listnotELS = ",listnotELS));
1635    sol = 1; loc = 0;
1636    for( i = 1, #listELS,
1637      for( j = 1, #listnotELS,
1638        if( bitxor(listELS[i],listnotELS[j]) == iwhile,
1639if( DEBUGLEVEL_ell >= 3, print(" Non ELS d'apres la structure de groupe"));
1640          listnotELS = concat(listnotELS,[iwhile]);
1641          sol = 0; break(2))));
1642    if( !sol, iwhile++; next);
1643
1644    for( i = 1, #listELS,
1645      for( j = i+1, #listELS,
1646        if( bitxor(listELS[i],listELS[j]) == iwhile,
1647if( DEBUGLEVEL_ell >= 3, print(" ELS d'aptres la structure de "));
1648          listELS = concat(listELS,[iwhile]);
1649          loc = 2;
1650          break(2))));
1651
1652    iaux = vector(#listgen,i,bittest(iwhile,i-1))~;
1653    iaux = (LS2coordtilda*iaux)%2;
1654
1655    zc = unnf*prod( i = 1, #LS2genunit,LS2genunit[i]^iaux[i]);
1656
1657if( DEBUGLEVEL_ell >= 2, print("zc = ",zc));
1658    liftzc = lift(zc);
1659
1660\\ Est-ce un point trivial ?
1661    ispointtriv = 0;
1662    for( i = 1, #listpointstriv,
1663      point = listpointstriv[i];
1664      if( #point == 2 || point[3] != 0,
1665        if( nfissquare(Lrnf.nf,subst((lift(point[1])-x)*lift(liftzc),y,lift(ext[2][2]))),
1666if( DEBUGLEVEL_ell >= 2, print(" vient du point trivial ",point));
1667          listpointsmwr = concat(listpointsmwr,[point]);
1668          m1++;
1669          listELS = concat(listELS,[iwhile]);
1670          if( degre(iwhile) > lastloc, m2++);
1671          sol = found = ispointtriv = 1; break
1672          )));
1673
1674\\ Ce n'est pas un point trivial
1675    if( !ispointtriv,
1676      c = polcoeff(liftzc,2);
1677      b = -polcoeff(liftzc,1);
1678      a = polcoeff(liftzc,0);
1679      sol = 0; found = 0;
1680\\ \\\\\\\\\\\\\
1681\\ On cherche a ecrire zc sous la forme a-b*theta
1682\\ \\\\\\\\\\\\\
1683      if( c == 0, sol = 1,
1684        alphac = (A*b+B*c-a)*c+b^2;
1685if( DEBUGLEVEL_ell >= 3, print("alphac = ",alphac));
1686        r = nfsqrt(nf,norm(zc))[1];
1687        if( alphac == 0,
1688\\ cas particulier
1689if( DEBUGLEVEL_ell >= 3, print(" on continue avec 1/zc"));
1690          sol = 1; zc = norm(zc)*(1/zc);
1691if( DEBUGLEVEL_ell >= 2, print(" zc = ",zc))
1692        ,
1693\\ Il faut resoudre une forme quadratique
1694\\ Existence (locale = globale) d'une solution :
1695          denc = deno(lift(c));
1696          if( denc != 1, cp = c*denc^2, cp = c);
1697          dena = deno(lift(alphac));
1698          if( dena != 1, alphacp = alphac*dena^2, alphacp = alphac);
1699if( DEBUGLEVEL_ell >= 2, print1(" symbole de Hilbert (",alphacp,",",cp,") = "));
1700          sol = loc || (mynfhilbert(nf, alphacp,cp)+1);
1701if( DEBUGLEVEL_ell >= 2, print(sol-1));
1702          if( sol,
1703            beta = A*(A*b*c+B*c^2+b^2)-C*c^2+a*b;
1704            mattr = matid(3);
1705            mattr[1,1] = c ;mattr[2,2] = alphac ;
1706            mattr[3,3] = r ;mattr[2,3] = -beta;
1707            mattr[1,2] = -(b +A*c) ;mattr[1,3] = a-B*c+A*(A*c+b);
1708if( DEBUGLEVEL_ell >= 2, print1(" sol de Legendre = "));
1709            vec = bnfqfsolve(bnf,alphacp,cp,flag3,auto);
1710if( DEBUGLEVEL_ell >= 2, print(lift(vec)));
1711            aux = vec[2]*dena;
1712            vec[2] = vec[1];vec[1] = aux;
1713            vec[3] = vec[3]*denc;
1714            vec = (mattr^(-1))*vec;
1715            vec /= content(lift(vec));
1716            z1 = (vec[3]*ttheta+vec[2])*ttheta+vec[1];         
1717if( DEBUGLEVEL_ell >= 3, print(" z1 = ",z1));
1718            zc *= z1^2;
1719if( DEBUGLEVEL_ell >= 2, print(" zc*z1^2 = ",zc));
1720          )
1721        )
1722      )
1723    );
1724   
1725    if( !sol,
1726      listnotELS = concat(listnotELS,[iwhile]);
1727      iwhile++;
1728      next
1729    );
1730
1731\\ \\\\\\\\\\
1732\\ Maintenant zc est de la forme a-b*theta
1733\\ \\\\\\\\\\
1734    if( !ispointtriv,
1735      liftzc = lift(zc);
1736if( DEBUGLEVEL_ell >= 3, print(" zc = ",liftzc));
1737if( DEBUGLEVEL_ell >= 4, print(" N(zc) = ",norm(zc)));
1738      if( poldegree(liftzc) >= 2, error(" bnfell2descent_gen : c <> 0"));
1739      b = -polcoeff(liftzc,1);
1740      a = polcoeff(liftzc,0);
1741if( DEBUGLEVEL_ell >= 4, print(" on factorise (a,b) = ",idealadd(nf,a,b)));
1742      ff = idealfactor(nf,idealadd(nf,a,b));
1743if( DEBUGLEVEL_ell >= 4, print(" ff = ",ff));
1744      cont = 1;
1745      for( i = 1, #ff[,1],
1746        cont = idealmul(nf,cont,idealpow(nf,ff[i,1],ff[i,2]\2)));
1747      cont = idealinv(nf,cont);
1748      cont = nfbasistoalg(nf,bnfisprincipal(bnf,cont,3)[2])^2;
1749      a *= cont; b *= cont; zc *= cont;
1750if( DEBUGLEVEL_ell >= 4, print(" [a,b] = ",[a,b]));
1751      if( nfissquare(nf,b),
1752if( DEBUGLEVEL_ell >= 3, print(" b est un carre"));
1753        point = [a/b,nfsqrt(nf,(a/b)^3+A*(a/b)^2+B*(a/b)+C)[1]];
1754if( DEBUGLEVEL_ell >= 2, print("point trouve = ",point));
1755        listpointsmwr = concat(listpointsmwr,[point]);
1756        m1++;
1757        if( degre(iwhile) > lastloc, m2++);
1758        found = ispointtriv = 1
1759      )
1760    );
1761
1762\\ \\\\\\\\\\\
1763\\ Construction de la quartique
1764\\ \\\\\\\\\\\
1765    if( !ispointtriv,
1766      r = nfsqrt(nf,norm(zc))[1];
1767if( DEBUGLEVEL_ell >= 4, print(" r = ",r));
1768      c = -2*(A*b+3*a);
1769if( DEBUGLEVEL_ell >= 4, print(" c = ",c));
1770      d = 8*r;
1771if( DEBUGLEVEL_ell >= 4, print(" d = ",d));
1772      e = (A^2*b^2 - 2*A*a*b-4*B*b^2-3*a^2);
1773if( DEBUGLEVEL_ell >= 4, print(" e = ",e));
1774      polorig = b*(x^4+c*x^2+d*x+e)*unnf;
1775if( DEBUGLEVEL_ell >= 2, print(" quartique : (",lift(b),")*Y^2 = ",lift(polorig/b)));
1776      list[3] = b;
1777      pol = polorig;
1778      if( bigflag,
1779        redq = bnfredquartique2(bnf,pol,r,a,b);
1780if( DEBUGLEVEL_ell >= 2, print(" reduite: Y^2 = ",lift(redq[1])));
1781        pol = redq[1]; transl = redq[2]; multip = redq[3]
1782      );
1783      point = nfratpoint(nf,pol,LIM1,1);
1784      if( point != [],
1785if( DEBUGLEVEL_ell >= 2, print("point = ",point));
1786        m1++;
1787        if( bigflag,
1788          point[1] = point[1]*multip+transl;
1789          point[2] = nfsqrt(nf,subst(polorig,x,point[1]/point[3]))[1]);
1790        mattr = matid(3);
1791        mattr[1,1] = -2*b^2; mattr[1,2] = (A*b+a)*b;
1792        mattr[1,3] = a^2+(2*B-A^2)*b^2; mattr[2,2] = -b;
1793        mattr[2,3] = a+A*b; mattr[3,3] =r;
1794        UVW = [point[1]^2,point[3]^2,point[1]*point[3]]~;
1795        vec = (mattr^(-1))*UVW;
1796        z1 = (vec[3]*ttheta+vec[2])*ttheta+vec[1];
1797        zc *= z1^2;
1798        zc /= -polcoeff(lift(zc),1);
1799if( DEBUGLEVEL_ell >= 3, print("zc*z1^2 = ",zc));
1800        pointxx = polcoeff(lift(zc),0);
1801        point2 = [ pointxx, nfsqrt(nf,subst(x^3+A*x^2+B*x+C,x,pointxx))[1]];
1802if( DEBUGLEVEL_ell >= 1, print(" point trouve = ",point2));
1803        listpointsmwr = concat(listpointsmwr,[point2]);
1804        if( degre(iwhile) > lastloc, m2++);
1805        found = 1; lastloc = -1
1806      , 
1807        if( loc
1808             || (!bigflag && nflocallysoluble(nf,pol,r,a,b))
1809             || (bigflag && nflocallysoluble(nf,pol,0,1,1)),
1810
1811          if( !loc, listlistELS=concat(listELS,[iwhile]));
1812          if( degre(iwhile) > lastloc, m2++; lastloc = degre(iwhile));
1813          point = nfratpoint(nf,pol,LIM3,1);
1814          if( point != [],
1815            m1++;
1816            if( bigflag,
1817              point[1] = point[1]*multip+transl;
1818              point[2] = nfsqrt(nf,subst(polorig,x,point[1]/point[3]))[1];
1819            );
1820            mattr = matid(3);
1821            mattr[1,1] = -2*b^2; mattr[1,2] = (A*b+a)*b;
1822            mattr[1,3] = a^2+(2*B-A^2)*b^2; mattr[2,2] = -b;
1823            mattr[2,3] = a+A*b; mattr[3,3] =r;
1824            UVW = [point[1]^2,point[3]^2,point[1]*point[3]]~;
1825            vec = (mattr^(-1))*UVW;
1826            z1 = (vec[3]*ttheta+vec[2])*ttheta+vec[1];
1827            zc *= z1^2;
1828            zc /= -polcoeff(lift(zc),1);
1829if( DEBUGLEVEL_ell >= 3, print(" zc*z1^2 = ",zc));
1830            pointxx = polcoeff(lift(zc),0);
1831            point2 = [ pointxx, nfsqrt(nf,subst(x^3+A*x^2+B*x+C,x,pointxx))[1]];
1832if( DEBUGLEVEL_ell >= 1, print(" point trouve = ",point2));
1833            listpointsmwr = concat(listpointsmwr,[point2]);
1834            found = 1; lastloc = -1;
1835          )
1836        , listnotELS = concat(listnotELS,[iwhile])
1837        )
1838      )
1839    );
1840    if( found,
1841      found = 0; lastloc = -1;
1842      v = degre(iwhile);
1843      iwhile = 1<<v;
1844      maskwhile >>= 1;
1845      LS2coordtilda = vecextract(LS2coordtilda,1<<#listgen-iwhile-1);
1846      listgen = vecextract(listgen,1<<#listgen-iwhile-1);
1847      while( listELS[#listELS] >= iwhile,
1848        listELS = vecextract(listELS,1<<(#listELS-1)-1));
1849      while( #listnotELS && listnotELS[#listnotELS] >= iwhile,
1850        listnotELS = vecextract(listnotELS,1<<(#listnotELS-1)-1))
1851    , iwhile ++
1852    )
1853  );
1854
1855if( DEBUGLEVEL_ell >= 2, 
1856  print("m1 = ",m1);
1857  print("m2 = ",m2));
1858if( DEBUGLEVEL_ell >= 1,
1859  print("#S(E/K)[2]    = ",1<<m2));
1860  if( m1 == m2,
1861if( DEBUGLEVEL_ell >= 1,
1862    print("#E(K)/2E(K)   = ",1<<m1);
1863    print("#III(E/K)[2]  = 1");
1864    print("rang(E/K)     = ",m1));
1865    rang = m1
1866  ,
1867if( DEBUGLEVEL_ell >= 1,
1868    print("#E(K)/2E(K)  >= ",1<<m1);
1869    print("#III(E/K)[2] <= ",1<<(m2-m1));
1870    print("rang(E/K)    >= ",m1));
1871    rang = m1;
1872    if( (m2-m1)%2,
1873if( DEBUGLEVEL_ell >= 1,
1874      print(" III devrait etre un carre, donc ");
1875      if( m2-m1 > 1,
1876        print("#E(K)/2E(K)  >= ",1<<(m1+1));
1877        print("#III(E/K)[2] <= ",1<<(m2-m1-1));
1878        print("rang(E/K)    >= ",m1+1)
1879      ,
1880        print("#E(K)/2E(K)  = ",1<<(m1+1));
1881        print("#III(E/K)[2] = 1");
1882        print("rang(E/K)    = ",m1+1)));
1883      rang = m1+1)
1884  );
1885if( DEBUGLEVEL_ell >= 1, print("listpointsmwr = ",listpointsmwr));
1886  for( i = 1, #listpointsmwr,
1887    if( #listpointsmwr[i] == 3,
1888      listpointsmwr[i] = vecextract(listpointsmwr[i],3));
1889    if( !ellisoncurve(ellnf,listpointsmwr[i]),
1890      error("bnfell2descent : MAUVAIS POINT ")));
1891if( DEBUGLEVEL_ell >= 4, print("fin de bnfell2descent_gen"));
1892  return([rang,m2,listpointsmwr]);
1893}
1894if( DEBUGLEVEL_ell >= 4, print("bnfellrank"));
1895{
1896bnfellrank(bnf,ell,help=[],bigflag=1,flag3=1) =
1897\\ Algorithme de la 2-descente sur la courbe elliptique ell.
1898\\ help est une liste de points connus sur ell.
1899
1900\\ attention bnf a un polynome en y.
1901\\ si bigflag !=0, on reduit les quartiques
1902\\ si flag3 != 0, on utilise bnfqfsolve2
1903local(urst,urst1,den,factden,eqtheta,rnfeq,bbnf,ext,rang);
1904
1905if( DEBUGLEVEL_ell >= 3, print("entree dans bnfellrank"));
1906  if( #ell <= 5, ell = ellinit(ell,1));
1907
1908\\ removes the coefficients a1 and a3
1909  urst = [1,0,0,0];
1910  if( ell.a1 != 0 || ell.a3 != 0,
1911    urst1 = [1,0,-ell.a1/2,-ell.a3/2];
1912    ell = ellchangecurve(ell,urst1);
1913    urst = ellcomposeurst(urst,urst1)
1914  );
1915
1916\\ removes denominators
1917  while( (den = idealinv(bnf,idealadd(bnf,idealadd(bnf,1,ell.a2),idealadd(bnf,ell.a4,ell.a6))))[1,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];
1923    urst1 = [1/den,0,0,0];
1924    ell = ellchangecurve(ell,urst1);
1925    urst = ellcomposeurst(urst,urst1);
1926  );
1927
1928  help = ellchangepoint(help,urst);
1929
1930\\ choix de l'algorithme suivant la 2-torsion
1931  ell *= Mod(1,bnf.pol);
1932  eqtheta = Pol([1,ell.a2,ell.a4,ell.a6]);
1933if( DEBUGLEVEL_ell >= 1, print("courbe elliptique : Y^2 = ",eqtheta));
1934  f = nfpolratroots(bnf,eqtheta);
1935
1936  if( #f == 0,                                  \\ cas 1: 2-torsion triviale
1937    rnfeq = rnfequation(bnf,eqtheta,1);
1938    urst1 = [1,-rnfeq[3]*Mod(y,bnf.pol),0,0];
1939    if( rnfeq[3] != 0,
1940      ell = ellchangecurve(ell,urst1);
1941      urst = ellcomposeurst(urst,urst1);
1942      eqtheta = subst(eqtheta,x,x-rnfeq[3]*Mod(y,bnf.pol));
1943      rnfeq = rnfequation(bnf,eqtheta,1);
1944if( DEBUGLEVEL_ell >= 2, print("translation : on travaille avec Y^2 = ",eqtheta));
1945    );
1946if( DEBUGLEVEL_ell >= 3, print1("bbnfinit "));
1947    bbnf = bnfinit(rnfeq[1],1);
1948if( DEBUGLEVEL_ell >= 3, print("done"));
1949    ext = [eqtheta, rnfeq, bbnf];
1950    rang = bnfell2descent_gen(bnf,ell,ext,help,bigflag,flag3)
1951  ,
1952  if( #f == 1,                                  \\ cas 2: 2-torsion = Z/2Z
1953    if( f[1] != 0,
1954      urst1 = [1,f[1],0,0];
1955      ell = ellchangecurve(ell,urst1);
1956      urst = ellcomposeurst(urst,urst1)
1957    );
1958    rang = bnfell2descent_viaisog(bnf,ell)
1959  ,                                             \\ cas 3: 2-torsion = Z/2Z*Z/2Z
1960    rang = bnfell2descent_complete(bnf,f[1],f[2],f[3],flag3)
1961  ));
1962 
1963  rang[3] = ellchangepointinverse(rang[3],urst);
1964if( DEBUGLEVEL_ell >= 3, print("fin de bnfellrank"));
1965
1966return(rang);
1967}
1968if( DEBUGLEVEL_ell >= 4, print("bnfell2descent_complete"));
1969{
1970bnfell2descent_complete(bnf,e1,e2,e3,flag3=1,auto=[y]) =
1971\\ calcul du rang d'une courbe elliptique
1972\\ par la methode de 2-descente complete.
1973\\ Y^2 = (x-e1)*(x-e2)*(x-e3);
1974\\ en suivant la methode decrite par J.Silverman
1975\\ si flag3 ==1 alors on utilise bnfqfsolve2 (equation aux normes)
1976\\ pour resoudre Legendre
1977\\ on pourra alors utiliser auto=nfgaloisconj(bnf.pol)
1978
1979\\ e1, e2 et e3 sont des entiers algebriques de bnf.
1980
1981local(KS2prod,oddclass,KS2gen,i,vect,selmer,rang,X,Y,b1,b2,vec,z1,z2,d31,quart0,quart,cont,fa,point,solx,soly,listepoints);
1982
1983if( DEBUGLEVEL_ell >= 2, print("Algorithme de la 2-descente complete"));
1984
1985\\ calcul de K(S,2)
1986
1987  KS2prod = (e1-e2)*(e2-e3)*(e3-e1)*2;
1988  oddclass = 0;
1989  while( !oddclass,
1990    KS2gen = bnfsunit(bnf,idealfactor(bnf,KS2prod)[,1]~);
1991    oddclass = (KS2gen[5][1]%2);
1992    if( !oddclass,
1993      KS2prod = idealmul(bnf,KS2prod,(KS2gen[5][3][1])));
1994  );
1995  KS2gen = KS2gen[1];
1996\\ A CHANGER : KS2gen = matbasistoalg(bnf,KS2gen);
1997  for( i = 1, #KS2gen,
1998    KS2gen[i] = nfbasistoalg(bnf, KS2gen[i]));
1999  KS2gen = concat(Mod(lift(bnf.tufu),bnf.pol),KS2gen);
2000if( DEBUGLEVEL_ell >= 2,
2001  print("#K(S,2)gen          = ",#KS2gen);
2002  print("K(S,2)gen = ",KS2gen));
2003
2004\\ parcours de K(S,2)*K(S,2)
2005
2006  vect = vector(#KS2gen,i,[0,1]);
2007  selmer = 0;
2008  rang = 0;
2009  listepoints = [];
2010
2011  forvec( X = vect,
2012    b1 = prod( i = 1, #KS2gen, KS2gen[i]^X[i]);
2013    forvec( Y = vect,
2014      b2 = prod( i = 1, #KS2gen, KS2gen[i]^Y[i]);
2015
2016if( DEBUGLEVEL_ell >= 3, print("[b1,b2] = ",lift([b1,b2])));
2017
2018\\ points triviaux provenant de la 2-torsion
2019
2020      if( b1==1 && b2==1,
2021if( DEBUGLEVEL_ell >= 2, print(" point trivial [0]"));
2022        selmer++; rang++; next);
2023      if( nfissquare(bnf.nf,(e2-e1)*b1)
2024        && nfissquare(bnf.nf,(e2-e3)*(e2-e1)*b2),
2025if( DEBUGLEVEL_ell >= 2, print(" point trivial [e2,0]"));
2026        selmer++; rang++; next);
2027      if( nfissquare(bnf.nf,(e1-e2)*b2)
2028        && nfissquare(bnf.nf,(e1-e3)*(e1-e2)*b1),
2029if( DEBUGLEVEL_ell >= 2, print(" point trivial [e1,0]"));
2030        selmer++; rang++; next);
2031      if( nfissquare(bnf.nf,(e3-e1)*b1)
2032        && nfissquare(bnf.nf,(e3-e2)*b2),
2033if( DEBUGLEVEL_ell >= 2, print(" point trivial [e3,0]"));
2034        selmer++; rang++; next);
2035
2036\\ premier critere local : sur les formes quadratiques
2037
2038      if( mynfhilbert(bnf.nf,b1*b2,b1*(e2-e1)) < 0
2039        || mynfhilbert(bnf.nf,b2,b1*(e3-e1)) < 0
2040        || mynfhilbert(bnf.nf,b1,b2*(e3-e2)) < 0
2041        ,
2042if( DEBUGLEVEL_ell >= 3, print("non ELS"));
2043        next);
2044
2045if( DEBUGLEVEL_ell >= 2, print("[b1,b2] = ",lift([b1,b2])));
2046if( DEBUGLEVEL_ell >= 2, print("qf loc soluble"));
2047
2048\\ solution de la premiere forme quadratique
2049
2050      if( b1 != b2,
2051        vec = bnfqfsolve(bnf,b1*b2,b1*(e2-e1),flag3);
2052if( DEBUGLEVEL_ell >= 3, print("sol part = ",vec));
2053        if( vec[3] == 0, error("bnfell2descent_complete : BUG !!! : vec[3]=0 "));
2054        z1 = vec[1]/vec[3]/b1; z2 = vec[2]/vec[3]
2055      ,
2056        z1 = (1+(e2-e1)/b1)/2; z2 = z1-1
2057      );
2058      d31 = e3-e1;
2059      quart0 = b2^2*(z1^2*b1 - d31)*x^4 - 4*z1*b2^2*z2*b1*x^3
2060           + 2*b1*b2*(z1^2*b1 + 2*b2*z2^2 + d31)*x^2 - 4*z1*b2*z2*b1^2*x
2061           + b1^2*(z1^2*b1 - d31);
2062      quart = quart0*b1*b2;
2063if( DEBUGLEVEL_ell >= 4, print("quart = ",quart));
2064      quart *= denominator(simplify(content(quart)))^2;
2065      cont = simplify(content(lift(quart)));
2066      fa = factor(cont);
2067      for( i = 1, #fa[,1], quart /= fa[i,1]^(2*(fa[i,2]\2)));
2068if( DEBUGLEVEL_ell >= 3, print("quart red = ",quart));
2069
2070\\ la quartique est-elle localement soluble ?
2071   
2072      if( !nflocallysoluble(bnf.nf,quart),
2073if( DEBUGLEVEL_ell >= 2, print(" quartique non ELS "));
2074        next);
2075      selmer++;
2076
2077\\ recherche de points sur la quartique.
2078   
2079      point = nfratpoint(bnf.nf,quart,LIM3,1);
2080      if( point != [],
2081if( DEBUGLEVEL_ell >= 2, print("point trouve sur la quartique !!"));
2082if( DEBUGLEVEL_ell >= 3, print(point));
2083        if( point[3],
2084          point /= point[3];
2085          z1 = (2*b2*point[1]*z2-z1*(b1+b2*point[1]^2))/(b1-b2*point[1]^2);
2086          solx = b1*z1^2+e1;
2087          soly = nfsqrt(bnf.nf,(solx-e1)*(solx-e2)*(solx-e3))[1];
2088          listepoints = concat(listepoints,[[solx,soly]]);
2089if( DEBUGLEVEL_ell >= 1, print("point sur la courbe elliptique =",[solx,soly]))
2090        );
2091        rang++
2092      ,
2093if( DEBUGLEVEL_ell >= 2, print("aucun point trouve sur la quartique"))
2094      )
2095    )
2096  );
2097
2098\\ fin
2099
2100if( DEBUGLEVEL_ell >= 1,
2101  print("#S^(2)      = ",selmer));
2102  if( rang > selmer/2, rang = selmer);
2103if( DEBUGLEVEL_ell >= 1,
2104  strange = rang != selmer;
2105  if( strange,
2106  print("#E[K]/2E[K]>= ",rang)
2107, print("#E[K]/2E[K] = ",rang));
2108  print("#E[2]       = 4");
2109);
2110  rang = ceil(log(rang)/log(2))-2;
2111  selmer = valuation(selmer,2);
2112if( DEBUGLEVEL_ell >= 1,
2113  if( strange,
2114  print(selmer-2," >= rang  >= ",rang)
2115, print("rang        = ",rang));
2116  if( rang, print("points = ",listepoints));
2117);
2118
2119return([rang,selmer,listepoints]);
2120}