Ticket #6955: ellQ.gp

File ellQ.gp, 49.8 KB (added by John Cremona, 13 years ago)
Line 
1\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
2\\       Copyright (C) 2008 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/ellQ.gp
21\\
22\\  *********************************************
23\\  *          VERSION 29/04/2008               *
24\\  *********************************************
25\\
26\\ Programme de calcul du rang des courbes elliptiques sur Q.
27\\ langage: GP
28\\ pour l'utiliser, lancer gp, puis taper
29\\ \r ellQ.gp
30\\
31\\ Ce programme utilise le module de resolution des formes quadratiques
32\\ situe a l'adresse
33\\ www.math.unicaen.fr/~simon/qfsolve.gp
34\\ Il faut donc aussi taper:
35\\ \r qfsolve.gp
36\\
37\\ Explications succintes :
38\\ La fonction ellrank() accepte toutes les courbes sous la forme
39\\ [a1,a2,a3,a4,a6]
40\\ Les coefficients peuvent etre entiers ou non.
41\\ L'algorithme utilise est celui de la 2-descente.
42\\ La 2-torsion peut etre quelconque.
43\\ Il suffit de taper :
44\\
45\\ gp > ell = [a1,a2,a3,a4,a6];
46\\ gp > ellrank(ell)
47\\
48\\ Retourne un vecteur [r,s,vec]
49\\ ou r est le rang probable (c'est toujours une minoration du rang),
50\\ s est le 2-rang du groupe de Selmer,
51\\ vec est une liste de points independants dans E(Q)/2E(Q).
52\\
53\\ Courbes de la forme: k*y^2 = x^3+A*x^2+B*x+C
54\\ sans 2-torsion, A,B,C entiers.
55\\ gp > bnf = bnfinit(x^3+A*x^2+B*x+C);
56\\ gp > ell = ellinit([0,A,0,B,C],1);
57\\ gp > rank = ell2descent_gen(ell,bnf,k);
58\\
59\\ Courbes avec #E[2](Q) >= 2:
60\\ ell doit etre sous la forme
61\\ y^2 = x^3 + A*^2 + B*x
62\\ avec A et B entiers.
63\\ gp > ell = [0,A,0,B,0]
64\\ gp > ell2descent_viaisog(ell)
65\\ = algorithme de la 2-descente par isogenies
66\\ Attention A et B doivent etre entiers
67\\
68\\ Courbes avec #E[2](Q) = 4: y^2 = (x-e1)*(x-e2)*(x-e3)
69\\ gp > ell2descent_complete(e1,e2,e3)
70\\ = algorithme de la 2-descente complete
71\\ Attention: les ei doivent etre entiers.
72\\
73\\
74\\ On peut avoir plus ou moins de details de calculs avec
75\\ DEBUGLEVEL_ell = 0;
76\\ DEBUGLEVEL_ell = 1; 2; 3;...
77\\
78\\
79\\
80
81{
82\\
83\\ Variables globales usuelles
84\\
85
86  DEBUGLEVEL_ell = 1; \\ pour avoir plus ou moins de details
87  LIM1 = 5;       \\ limite des points triviaux sur les quartiques
88  LIM3 = 50;      \\ limite des points sur les quartiques ELS
89  LIMTRIV = 50;   \\ limite des points triviaux sur la courbe elliptique
90
91\\
92\\  Variables globales techniques
93\\
94
95  MAXPROB = 20;
96  LIMBIGPRIME = 30; \\ pour distinguer un petit nombre premier d'un grand
97                    \\ utilise un test probabiliste pour les grands
98                    \\ si LIMBIGPRIME = 0, n'utilise aucun test probabiliste
99  ELLREDGENFLAG = 1;\\ pour reduire les genereteurs a la fin de l'algorithme
100}
101
102\\
103\\  Programmes
104\\
105
106\\
107\\ Fonctions communes ell.gp et ellQ.gp
108\\
109{
110ellinverturst(urst) =
111local(u = urst[1], r = urst[2], s = urst[3], t = urst[4]);
112  [1/u,-r/u^2,-s/u,(r*s-t)/u^3];
113}
114{
115ellchangecurveinverse(ell,v) = ellchangecurve(ell,ellinverturst(v));
116}
117{
118ellchangepointinverse(pt,v) = ellchangepoint(pt,ellinverturst(v));
119}
120{
121ellcomposeurst(urst1,urst2) =
122local(u1 = urst1[1], r1 = urst1[2], s1 = urst1[3], t1 = urst1[4],
123      u2 = urst2[1], r2 = urst2[2], s2 = urst2[3], t2 = urst2[4]);
124  [u1*u2,u1^2*r2+r1,u1*s2+s1,u1^3*t2+s1*u1^2*r2+t1];
125}
126{
127ellinverturst(urst) =
128local(u = urst[1], r = urst[2], s = urst[3], t = urst[4]);
129  [1/u,-r/u^2,-s/u,(r*s-t)/u^3];
130}
131{
132ellchangecurveinverse(ell,v) = ellchangecurve(ell,ellinverturst(v));
133}
134{
135ellchangepointinverse(pt,v) = ellchangepoint(pt,ellinverturst(v));
136}
137{
138ellcomposeurst(urst1,urst2) =
139local(u1 = urst1[1], r1 = urst1[2], s1 = urst1[3], t1 = urst1[4],
140      u2 = urst2[1], r2 = urst2[2], s2 = urst2[3], t2 = urst2[4]);
141  [u1*u2,u1^2*r2+r1,u1*s2+s1,u1^3*t2+s1*u1^2*r2+t1];
142}
143{
144polratroots(pol) =
145local(f,ans);
146  f = factor(pol)[,1];
147  ans = [];
148  for( j = 1, #f,
149    if( poldegree(f[j]) == 1,
150      ans = concat(ans,[-polcoeff(f[j],0)/polcoeff(f[j],1)])));
151  return(ans);
152}
153if( DEBUGLEVEL_ell >= 4, print("mysubst"));
154{
155mysubst(polsu,subsx) =
156  if( type(lift(polsu)) == "t_POL",
157    return(simplify(subst(lift(polsu),variable(lift(polsu)),subsx)))
158  , return(simplify(lift(polsu))));
159}
160if( DEBUGLEVEL_ell >= 4, print("degre"));
161{
162degre(idegre) =
163local(ideg,jdeg);
164
165  ideg = idegre; jdeg = 0;
166  while( ideg >>= 1, jdeg++);
167  return(jdeg);
168}
169if( DEBUGLEVEL_ell >= 4, print("nfissquare"));
170{
171nfissquare(nf, a) = #nfsqrt(nf,a) > 0;
172}
173if( DEBUGLEVEL_ell >= 4, print("nfsqrt"));
174{
175nfsqrt( nf, a) =
176\\ si a est un carre, renvoie [sqrt(a)], sinon [].
177local(alift,ta,res,pfact,r1,py);
178
179if( DEBUGLEVEL_ell >= 5, print("entree dans nfsqrt ",a));
180  if( a==0 || a==1,
181if( DEBUGLEVEL_ell >= 5, print("fin de nfsqrt"));
182    return([a]));
183
184  alift = lift(a);
185  ta = type(a);
186  if( !poldegree(alift), alift = polcoeff(alift,0));
187
188  if( type(alift) != "t_POL",
189    if( issquare(alift),
190if( DEBUGLEVEL_ell >= 5, print("fin de nfsqrt"));
191      return([sqrtrat(alift)])));
192
193  if( poldegree(nf.pol) <= 1,
194if( DEBUGLEVEL_ell >= 5, print("fin de nfsqrt"));
195    return([]));
196  if( ta == "t_POL", a = Mod(a,nf.pol));
197
198\\ tous les plgements reels doivent etre >0
199\\   
200  r1 = nf.sign[1];
201  for( i = 1, r1,
202    py = mysubst(alift,nf.roots[i]);
203    if( sign(py) < 0,
204if( DEBUGLEVEL_ell >= 5, print("fin de nfsqrt"));
205      return([])));
206\\ factorisation sur K du polynome X^2-a :
207  if( variable(nf.pol) == x,
208    py = subst(nf.pol,x,y);
209    pfact = lift(factornf(x^2-mysubst(alift,Mod(y,py)),py)[1,1])
210  ,
211    pfact = lift(factornf(x^2-a,nf.pol)[1,1]));
212  if( poldegree(pfact) == 2,
213if( DEBUGLEVEL_ell >= 5, print("fin de nfsqrt"));
214    return([]));
215if( DEBUGLEVEL_ell >= 5, print("fin de nfsqrt"));
216  return([subst(polcoeff(pfact,0),y,Mod(variable(nf.pol),nf.pol))]);
217}
218if( DEBUGLEVEL_ell >= 4, print("sqrtrat"));
219{
220sqrtrat(a) =
221  sqrtint(numerator(a))/sqrtint(denominator(a));
222}
223\\
224\\ Fonctions propres a ellQ.gp
225\\
226
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}
316if( DEBUGLEVEL_ell >= 4, print("polratroots"));
317{
318polratroots(pol) =
319local(f,ans);
320  f = factor(pol)[,1];
321  ans=[];
322  for( j = 1, #f,
323    if( poldegree(f[j]) == 1,
324      ans = concat(ans,[-polcoeff(f[j],0)/polcoeff(f[j],1)])));
325  return(ans);
326}
327if( DEBUGLEVEL_ell >= 4, print("ratpoint"));
328{
329ratpoint(pol,lim=1,singlepoint=1,tryhard=0) =
330\\ Recherche de points sur y^2=pol(x).
331\\ Les coeff de pol sont entiers.
332\\ Si singlepoint >= 1, cherche un seul point, sinon plusieurs.
333\\ Si tryhard == 1, on essaye une autre strategie quand pol est imprimitif.
334
335local(listpoints,point1,odd,deg4,pol16,tab16,pol9,tab9,pol5,tab5,pol0,vecz,vecx,lead,zz,xx,x16,x9,x5,evpol,ix,iz,K,factK,e,cont,ind,p,sol,factpol,r,M,U);
336
337if( DEBUGLEVEL_ell >= 4, print("entree dans ratpoint avec pol = ",pol); print("lim = ",lim););
338  if( !singlepoint, listpoints = []);
339  point1 = [];
340
341\\          cas triviaux
342  if( issquare(pollead(pol)),
343    point1 = [ 1, sqrtrat(pollead(pol)), 0];
344if( DEBUGLEVEL_ell >= 3, print("solution triviale: a est un carre"));
345    if( singlepoint,
346if( DEBUGLEVEL_ell >= 4, print("fin de ratpoint"));
347      return(point1));
348    listpoints = concat(listpoints,[point1]));
349  if( issquare(polcoeff(pol,0)),
350    point1 = [ 0, sqrtrat(polcoeff(pol,0)) ];
351if( DEBUGLEVEL_ell >= 3, print("solution triviale: e est un carre"));
352    if( singlepoint,
353if( DEBUGLEVEL_ell >= 4, print("fin de ratpoint"));
354      return(point1));
355    listpoints = concat(listpoints,[point1]));
356  odd = poldegree(pol)%2;
357  deg4 = poldegree(pol) == 4;
358
359\\ initialisation du crible modulo 16, 9 et 5
360  if( deg4,
361    pol16 = (Vec(pol)*Mod(1,16))~;
362    tab16 = matrix(16,16);
363    for(xx = 0, 16-1,
364      for(zz = 0, 16-1,
365        tab16[xx+1,zz+1] = !issquare([xx^4,xx^3*zz,xx^2*zz^2,xx*zz^3,zz^4]*pol16)));
366    pol9 = (Vec(pol)~)*Mod(1,9);
367    tab9 = matrix(9,9);
368    for(xx = 0, 9-1,
369      for(zz = 0, 9-1,
370        tab9[xx+1,zz+1] = !issquare([xx^4,xx^3*zz,xx^2*zz^2,xx*zz^3,zz^4]*pol9)));
371    pol5 = (Vec(pol)~)*Mod(1,5);
372    tab5 = matrix(5,5);
373    for(xx = 0, 5-1,
374      for(zz = 0, 5-1,
375        tab5[xx+1,zz+1] = !issquare([xx^4,xx^3*zz,xx^2*zz^2,xx*zz^3,zz^4]*pol5)))
376  );
377 
378  lead = pollead(pol);
379  pol0 = polcoeff(pol,0);
380
381  if( odd,
382    vecz = vector(lim,i,i^2);
383  ,
384\\ si le degre de pol est pair, il faut que le coeff dominant soit
385\\ un carre mod zz.
386    vecz = vector(lim);
387    zz = 0;
388    for( i = 1, lim,
389      zz++; while( !issquare(Mod(lead,zz)),zz++); vecz[i] = zz
390  ));
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);
396
397if( DEBUGLEVEL_ell >= 4, print("xmax = ",vecx[lim]));
398if( DEBUGLEVEL_ell >= 4, print("zmax = ",vecz[lim]));
399
400if( DEBUGLEVEL_ell >= 5, print("vecx = ",vecx));
401if( DEBUGLEVEL_ell >= 5, print("vecz = ",vecz));
402
403\\ boucle sur x = xx/zz
404  for( somme = 2, 2*lim,
405    for( ix = max(1,somme-lim), min(lim,somme-1),
406      xx = vecx[ix]; iz = somme-ix; zz = vecz[iz];
407      if( gcd(zz,xx) > 1, next);
408      if( odd && !issquare(lead*Mod(xx,zz)), next);
409      for( eps = 1, 2, if( eps == 2, zz = -zz);
410      if( deg4 &&
411        (tab16[xx%16+1,zz%16+1] || tab9[xx%9+1,zz%9+1] || tab5[xx%5+1,zz%5+1])
412      , next);
413      evpol = subst(pol,variable(pol),xx/zz);
414      if( issquare(evpol),
415        point1 = [xx/zz,sqrtrat(evpol)];
416        if( singlepoint, break(3));
417        listpoints = concat(listpoints,[point1])
418  ))));
419
420  if( point1 != [],
421if( DEBUGLEVEL_ell >= 3, print("point trouve par ratpoint = ",point1));
422if( DEBUGLEVEL_ell >= 4, print("sortie de ratpoint "));
423    if( singlepoint, return(point1), return(listpoints))
424  );
425
426\\
427\\ Essaye une autre strategie quand pol a un content non trivial
428\\
429
430  if( !odd && tryhard,
431if( DEBUGLEVEL_ell >= 4, print(" Autre strategie dans ratpoint **********"));
432    K = content(pol);
433    if( K != 1,
434      pol /= K;
435      factK = factor(K);
436      e = factK[,2]\2;
437      cont = factorback(factK[,1],e);
438      K /= cont^2;
439      if(K != 1,
440        e = factK[,2]%2;
441        ind = #e; while( !e[ind], ind--);
442        p = factK[ind,1];
443        if( valuation( pollead(pol), p) == 1 ||
444            ( valuation( pollead(pol), p) >= 2 && valuation( polcoeff(pol,poldegree(pol)-1), p) == 0),
445if( DEBUGLEVEL_ell >= 4, print(" utilise une racine de pol mod p = ",p));
446          sol = ratpoint(K/p^2*subst(polrecip(pol),variable(pol),p*variable(pol)),lim,singlepoint,1);
447          if( #sol > 0,
448            point1 = [ 1/(sol[1]*p),
449              sol[2]*cont*p/(p*sol[1])^(poldegree(pol)/2) ];
450if( DEBUGLEVEL_ell >= 4, print("sortie de ratpoint ",point1));
451            return(point1)
452          )
453        );
454        factpol = factormod(pol,p)[,1];
455        for( i = 1, #factpol,
456          if( poldegree(factpol[i]) !=1, next);
457if( DEBUGLEVEL_ell >= 4, print(" utilise une racine de pol mod p = ",p));
458          r = -centerlift(polcoeff(factpol[i],0));
459          if( valuation(subst(pol,variable(pol),r),p) > 2, next);
460          M = [p,r;0,1];
461          U = redquartique(subst(K*pol,variable(pol),p*variable(pol)+r));
462          if( content(U[1]) != p, next);
463          sol = ratpoint(K/p^2*U[1],lim,singlepoint,1);
464          if( #sol > 0,
465            M = (M*U[2])*[sol[1], #sol == 2]~;
466            point1 = [ M[1]/M[2], sol[2]*cont*p/M[2]^(poldegree(pol)/2) ];
467if( DEBUGLEVEL_ell >= 4, print("sortie de ratpoint ",point1));
468            return(point1)
469          )
470        )
471      )
472    )
473  );
474 
475  return([]);
476}
477if( DEBUGLEVEL_ell >= 5, print("psquare"));
478{
479psquare( a, p) =
480local(ap,v);
481
482if( DEBUGLEVEL_ell >= 5, print("entree dans psquare ",[a,p]));
483\\ a est un entier
484\\ renvoie 1 si a est un carre dans Zp 0 sinon
485  if( a == 0,
486if( DEBUGLEVEL_ell >= 5, print("fin de psquare 1"));
487    return(1));
488\\
489  v = valuation(a,p);
490  if( v%2,
491if( DEBUGLEVEL_ell >= 5, print("fin de psquare 0"));
492    return(0));
493  if( p == 2,
494    ap = (a>>v)%8-1,
495    ap = kronecker(a/p^v,p)-1
496  );
497if( DEBUGLEVEL_ell >= 5, print("fin de psquare ", !ap));
498  return(!ap);
499}
500if( DEBUGLEVEL_ell >= 4, print("lemma6"));
501{
502lemma6(pol, p, nu, xx) =
503local(gx,gpx,lambda,mu);
504
505\\ pour les p <> 2
506  gx = subst( pol, variable(pol), xx);
507  if( psquare(gx,p), return(1));
508  gpx = subst( pol', variable(pol), xx);
509  lambda = valuation(gx,p); mu = valuation(gpx,p);
510
511  if( lambda > 2*mu, return(1));
512\\  if( (lambda >= mu+nu) && (nu > mu), return(1));
513  if( (lambda >= 2*nu) && (mu >= nu), return(0));
514  return(-1);
515}
516if( DEBUGLEVEL_ell >= 4, print("lemma7"));
517{
518lemma7( pol, nu, xx) =
519local(gx,gpx,lambda,mu,q);
520
521\\ pour p = 2
522  gx = subst( pol, variable(pol), xx);
523  if( psquare(gx,2), return(1));
524  gpx = subst( pol', variable(pol), xx);
525  lambda = valuation(gx,2); mu = valuation(gpx,2);
526  if( lambda > 2*mu, return(1));
527  if( nu > mu,
528    if( lambda%2, return(-1));
529    q = mu+nu-lambda;
530    if( q == 1, return(1));
531    if( q == 2 && (gx>>lambda)%4 == 1, return(1));
532    return(-1));
533  q = lambda-2*nu;
534  if( q >= 0, return(0));
535  if( q == -2 && (gx>>lambda)%4 == 1, return(0));
536  return(-1);
537}
538if( DEBUGLEVEL_ell >= 4, print("zpsoluble"));
539{
540zpsoluble(pol, p, nu, pnu, x0, pnup) =
541local(result,pol2,fact,x1);
542
543if( DEBUGLEVEL_ell >= 5, print("entree dans zpsoluble ",[pol,p,x0,nu]));
544  if( p == 2,
545    result = lemma7(pol,nu,x0),
546    result = lemma6(pol,p,nu,x0));
547  if( result == +1,
548if( DEBUGLEVEL_ell >= 5, print("fin de zpsoluble 1 lemma"));
549    return(1));
550  if( result == -1,
551if( DEBUGLEVEL_ell >= 5, print("fin de zpsoluble 0 lemma"));
552    return(0));
553  pnup = pnu*p;
554  nu++;
555  if( p< LIMBIGPRIME || !LIMBIGPRIME,
556    for( i = 0, p-1,
557      if( zpsoluble(pol,p,nu,pnup,x0+pnu*i),
558if( DEBUGLEVEL_ell >= 5, print("fin de zpsoluble"));
559        return(1)))
560  ,
561    pol2 = subst(pol,variable(pol),x0+pnu*variable(pol));
562    pol2 /= content(pol2);
563    pol2 = pol2*Mod(1,p);
564    if( !poldegree(pol2), return(0));
565    fact = factormod(pol2,p)[,1];
566    for( i = 1, #fact,
567      x1 = -centerlift(polcoeff(fact[i],0));
568      if( zpsoluble(pol,p,nu,pnup,x0+pnu*x1),
569if( DEBUGLEVEL_ell >= 5, print("fin de zpsoluble"));
570        return(1)));
571    for( i = 1, MAXPROB,
572      x1 = random(p);
573      if( zpsoluble(pol,p,nu,pnup,x0+pnu*x1),
574if( DEBUGLEVEL_ell >= 5, print("fin de zpsoluble"));
575        return(1)))
576  );
577if( DEBUGLEVEL_ell >= 2,
578  if( p >= LIMBIGPRIME,
579    print("******* test probabiliste en p = ",p,"*******")));
580if( DEBUGLEVEL_ell >= 5, print("fin de zpsoluble"));
581  return(0);
582}
583if( DEBUGLEVEL_ell >= 4, print("qpsoluble"));
584{
585qpsoluble(pol, p) =
586if( DEBUGLEVEL_ell >= 5, print("entree dans qpsoluble ",p); print("pol = ",pol));
587  if( psquare(pollead(pol),p),
588if( DEBUGLEVEL_ell >= 5, print("fin de qpsoluble 1"));
589    return(1));
590  if( psquare(polcoeff(pol,0),p),
591if( DEBUGLEVEL_ell >= 5, print("fin de qpsoluble 1"));
592    return(1));
593  if( zpsoluble(pol,p,0,1,0),
594if( DEBUGLEVEL_ell >= 5, print("fin de qpsoluble 1"));
595    return(1));
596  if( zpsoluble(polrecip(pol),p,1,p,0),
597if( DEBUGLEVEL_ell >= 5, print("fin de qpsoluble 1"));
598    return(1));
599if( DEBUGLEVEL_ell >= 5, print("fin de qpsoluble 0"));
600  return(0);
601}
602if( DEBUGLEVEL_ell >= 4, print("locallysoluble"));
603{
604locallysoluble(pol) =
605\\ teste l'existence locale de solutions de y^2 = pol(x,z)
606local(plist,disc0,p,c,vc);
607
608if( DEBUGLEVEL_ell >= 4, print("entree dans locallysoluble :",pol));
609
610\\ place reelle
611  if( !(poldegree(pol)%2) && sign(pollead(pol)) < 0
612         && sign(polcoeff(pol,0)) < 0 && polsturm(pol) == 0,
613if( DEBUGLEVEL_ell >= 3, print(" non ELS a l'infini"));
614if( DEBUGLEVEL_ell >= 4, print("fin de locallysoluble"));
615    return(0));
616
617\\
618\\ places finies de plist */
619\\
620  pol *= denominator(content(pol))^2;
621  c = content(pol);
622
623  disc0 = poldisc(pol);
624  plist = factor (abs(2*disc0));
625if( DEBUGLEVEL_ell >= 4, print("liste de premiers = ",plist));
626  for( i = 1, #plist[,1],
627    p = plist[i,1];
628if( DEBUGLEVEL_ell >= 4, print("p = ",p));
629    vc = valuation(c,p);
630    if( vc >= 2,
631      pol /= p^(2*(vc\2));
632      plist[i,2] -= 2*(vc\2)*(2*poldegree(pol)-2)
633    );
634    if( poldegree(pol) == 4 && p != 2 && plist[i,2] < 2, next);
635    if( !qpsoluble(pol,p),
636if( DEBUGLEVEL_ell >= 3, print(" non ELS en ",p));
637if( DEBUGLEVEL_ell >= 4, print("fin de locallysoluble"));
638      return(0)));
639
640if( DEBUGLEVEL_ell >= 2, print(" quartique ELS"));
641if( DEBUGLEVEL_ell >= 4, print("fin de locallysoluble"));
642  return(1);
643}
644if( DEBUGLEVEL_ell >= 4, print("redquartique"));
645{
646redquartique(pol) =
647\\ reduction d'une quartique.
648\\ ou plus generalement d'un polynome de degre deg.
649local(prec,prec0,d,disc2,test,normderiv,disc2v,r,q,M);
650
651if( DEBUGLEVEL_ell >= 4, print("entree dans redquartique"));
652if( DEBUGLEVEL_ell >= 3, print(" reduction de la quartique ",pol));
653
654\\ choix de la precision des calculs.
655  prec = prec0 = default(realprecision);
656  d = poldegree(pol);
657  disc2 = poldisc(pol)^2;
658  test = 0;
659  while( test == 0,
660if( DEBUGLEVEL_ell >= 4, print(" precision = ",prec));
661    r = polroots(pol);
662    normderiv = vector( d, i, norm(subst(pol',variable(pol),r[i])));
663    disc2v = prod( i = 1, d, normderiv[i]) * pollead(pol)^(2*d-4);
664    test = abs(disc2v-disc2) < 10^(-prec\2);
665    if( !test, default(realprecision, prec *= 2))
666  );
667
668\\ On n'utilise plus
669\\  q = Vec(sum( i = 1, d, norm(x-r[i])));
670\\ mais la normalisation de Cremona-Stoll
671  q = Vec(sum( i = 1, d, norm(x-r[i]) / normderiv[i]^(1/(d-2))));
672  M = QfbReduce([q[1],q[2]/2;q[2]/2,q[3]]);
673  pol = subst(pol,variable(pol),Pol(M[1,])/Pol(M[2,]))*Pol(M[2,])^poldegree(pol);
674
675  if( prec != prec0, default(realprecision,prec0));
676
677if( DEBUGLEVEL_ell >= 3, print(" quartique reduite = ",pol));
678if( DEBUGLEVEL_ell >= 4, print("sortie de redquartique"));
679
680  return([pol,M]);
681}
682if( DEBUGLEVEL_ell >= 4, print("reducemodsquares"));
683{
684reducemodsquares(delta,d) =
685\\ reduction du coefficient de x^d dans ( delta modulo les carres )
686local(deg,xx,z,qd,Qd,reduc);
687
688  deg = poldegree(delta.mod);
689  xx = Mod(x,delta.mod);
690  z = subst(Pol(vector(deg,i,eval(Str("a"i)))),x,xx);
691  qd = polcoeff(lift(delta*z^2),d,x);
692  Qd = simplify(matrix(deg,deg,i,j,deriv(deriv(qd,eval(Str("a"i))),eval(Str("a"j)))/2));
693
694  reduc = IndefiniteLLL(Qd);
695  if( #reduc == 2, reduc = reduc[2][,1]);
696
697  return(delta*subst(Pol(reduc),x,xx)^2);
698}
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}
715if( DEBUGLEVEL_ell >= 4, print("ellredgen"));
716{
717ellredgen(ell,listgen,K=1) =
718\\ reduction des generateurs de listgen
719\\ sur la courbe ell = [a1,a2,a3,a4,a6]
720\\ ou K*y^2 = x^3 + a2*x^2 + a4*x + a6 (lorsque a1 = a3 = 0);
721local(d,sqrtK,urst,M,U,limgoodrelations,listgen2);
722
723if( DEBUGLEVEL_ell >= 3, print("Reduction des generateurs ",listgen));
724if( DEBUGLEVEL_ell >= 5, print("ell=",ell));         
725  d = #listgen;
726  if( d == 0, return([]));
727
728  if( #ell < 19, ell = ellinit(ell));
729  if( K != 1,
730    if( ell.a1 != 0 || ell.a3 != 0, error(" ellredgen : a1*a3 != 0"));
731    ell[2] *= K; ell[4] *= K^2; ell[5] *= K^3;
732    ell[6] *= K; ell[7] *= K^2; ell[8] *= K^3; ell[9] *= K^4;
733    ell[10] *= K^2; ell[11] *= K^3; ell[12] *= K^6;
734    sqrtK = sqrt(K);
735    ell[14] *= K;
736    ell[15] /= sqrtK; ell[16] /= sqrtK;
737    ell[17] *= sqrtK; ell[18] *= sqrtK;
738    ell[19] /= K;
739
740    for( i = 1, d,
741      for( j = 1, #listgen[i],
742        listgen[i][j] *= K^j))
743  );
744
745  ell = ellminimalmodel(ell,&urst);
746  listgen = ellchangepoint(listgen,urst);
747
748\\ Recherche des relations entre les points de listgen
749\\ par recherche du noyau de la matrice des hauteurs
750
751  M = ellheightmatrix(ell,listgen);
752if( DEBUGLEVEL_ell >= 4, print("matrice des hauteurs = ",M));
753  M = round( M*10^(default(realprecision)-10) );
754  U = qflll(M,4);
755  U = concat(U[1],U[2]);
756  limgoodrelations = 0;
757  while( limgoodrelations+1 <= d
758    && vecmax(abs(U[,limgoodrelations+1])) < 20, limgoodrelations++);
759  U = vecextract(U,1<<limgoodrelations-1);
760  U = completebasis(U);
761if( DEBUGLEVEL_ell >= 4, print("changement de base = ",U));
762 
763  listgen2 = vector(d);
764  for( i = 1, d,
765    listgen2[i] = [0];
766    for( j = 1, d,
767      listgen2[i] = elladd(ell,listgen2[i],ellpow(ell,listgen[j],U[j,i]))));
768  listgen = listgen2;
769
770\\ Tri des points d'ordre infini
771
772  listgen2 = [];
773  for( i = 1, d,
774    if( !ellorder(ell,listgen[i]),
775      listgen2 = concat(listgen2,[listgen[i]])));
776  listgen = listgen2;
777if( DEBUGLEVEL_ell >= 3, print("points d'ordre infini = ",listgen));
778  d = #listgen;
779  if( d == 0, return([]));
780
781\\ Reduction des points d'ordre infini
782
783  if( d > 1,
784    M = ellheightmatrix(ell,listgen);
785if( DEBUGLEVEL_ell >= 4, print("matrice des hauteurs = ",M));
786    U = qflllgram(M);
787if( DEBUGLEVEL_ell >= 4, print("changement de base = ",U));
788
789    listgen2 = vector(d);
790    for( i = 1, d,
791      listgen2[i] = [0];
792      for( j = 1, d,
793        listgen2[i] = elladd(ell,listgen2[i],ellpow(ell,listgen[j],U[j,i]))));
794    listgen = listgen2
795  );
796
797  listgen = ellsort(listgen);
798if( DEBUGLEVEL_ell >= 4, print("generateurs tries = ",listgen));
799
800  listgen = ellchangepointinverse(listgen,urst);
801  if( K != 1,
802    for( i = 1, d,
803      for( j = 1, 2,
804        listgen[i][j] /= K^j)));
805
806\\ on ne garde que les points (x,y) avec y >= 0
807
808  if( ell.a1 == 0 && ell.a3 == 0,
809    for( i = 1, d,
810      if( #listgen[i] == 2,
811        listgen[i][2] = abs(listgen[i][2]))));
812
813if( DEBUGLEVEL_ell >= 2, print("generateurs reduits = ",listgen));
814  return(listgen);
815}
816if( DEBUGLEVEL_ell >= 4, print("ell2descent_gen"));
817{
818ell2descent_gen(ell,ext,K=1,help=[],redflag=0) =
819\\ si ell= K*y^2=P(x), alors ext est le buchinitfu de l'extension.
820\\ theta est une racine de P.
821\\ dans la suite ext est note L = Q(theta).
822\\ help est une liste de points deja connus sur ell.
823\\ ell est de la forme K*y^2=x^3+A*x^2+B*x+C */
824\\ ie ell=[0,A,0,B,C], avec A,B et C entiers */
825\\
826\\ si redflag != 0, on utilise la reduction modulo les carres.
827\\
828
829local(A,B,C,S,SLprod,SLlist,aux,oddclass,LS2gen,fact,trouve,i,polrel,ttheta,polprime,KS2gen,LS2genunit,normLS2gen,normcoord,LS2coordtilda,LS2,listgen,listpoints,listpointstriv,listpointsmwr,list,m1,m2,lastloc,maskwhile,iwhile,zc,j,iaux,liftzc,den,ispointtriv,point,found,idealfactorzc,idealzc,baseidealzc,q2,sol,q1,param,pol,redq,q0,pointxx,point2,v,rang);
830
831if( DEBUGLEVEL_ell >= 4, print("entree dans ell2descent_gen"));
832
833\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
834\\      construction de L(S,2)         \\
835\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
836
837  if( #ell < 13, ell = ellinit(ell,1));
838
839  if( ell.a1 != 0 || ell.a3 != 0,
840    error(" ell2descent_gen : la courbe n'est pas sous la forme [0,a,0,b,c]"));
841  if( denominator(ell.a2) > 1 || denominator(ell.a4) > 1 || denominator(ell.a6) >1,
842    error(" ell2descent_gen : coefficients non entiers"));
843
844  A = ell.a2; if( DEBUGLEVEL_ell >= 2, print("A = ",A));
845  B = ell.a4; if( DEBUGLEVEL_ell >= 2, print("B = ",B));
846  C = ell.a6; if( DEBUGLEVEL_ell >= 2, print("C = ",C));
847
848  polrel = Pol([1,A,B,C]);
849  if( !ext,
850if( DEBUGLEVEL_ell >= 2, print("bnfinit(",polrel,")"));
851    ext = bnfinit(polrel,1));
852
853  S = -abs(K*ext.index);
854  SLprod = idealmul(ext,K,idealadd(ext,ext.pol',ext.index));
855  SLlist = idealfactor(ext,SLprod)[,1]~;
856  aux = []; SLprod = 1;
857   for( i = 1, #SLlist,
858     if( !(K%SLlist[i][1]) || valuation(ext.index,SLlist[i][1]),
859       SLprod = idealmul(ext,SLprod,SLlist[i]);
860       aux = concat(aux,[SLlist[i]])));
861  SLlist = aux;
862  oddclass = 0;
863  while( !oddclass,
864\\ Constructoin de S:
865if( DEBUGLEVEL_ell >= 4, print("SLlist = ",SLlist));
866\\ Construction des S-unites
867    LS2gen = bnfsunit(ext,SLlist);
868if( DEBUGLEVEL_ell >= 4, print("LS2gen = ",LS2gen));
869\\ on ajoute la partie paire du groupe de classes.
870    oddclass = LS2gen[5][1]%2;
871    if( !oddclass,
872if( DEBUGLEVEL_ell >= 3, print("Groupe de classes pair"));
873if( DEBUGLEVEL_ell >= 4, print(LS2gen[5]));
874      S *= LS2gen[5][3][1][1,1];
875      SLprod = idealmul(ext,SLprod,LS2gen[5][3][1]);
876      fact = idealfactor(ext,LS2gen[5][3][1])[,1];
877      trouve = 0; i = 0;
878      while( !trouve,
879        i++; trouve = 1;
880        for( j = 1, #SLlist,
881          if( SLlist[j] == fact[i], trouve = 0; break)));
882      SLlist = concat(SLlist,[fact[i]]))
883  );
884
885if( DEBUGLEVEL_ell >= 4, print("S = ",S));
886
887  ttheta = Mod(x,polrel);            \\ ttheta est la racine de P(x)
888  polprime = Mod(polrel',polrel);
889
890  KS2gen = factor(S)[,1]~;
891 
892if( DEBUGLEVEL_ell >= 3, print("#KS2gen = ",#KS2gen));
893if( DEBUGLEVEL_ell >= 3, print("KS2gen = ",KS2gen));
894
895  LS2genunit = ext.tufu;
896  LS2genunit = concat(LS2gen[1],LS2genunit);
897
898  LS2genunit = subst(LS2genunit,x,ttheta);
899  LS2genunit = LS2genunit*Mod(1,polrel);
900if( DEBUGLEVEL_ell >= 3, print("#LS2genunit = ",#LS2genunit));
901if( DEBUGLEVEL_ell >= 3, print("LS2genunit = ",LS2genunit));
902   
903\\ dans LS2gen, on ne garde que ceux dont la norme
904\\ est un carre.
905
906  normLS2gen = norm(LS2genunit);
907if( DEBUGLEVEL_ell >= 4, print("normLS2gen = ",normLS2gen));
908
909\\ matrice de l'application norme
910
911  normcoord = matrix(#KS2gen,#normLS2gen);
912  for( j = 1, #normLS2gen,
913    normcoord[1,j] = (sign(normLS2gen[j]) < 0);
914    for( i = 2, #KS2gen,
915      normcoord[i,j] = valuation(normLS2gen[j],KS2gen[i])));
916if( DEBUGLEVEL_ell >= 4, print("normcoord = ",normcoord));
917
918\\ construction du noyau de la norme
919
920  LS2coordtilda = lift(matker(normcoord*Mod(1,2)));
921if( DEBUGLEVEL_ell >= 4, print("LS2coordtilda = ",LS2coordtilda));
922  LS2 = vector(#LS2coordtilda[1,],i,0);
923  for( i = 1, #LS2coordtilda[1,],
924    aux = 1;
925    for( j = 1, #LS2coordtilda[,i],
926      if( sign(LS2coordtilda[j,i]),
927        aux *= LS2genunit[j]));
928    LS2[i] = aux
929  );
930if( DEBUGLEVEL_ell >= 3, print("LS2 = ",LS2));
931if( DEBUGLEVEL_ell >= 3, print("norm(LS2) = ",norm(LS2)));
932
933\\ Reduction des generateurs de LS2
934 
935  if( redflag,
936    for( i = 1, #LS2,
937      LS2[i] = reducemodsquares(LS2[i],2)));
938
939\\ Fin de la construction de LS2
940
941  listgen = LS2;
942if( DEBUGLEVEL_ell >= 2, print("LS2gen = ",listgen));
943if( DEBUGLEVEL_ell >= 2, print("#LS2gen = ",#listgen));
944  listpoints = [];
945
946if( DEBUGLEVEL_ell >= 3, print("(A,B,C) = ",[A,B,C]));
947
948\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
949\\   Recherche de points triviaux.  \\
950\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
951
952if( DEBUGLEVEL_ell >= 2, print("Recherche de points triviaux sur la courbe"));
953  listpointstriv = ratpoint(K^3*subst(polrel,x,x/K),LIMTRIV,0,0);
954  for( i = 1, #listpointstriv,
955    if( #listpointstriv[i] == 3,
956      listpointstriv[i] = [0]
957    , for( j = 1, 2, listpointstriv[i][j] /= K^j))
958   );
959  listpointstriv = concat(help,listpointstriv);
960if( DEBUGLEVEL_ell >= 1, print("Points triviaux sur la courbe = ",listpointstriv));
961
962
963\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
964\\          parcours de L(S,2)         \\
965\\ \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
966
967  listpointsmwr = [];
968  list = [ 6, ell.disc, 0 ];
969  m1 = 0; m2 = 0; lastloc = -1;
970  maskwhile = 1<<#listgen;
971  iwhile = 1;
972  while( iwhile < maskwhile,
973if( DEBUGLEVEL_ell >= 4, print("iwhile = ",iwhile); print("listgen = ",listgen));
974    zc = Mod(1,polrel); j = 1; iaux = iwhile;
975    while( iaux,
976      if( iaux%2, zc *= listgen[j]);
977      iaux >>= 1; j++);
978if( DEBUGLEVEL_ell >= 2, print(); print("zc = ",zc));
979    liftzc = lift(zc);
980    if( redflag,
981      zc = reducemodsquares(zc,2);
982      liftzc = lift(zc);
983      den = denominator(content(liftzc))^2;
984      zc *= den; liftzc *= den;
985if( DEBUGLEVEL_ell >= 2, print("zcred = ",zc))
986    );
987
988\\ Est-ce un point trivial ?
989    ispointtriv = 0;
990    for( i = 1, #listpointstriv,
991      point = listpointstriv[i];
992      if( #point == 2,
993        if( nfissquare(ext.nf,K*(point[1]-x)*liftzc),
994if( DEBUGLEVEL_ell >= 2, print(" vient du point trivial ",point));
995          listpointsmwr = concat(listpointsmwr,[point]);
996          m1++;
997          if( degre(iwhile) > lastloc, m2++);
998          found = (ispointtriv = 1);
999          break
1000    )));
1001   
1002\\ Ce n'est pas un point trivial
1003    if( !ispointtriv,
1004\\ Il faut resoudre une forme quadratique
1005\\      q2 = matrix(3,3,i,j,trace(zc*ttheta^(i+j-2)/polprime));
1006\\if( DEBUGLEVEL_ell >= 4, print("q2 = ",q2));
1007      idealfactorzc = idealfactor(ext,zc);
1008      idealfactorzc[,2] *= -1;
1009      idealfactorzc[,2] \= 2;
1010      idealzc = factorback(idealfactorzc,,ext);
1011      if( idealzc == 1, idealzc = matid(3));
1012      baseidealzc = vector(3,i,nfbasistoalg(ext,idealzc[,i]));
1013      q2 = matrix(3,3,i,j,trace(zc*baseidealzc[i]*baseidealzc[j]/polprime));
1014if( DEBUGLEVEL_ell >= 4, print("q2 = ",q2));
1015\\      q2 *= ext.index;
1016\\ if( DEBUGLEVEL_ell >= 4, print("q2 = ",q2));
1017if( DEBUGLEVEL_ell >= 4, print("q2/content(q2) = ",q2/content(q2)));
1018      sol = Qfsolve(q2/content(q2));
1019if( DEBUGLEVEL_ell >= 4,print("sol = ",sol));
1020      if( type(sol) == "t_INT",
1021if( DEBUGLEVEL_ell >= 3, print("non ELS en ",sol));
1022        iwhile++; next
1023      );
1024
1025\\ \\\\\\\\\\\
1026\\ Construction de la quartique
1027\\ \\\\\\\\\\\
1028      q1 = -matrix(3,3,i,j,trace(zc*baseidealzc[i]*baseidealzc[j]*(ttheta+A)/polprime));
1029      param = Qfparam(q2,sol)*[x^2,x,1]~;
1030      param /= content(param);
1031      pol = param~*q1*param;
1032if( DEBUGLEVEL_ell >= 2, print(" quartique: ",K,"*Y^2 = ",pol));
1033      redq = redquartique(pol);
1034if( DEBUGLEVEL_ell >= 2, print(" reduite: ",K,"*Y^2 = ",redq[1]));
1035      pol = redq[1];
1036      den = denominator(content(K*pol));
1037      pol *= den^2;
1038
1039\\ \\\\\\\\\\\
1040\\ Recherche de points sur la quartique
1041\\ \\\\\\\\\\\
1042
1043      point = ratpoint(K*pol,LIM1,1,0);
1044      if( point != [],
1045        m1++;
1046        if( #point == 2, point = concat(point,[1]));
1047        point = concat(redq[2]*[point[1],point[3]]~,[point[2]/den]~);
1048        param = subst(param,x,x/y)*y^2;
1049        param = subst(subst(param,x,point[1]),y,point[2]);
1050if( DEBUGLEVEL_ell >= 2, print(" point sur la quartique = ",point));
1051        param *= K/point[3];
1052if( DEBUGLEVEL_ell >= 3, print("reconstruction du point sur la courbe"));
1053        q0 = matrix(3,3,i,j,trace(zc*baseidealzc[i]*baseidealzc[j]*(ttheta^2+A*ttheta+B)/polprime));
1054        pointxx = param~*q0*param/K;
1055        point2 = [ pointxx, sqrtrat(subst(x^3+A*x^2+B*x+C,x,pointxx)/K)];
1056if( DEBUGLEVEL_ell >= 1, print("point sur la courbe = ",point2));
1057        listpointsmwr = concat(listpointsmwr,[point2]);
1058        if( degre(iwhile) > lastloc, m2++);
1059        found = 1
1060      ,
1061        if( locallysoluble(K*pol),
1062          if( degre(iwhile) > lastloc, m2++; lastloc = degre(iwhile));
1063          point = ratpoint(K*pol,LIM3,1,1);
1064          if( #point > 0,
1065            m1++;
1066            if( #point == 2, point = concat(point,[1]));
1067            point = concat(redq[2]*[point[1],point[3]]~,[point[2]/den]~);
1068            param = subst(param,x,x/y)*y^2;
1069            param = subst(subst(param,x,point[1]),y,point[2]);
1070if( DEBUGLEVEL_ell >= 3, print(" point sur la quartique = ",point));
1071            param *= K/point[3];
1072if( DEBUGLEVEL_ell >= 3, print("reconstruction du point sur la courbe"));
1073            q0 = matrix(3,3,i,j,trace(zc*baseidealzc[i]*baseidealzc[j]*(ttheta^2+A*ttheta+B)/polprime));
1074            pointxx = param~*q0*param/K;
1075            point2 = [ pointxx, sqrtrat(subst(x^3+A*x^2+B*x+C,x,pointxx)/K)];
1076if( DEBUGLEVEL_ell >= 1, print("point sur la courbe = ",point2));
1077            listpointsmwr = concat(listpointsmwr,[point2]);
1078            found = 1
1079          )
1080        )
1081      )
1082    );
1083    if( found || ispointtriv,
1084      found = 0; lastloc = -1;
1085      v = 0; iaux = (iwhile>>1);
1086      while( iaux, iaux >>= 1; v++);
1087      maskwhile >>= 1;
1088      listgen = vecextract(listgen,1<<#listgen-1<<v-1);
1089      iwhile = 1<<v
1090    , iwhile ++
1091    )
1092  );
1093
1094if( DEBUGLEVEL_ell >= 2,
1095  print();
1096  print("rang des points trouves = ",m1);
1097  print("rang du groupe de Selmer = ",m2));
1098if( DEBUGLEVEL_ell >= 1,
1099  print("#S(E/Q)[2]    = ",1<<m2));
1100  if( m1 == m2,
1101if( DEBUGLEVEL_ell >= 1,
1102    print("#E(Q)/2E(Q)   = ",1<<m1);
1103    print("#III(E/Q)[2]  = 1");
1104    print("rang(E/Q)     = ",m1));
1105    rang = m1
1106  ,
1107if( DEBUGLEVEL_ell >= 1,
1108    print("#E(Q)/2E(Q)  >= ",1<<m1);
1109    print("#III(E/Q)[2] <= ",1<<(m2-m1));
1110    print("rang(E/Q)    >= ",m1));
1111    rang = m1;
1112    if( (m2-m1)%2,
1113if( DEBUGLEVEL_ell >= 1,
1114      print(" III devrait etre un carre, donc ");
1115      if( m2-m1 > 1,
1116        print("#E(Q)/2E(Q)  >= ",1<<(m1+1));
1117        print("#III(E/Q)[2] <= ",1<<(m2-m1-1));
1118        print("rang(E/Q)    >= ",m1+1)
1119      ,
1120        print("#E(Q)/2E(Q)  = ",1<<(m1+1));
1121        print("#III(E/Q)[2] = 1");
1122        print("rang(E/Q)    = ",m1+1)));
1123        rang = m1+1
1124    )
1125  );
1126if( DEBUGLEVEL_ell >= 1, print("listpointsmwr = ",listpointsmwr));
1127  for( i = 1, #listpointsmwr,
1128    if( subst(polrel,x,listpointsmwr[i][1])-K*listpointsmwr[i][2]^2,
1129      error(" ell2descent_gen : MAUVAIS POINT = ",listpointsmwr[i])));
1130  if( #listpointsmwr >= 2,
1131    listpointsmwr = ellredgen(ell,listpointsmwr,K));
1132if( DEBUGLEVEL_ell >= 4, print("fin de ell2descent_gen"));
1133  return([rang,m2,listpointsmwr]);
1134}
1135if( DEBUGLEVEL_ell >= 4, print("ellrank"));
1136{
1137ellrank(ell,help=[]) =
1138\\ Algorithme de la 2-descente sur la courbe elliptique ell.
1139\\ help est une liste de points connus sur ell.
1140local(urst,urst1,den,eqell,tors2,bnf,rang,time1);
1141
1142if( DEBUGLEVEL_ell >= 3, print("entree dans ellrank"));
1143  if( #ell < 13, ell = ellinit(ell,1));
1144
1145\\ supprime les coefficients a1 et a3
1146  urst = [1,0,0,0];
1147  if( ell.a1 != 0 || ell.a3 != 0,
1148    urst1 = [1,0,-ell.a1/2,-ell.a3/2];
1149    ell = ellchangecurve(ell,urst1);
1150    urst = ellcomposeurst(urst,urst1)
1151  );
1152
1153\\ supprime les denominateurs
1154  while( (den = denominator([ell.a2,ell.a4,ell.a6])) > 1,
1155    den = factor(den); den[,2] = vectorv(#den[,2],i,1);
1156    den = factorback(den);
1157    urst1 = [1/den,0,0,0];
1158    ell = ellchangecurve(ell,urst1);
1159    urst = ellcomposeurst(urst,urst1)
1160  );
1161
1162  help = ellchangepoint(help,urst);
1163  eqell = Pol([1,ell.a2,ell.a4,ell.a6]);
1164if( DEBUGLEVEL_ell >= 1, print("courbe elliptique : Y^2 = ",eqell));
1165
1166\\ choix de l'algorithme suivant la 2-torsion
1167
1168  tors2 = ellhalf(ell,[0]);
1169if( DEBUGLEVEL_ell >= 1, print("E[2] = ",tors2)); 
1170
1171  if( #tors2 == 1,                              \\ cas 1: 2-torsion triviale
1172if( DEBUGLEVEL_ell >= 3, print1("bnfinit "));
1173if( DEBUGLEVEL_ell >= 4, gettime());
1174    bnf = bnfinit(eqell,1);
1175if( DEBUGLEVEL_ell >= 4, time1 = gettime());
1176if( DEBUGLEVEL_ell >= 3, print("ok"));
1177    rang = ell2descent_gen(ell,bnf,1,help);
1178if( DEBUGLEVEL_ell >= 4, print("temps dans bnfinit = ",time1));
1179if( DEBUGLEVEL_ell >= 4, print("temps pour le reste = ",gettime()));
1180  ,
1181  if( #tors2 >= 2,                              \\ cas 2: 2-torsion >= Z/2Z
1182    if( ell.a6 != 0,
1183      urst1 = [1,tors2[2][1],0,0];
1184      ell = ellchangecurve(ell,urst1);
1185      urst = ellcomposeurst(urst,urst1)
1186    );
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)
1191  ,                                             \\ cas 3: 2-torsion = Z/2Z*Z/2Z
1192\\    rang = ell2descent_complete(tors2[2][1],tors2[3][2],tors2[4][3])
1193  ));
1194
1195  rang[3] = ellchangepointinverse(rang[3],urst);
1196if( DEBUGLEVEL_ell >= 3, print("fin de ellrank"));
1197
1198return(rang);
1199}
1200if( DEBUGLEVEL_ell >= 4, print("ell2descent_complete"));
1201{
1202ell2descent_complete(e1,e2,e3) =
1203\\ calcul du rang d'une courbe elliptique
1204\\ par la methode de 2-descente complete.
1205\\ Y^2 = (x-e1)*(x-e2)*(x-e3);
1206\\ en suivant la methode decrite par J.Silverman
1207\\ renvoie [r,s,v] avec
1208\\   r est une borne inferieure du rang de E(Q)
1209\\   s est le rang du 2-groupe de Selmer
1210\\   v est un systeme de points independants dans E(Q)/2E(Q)
1211
1212\\ e1, e2 et e3 sont des entiers.
1213
1214local(d32,d31,d21,G1,G2,vect1,vect2,selmer,rang,listepoints,b1,b2,q1,q2,sol1,param1,param1x,sol2,quart,point,z1,solx,soly,strange);
1215
1216if( DEBUGLEVEL_ell >= 2, print("Algorithme de la 2-descente complete"));
1217
1218\\ calcul des groupes G1 et G2
1219
1220  d32 = e3-e2; d31 = e3-e1; d21 = e2-e1;
1221  G1 = factor(-abs(d31*d21))[,1];
1222  G2 = factor(-abs(d32*d21))[,1];
1223
1224if( DEBUGLEVEL_ell >= 3, print("G1 = ",G1));
1225if( DEBUGLEVEL_ell >= 3, print("G2 = ",G2));
1226
1227\\ parcours de G1*G2
1228
1229  vect1 = vector(#G1,i,[0,1]);
1230  vect2 = vector(#G2,i,[0,1]);
1231  selmer = 0;
1232  rang = 0;
1233  listepoints = [];
1234
1235  forvec( X = vect1,
1236    b1 = prod( i = 1, #G1, G1[i]^X[i]);
1237    forvec( Y = vect2,
1238      b2 = prod( i = 1, #G2, G2[i]^Y[i]);
1239
1240if( DEBUGLEVEL_ell >= 3, print("[b1,b2] = ",lift([b1,b2])));
1241
1242\\ points triviaux provenant de la 2-torsion
1243
1244      if( b1==1 && b2==1,
1245if( DEBUGLEVEL_ell >= 4, print(" point trivial [0]"));
1246        selmer++; rang++; next);
1247      if( issquare(-d21*b2) && issquare(d31*d21*b1),
1248if( DEBUGLEVEL_ell >= 3, print(" point trivial [e1,0]"));
1249        selmer++; rang++; listepoints = concat(listepoints,[[e1,0]]); next);
1250      if( issquare(d21*b1) && issquare(-d32*d21*b2),
1251if( DEBUGLEVEL_ell >= 3, print(" point trivial [e2,0]"));
1252        selmer++; rang++; listepoints = concat(listepoints,[[e2,0]]); next);
1253      if( issquare(d31*b1) && issquare(d32*b2),
1254if( DEBUGLEVEL_ell >= 3, print(" point trivial [e3,0]"));
1255        selmer++; rang++; listepoints = concat(listepoints,[[e3,0]]); next);
1256
1257\\ il faut resoudre 2 equations quadratiques:
1258\\ (1) b1*z1^2-b2*z2^2 = e2-e1
1259\\ (2) b1*z1^2-b1*b2*z3^2 = e3-e1
1260\\ on aura alors (x,y) = (b1*z1^2+e1,b1*b2*z1*z2*z3)
1261
1262      q1 = matdiagonal([b1,-b2,-d21]);
1263if( DEBUGLEVEL_ell >= 3, print(" q1 = ",q1));
1264      q2 = matdiagonal([b1,-b1*b2,-d31]);
1265if( DEBUGLEVEL_ell >= 3, print(" q2 = ",q2));
1266
1267\\ solution de la premiere forme quadratique
1268
1269      sol1 = Qfsolve(q1);
1270      if( type(sol1) == "t_INT",
1271if( DEBUGLEVEL_ell >= 3, print(" q1 non ELS en ",sol1));
1272        next);
1273if( DEBUGLEVEL_ell >= 3, print(" sol part de q1 = ",sol1));
1274      param1 = Qfparam(q1,sol1,1);
1275if( DEBUGLEVEL_ell >= 3, print(" param de q1 = ",param1));
1276      param1x = param1*[x^2,x,1]~;
1277
1278      sol2 = Qfsolve(q2);
1279      if( type(sol2) == "t_INT",
1280if( DEBUGLEVEL_ell >= 3, print(" q2 non ELS en ",sol2));
1281        next);
1282
1283      quart = b1*b2*(b1*param1x[1]^2-d31*param1x[3]^2);
1284if( DEBUGLEVEL_ell >= 3, print(" quart = ",quart));
1285
1286\\ la quartique est-elle localement soluble ?
1287   
1288      if( !locallysoluble(quart),
1289if( DEBUGLEVEL_ell >= 3, print(" quartique non ELS "));
1290        next);
1291if( DEBUGLEVEL_ell >= 2, print(" y^2 = ",quart));
1292      selmer++;
1293
1294\\ recherche de points sur la quartique.
1295   
1296      point = ratpoint(quart,LIM3,1);
1297      if( point != [],
1298if( DEBUGLEVEL_ell >= 2, print("point trouve sur la quartique !!"));
1299if( DEBUGLEVEL_ell >= 3, print(point));
1300        if( #point == 2,
1301          z1 = subst(param1x[1],x,point[1])/subst(param1x[3],x,point[1])
1302        , z1 = param1[1,1]/param1[3,1]);
1303        solx = b1*z1^2+e1;
1304        soly = sqrtrat((solx-e1)*(solx-e2)*(solx-e3));
1305        listepoints = concat(listepoints,[[solx,soly]]);
1306if( DEBUGLEVEL_ell >= 1, print("point sur la courbe elliptique = ",[solx,soly]));
1307        rang++
1308      ,
1309if( DEBUGLEVEL_ell >= 2, print("aucun point trouve sur la quartique"))
1310      )
1311    )
1312  );
1313
1314\\ fin
1315
1316if( DEBUGLEVEL_ell >= 1,
1317  print("#S^(2)      = ",selmer));
1318  if( rang > selmer/2, rang = selmer);
1319if( DEBUGLEVEL_ell >= 1,
1320  strange = rang != selmer;
1321  if( strange,
1322  print("#E[K]/2E[K]>= ",rang)
1323, print("#E[K]/2E[K] = ",rang));
1324  print("#E[2]       = 4");
1325);
1326  rang = ceil(log(rang)/log(2))-2;
1327  selmer = valuation(selmer,2);
1328if( DEBUGLEVEL_ell >= 1,
1329  if( strange,
1330  print(selmer-2," >= rang  >= ",rang)
1331, print("rang        = ",rang));
1332  if( rang, print("points = ",listepoints));
1333);
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);
1337
1338return([rang,selmer,listepoints]);
1339}
1340if( DEBUGLEVEL_ell >= 4, print("ellcount"));
1341{
1342ellcount( c, d, KS2gen, listpointstriv=[]) =
1343local(found,listgen,listpointscount,m1,m2,lastloc,mask,i,d1,iaux,j,triv,pol,point,deuxpoints,v);
1344
1345if( DEBUGLEVEL_ell >= 4, print("entree dans ellcount ",[c,d]));
1346if( DEBUGLEVEL_ell >= 4, print("KS2gen = ",KS2gen));
1347if( DEBUGLEVEL_ell >= 4, print("listpointstriv = ",listpointstriv));
1348
1349  found = 0;
1350  listgen = KS2gen;
1351  listpointscount = [];
1352
1353  m1 = m2 = 0; lastloc = -1;
1354
1355  mask = 1 << #KS2gen;
1356  i = 1;
1357  while( i < mask,
1358    d1 = 1; iaux = i; j = 1;
1359    while( iaux,
1360      if( iaux%2, d1 *= listgen[j]);
1361      iaux >>= 1; j++);
1362if( DEBUGLEVEL_ell >= 2, print("d1 = ",d1));
1363    triv = 0;
1364    for( j = 1, #listpointstriv,
1365      if( listpointstriv[j][1] && issquare(d1*listpointstriv[j][1]),
1366        listpointscount = concat(listpointscount,[listpointstriv[j]]);
1367if( DEBUGLEVEL_ell >= 2, print("point trivial"));
1368        triv = 1; m1++;
1369        if( degre(i) > lastloc, m2++);
1370        found = 1; lastloc = -1; break));
1371    if( !triv,
1372    pol = Pol([d1,0,c,0,d/d1]);
1373if( DEBUGLEVEL_ell >= 3, print("quartique = y^2 = ",pol));
1374    point = ratpoint(pol,LIM1,1);
1375    if( point != [],
1376if( DEBUGLEVEL_ell >= 2, print("point sur la quartique"));
1377if( DEBUGLEVEL_ell >= 3, print(point));
1378      m1++;
1379      listpointscount = concat(listpointscount,[d1*point[1]*point]);
1380      if( degre(i) > lastloc, m2++);
1381      found = 1; lastloc = -1
1382    ,
1383      if( locallysoluble(pol),
1384        if( degre(i) > lastloc, m2++; lastloc = degre(i));
1385        point = ratpoint(pol,LIM3,1);
1386        if( point != [],
1387if( DEBUGLEVEL_ell >= 2, print("point sur la quartique"));
1388if( DEBUGLEVEL_ell >= 3, print(point));
1389          m1++;
1390          listpointscount = concat(listpointscount,[d1*point[1]*point]);
1391          if( degre(i) > lastloc, m2++);
1392          found = 1; lastloc = -1
1393        ,
1394if( DEBUGLEVEL_ell >= 2, print(" pas de point trouve sur la quartique"))
1395          ))));
1396    if( found,
1397      found = 0;
1398      v = 0; iaux = (i>>1);
1399      while( iaux, iaux >>= 1; v++);
1400      mask >>= 1;
1401      listgen = vecextract(listgen,(1<<#listgen)-(1<<v)-1);
1402      i = (1<<v)
1403    , i++)
1404  );
1405
1406  for( i = 1, #listpointscount,
1407   if( #listpointscount[i] > 1,
1408    if( subst(x^3+c*x^2+d*x,x,listpointscount[i][1])-listpointscount[i][2]^2 != 0,
1409      error(" ellcount : MAUVAIS POINT "))));
1410if( DEBUGLEVEL_ell >= 4, print("fin de ellcount"));
1411  return([listpointscount,[m1,m2]]);
1412}
1413if( DEBUGLEVEL_ell >= 4, print("ell2descent_viaisog"));
1414{
1415ell2descent_viaisog(ell,help=[]) =
1416\\ Calcul du rang des courbes elliptiques avec 2-torsion
1417\\ par la methode des 2-isogenies.
1418\\
1419\\ ell = [a1,a2,a3,a4,a6]
1420\\ y^2+a1xy+a3y=x^3+a2x^2+a4x+a6
1421\\
1422\\ ell doit etre sous la forme
1423\\ y^2=x^3+ax^2+bx -> ell = [0,a,0,b,0]
1424\\ avec a et b entiers.
1425local(i,P,Pfact,tors,listpointstriv,apinit,bpinit,plist,KS2prod,oddclass,KS2gen,listpoints,pointgen,n1,n2,certain,np1,np2,listpoints2,aux1,aux2,certainp,rang,strange);
1426
1427if( DEBUGLEVEL_ell >= 2, print("Algorithme de la 2-descente par isogenies"));
1428  if( #ell < 13, ell = ellinit(ell,1));
1429
1430  if( ell.disc == 0,
1431    error(" ell2descent_viaisog : courbe singuliere !!"));
1432  if( ell.a1 != 0 || ell.a3 != 0 || ell.a6 != 0,
1433    error(" ell2descent_viaisog : la courbe n'est pas sous la forme [0,a,0,b,0]"));
1434  if( denominator(ell.a2) > 1 || denominator(ell.a4) > 1,
1435    error(" ell2descent_viaisog : coefficients non entiers"));
1436
1437  P = Pol([1,ell.a2,ell.a4]);
1438  Pfact = factor(P)[,1];
1439  tors = #Pfact;
1440  listpointstriv = concat(help,elltorseven(ell)[3]);
1441 
1442  apinit = -2*ell.a2; bpinit = ell.a2^2-4*ell.a4;
1443
1444  plist = factor(6*ell.disc)[,1];
1445
1446if( DEBUGLEVEL_ell >= 3, print(" Recherche de points triviaux sur la courbe"));
1447  P *= x;
1448if( DEBUGLEVEL_ell >= 3, print("Y^2 = ",P));
1449  listpointstriv = concat( listpointstriv, ratpoint(P,LIMTRIV,0));
1450if( DEBUGLEVEL_ell >= 1, print("points triviaux sur E(Q) = ",listpointstriv); print());
1451
1452  KS2prod = -abs(ell.a4);
1453  if( bpinit < 0, KS2prod = -KS2prod);
1454  KS2gen = factor(KS2prod)[,1];
1455
1456if( DEBUGLEVEL_ell >= 2,
1457  print("#K(b,2)gen          = ",#KS2gen);
1458  print("K(b,2)gen = ",KS2gen));
1459
1460  listpoints = ellcount(ell.a2,ell.a4,KS2gen,listpointstriv);
1461  pointgen = listpoints[1];
1462if( DEBUGLEVEL_ell >= 1, print("points sur E(Q) = ",pointgen); print());
1463  n1 = listpoints[2][1]; n2 = listpoints[2][2];
1464
1465  certain = (n1 == n2);
1466if( DEBUGLEVEL_ell >= 1,
1467  if( certain,
1468    print("[E(Q):phi'(E'(Q))]  = ",1<<n1);
1469    print("#S^(phi')(E'/Q)     = ",1<<n2);
1470    print("#III(E'/Q)[phi']    = 1"); print()
1471  ,
1472    print("[E(Q):phi'(E'(Q))] >= ",1<<n1);
1473    print("#S^(phi')(E'/Q)     = ",1<<n2);
1474    print("#III(E'/Q)[phi']   <= ",1<<(n2-n1)); print())
1475);
1476
1477  KS2prod = -abs(bpinit);
1478  if( ell.a4 < 0, KS2prod = -KS2prod);
1479  KS2gen = factor(KS2prod)[,1];
1480
1481if( DEBUGLEVEL_ell >= 2, 
1482  print("#K(a^2-4b,2)gen     = ",#KS2gen);
1483  print("K(a^2-4b,2)gen     = ",KS2gen));
1484
1485  P = Pol([1,apinit,bpinit]);
1486  listpointstriv = elltorseven([0,apinit,0,bpinit,0])[3];
1487   
1488if( DEBUGLEVEL_ell >= 3, print(" Recherche de points triviaux sur la courbe"));
1489  P *= x;
1490if( DEBUGLEVEL_ell >= 3, print("Y^2 = ",P));
1491  listpointstriv = concat( listpointstriv, ratpoint(P,LIMTRIV,0));
1492if( DEBUGLEVEL_ell >= 1, print("points triviaux sur E'(Q) = ",listpointstriv); print());
1493
1494  listpoints = ellcount(apinit,bpinit,KS2gen,listpointstriv);
1495if( DEBUGLEVEL_ell >= 1, print("points sur E'(Q) = ",listpoints[1]));
1496  np1 = listpoints[2][1]; np2 = listpoints[2][2];
1497  listpoints2 = vector(#listpoints[1],i,0);
1498  for( i = 1, #listpoints[1],
1499    listpoints2[i] = [0,0];
1500    aux1 = listpoints[1][i][1]^2;
1501    if( aux1 != 0,
1502      aux2 = listpoints[1][i][2];
1503      listpoints2[i][1] = aux2^2/aux1/4;
1504      listpoints2[i][2] = aux2*(bpinit-aux1)/aux1/8
1505    , listpoints2[i] = listpoints[1][i]));
1506if( DEBUGLEVEL_ell >= 1, print("points sur E(Q) = ",listpoints2); print());
1507  pointgen = concat(pointgen,listpoints2);
1508
1509  certainp = (np1 == np2);
1510if( DEBUGLEVEL_ell >= 1,
1511  if( certainp,
1512    print("[E'(Q):phi(E(Q))]   = ",1<<np1);
1513    print("#S^(phi)(E/Q)       = ",1<<np2);
1514    print("#III(E/Q)[phi]      = 1"); print()
1515  ,
1516    print("[E'(Q):phi(E(Q))]  >= ",1<<np1);
1517    print("#S^(phi)(E/Q)       = ",1<<np2);
1518    print("#III(E/Q)[phi]     <= ",1<<(np2-np1)); print());
1519
1520  if( !certain && (np2 > np1), print1(1<<(np2-np1)," <= "));
1521  print1("#III(E/Q)[2]       ");
1522  if( certain && certainp, print1(" "), print1("<"));
1523  print("= ",1<<(n2+np2-n1-np1));
1524
1525  print("#E(Q)[2]            = ",1<<tors);
1526);
1527  rang = n1+np1-2;
1528if( DEBUGLEVEL_ell >= 1,
1529  if( certain && certainp,
1530    print("#E(Q)/2E(Q)         = ",(1<<(rang+tors)));
1531    print("rang                = ",rang); print()
1532  ,
1533    print("#E(Q)/2E(Q)        >= ",(1<<(rang+tors))); print();
1534    print(rang," <= rang          <= ",n2+np2-2); print()
1535  ));
1536
1537  strange = (n2+np2-n1-np1)%2;
1538  if( strange,
1539if( DEBUGLEVEL_ell >= 1,
1540      print(" !!! III doit etre un carre !!!"); print("donc"));
1541    if( certain,
1542      np1++;
1543      certainp = (np1 == np2);
1544if( DEBUGLEVEL_ell >= 1,
1545        if( certainp,
1546          print("[E'(Q):phi(E(Q))]   = ",1<<np1);
1547          print("#S^(phi)(E/Q)       = ",1<<np2);
1548          print("#III(E/Q)[phi]      = 1"); print()
1549        ,
1550          print("[E'(Q):phi(E(Q))]  >= ",1<<np1);
1551          print("#S^(phi)(E/Q)       = ",1<<np2);
1552          print("#III(E/Q)[phi]     <= ",1<<(np2-np1)); print())
1553      )
1554    ,
1555    if( certainp,
1556      n1++;
1557      certain = (n1 == n2);
1558if( DEBUGLEVEL_ell >= 1,
1559        if( certain,
1560          print("[E(Q):phi'(E'(Q))]   = ",1<<n1);
1561          print("#S^(phi')(E'/Q)       = ",1<<n2);
1562          print("#III(E'/Q)[phi']      = 1"); print()
1563        ,
1564          print("[E(Q):phi'(E'(Q))]  >= ",1<<n1);
1565          print("#S^(phi')(E'/Q)      = ",1<<n2);
1566          print("#III(E'/Q)[phi']    <= ",1<<(n2-n1)); print())
1567      )
1568    , n1++)
1569  );
1570
1571if( DEBUGLEVEL_ell >= 1,
1572  print("#S^(2)(E/Q)           = ",1<<(n2+np2-1));
1573  if( !certain && (np2 > np1), print1(1<<(np2-np1)," <= "));
1574  print1("#III(E/Q)[2]       ");
1575  if( certain && certainp, print1(" "), print1("<"));
1576  print("= ",1<<(n2+np2-n1-np1));
1577  print("#E(Q)[2]            = ",1<<tors);
1578);
1579  rang = n1+np1-2;
1580if( DEBUGLEVEL_ell >= 1,
1581  if( certain && certainp,
1582    print("#E(Q)/2E(Q)         = ",(1<<(rang+tors))); print();
1583    print("rang                = ",rang); print()
1584  ,
1585    print("#E(Q)/2E(Q)        >= ",(1<<(rang+tors))); print();
1586    print(rang," <= rang          <= ",n2+np2-2); print())
1587  ));
1588
1589\\ fin de strange
1590 
1591  if( ELLREDGENFLAG, pointgen = ellredgen(ell,pointgen));
1592  pointgen = concat(ellsort(elltorseven(ell)[3]),pointgen);
1593if( DEBUGLEVEL_ell >= 1, print("points = ",pointgen));
1594if( DEBUGLEVEL_ell >= 3, print("fin de ell2descent_viaisog"));
1595  return([rang,n2+np2-2+tors,pointgen]);
1596}
1597
1598