1 // mrank2.cc: implementation of class rank2 for descent via 2-isogeny
2 //////////////////////////////////////////////////////////////////////////
3 //
4 // Copyright 1990-2012 John Cremona
5 //
6 // This file is part of the eclib package.
7 //
8 // eclib is free software; you can redistribute it and/or modify it
9 // under the terms of the GNU General Public License as published by the
10 // Free Software Foundation; either version 2 of the License, or (at your
11 // option) any later version.
12 //
13 // eclib is distributed in the hope that it will be useful, but WITHOUT
14 // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
16 // for more details.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with eclib; if not, write to the Free Software Foundation,
20 // Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA
21 //
22 //////////////////////////////////////////////////////////////////////////
23
24 #include <iomanip> // for setbase(), used for hex output of codes
25
26 #include <eclib/bitspace.h>
27 #include <eclib/mlocsol.h>
28 #include <eclib/mglobsol.h>
29 #include <eclib/mrank2.h>
30 #include <eclib/sqfdiv.h>
31 #include <eclib/desc2.h>
32
33 #ifndef QSIEVE_OPT
34 #define QSIEVE_OPT 0 // uses Stoll's sieve
35 #endif
36
37 //#define DEBUG_ROOTS
38 //#define DEBUG_SQF
39 //#define DEBUG_ELS
40 //#define DEBUG_GLS
41
42 // bound on first descent search when a second descent is going to be done
43 #define FIRST_DESCENT_BOUND 8
44
45 #define MAX_R 5 // will not attempt to list all coset reps
46 // for 2E(Q) in E(Q) if rank is more than this
47
makepoint(const bigint & c,const bigint & d1,const bigint & d2,const bigint & x,const bigint & y,const bigint & z,int which)48 void rank2::makepoint(const bigint& c,const bigint& d1,const bigint& d2,
49 const bigint& x, const bigint& y, const bigint& z,
50 int which)
51 {
52 Point P(ee);
53 if (verbose) cout<<" (x:y:z) = ("<<x<<":"<<y<<":"<<z<<")\n";
54 if(which==0)
55 {
56 if(sign(z)!=0) P.init(ee, d1*x*x*z, d1*x*y, pow(z,3));
57 if(verbose)
58 {
59 cout<<"\tCurve E " /* <<(Curve)ee */ <<"\tPoint "<< P;
60 bigfloat h = height(P);
61 cout << ", height = " << h;
62 if(!P.isvalid()) {cout << " --warning: NOT on curve! " ;}
63 cout << endl;
64 }
65 }
66 else
67 {
68 if(verbose) // display point on isogenous curve too
69 {
70 Point Q(eedash);
71 if(sign(z)!=0) Q.init(eedash, d1*x*x*z, d1*x*y, pow(z,3));
72 cout<<"\tCurve E' " /* <<(Curve)eedash */ <<"\tPoint "<< Q;
73 bigfloat h = height(Q);
74 cout << ", height = " << h;
75 if(!Q.isvalid()) {cout << " --warning: NOT on curve! " ;}
76 cout << endl;
77 }
78 const bigint& xz=x*z, x2=x*x, z2=z*z;
79 if(sign(xz)!=0) P.init(ee,2*y*y*xz,y*(d1*x2*x2-d2*z2*z2),pow(2*xz,3));
80 if(verbose)
81 {
82 cout<<"\tCurve E " /* <<(Curve)ee */ <<"\tPoint "<< P;
83 bigfloat h = height(P);
84 cout << ", height = " << h;
85 if(!P.isvalid()) {cout << " --warning: NOT on curve! " ;}
86 cout << endl;
87 }
88 }
89 if(order(P)<0) {pointlist.push_back(P); npoints++;} // else torsion so ignore
90 }
91
testquartic(const bigint & c,const bigint & d1,const bigint & d2,int which)92 int rank2::testquartic(const bigint& c,const bigint& d1,const bigint& d2,int which)
93 {
94 // creates quartic (d1,0,c,0,d2), assumed els, and tries to find a rational point
95 // returns +1 if rational point found (handled my makepoint())
96 // 0 if undecided
97 static const bigint zero = BIGINT(0);
98 static const bigint one = BIGINT(1);
99 quartic q(d1, zero, c, zero, d2);
100 if (verbose) cout<<q<<": ";
101
102 bigint x,y,z;
103
104 // First a quick search for a small point:
105 if (ratpoint(q,one,BIGINT(lim1),x,y,z))
106 {
107 makepoint(c,d1,d2,x,y,z,which);
108 return 1;
109 }
110
111 // Next a more serious search using a sieve,
112 // but not too far if we are going to use second descent
113 quartic_sieve qs(&q,QSIEVE_OPT,0);
114 long lim3=lim2;
115 if(do_second_descent) // save the long search till later...
116 if(lim3>FIRST_DESCENT_BOUND) lim3=FIRST_DESCENT_BOUND;
117 if(qs.search(lim3, 1, 1)) // maxnpoints, pos-x-only
118 {
119 qs.getpoint(x,y,z);
120 makepoint(c,d1,d2,x,y,z,which);
121 return 1;
122 }
123 if (verbose)
124 cout<<" no rational point found (hlim="<<lim3<<")\n";
125 return 0;
126 } /* end of testquartic*/
127
second_descent(const bigint & c,const bigint & d1,const bigint & d2,int which)128 int rank2::second_descent(const bigint& c, const bigint& d1, const bigint& d2, int which)
129 // creates quartic (d1,0,c,0,d2), assumed els, and tries to find a rational point
130 // by second descent (after testquartic() has failed using direct search)
131 // returns +1 if rational point found (handled by makepoint())
132 // 0 if undecided
133 // -1 if no rational point exists (proved by second descent)
134 {
135 int res, verb=verbose; if(verb) verb--; // reduced verbosity level within desc2()
136 bigint x,y,z;
137 if(verbose) cout<<"d1="<<d1<<": "<<flush;
138 if(which)
139 res = desc2(c,d1,d2,badprimes, supp0, elsgens0, mask0, lim2,x,y,z,verb);
140 else
141 res = desc2(c,d1,d2,badprimes, supp1, elsgens1, mask1, lim2,x,y,z,verb);
142
143 if(verbose==1) cout<<endl;
144
145 switch(res)
146 {
147 case 1: // Positive success, point found
148 makepoint(c,d1,d2,x,y,z,which);
149 if(verbose)
150 cout<<"Second descent successfully found rational point for d1="<<d1<<"\n"<<endl;
151 break;
152 case -1: // Negative success, no point exists
153 if(verbose)
154 cout<<"Second descent shows that no rational point exists for d1="<<d1<<"\n"<<endl;
155 break;
156 case 0: // Undecided
157 if (verbose)
158 cout<<"Second descent inconclusive for d1="<<d1<<": ELS descendents exist but no rational point found\n"<<endl;
159 }
160 return res;
161 }
162
163 // First local descent: determine
164 // S^(phi)(E') if which==0
165 // S^(phi')(E) if which==1
166
find_elsgens(int which,const bigint & c,const bigint & d)167 void rank2::find_elsgens(int which, const bigint& c, const bigint& d)
168 {
169 static const bigint zero = BIGINT(0);
170 if (verbose>1)
171 {
172 if(which) cout<<"\n";
173 cout<<"Finding els gens for E"; if(which) cout<<"'";
174 cout<<" (c"; if(which) cout<<"'";
175 cout<<"= "<<c<<", d"; if(which) cout<<"'";
176 cout<<"= "<<d<<")"<<endl;
177 }
178 bigint ddash = c*c-4*d;
179 int posd1, istep=1;
180 if (sign(ddash)< 0)
181 posd1=1;
182 else
183 posd1= (c+sqrt(ddash)< 0);
184 if(posd1) istep=2; // this will skip negative d1
185
186 int extra2torsion=0;
187 bigint ee2, ee3;
188 if(which)
189 {
190 if(d_is_sq)
191 {
192 extra2torsion=1;
193 ee2 = e2dash;
194 ee3 = e3dash;
195 }
196 }
197 else
198 {
199 if(ddash_is_sq)
200 {
201 extra2torsion=1;
202 ee2 = e2;
203 ee3 = e3;
204 }
205 }
206
207 bigint d1, d2, badp;
208 const vector<bigint>& supp = (which? supp1 : supp0);
209 long ns = supp.size();
210 if(verbose>1)
211 {
212 cout<<"Support (length "<<ns<<"): "<<supp<<endl;
213 }
214 if(ns>=NTL_BITS_PER_LONG)
215 {
216 cerr<<"Too many primes dividing discriminant!\n";
217 cerr<<"mwrank cannot handle supports of more than "<<NTL_BITS_PER_LONG<< "primes"<<endl;
218 cerr<<"Failure in find_elsgens"<<endl;
219 success = 0;
220 return;
221 }
222
223 long mask=0, maxn=1<<ns, index, j, nelsgens=0;
224 vector<bigint> elsgens;
225
226 // use all torsion: added 24/6/02
227 // Find and process torsion
228
229 vector<Point> torsion;
230 if(which)
231 torsion=torsion_points(eedash);
232 else
233 torsion=torsion_points(ee);
234 long it,ntorsion=torsion.size();
235 #ifdef DEBUG_ELS
236 cout<<"Number of torsion points = "<<ntorsion<<endl;
237 cout<<"Torsion subgroup is ";
238 if(extra2torsion) cout<<"not ";
239 cout<<"cyclic"<<endl;
240 #endif
241
242 // find generators of torsion mod 2*torsion
243
244 Point T1,T2; long orderT1=1;
245 for(it=0; it<ntorsion; it++)
246 {
247 Point T = torsion[it];
248 long orderT = order(T);
249 if(orderT>orderT1) {T1=T; orderT1=orderT;}
250 }
251 #ifdef DEBUG_ELS
252 cout<<"Generator of maximal order = "<<T1<<" of order "<<orderT1<<endl;
253 #endif
254 if(extra2torsion)
255 // then we need a second generator, for which we can take either of
256 // the 2-torsion points which is not the multiple of T1
257 {
258 T2=(ntorsion/4)*T1;
259 bigint xT2=T2.getX();
260 if(xT2==zero)
261 T2.init(ee,ee2,zero);
262 else
263 T2.init(ee,zero,zero);
264 #ifdef DEBUG_ELS
265 cout<<"Second generator = "<<T2<<" of order 2"<<endl;
266 #endif
267 ntorsion=2;
268 torsion[0]=T1;
269 torsion[1]=T2;
270 }
271 else
272 {
273 ntorsion=1;
274 torsion[0]=T1;
275 }
276 #ifdef DEBUG_ELS
277 cout<<"Processing generating torsion points"<<endl;
278 #endif
279
280 bigint d1x;
281 for(it=0; it<ntorsion; it++)
282 {
283 Point T = torsion[it];
284 if(T.is_zero()) continue;
285 d1=T.getX();
286 if(d1==zero) d1=d;
287 #ifdef DEBUG_ELS
288 if(verbose>1) cout<<"Processing torsion d1 = " << d1 << endl;
289 #endif
290 index=makeindex(supp, d1, d1x);
291 d1=d1x; // the square-free part, = d1/square
292 #ifdef DEBUG_ELS
293 if(verbose>1) cout<<"Squarefree part of d1 = " << d1 << endl;
294 #endif
295 if(mask&index) continue;
296 if(index)
297 {
298 elsgens.push_back(d1); nelsgens++;
299 if(verbose>1)
300 cout<<"Adding (torsion) els generator #"<<nelsgens<<": d1 = " << d1;
301 for(j=ns-1; j>=0; j--)
302 // if(testbit(index,j)) {setbit(mask,j); break;}
303 if(index&(1<<j)) {mask|=(1<<j); break;}
304 #ifdef DEBUG_ELS
305 if(verbose>1) cout<<" (pivot = "<<j<<": "<<supp[j]<<")";
306 #endif
307 if(verbose>1) cout<<endl;
308 }
309 }
310
311 // Record that the first 0, 1 or 2 elsgens come from known torsion points,
312 // so can be ingored in the gls step later:
313 if(which)
314 {
315 nt2gens1=nelsgens;
316 }
317 else
318 {
319 nt2gens0=nelsgens;
320 }
321
322 // Now systematically go through square-free divisors d1 of d, always modulo
323 // the subgroup of those which we know are els
324
325 int res;
326 for(index=istep; index<maxn; index+=istep)
327 {
328 if(mask&index) continue;
329 d1 = makenum(supp,index);
330 d2 = d/d1;
331 #ifdef DEBUG_ELS
332 if(verbose>1) cout<<"Testing d1 = "<<d1<<":\t";
333 #endif
334 res = locallysoluble(d1,c,d2,badprimes,badp);
335 if(res)
336 {
337 elsgens.push_back(d1); nelsgens++;
338 if(verbose>1)
339 cout<<"Adding els generator #"<<nelsgens<<": d1 = " << d1;
340 for(j=ns-1; j>=0; j--)
341 // if(testbit(index,j)) {setbit(mask,j); break;}
342 if(index&(1<<j)) {mask|=(1<<j); break;}
343 #ifdef DEBUG_ELS
344 if(verbose>1) cout<<" (pivot = "<<j<<": "<<supp[j]<<")";
345 #endif
346 if(verbose>1) cout<<endl;
347 }
348 else
349 {
350 #ifdef DEBUG_ELS
351 if(verbose>1) cout<<"not locally soluble at p = "<<badp<<endl;
352 #endif
353 }
354 }
355 if(verbose>1)
356 {
357 cout<<"After els sieving, nelsgens = " << nelsgens;
358 cout << endl;
359
360 cout<<"2-rank of S^{phi"; if(which) cout<<"'";
361 cout<<"}(E"; if(!which) cout<<"'";
362 cout<<") = "<<nelsgens<<endl;
363
364 if(nelsgens>0) cout<<"(els)gens: "<<elsgens<<endl;
365 // cout<<"mask = "<<mask<<endl;
366 }
367 if(which)
368 {
369 elsgens1=elsgens; mask1=mask; els1=nelsgens;
370 }
371 else
372 {
373 elsgens0=elsgens; mask0=mask; els0=nelsgens;
374 }
375 }
376
377 // Second local descent: determine
378 // phi'(S^2(E)) in S^(phi)(E') if which==0
379 // phi(S^2(E')) in S^(phi')(E) if which==1
380
find_els2gens(int which,const bigint & c,const bigint & d)381 void rank2::find_els2gens(int which, const bigint& c, const bigint& d)
382 {
383 if (verbose>1)
384 {
385 if(which) cout<<"\n";
386 cout<<"Finding els2 gens for E"; if(which) cout<<"'";
387 cout<<" (c"; if(which) cout<<"'";
388 cout<<"= "<<c<<", d"; if(which) cout<<"'";
389 cout<<"= "<<d<<") which lift to S^2(E";
390 if(which) cout<<"'"; cout<<")"<<endl;
391 }
392
393 vector<bigint>& elsgens = (which? elsgens1: elsgens0);
394 long nelsgens = (which? els1: els0);
395 long nt2gens = (which? nt2gens1: nt2gens0);
396
397 bigint d1, d2, badp, x,y,z;
398 unsigned long els2mask; long index;
399 long maxn=1<<nelsgens, nels2gens=0, els2piv;
400 vector<bigint> els2gens;
401 bitspace els2_space(nelsgens);
402
403 // first record the torsion contribution:
404
405 for(index=0; index<nt2gens; index++)
406 {
407 #ifdef DEBUG_ELS
408 d1=elsgens[index];
409 if(verbose>1) cout<<"Processing torsion d1 = " << d1 << ":"<<endl;
410 #endif
411 els2mask=(1<<index);
412 if(els2_space.mask(els2mask)) continue; // we work mod the els2 subgp
413 els2piv=els2_space.reduce(els2mask);
414 if(els2piv<0) continue; // means we're in the subgp; won't happen
415 #ifndef DEBUG_ELS // else done earlier
416 d1=elsgens[index];
417 #endif
418 els2gens.push_back(d1); nels2gens++;
419 els2_space.augment(els2mask,els2piv);
420 if(verbose>1) cout<<"Adding (torsion) els2 generator #"<<(nels2gens)
421 <<": d1 = " << d1 <<endl;
422 #ifdef DEBUG_ELS
423 else
424 cout<<"Just incremented nels2gens to "<<nels2gens<<endl;
425 cout<<"now bitmask = "<<els2_space.getbitmask()<<endl;
426 #endif
427 } // end of torsion loop
428
429 int res, verb=0; if(verbose>2)verb=verbose-2;
430 // reduced verbosity level within desc2()
431
432 for(index=1; (index<maxn)&&(nels2gens<nelsgens); index++)
433 {
434 // First mask against known els2 subgroup:
435 if(els2_space.mask(index)) continue; // we work mod the els2 subgp
436 els2mask=index;
437 els2piv=els2_space.reduce(els2mask);
438 #ifdef DEBUG_ELS
439 cout<<"index = "<<index<<endl;
440 cout<<"els2mask = "<<els2mask<<endl;
441 cout<<"els2piv = "<<els2piv <<endl;
442 #endif
443 if(els2piv<0) continue; // means w're in the els2 subgp; won't happen
444 d1 = makenum(elsgens,index);
445 d2 = d/d1;
446 #ifdef DEBUG_ELS
447 if(verbose>1) cout<<"Processing d1 = " << d1 << ":\t";
448 #endif
449 if(which)
450 res = 1+desc2(c,d1,d2,badprimes,supp0,elsgens0,mask0,lim2,x,y,z,verb,1);
451 else
452 res = 1+desc2(c,d1,d2,badprimes,supp1,elsgens1,mask1,lim2,x,y,z,verb,1);
453 #ifdef DEBUG_ELS
454 if(verbose>1) cout<<"res = " << res << endl;
455 #endif
456 if(res)
457 {
458 els2gens.push_back(d1); nels2gens++;
459 els2_space.augment(els2mask,els2piv);
460 #ifdef DEBUG_ELS
461 cout<<"now bitmask = "<<els2_space.getbitmask()<<endl;
462 #endif
463 if(verbose>1)
464 cout<<"Adding els2 generator #"<<nels2gens<<": d1 = " << d1<<endl;
465 }
466 }
467
468 if(verbose>1)
469 {
470 cout<<"After els2 sieving, nels2gens = " << nels2gens;
471 cout << endl;
472
473 cout<<"2-rank of phi";if(!which) cout<<"'";
474 cout<<"(S^2(E"; if(which) cout<<"'";
475 cout<<")) = "<<nels2gens<<endl;
476
477 if(nels2gens>0) cout<<"(els2)gens: "<<els2gens<<endl;
478 }
479 if(which)
480 {
481 els2gens1=els2gens; els21=nels2gens;
482 }
483 else
484 {
485 els2gens0=els2gens; els20=nels2gens;
486 }
487 }
488
find_glsgens(int which,const bigint & c,const bigint & d)489 void rank2::find_glsgens(int which, const bigint& c, const bigint& d)
490 {
491 vector<bigint>& elsgens = (which? els2gens1: els2gens0);
492 long nelsgens = (which? els21: els20);
493 long nt2gens = (which? nt2gens1: nt2gens0);
494 vector<bigint> gls_gens;
495 bitspace gls_space(nelsgens);
496 long glspiv, maxn = 1<<nelsgens;
497 long nglsgens=0;
498 unsigned long glsmask;
499 long index, stage, nstages=1; if(do_second_descent) nstages=2;
500 long shortfall1, shortfall2;
501 bigint d1, d2;
502 int res;
503
504 // first record the torsion contribution:
505
506 for(index=0; index<nt2gens; index++)
507 {
508 glsmask=(1<<index);
509 d1=elsgens[index];
510 glspiv=gls_space.reduce(glsmask);
511 if(glspiv<0)
512 {
513 #ifdef DEBUG_GLS
514 cout<<"d1="<<d1<<": known gls (torsion)"<<endl;
515 #endif
516 continue; // as we are certainly in the gls subgroup
517 }
518 if(verbose>1)
519 {
520 cout<<"Adding (torsion) gls generator #"<<(nglsgens+1)
521 <<": d1 = " << d1 <<endl;
522 }
523 gls_gens.push_back(d1); nglsgens++;
524 gls_space.augment(glsmask,glspiv);
525 #ifdef DEBUG_GLS
526 cout<<"Just incremented nglsgens to "<<nglsgens<<endl;
527 #endif
528 } // end of torsion loop
529
530 // Next test all first descent curves for global solubility (stage 1)
531 // and (optionally) do a second descent (stage 2) on unsuccessful ones.
532
533 // The next lines are in case torsion accounts for everything,
534 // as then the following loop is not executed at all:
535 shortfall1 = 0; shortfall2 = 0;
536 if(which)
537 {gls1=gls21=nglsgens;}
538 else
539 {gls0=gls20=nglsgens;}
540
541 if((nglsgens==nelsgens)&&verbose)
542 cout<<"This component of the rank is 0\n";
543
544 for(stage=1; (stage<=nstages)&&(nglsgens<nelsgens); stage++)
545 {
546 if(verbose&&do_second_descent)
547 {
548 if(stage==1)
549 cout << "First stage (no second descent yet)..."<<endl;
550 else
551 cout << "Second stage (using second descent)..."<<endl;
552 }
553 for(index=1; (index<maxn)&&(nglsgens<nelsgens); index++)
554 {
555 #ifdef DEBUG_GLS
556 d1 = makenum(elsgens,index); // do it now for output purposes
557 cout<<"index="<<setbase(16)<<index<<setbase(10);
558 cout<<"\td1 = "<<d1<<":\t";
559 #endif
560 // First mask against gls subgroup:
561 glsmask=index;
562 glspiv=gls_space.reduce(glsmask);
563 if(glspiv<0)
564 {
565 #ifdef DEBUG_GLS
566 cout<<"known to be gls"<<endl;
567 #endif
568 continue; // as we are certainly in the gls subgroup
569 }
570 #ifdef DEBUG_GLS
571 cout<<"not known gls"
572 <<", mask="<<setbase(16)<<glsmask<<setbase(10)
573 <<endl;
574 cout<<"Keeping j="<<glspiv
575 <<" with mask "<<setbase(16)<<glsmask<<setbase(10)
576 <<" for possible use as pivot"<<endl;
577 #else
578 d1 = makenum(elsgens,index); // else we did it earlier
579 #endif
580 d2 = d/d1;
581 if(stage==1) res = testquartic(c,d1,d2,which);
582 else res = second_descent(c,d1,d2,which);
583 #ifdef DEBUG_GLS
584 cout<<"result = "<<res<<endl;
585 #endif
586 switch(res){
587 case -1: // should not happen
588 {
589 cout<<"Problem in 2nd descent!" <<endl;
590 success=0;
591 return;
592 }
593 case 0: // no point found, nothing to do
594 {
595 break;
596 }
597 case 1:
598 {
599 gls_space.augment(glsmask, glspiv);
600 gls_gens.push_back(d1);
601 nglsgens++;
602 #ifdef DEBUG_GLS
603 cout<<"Just incremented nglsgens to "<<nglsgens<<endl;
604 #endif
605 if(verbose>1)
606 {
607 cout<<"Adding gls generator #"<<(nglsgens)
608 <<": d1 = " << d1;
609 #ifdef DEBUG_GLS
610 cout<<" (g_pos = "<<glspiv<<", g_gen = "
611 <<setbase(16)<<glsmask<<setbase(10)<<")";
612 #endif
613 cout<<endl;
614 }
615 break;
616 }
617 } // end of switch(res)
618 } // end of index loop
619 #ifdef DEBUG_GLS
620 cout<<"At bottom of index loop\n";
621 cout<<"nglsgens = "<<nglsgens<<"\tgls_gens = ";
622 cout << gls_gens << endl;
623 #endif
624 if(stage==1) // finished stage 1, i.e. first descent
625 {
626 shortfall1 = nelsgens-nglsgens;
627 if(verbose)
628 {
629 cout << "After first global descent, this component of the rank";
630 if(shortfall1)
631 {
632 cout << "\n\thas lower bound " << nglsgens-nt2gens
633 << "\n\tand upper bound " << nelsgens-nt2gens
634 << "\n\t(difference = " << shortfall1 << ")\n";
635 if(nstages==2)
636 cout<<"Second descent will attempt to reduce this"<<endl;
637 }
638 else
639 cout << " = " << nelsgens <<endl;
640 if(verbose>1) cout<<endl;
641 }
642 if(which)
643 gls1=gls21=nglsgens;
644 else
645 gls0=gls20=nglsgens;
646 }
647 else // finished stage==2, i.e. 2nd descent
648 {
649 shortfall2 = nelsgens-nglsgens;
650 if(verbose)
651 {
652 cout << "After second global descent, this component of the rank";
653 if(shortfall2)
654 {
655 cout << "\n\thas lower bound " << nglsgens-nt2gens
656 << "\n\tand upper bound " << nelsgens-nt2gens
657 << "\n\t(difference = " << shortfall2 << ")\n";
658 }
659 else
660 cout << " = " << nelsgens <<endl;
661 if(verbose>1) cout<<endl;
662 }
663 if(which)
664 gls21=nglsgens;
665 else
666 gls20=nglsgens;
667 }
668
669 } // end of stage loop
670
671 if(which)
672 glsgens1=gls_gens;
673 else
674 glsgens0=gls_gens;
675
676 if(verbose>1)
677 {
678 cout<<"After gls sieving, nglsgens = " << nglsgens << endl;
679 cout<< "shortfall in rank from first descent = "<<shortfall1<<endl;
680 if(do_second_descent)
681 cout<< "shortfall in rank from second descent = "<<shortfall2<<endl;
682 if(nglsgens>0) cout<<"(gls)gens: "<<gls_gens<<endl;
683 }
684 }
685
local_descent(const bigint & x0)686 void rank2::local_descent(const bigint& x0)
687 {
688 bigint c,d,cdash,ddash,disc,rootd;
689 const bigint zero = BIGINT(0), two=BIGINT(2), minusone=BIGINT(-1);
690
691 c = 3 * x0 + s2;
692 d = x0*(c+s2) + s4;
693 cdash = - 2 * c;
694 ddash = c*c - 4*d;
695 disc = 2*d*ddash;
696 if (is_zero(disc)) // this should have been caught by the calling program!
697 {
698 cerr << "Curve is singular"<<endl;
699 success = 0;
700 return;
701 }
702
703 if(verbose>1)
704 {
705 cout<<"(c,d) =("<<c<<","<<d<<")"<<endl;
706 cout<<"(c',d')=("<<cdash<<","<<ddash<<")"<<endl;
707 }
708
709 ee = Curvedata(zero, c, zero, d, zero);
710 eedash = Curvedata(zero,cdash,zero,ddash,zero);
711 Eprime = Curvedata(eedash,1); // minimal model
712 if(verbose)
713 {
714 cout<<"Using 2-isogenous curve "<<(Curve)(eedash);
715 if (eedash!=Eprime)
716 cout<<" (minimal model "<<(Curve)(Eprime)<<")";
717 cout<<endl;
718 }
719 supp0 = pdivs(d);
720 supp1 = pdivs(ddash);
721 badprimes.clear();
722 set_union(supp0.begin(),supp0.end(),supp1.begin(),supp1.end(),back_inserter(badprimes));
723 if(find(badprimes.begin(),badprimes.end(),two)==badprimes.end())
724 badprimes.insert(badprimes.begin(),two);
725 // cout<<"supp0="<<supp0<<", supp1="<<supp1<<endl;
726 // cout<<"badprimes="<<badprimes<<endl;
727 // NB -1 MUST be the first entry in the supports!
728 supp0.insert(supp0.begin(),minusone);
729 supp1.insert(supp1.begin(),minusone);
730 // cout<<"supp0="<<supp0<<", supp1="<<supp1<<endl;
731
732 d_is_sq = ddash_is_sq = 0;
733 if(isqrt(ddash,rootd))
734 {
735 ddash_is_sq=1;
736 e2 = (-c-rootd)/2;
737 e3 = (-c+rootd)/2;
738 }
739 if(isqrt(d,rootd))
740 {
741 d_is_sq=1;
742 e2dash = c-2*rootd;
743 e3dash = c+2*rootd;
744 }
745
746 if(verbose)
747 {
748 cout<<"-------------------------------------------------------\n";
749 cout<<"First step, determining 1st descent Selmer groups\n";
750 cout<<"-------------------------------------------------------\n";
751 }
752 // cout<<"Calling find_elsgens(0,...)"<<endl;
753 find_elsgens(0,c,d);
754 if (!success) return;
755 // cout<<"Calling find_elsgens(1,...)"<<endl;
756 find_elsgens(1,cdash,ddash);
757 if (!success) return;
758 selmer_rank_phi_Eprime = els0;
759 selmer_rank_phiprime_E = els1;
760 rank_bound = els0+els1-2;
761
762 if(verbose)
763 {
764 cout<<"After first local descent, rank bound = "<<rank_bound<<"\n";
765 cout<<"rk(S^{phi}(E'))=\t"<<selmer_rank_phi_Eprime<<endl;
766 cout<<"rk(S^{phi'}(E))=\t"<<selmer_rank_phiprime_E<<endl;
767 cout<<endl;
768 cout<<"-------------------------------------------------------\n";
769 cout<<"Second step, determining 2nd descent Selmer groups\n";
770 cout<<"-------------------------------------------------------\n";
771 }
772
773 if((rank_bound==0)||(!do_second_descent))
774 {
775 els20=els0; els21=els1;
776 els2gens0=elsgens0; els2gens1=elsgens1;
777 if(verbose)
778 {
779 if(do_second_descent)
780 cout<<"...skipping since we already know rank=0"<<endl;
781 else
782 cout<<"...skipping -- no second descent requested"<<endl;
783 }
784 }
785 else
786 {
787 find_els2gens(0,c,d);
788 find_els2gens(1,cdash,ddash);
789 rank_bound = els20+els21-2;
790 }
791 // NB these are only upper bounds if !do_second_descent
792 selmer_rank = els1+els20-1+ddash_is_sq;
793 selmer_rank_Eprime = els0+els21-1+d_is_sq;
794 if(verbose)
795 {
796 if(do_second_descent)
797 {
798 cout<<"After second local descent, rank bound = "<<rank_bound<<"\n";
799 cout<<"rk(phi'(S^{2}(E)))=\t"<<els20<<endl;
800 cout<<"rk(phi(S^{2}(E')))=\t"<<els21<<endl;
801 cout<<"rk(S^{2}(E))=\t"<<selmer_rank<<endl;
802 cout<<"rk(S^{2}(E'))=\t"<<selmer_rank_Eprime<<endl;
803 cout<<endl;
804 }
805 else
806 {
807 cout<<"No second local descent done, rank bound = "<<rank_bound<<"\n";
808 cout<<"rk(phi'(S^{2}(E)))<=\t"<<els20<<endl;
809 cout<<"rk(phi(S^{2}(E')))<=\t"<<els21<<endl;
810 cout<<"rk(S^{2}(E))<=\t"<<selmer_rank<<endl;
811 cout<<"rk(S^{2}(E'))<=\t"<<selmer_rank_Eprime<<endl;
812 cout<<endl;
813 }
814 }
815 }
816
rank2(Curvedata * ec,int verb,int sel,long l1,long l2,int second)817 rank2::rank2(Curvedata* ec, int verb, int sel, long l1, long l2, int second)
818 : rank12(ec,verb,sel,l1,l2,0,second)
819 {
820 bigint a1, a2, a3, a4, a6;
821 ec->getai(a1,a2,a3,a4,a6);
822 fullnpoints = npoints = 0;
823 rank = 0; // default value if failure occurs
824 int best_isogeny=0, best_rank_bound=999999;
825
826 success = 1;
827 int n, scaled=0;
828 if (odd(a1) || odd(a3))
829 { s2= a1*a1+ 4*a2;
830 s4= 8*(a1*a3+ 2*a4);
831 s6= 16*(a3*a3+ 4*a6);
832 scaled=1;
833 }
834 else
835 { s2=a2; s4=a4; s6=a6;
836 }
837
838 vector<bigint> xlist = Introotscubic(s2,s4,s6);
839 ntwo_torsion = xlist.size();
840 if (ntwo_torsion==0)
841 {
842 success=0;
843 if (verbose) cerr << "No points of order 2, cannot create rank2 class"<<endl;
844 return;
845 }
846
847 if(verbose) cout << "\n" << ntwo_torsion << " points of order 2:\n";
848
849 // If there are 3 points of order 2, we order them (for consistency:
850 // otherwise the order can be machine dependent)
851 if(ntwo_torsion==3) sort(xlist.begin(),xlist.end());
852
853 two_torsion.resize(ntwo_torsion);
854 for(n=0; n<ntwo_torsion; n++)
855 {
856 bigint ei = xlist[n];
857 if(scaled) two_torsion[n].init(the_curve,2*ei,-a1*ei-4*a3,BIGINT(8));
858 else two_torsion[n].init(the_curve,ei,BIGINT(0),BIGINT(1));
859 if(verbose)
860 { if(n>0) cout<<", "; cout<<two_torsion[n];}
861 }
862 if(verbose)cout<<endl<<endl;
863
864 long mw0,mw1,sel0,sel1,sha0,sha1,sel20,sel21,sha20,sha21,lindex;
865
866 // we loop over the different 2-isogenies (if there are 3), doing the
867 // local descent 3 times to get the best rank bound (or stopping if
868 // that bound is 0).
869
870 // Fudge: we have to repeat the best one if it was not the last in
871 // order to get the class-global variables right for the global
872 // descent. This should obviously be fixed.
873
874 // The relevant variables are:
875 // ee,eedash,supp0,supp1,badprimes,d_is_sq,ddash_is_sq,
876 // e2,e3,e2dash,e3dash,els20,els21,elsgens0,elsgens1
877
878 for(n=0; n<ntwo_torsion; n++)
879 {
880 if(verbose&&(ntwo_torsion>1))
881 {
882 cout<<"****************************"<<endl;
883 cout<<"* Using 2-isogeny number "<<(n+1)<<" *"<<endl;
884 cout<<"****************************\n"<<endl;
885 }
886 local_descent(xlist[n]);
887 if (!success)
888 {
889 if (verbose) cout << "Failure in local descent"<<endl;
890 return;
891 }
892 if((n==0)||(rank_bound<=best_rank_bound))
893 {
894 best_rank_bound=rank_bound;
895 best_isogeny=n;
896 }
897 if(best_rank_bound==0) break;
898 }
899
900 int rerun_needed=(rank_bound>best_rank_bound);
901 rank_bound=best_rank_bound;
902
903 if(verbose&&best_isogeny>0)
904 {
905 cout<<"After second local descent, combined upper bound on rank = "
906 <<rank_bound<<"\n";
907 }
908
909 if(selmer_only) // Nothing more to do in this case
910 {
911 rank=0;
912 certain=(rank_bound==0);
913 }
914 else // do the global stages...
915 {
916
917 // Redo the best local descent if necessary
918 bigint x0=xlist[best_isogeny];
919 if(rerun_needed)
920 {
921 if(verbose) cout<<"Rerunning the local descent for isogeny number "
922 <<(best_isogeny+1)
923 << ", which gave the best upper bound on the rank"
924 <<endl;
925 local_descent(x0);
926 }
927
928 bigint c = 3 * x0 + s2;
929 bigint d = x0*(c+s2) + s4;
930 bigint cdash = - 2 * c;
931 bigint ddash = c*c - 4*d;
932
933 if(verbose)
934 {
935 cout<<"Third step, determining E(Q)/phi(E'(Q)) and E'(Q)/phi'(E(Q))\n";
936 cout<<"-------------------------------------------------------\n";
937 cout<<"1. E(Q)/phi(E'(Q))\n";
938 cout<<"-------------------------------------------------------\n";
939 }
940
941 if(verbose)
942 {
943 cout<<"(c,d) =("<<c<<","<<d<<")"<<endl;
944 cout<<"(c',d')=("<<cdash<<","<<ddash<<")"<<endl;
945 }
946
947 find_glsgens(0,c,d);
948 if (!success)
949 {
950 if (verbose) cout << "Failure in global descent with c,d="<<c<<","<<d<<endl;
951 return;
952 }
953
954 npoints1=npoints;
955 if(verbose)
956 {
957 cout<<"-------------------------------------------------------\n";
958 cout<<"2. E'(Q)/phi'(E(Q))\n";
959 cout<<"-------------------------------------------------------\n";
960 }
961 find_glsgens(1,cdash,ddash);
962 if (!success)
963 {
964 if (verbose) cout << "Failure in global descent with c',d'="<<cdash<<","<<ddash<<endl;
965 return;
966 }
967
968 //Debug only:
969 /*
970 cout<<"gls20 = "<<gls20<<endl;
971 cout<<"gls21 = "<<gls21<<endl;
972 cout<<"els0 = "<<els0<<endl;
973 cout<<"els1 = "<<els1<<endl;
974 cout<<"els20 = "<<els20<<endl;
975 cout<<"els21 = "<<els21<<endl;
976 */
977
978 rank = gls20+gls21-2; // lower bound for rank
979 //rank_bound = els20+els21-2; // upper bound for rank, from above
980
981 mw0 = 1<<gls20;
982 sel0 = 1<<els0; sha0=sel0/mw0;
983 sel20 = 1<<els20; sha20=sel20/mw0;
984
985 // gls20 is the 2-rank of E/phi(E') (lower bound)
986 // els0 is the 2-rank of S^(phi)(E') (exact)
987 // els20 is the 2-rank of its subgroup phi'(S^(2)(E)) (exact if 2nd descent)
988 // els0-gls20 is the 2-rank of III(E')[phi) (upper bound)
989 // els20-gls20 is the 2-rank of its subgroup phi'(III(E)[2]) (upper bound)
990
991 mw1 = 1<<gls21;
992 sel1 = 1<<els1; sha1=sel1/mw1;
993 sel21 = 1<<els21; sha21=sel21/mw1;
994
995 // gls21 is the 2-rank of E'/phi'(E) (lower bound)
996 // els1 is the 2-rank of S^(phi')(E) (exact)
997 // els21 is the 2-rank of its subgroup phi(S^(2)(E')) (exact if 2nd descent)
998 // els1-gls21 is the 2-rank of III(E)[phi'] (upper bound)
999 // els21-gls21 is the 2-rank of its subgroup phi(III(E')[2]) (upper bound)
1000
1001 index2=mw0*mw1;
1002 lindex=gls20+gls21;
1003 if(!ddash_is_sq) {index2/=2; lindex-=1;}
1004
1005 int certain0 = (els20==gls20); // i.e. descent was conclusive
1006 int certain1 = (els21==gls21); //
1007 certain = certain0&&certain1;
1008
1009 if (verbose)
1010 {
1011 cout<<"\n";
1012 cout<<"-------------------------------------------------------\n";
1013 cout<<"Summary of results:\n";
1014 cout<<"-------------------------------------------------------\n";
1015
1016 if(certain) cout << "\trank(E) " << "= " << rank;
1017 else cout << "\t"<<rank<< " <= rank(E) <= " << rank_bound;
1018 cout<<"\n";
1019 cout << "\t#E(Q)/2E(Q) ";
1020 if(!certain) cout<<">";
1021 cout<<"= " << index2 << "\n\n";
1022
1023 long lb, ub; int lb1, eq;
1024
1025 if(do_second_descent)
1026 {
1027 cout<<"Information on III(E/Q):\n";
1028
1029 lb = 1; lb1=1; ub = sha1; eq=(lb==ub);
1030 cout<<"\t#III(E/Q)[phi'] ";
1031 if(eq) cout<<"= "<<ub<<"\n";
1032 else
1033 {
1034 if(lb1) cout<<"<= "<<ub<<"\n";
1035 else cout<<"is between "<<lb<<" and "<<ub<<"\n";
1036 }
1037
1038 lb = sel1/sel21; lb1=(lb==1); ub = sha1*sha20; eq=(lb==ub);
1039 cout<<"\t#III(E/Q)[2] ";
1040 if(eq) cout<<"= "<<ub<<"\n";
1041 else
1042 {
1043 if(lb1) cout<<"<= "<<ub<<"\n";
1044 else cout<<"is between "<<lb<<" and "<<ub<<"\n";
1045 }
1046 cout<<endl;
1047
1048 cout<<"Information on III(E'/Q):\n";
1049
1050 lb = 1; lb1=1; ub = sha20; eq = (lb==ub);
1051 cout<<"\t#phi'(III(E/Q)[2]) ";
1052 if(!eq) cout<<"<";
1053 cout<<"= " << ub << "\n";
1054
1055 lb = sel0/sel20; lb1=(lb==1); ub = sha0; eq=(lb==ub);
1056 cout<<"\t#III(E'/Q)[phi] ";
1057 if(eq) cout<<"= "<<ub<<"\n";
1058 else
1059 {
1060 if(lb1) cout<<"<= "<<ub<<"\n";
1061 else cout<<"is between "<<lb<<" and "<<ub<<"\n";
1062 }
1063
1064 lb = sel0/sel20; lb1=(lb==1); ub = sha0*sha21; eq=(lb==ub);
1065 cout<<"\t#III(E'/Q)[2] ";
1066 if(eq) cout<<"= "<<ub<<"\n";
1067 else
1068 {
1069 if(lb1) cout<<"<= "<<ub<<"\n";
1070 else cout<<"is between "<<lb<<" and "<<ub<<"\n";
1071 }
1072 cout<<endl;
1073 } // end of do_second_descent branch
1074 else
1075 {
1076 cout<<"Information on III(E/Q):\n";
1077
1078 ub = sha1;
1079 if(ub==1) cout<<"\t III(E/Q)[phi'] is trivial\n";
1080 else cout<<"\t#III(E/Q)[phi'] <= " << ub << "\n";
1081
1082 ub = sha0*sha1;
1083 if(ub==1) cout<<"\t III(E/Q)[2] is trivial\n";
1084 else cout<<"\t#III(E/Q)[2] <= " << ub << "\n";
1085
1086 cout<<"Information on III(E'/Q):\n";
1087
1088 ub = sha0;
1089 if(ub==1) cout<<"\t III(E'/Q)[phi] is trivial\n";
1090 else cout<<"\t#III(E'/Q)[phi] <= " << ub << "\n";
1091
1092 ub = sha0*sha1;
1093 if(ub==1) cout<<"\t III(E'/Q)[2] is trivial\n";
1094 else cout<<"\t#III(E'/Q)[2] <= " << ub << "\n";
1095
1096 cout<<endl;
1097 } // end of no_second_descent branch
1098
1099 }
1100
1101 if(npoints>0) makegens();
1102 if(rank<=MAX_R)makepoints();
1103 } // end of "else" clause after if(selmer_only)
1104 }
1105
makegens()1106 void rank2::makegens()
1107 {
1108 Curvedata ee_min;
1109 bigint u, r, s, t, x, y, z; int i;
1110 ee_min=ee.minimalize(u,r,s,t);
1111 if(verbose)
1112 {
1113 cout<<"-------------------------------------------------------\n";
1114 cout << "\nList of points on E = " << (Curve)ee_min << ":\n";
1115 cout<<"\nI. Points on E mod phi(E')\n";
1116 if(npoints1==0) cout << "--none (modulo torsion).\n\n";
1117 }
1118 for(i=0; i<npoints; i++)
1119 {
1120 if(verbose&&(i==npoints1)) {cout<<"\nII. Points on phi(E') mod 2E\n";}
1121 Point q = transform(pointlist[i],the_curve,u,r,s,t);
1122 bigfloat h = height(q);
1123 int valid = q.isvalid();
1124 if(verbose||!valid) cout << "Point " << q << ", height = " << h;
1125 if(!valid) cout << " --warning: NOT on curve!";
1126 if(verbose||!valid) cout << "\n";
1127 pointlist[i]=q;
1128 }
1129 if(verbose&&(npoints1==npoints))
1130 {
1131 cout<<"\nII. Points on phi(E') mod 2E\n";
1132 cout << "--none (modulo torsion).\n\n";
1133 }
1134 }
1135
listgens()1136 void rank2::listgens()
1137 {
1138 long i;
1139 cout << "List of generators of E(Q)/2E(Q) for E = "
1140 << (Curve)(*the_curve) << ": \n";
1141 for(i=0; i<npoints; i++)
1142 {
1143 Point p = pointlist[i];
1144 cout << "Point " <<
1145 // "on " << (Curve)(*CD_orig) << ": " <<
1146 p;
1147 bigfloat h = height(p);
1148 cout << ", height = " << h;
1149 if(!p.isvalid()) cout << " --warning: NOT on curve!";
1150 cout << "\n";
1151 }
1152 }
1153
listgens(Curvedata * CD_orig,const bigint & u,const bigint & r,const bigint & s,const bigint & t)1154 void rank2::listgens(Curvedata* CD_orig, const bigint& u, const bigint& r,
1155 const bigint& s, const bigint& t)
1156 {
1157 long i;
1158 cout << "List of generators of E(Q)/2E(Q) (mod torsion) for E = "
1159 << (Curve)(*CD_orig) << ": \n";
1160 for(i=0; i<npoints; i++)
1161 {
1162 Point p = transform(pointlist[i],CD_orig,u,r,s,t,1);
1163 cout << "Point " << (i+1) <<
1164 // "on " << (Curve)(*CD_orig) <<
1165 ": " << p;
1166 bigfloat h = height(p);
1167 cout << ", height = " << h;
1168 if(!p.isvalid()) cout << " --warning: NOT on curve!";
1169 cout << "\n";
1170 }
1171 }
1172
makepoints()1173 void rank2::makepoints()
1174 {
1175 if(fullnpoints>0) return; // avoids calling this twice
1176 int i, j;
1177 long smallindex = index2/(1+ntwo_torsion);
1178 fullnpoints=1; // will be smallindex
1179 fullpointlist.resize(smallindex);
1180 fullpointlist[0]=Point(the_curve);
1181
1182 if(verbose&&(rank>0))
1183 {
1184 cout<<"-------------------------------------------------------\n";
1185 cout << "Computing full set of "<<smallindex<<" coset representatives for\n";
1186 cout << "2E(Q) in E(Q) (modulo torsion), and sorting into height order...." << flush;
1187 }
1188 for(i=0; i<rank; i++)
1189 {
1190 for(j=0; j<fullnpoints; j++)
1191 fullpointlist[j+fullnpoints] = fullpointlist[j]+pointlist[i];
1192 fullnpoints*=2;
1193 }
1194 if(fullnpoints!=smallindex)
1195 {
1196 cout << "Problem: index = " << index2 << " but " <<fullnpoints<<" cosets\n";
1197 }
1198 //
1199 // Now reorder points into increasing height order:
1200 //
1201 for(i=0; i<fullnpoints; i++)
1202 for(j=i+1; j<fullnpoints; j++)
1203 if(height(fullpointlist[j])<height(fullpointlist[i]))
1204 {
1205 Point temp = fullpointlist[i];
1206 fullpointlist[i]=fullpointlist[j];
1207 fullpointlist[j]=temp;
1208 }
1209 if(verbose&&(rank>0)) cout << "done.\n" << endl;
1210 }
1211
listpoints()1212 void rank2::listpoints()
1213 {
1214 makepoints();
1215 cout << "Points on curve E = " << (Curve)(*the_curve)
1216 << " covering E(Q)/2E(Q), modulo torsion:";
1217 if(rank==0) cout<<" none.";
1218 else
1219 if(rank>MAX_R)
1220 cout << "Too many to list ("<<fullnpoints<<" mod torsion)\n";
1221 else
1222 {
1223 cout<<"\n"<<fullnpoints<<" points, [0:1:0] and:\n";
1224 for (long i=1; i<fullnpoints; i++)
1225 {
1226 Point p = fullpointlist[i];
1227 cout << "Point " << p;
1228 bigfloat h = height(p);
1229 cout << ", height = " << h;
1230 if(!p.isvalid()) {cout << " --warning: NOT on curve! " ;}
1231 cout << "\n";
1232 }
1233 }
1234 cout<<"\n\n";
1235 }
1236
listpoints(Curvedata * CD_orig,const bigint & u,const bigint & r,const bigint & s,const bigint & t)1237 void rank2::listpoints(Curvedata* CD_orig, const bigint& u, const bigint& r,
1238 const bigint& s, const bigint& t)
1239 {
1240 makepoints();
1241 cout << "Points on original curve E = " << (Curve)(*CD_orig)
1242 << " covering E(Q)/2E(Q), modulo torsion:";
1243 if(rank==0) cout<<" none.";
1244 else if(rank>MAX_R)
1245 cout << "Too many to list ("<<fullnpoints<<" mod torsion)\n";
1246 else
1247 {
1248 cout<<"\n"<<fullnpoints<<" points:"<<endl;
1249 for (long i=1; i<fullnpoints; i++)
1250 {
1251 Point p0 = fullpointlist[i];
1252 // cout << "Point on " << (Curve)(*the_curve) << ": " << p0;
1253 // if(!p0.isvalid()) cout << " --warning: NOT on curve!";
1254 // cout << "\n";
1255 Point p = transform(p0,CD_orig,u,r,s,t,1);
1256 cout << "Point " <<
1257 // "on " << (Curve)(*CD_orig) << ": " <<
1258 p;
1259 bigfloat h = height(p);
1260 cout << ", height = " << h;
1261 if(!p.isvalid()) cout << " --warning: NOT on curve!";
1262 cout << "\n";
1263 }
1264 }
1265 cout<<"\n\n";
1266 }
1267