1 // FILE NEWFORMS.CC: implementation of newforms class
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>
25 #include <eclib/newforms.h>
26 #include <eclib/periods.h>
27 #include <eclib/curvesort.h>
28
29 // Functions for ordering newforms
30 // (1) Old ordering (first aq, then ap for good p);
31 // with the order for eigenvalues being
32 // 1,-1 or 0,1,-1,2,-2,...
33 // (2) New ordering (ap for all p in natural order)
34 // with the order for eigenvalues being
35 // -1,1 or ...,-2,-1,0,1,2,... (plain numerical order)
36
37 // less_ap(a,b) returns +1,0,-1 according to whether a is
38 // before/equal/after b in the above ordering
less_ap(long a,long b,int old=0)39 int less_ap(long a, long b, int old=0)
40 {
41 if(!old) return sign(b-a); // the simple new ordering!
42 if(a==b) return 0;
43 int s = sign(abs(b)-abs(a));
44 if(s) return s;
45 else return sign(a-b); // this way round! +1 before -1 etc
46 }
47
48 // Compare two ap-vectors lexicographically, using less_ap(.,.,old):
49 int less_apvec(const vector<long>& v, const vector<long>& w, int old=0);
less_apvec(const vector<long> & v,const vector<long> & w,int old)50 int less_apvec(const vector<long>& v, const vector<long>& w, int old)
51 {
52 vector<long>::const_iterator vi=v.begin(), wi=w.begin();
53 while(vi!=v.end())
54 {
55 int s = less_ap(*vi++,*wi++,old);
56 if(s) return s;
57 }
58 return 0;
59 }
60
61 struct less_newform_old : public binary_function<newform, newform, bool> {
operator ()less_newform_old62 bool operator()(const newform& f, const newform& g)
63 {
64 int s = less_apvec(f.aqlist,g.aqlist,1);
65 if(s==0) s = less_apvec(f.aplist,g.aplist,1);
66 return (s==1);
67 }
68 };
69
70 struct less_apvec_function : public binary_function<const vector<long>&, const vector<long>&, bool> {
operator ()less_apvec_function71 bool operator()(const vector<long>& f, const vector<long>& g)
72 {
73 return 1==less_apvec(f,g);
74 }
75 };
76
eiglist(const newform & f,int oldorder)77 vector<long> eiglist(const newform& f, int oldorder)
78 {
79 /*
80 cout<<"Entering eiglist with f.aqlist="<<f.aqlist<<"\nand f.aplist=";
81 vec_out(cout,f.aplist,10);
82 cout<<endl;
83 */
84 vector<long> eigs;
85 primevar pr;
86 long N = (f.nf)->modulus;
87 vector<long>::const_iterator aqi=f.aqlist.begin();
88 vector<long>::const_iterator api=f.aplist.begin();
89 vector<long>::iterator eigsi;
90 if(oldorder)
91 {
92 eigs.resize(f.aplist.size());
93 eigsi=eigs.begin();
94 while(aqi!=f.aqlist.end())
95 *eigsi++ = *aqi++;
96 while(api!=f.aplist.end())
97 {
98 if(ndivides(pr,N)) *eigsi++ = *api;
99 api++; pr++;
100 }
101 }
102 else
103 {
104 eigs=f.aplist; // copy; now adjust the aq:
105 eigsi=eigs.begin();
106 while((aqi!=f.aqlist.end())&&(eigsi!=eigs.end()))
107 {
108 if(::divides(pr.value(),N)) *eigsi = (*aqi++);
109 eigsi++; pr++;
110 }
111 }
112 /*
113 cout<<"Leaving eiglist with eigs=";
114 vec_out(cout,eigs,10);
115 cout<<endl;
116 */
117 return eigs;
118 }
119
120 struct less_newform_new : public binary_function<newform, newform, bool> {
operator ()less_newform_new121 bool operator()(const newform& f, const newform& g)
122 {
123 // return less_apvec(eiglist(f),eiglist(g),0)==1;
124 return less_apvec(f.aplist,g.aplist,0)==1;
125 }
126 };
127
128 // Newform constructor given the ap and aq lists and extra data (but
129 // no homology basis), e.g. after reading from newforms file
130
newform(const vector<int> & data,const vector<long> & aq,const vector<long> & ap,newforms * nfs)131 newform::newform(const vector<int>& data, const vector<long>& aq, const vector<long>& ap, newforms* nfs) : nf(nfs)
132 {
133 sfe=data[0];
134 ap0=data[1];
135 np0=data[2];
136 dp0=data[3];
137 loverp=rational(dp0,np0);
138 lplus=data[4];
139 mplus=data[5];
140 lminus=data[6];
141 mminus=data[7];
142 a=data[8];
143 b=data[9];
144 c=data[10];
145 d=data[11];
146 dotplus=data[12];
147 dotminus=data[13];
148 type=data[14];
149 degphi=data[15];
150 aqlist=aq;
151 aplist=ap;
152 index=-1;
153 pdot=0;
154 rk=-1;
155 }
156
157 // Newform constructor, given the homology basis vector(s) and
158 // Hecke eigenvalues
159
newform(const vec & vplus,const vec & vminus,const vector<long> & ap,newforms * nfs,long ind)160 newform::newform(const vec& vplus, const vec& vminus, const vector<long>& ap, newforms* nfs,long ind)
161 :nf(nfs), sign(nfs->sign), bplus(vplus),bminus(vminus),index(ind),aplist(ap),rk(-1)
162 {
163 int verbose=(nf->verbose);
164
165 if(verbose)
166 {
167 cout<<"Creating H1";
168 if(sign==+1) cout<<"+";
169 if(sign==-1) cout<<"-";
170 cout<<" newform from aplist..."<<endl;
171 if(verbose>2)
172 {
173 if(sign!=-1) cout<<"bplus = "<<bplus<<endl;
174 if(sign!=+1) cout<<"bminus = "<<bminus<<endl;
175 }
176 }
177
178 // check_expand_contract();
179
180 // Fixing the eigenvalue lists: ap is indexed by primes in natural
181 // order we need to extract aq (computing any not yet there).
182
183 // At the same time we change the entries in aplist for bad primes q
184 // from the Wq-eigenvalue to the newform coefficient.
185 fixup_eigs();
186
187 // Compute cuspidalfactors and type (only if sign=0):
188
189 type = 0;
190 find_cuspidal_factors();
191
192 // Compute coordsplus/minus and denomplus/minus
193
194 find_coords_plus_minus();
195
196 // Compute pdot, dp0, loverp (unless sign is -1)
197
198 find_bsd_ratio();
199
200 // Find deg(phi) (only if sign is 0)
201 degphi = 0;
202 find_degphi();
203
204 // Find twisting primes if N non-square
205
206 lplus=mplus=0;
207 lminus=mminus=0;
208 find_twisting_primes();
209
210 // find a,b,c,d,dotplus,dotminus
211
212 a=b=c=d=0;
213 dotplus=dotminus=0;
214 find_matrix();
215
216 // Set default values for optimalityfactrplus/minus (will be reset if constructing from curves)
217 optimalityfactorplus = 1;
218 optimalityfactorminus = 1;
219 }
220
check_expand_contract()221 int newform::check_expand_contract()
222 {
223 int success=1;
224 long denom = nf->h1->h1denom();
225 vec bplusx, bminusx, tvec;
226 if (sign!=-1)
227 {
228 bplusx= nf->h1->extend_coords(bplus);
229 tvec = nf->h1->contract_coords(bplusx);
230 tvec /= denom;
231 if (tvec!=bplus)
232 {
233 success=0;
234 cout<<"! bplus ="<<bplus<<" extends to "<<bplusx<<" which contracts to "<<tvec<<endl;
235 }
236 }
237 if (sign!=+1)
238 {
239 bminusx= nf->h1->extend_coords(bminus);
240 tvec = nf->h1->contract_coords(bminusx);
241 tvec /= denom;
242 if (tvec!=bminus)
243 {
244 success=0;
245 cout<<"! bminus="<<bminus<<" extends to "<<bminusx<<" which contracts to "<<tvec<<endl;
246 }
247 }
248 return success;
249 }
250
fixup_eigs()251 void newform::fixup_eigs()
252 {
253 long denom = nf->h1->h1denom();
254 aqlist.resize(nf->npdivs);
255 vector<long>::iterator api=aplist.begin(), pi=nf->plist.begin();
256 vector<long>::iterator aqi=aqlist.begin();
257 primevar pr; long q, i;
258 long n = nf->modulus;
259 while((api!=aplist.end())&&(aqi!=aqlist.end()))
260 {
261 q=pr.value(); pr++;
262 if(::divides(q,n))
263 {
264 *aqi++=*api;
265 *api=(::divides(q*q,n)? 0: -*api);
266 pi++;
267 }
268 api++;
269 }
270 if(aqi!=aqlist.end()) // compute missing aq
271 {
272 long piv;
273 ssubspace espace;
274 if(sign==-1)
275 espace=make1d(bminus,piv);
276 else
277 espace=make1d(bplus,piv);
278 piv*=denom;
279 while(aqi!=aqlist.end()) // compute missing aq
280 {
281 q=*pi++;
282 if(nf->verbose) cout<<"Computing Wq for q="<<q<<"..."<<flush;
283 smat Wq = nf->h1->s_heckeop_restricted(q,espace,1,0);
284 long aq = Wq.elem(1,1) / piv;
285 if(nf->verbose) cout<<"aq ="<<aq<<endl;
286 *aqi++=aq;
287 }
288 }
289 if(nf->verbose) cout<<"aqlist = "<<aqlist<<endl;
290
291 //Compute sfe:
292
293 sfe=-1;
294 for(i=0; i<(nf->npdivs); i++) sfe*=aqlist[i];
295 if(nf->verbose) cout<<"sfe = "<<sfe<<endl;
296 }
297
298 // Before recovering eigenbases, we need to put back the aq into the
299 // aplist (and resort, for efficiency).
unfix_eigs()300 void newform::unfix_eigs()
301 {
302 vector<long>::iterator api=aplist.begin();
303 vector<long>::iterator aqi=aqlist.begin();
304 primevar pr;
305 long n = nf->modulus;
306 while((api!=aplist.end())&&(aqi!=aqlist.end()))
307 {
308 if(::divides(pr.value(),n)) *api=*aqi++;
309 api++;
310 pr++;
311 }
312 }
313
314 // After recovering eigenbases, we need to replace the ap for bad p
refix_eigs()315 void newform::refix_eigs()
316 {
317 vector<long>::iterator api=aplist.begin();
318 primevar pr;
319 long n = nf->modulus, np = nf->npdivs, ip=0, q;
320 while((api!=aplist.end())&&(ip<np))
321 {
322 q=pr.value();
323 if(::divides(q,n))
324 {
325 *api=(::divides(q*q,n)? 0: -*api);
326 ip++;
327 }
328 api++;
329 pr++;
330 }
331 }
332
find_bsd_ratio()333 void newform::find_bsd_ratio()
334 {
335 // get ap for p=p0:
336
337 primevar pr;
338 vector<long>::const_iterator api=aplist.begin();
339 while(pr.value()!=nf->p0) {pr++; api++;}
340 ap0=*api;
341 np0 = 1 + (nf->p0) - ap0;
342 if(nf->verbose) cout<<"ap0 = "<<ap0<<"\tnp0 = "<<np0<<endl;
343
344 if(sign==-1) return;
345
346 pdot = (nf->mvp)*bplus; // should be negative since L(f,1)>=0
347 if (pdot>0)
348 // NB This will ensure that plus modular symbols have the right
349 // sign for curves where L(E,1) is nonzero, but more work is
350 // necessary for the plus symbols when L(Em1)=0, and for minus
351 // symbols. The additional work is done in find_matrix().
352 {
353 coordsplus *= -1;
354 bplus *= -1;
355 pdot *= -1;
356 }
357 dp0=abs(pdot);
358 // DO NOT scale pdot by denom: factor will cancel when used to compute ap
359 // DO scale dp0 since it is used to compute L/P
360 if(dp0!=0)
361 {
362 if(denomplus>1)
363 {
364 if(::divides(denomplus,dp0)) dp0/=denomplus;
365 else
366 cout<<"newform constructor error: dp0 not divisible by denomplus!"
367 <<endl;
368 }
369 }
370 loverp = rational(dp0,np0);
371 if(nf->verbose)
372 {
373 cout<<"pdot = "<<pdot<<endl;
374 cout<<"dp0 = "<<dp0<<endl;
375 cout<<"np0 = "<<np0<<endl;
376 cout<<"loverp = "<<loverp<<endl;
377 }
378 }
379
find_coords_plus_minus()380 void newform::find_coords_plus_minus()
381 {
382 int verbose = nf->verbose;
383 int i, ncoords=nf->h1->coord_vecs.size()-1;
384 // ncoords is the same as ngens in homspace, i.e. the number of symbols aftre 2-term relations
385 svec cvi;
386 if(sign!=-1)
387 coordsplus=vec(ncoords);
388 if(sign!=+1)
389 coordsminus=vec(ncoords);
390 // if(verbose) cout<<"About to compute coordsplus/minus"<<endl;
391 for(i=1; i<=ncoords; i++)
392 {
393 cvi = nf->h1->coord_vecs[i];
394 if(sign!=-1)
395 coordsplus[i]=dotmodp(cvi,bplus,MODULUS);
396 if(sign!=+1)
397 coordsminus[i]=dotmodp(cvi,bminus,MODULUS);
398 }
399 contplus=vecgcd(coordsplus);
400 if (contplus>1) coordsplus/=contplus;
401 contminus=vecgcd(coordsminus);
402 if (contminus>1) coordsminus/=contminus;
403
404 if(sign!=+1)
405 {
406 denomminus=contminus*cuspidalfactorminus;
407 if(verbose>1) cout<<"coordsminus = "<<coordsminus<<endl;
408 if(verbose) cout<<"denomminus = "<<denomminus<<endl;
409 }
410 if(sign!=-1)
411 {
412 denomplus=contplus*cuspidalfactorplus;
413 if(verbose>1) cout<<"coordsplus = "<<coordsplus<<endl;
414 if(verbose) cout<<"denomplus = "<<denomplus<<endl;
415 }
416 }
417
find_cuspidal_factors()418 void newform::find_cuspidal_factors()
419 {
420 vec bplusc, bminusc;
421 int verbose = nf->verbose;
422
423 cuspidalfactorplus=1;
424 cuspidalfactorminus=1;
425
426 if(!(nf->h1->cuspidal))
427 {
428 if(sign!=-1) // do this if sign = 0,1
429 {
430 bplusc=(nf->h1->tkernbas)*bplus;
431 cuspidalfactorplus = vecgcd(bplusc);
432 bplusc /= cuspidalfactorplus;
433 }
434 if(sign!=+1) // do this if sign = 0,-1
435 {
436 bminusc=(nf->h1->tkernbas)*bminus;
437 cuspidalfactorminus = vecgcd(bminusc);
438 bminusc/= cuspidalfactorminus;
439 }
440 if(sign==0) // do this only if sign = 0
441 {
442 type=3-vecgcd(bplusc-bminusc);
443 if(verbose) cout<<"Lattice type = "<<type<<endl;
444 if((type!=1)&&(type!=2))
445 {
446 cerr<<"Error: lattice type computed to be "<<type<<", should be 1 or 2!"<<endl;
447 }
448 }
449
450 if(verbose&&(cuspidalfactorplus*cuspidalfactorminus>1))
451 {
452 if(sign!=-1)
453 {
454 cout<<"cuspidalfactorplus = "<<cuspidalfactorplus<<endl;
455 if(verbose>2) cout<<"bplusc = "<<bplusc<<endl;
456 }
457 if(sign!=+1)
458 {
459 cout<<"cuspidalfactorminus = "<<cuspidalfactorminus<<endl;
460 if(verbose>2) cout<<"bminusc = "<<bminusc<<endl;
461 }
462 }
463 }
464 }
465
find_degphi()466 void newform::find_degphi()
467 {
468 if(sign!=0) return;
469 #ifdef DEG_PHI
470 if(nf->verbose) cout<<"computing deg(phi)..."<<flush;
471 degphi=jumpinfo->degphi(bplusc,bminusc,type);
472 if(nf->verbose) cout<<"done..."<<flush;
473 #else
474 degphi=0;
475 #endif
476 }
477
find_twisting_primes()478 void newform::find_twisting_primes()
479 {
480 int verbose=(nf->verbose);
481 if(verbose) cout<<"computing twisting primes (sign="<<sign<<")..."<<flush;
482 if(sign!=-1)
483 {
484 if(dp0!=0)
485 {
486 lplus=1; // so we need not search for a prime 1(mod 4) below
487 mplus=1; // dummy value, not used
488 }
489 else
490 {
491 lplus=0;
492 mplus =0;
493 }
494 }
495 if(sign!=+1)
496 {
497 lminus=0;
498 mminus=0;
499 }
500 if(nf->squarelevel) return;
501
502 long n = nf->modulus;
503 for (primevar lvar; lvar.ok() &&
504 (((sign!=-1)&&(mplus==0)) ||
505 ((sign!=+1)&&(mminus==0))); lvar++)
506 {
507 //cout << "Trying l = " << lvar << endl;
508 while (n%lvar==0) {lvar++;}
509 long l = lvar;
510 //cout << "Trying l = " << l << endl;
511 if (legendre(-n,l)!=sfe) continue;
512 //cout << "Legendre condition passed... " << endl;
513
514 if((sign!=-1)&&(mplus==0)&&(l%4==1))
515 {
516 lplus = l;
517 //cout << "Trying lplus = " << l << "\n";
518 map<long,vec>::const_iterator vi = nf->mvlplusvecs.find(l);
519 if(vi==nf->mvlplusvecs.end())
520 mplus = (nf->mvlplusvecs[l]=nf->h1->manintwist(l))*bplus;
521 else
522 mplus = (vi->second)*bplus;
523 // We force mplus>0 to fix the sign of the modular symbol to
524 // agree with L(f*chi,1)>0, since L(f*chi,1) is real and a
525 // positive multiple of mplus. This uses the fact that the
526 // Gauus sum is +sqrt(l).
527 if (mplus<0)
528 {
529 mplus *= -1;
530 bplus *= -1;
531 coordsplus *= -1;
532 }
533 if((denomplus>1)&&(mplus!=0))
534 {
535 if(::divides(denomplus,mplus)) mplus/=denomplus;
536 else
537 cout<<"Warning in newform constructor: mplus not divisible by denomplus!"
538 <<endl;
539 }
540 }
541 if((sign!=+1)&&(mminus==0)&&(l%4==3))
542 {
543 lminus = l;
544 //cout << "Trying lminus = " << l << "\n";
545 map<long,vec>::const_iterator vi = nf->mvlminusvecs.find(l);
546 if(vi==nf->mvlminusvecs.end())
547 mminus = (nf->mvlminusvecs[l]=nf->h1->manintwist(l))*bminus;
548 else
549 mminus = (vi->second)*bminus;
550 // We force mminus<0 to fix the sign of the modular symbol
551 // to agree with L(f*chi,1)>0, since L(f*chi,1) is real and
552 // a negative multiple of mminus. This uses the fact that
553 // the Gauus sum is +i*sqrt(l).
554 if (mminus>0)
555 {
556 mminus *= -1;
557 bminus *= -1;
558 coordsminus *= -1;
559 }
560 if((denomminus>1)&&(mminus!=0))
561 {
562 if(::divides(denomminus,mminus)) mminus/=denomminus;
563 else
564 cout<<"Warning in newform constructor: mminus="<<mminus<<" is not divisible by denomminus="<<denomminus<<"!"
565 <<endl;
566 }
567 }
568 }
569
570 if(verbose)
571 {
572 cout<<"done..."<<flush;
573 cout<<"lplus = "<<lplus<<endl;
574 cout<<"mplus = "<<mplus<<endl;
575 cout<<"lminus = "<<lminus<<endl;
576 cout<<"mminus = "<<mminus<<endl;
577 }
578 }
579
find_matrix()580 void newform::find_matrix()
581 {
582 int verbose=(nf->verbose);
583 if(verbose) cout<<"computing a,b,c,d..."<<flush;
584 long n = nf->modulus;
585 int found=0;
586 vec v;
587 for(d=2; !found; d++)
588 {
589 if(1==gcd(d,n))
590 {
591 for(b=1; (b<d) && !found; b++)
592 {
593 if(1==bezout(d,-n*b,a,c))
594 {
595 // cout<<"b/d = "<<b<<"/"<<d<<": ";
596 v = nf->h1->coords(b,d).as_vec();
597 // cout<<"v="<<v<<endl;
598 if(sign!=-1)
599 {
600 dotplus=v*bplus;
601 if(::divides(denomplus,dotplus))
602 dotplus/=denomplus;
603 else
604 cout<<"Warning in find_matrix: dotplus not divisible by denomplus!"<<endl;
605 }
606 if(sign!=+1)
607 {
608 dotminus=v*bminus;
609 if(::divides(denomminus,dotminus))
610 dotminus/=denomminus;
611 else
612 cout<<"Warning in find_matrix: dotminus not divisible by denomminus!"<<endl;
613 }
614 found=(((dotplus!=0)||(sign==-1))&&
615 ((dotminus!=0)||(sign==+1)));
616 }
617 }
618 }
619 }
620 b--; d--; //because they get incremented BEFORE the loop end-test
621 if(d<0) {a=-a; b=-b; c=-c; d=-d;} // because we need d>0 for integration
622
623 if(verbose)
624 {
625 cout<<"done: ";
626 cout << "[(" <<a<<","<<b<<";"<<c
627 <<","<<d<<"),"<<dotplus<<","<<dotminus
628 <<";"<<type<<"]"<<endl;
629 }
630 }
631
add_more_ap(int nap)632 void newform::add_more_ap(int nap)
633 {
634 if((int)aplist.size()>=nap) return;
635 int verbose=(nf->verbose);
636 long piv, p, ap;
637 // Do not make the espace right away, as it is possible that the
638 // only ap we are missing are aq which we already have...
639 ssubspace espace;
640 int have_espace=0;
641
642 primevar pr(nap,aplist.size()+1);
643 while((int)aplist.size()<nap)
644 {
645 p=pr;
646 if(::divides(p,nf->modulus))
647 {
648 if(::divides(p*p,nf->modulus))
649 ap=0;
650 else
651 ap=-aqlist[find(nf->plist.begin(),nf->plist.end(),p)-nf->plist.begin()];
652 }
653 else
654 {
655 if(verbose>1) cout<<"Computing Tp for p="<<p<<endl;
656 if(!have_espace)
657 {
658 if(sign==-1)
659 espace=make1d(bminus,piv);
660 else
661 espace=make1d(bplus,piv);
662 piv*=nf->h1->h1denom();
663 have_espace=1;
664 }
665 ap = (nf->h1->s_heckeop_restricted(p,espace,1,0)).elem(1,1) / piv;
666 }
667 aplist.push_back(ap);
668 pr++;
669 }
670 if(verbose>1) cout<<"aplist = "<<aplist<<endl;
671 }
672
673 // Compute analytic rank and special value if not set
compute_rank()674 void newform::compute_rank()
675 {
676 if (rk==-1) // not yet computed
677 {
678 ldash1 x(nf, this);
679 Lvalue = abs(x.value()); // = L^{(r)}(f,1) -- note the r! factor!
680 rk = 0;
681 if (num(loverp)==0) // else trivially 0
682 rk = x.rank();
683 }
684 }
685
rank()686 long newform::rank()
687 {
688 compute_rank(); return rk;
689 }
690
special_value()691 bigfloat newform::special_value()
692 {
693 compute_rank(); return Lvalue;
694 }
695
696 // To find optimality factors when created from a curve E. Here
697 // CP_opt is the periods of this newform, and of the optimal curve E0,
698 // which we compute earlier in the newforms class since that is where
699 // getperiods() is implemented.
700
find_optimality_factors(const CurveRed & E,int i)701 void newform::find_optimality_factors(const CurveRed& E, int i)
702 {
703 int verbose=(nf->verbose);
704 bigcomplex w1,w2;
705 bigfloat x0, y0, xE, yE;
706
707 // Definitions: for the period lattice of the optimal curve or
708 // newform, x0 = (type)*(least real period) and y0 = (type)*(least
709 // imag period)/i. Here type = #components = 1 (Delta<0) or 2
710 // (Deta>0). Note that if we are in the plus or minus quotient then
711 // we can still compute these even though we do not now the type
712 // since they generate the projections of the period lattice to the
713 // real (resp. imaginary) axis.
714
715 // Similarly xE, yE for the input curve E, though here we do know
716 // the type.
717
718 // Compute the real and/or imginary periods of the newform, which
719 // are those of the optimal curve in the isogeny class:
720 if (sign==+1)
721 {
722 nf->get_real_period(i,x0);
723 x0 = 2*abs(x0);
724 }
725 else if (sign==-1)
726 {
727 // NB it is impossible to get the scaling right in this case
728 nf->get_imag_period(i,y0);
729 y0 = abs(y0);
730 }
731 else
732 {
733 Cperiods CP_opt = nf->getperiods(i);
734 int opt_type = CP_opt.getwRI(w1,w2);
735 x0 = opt_type * abs(w1.real());
736 y0 = abs(w2.imag());
737 }
738
739 // Compute the real and/or imginary periods of the input curve, which may not be optimal:
740 Cperiods CP(E);
741 int Etype = CP.getwRI(w1,w2);
742 xE = Etype * abs(w1.real()); // least real period
743 yE = abs(w2.imag()); // least imag period
744 // now xE, yE are twice the least real/imag part of a period of E in both cases
745
746 // Now we find rational approximations to the rations x0/xE and
747 // y0/yE. These will have very small numerators and denominators.
748 long n,d;
749 if (sign!=-1)
750 {
751 ratapprox(x0/xE,n,d,163);
752 optimalityfactorplus = rational(n,d);
753 if (verbose) cout << "x-ratio (optimalityfactorplus) = " << (x0/xE) << " = " <<n<<"/"<<d<<" = "<<optimalityfactorplus << endl;
754 }
755 if (sign!=+1)
756 {
757 ratapprox(y0/yE,n,d,163);
758 optimalityfactorminus = rational(n,d);
759 if (verbose) cout << "y-ratio (optimalityfactorminus) = " << (y0/yE) << " = " <<n<<"/"<<d<<" = "<< optimalityfactorminus << endl;
760 }
761 }
762
763 // Adjust the sign of dotplus/dotminus, mplus/mminus, pdot and the
764 // associated primitive coordinate vectors by computing a period
765 // numerically. Here we only need sufficient precision to determine
766 // the sign, knowing the values to be nonzero rationals with small
767 // denominator.
768
sign_normalize()769 void newform::sign_normalize()
770 {
771 int verbose = (nf->verbose);
772 int sign = (nf->sign);
773
774 periods_direct integrator(nf, this);
775 integrator.compute();
776 bigfloat x0 = integrator.rper();
777 bigfloat y0 = integrator.iper();
778 if(verbose>1)
779 cout<<"integral over {0,"<<b<<"/"<<d<<"} = ("<<x0<<")+i*("<<y0<<")"<<endl;
780 if (sign!=-1)
781 {
782 if (dotplus*x0<0)
783 {
784 if (verbose)
785 cout<<"flipping sign for plus symbols"<<endl;
786 coordsplus *= -1;
787 bplus *= -1;
788 dotplus *= -1;
789 pdot *= -1;
790 mplus *= -1;
791 }
792 }
793 if (sign!=+1)
794 {
795 if (dotminus*y0<0)
796 {
797 if (verbose)
798 cout<<"flipping sign for minus symbols"<<endl;
799 coordsminus *= -1;
800 dotminus *= -1;
801 bminus *= -1;
802 mminus *= -1;
803 }
804 }
805 if (verbose>1)
806 {
807 cout<<"Approximate numerical values:"<<endl;
808 if (sign==0)
809 {
810 cout<<"x = "<<(x0/dotplus)<<endl;
811 cout<<"y = "<<(y0/dotminus)<<endl;
812 cout<<"integral over {0,"<<b<<"/"<<d<<"} = ("<<x0<<")+i*("<<y0<<")"<<endl;
813 }
814 if (sign==+1)
815 {
816 cout<<"x = "<<(x0/dotplus)<<endl;
817 cout<<"integral over {0,"<<b<<"/"<<d<<"} has real part "<<x0<<endl;
818 }
819 if (sign==-1)
820 {
821 cout<<"y = "<<(y0/dotminus)<<endl;
822 cout<<"integral over {0,"<<b<<"/"<<d<<"} has imaginary part "<<x0<<endl;
823 }
824 }
825 }
826
~newforms(void)827 newforms::~newforms(void)
828 {
829 delete of;
830 delete h1plus;
831 delete h1minus;
832 delete h1full;
833 }
834
makeh1(int s)835 void newforms::makeh1(int s)
836 {
837 if(s==1)
838 {
839 if(!h1plus)
840 {
841 if(verbose) cout<<"Constructing H1 (with sign=+1) ..."<<flush;
842 h1plus = new homspace(modulus,1,0,0 /*verbose*/);
843 if(verbose) cout<<"done"<<endl;
844 }
845 h1 = h1plus;
846 return;
847 }
848 if(s==-1)
849 {
850 if(!h1minus)
851 {
852 if(verbose) cout<<"Constructing H1 (with sign=-1) ..."<<flush;
853 h1minus = new homspace(modulus,-1,0,0 /*verbose*/);
854 if(verbose) cout<<"done"<<endl;
855 }
856 h1 = h1minus;
857 return;
858 }
859 if(s==0)
860 {
861 if(!h1full)
862 {
863 if(verbose) cout<<"Constructing H1 (with sign=0) ..."<<flush;
864 h1full = new homspace(modulus,0,0,0 /*verbose*/);
865 if(verbose) cout<<"done"<<endl;
866 }
867 h1 = h1full;
868 return;
869 }
870 cout<<"Error in makeh1(s): s = "<<s<<" should be one of 0,1,-1"<<endl;
871 return;
872 }
873
createfromscratch(int s,long ntp)874 void newforms::createfromscratch(int s, long ntp)
875 {
876 sign = s;
877 makeh1(s);
878 // cout<<"Constructing oldforms with sign="<<sign<<endl;
879 of = new oldforms(ntp,h1,(verbose>1),sign); // h1 provides the level*
880 if(verbose>1) of->display();
881 maxdepth = of->nap;
882 long mindepth = npdivs; // must include at least one good p, and it
883 // is cheap to continue recursing after
884 // reaching dimension 1.
885 n1ds = 0;
886 int upperbound = h1->dimension-(of->totalolddim);
887 if(upperbound>0) // Else no newforms certainly so do no work!
888 {
889 mvp=h1->maninvector(p0);
890 // cout<<"mvp = "<<mvp<<endl;
891 if(verbose>1) cout<<"h1 denom = "<<h1->h1denom()<<endl;
892 form_finder ff(this,(sign!=0),maxdepth,mindepth,1,0,verbose);
893 basisflag=0;
894 ff.find();
895 }
896 if(verbose)
897 {
898 cout << "Total dimension = " << h1->dimension << endl;
899 cout << "Number of rational newforms = " << n1ds <<endl;
900 if(h1->dimension==of->totalolddim+n1ds)
901 cout<<"The whole space splits over Q" << endl;
902 }
903
904 if(n1ds==0) return;
905
906 int i,nap,maxnap=0;
907
908 if((n1ds>1)&&(modulus<130000)) // reorder into old order
909 {
910 if(verbose) cout<<"Reordering newforms into old order as N<130000"<<endl;
911 // if(verbose) cout<<"Before sorting:\n"; display();
912 sort(1);
913 // if(verbose) cout<<"After sorting:\n"; display();
914 }
915 // At this point the newforms may contain different numbers of ap,
916 // so we need to even these up, which we do by computing more ap for
917 // those which need it.
918 if(n1ds>0)
919 {
920 for(i=0; i<n1ds; i++)
921 if((nap=nflist[i].aplist.size())>maxnap) maxnap=nap;
922 if(verbose)
923 cout<<"Max number of ap in newforms so far is "<<maxnap
924 <<", increasing to " << DEFAULT_SMALL_NAP << endl;
925 if (maxnap < DEFAULT_SMALL_NAP) maxnap = DEFAULT_SMALL_NAP;
926 for(i=0; i<n1ds; i++)
927 {
928 if((nap=nflist[i].aplist.size())<maxnap)
929 {
930 if(verbose)
931 cout<<"Newform #"<<(i+1)<<" has only "<<nap
932 <<" ap so we need to compute more..."<<endl;
933 nflist[i].add_more_ap(maxnap);
934 }
935
936 // Now if necessary we adjust the sign of dotplus/dotminus and the
937 // associated primitive coordinate vectors by computing a period
938 // numerically. Here we only need sufficient precision to determine
939 // the sign, knowing the values to be nonzero rationals with small
940 // denominator.
941
942 if(verbose)
943 cout<<"Newform #"<<(i+1)<<": fixing sign normalization using approximate periods"<<endl;
944 nflist[i].sign_normalize();
945 }
946 }
947 // Compute homspace::projcoord, so proj_coords can be used
948 // Replaces coord_vecs of homspace with projections onto eigenspaces
949 // NB if #newforms>1 this MUST be re-called after any sorting of newforms
950 make_projcoord();
951
952 // Look for a j0 such that nflist[i].bplus/bminus[j0]!=0 for all i, or a set of such j
953 find_jlist();
954 }
955
956
find_jlist()957 void newforms::find_jlist()
958 {
959 int i, j, ok=0; j0=0;
960 for(j=1; (!ok)&&(j<=h1->h1dim()); j++)
961 {
962 ok=1;
963 for (i=0; (i<n1ds)&&ok; i++)
964 if(sign==-1)
965 ok=(nflist[i].bminus[j]!=0);
966 else
967 ok=(nflist[i].bplus[j]!=0);
968 if(ok) j0=j;
969 }
970 if(ok)
971 {
972 if(verbose>1) cout<<"j0="<<j0<<endl;
973 jlist.insert(j0);
974 for (i=0; i<n1ds; i++)
975 {
976 nflist[i].j0 = j0;
977 vec& bas = (sign==-1? nflist[i].bminus: nflist[i].bplus);
978 nflist[i].fac = bas[j0];
979 if (verbose>1)
980 {
981 cout<<"Newform #"<<(i+1)<<": bplus = "<<bas<<endl;
982 cout<<" fac = "<<nflist[i].fac<<endl;
983 }
984 }
985 }
986 else
987 {
988 if(verbose)
989 cout<<"Failed to find j0 such that nflist[i].bplus/bminus[j0]!=0 for all i"
990 <<endl;
991 // Find out which pivots we'll be using:
992 for (i=0; i<n1ds; i++)
993 {
994 vec& bas = nflist[i].bplus;
995 j=1; while(bas[j]==0) j++;
996 jlist.insert(j);
997 nflist[i].j0 = j;
998 nflist[i].fac = nflist[i].bplus[j];
999 if (verbose>1)
1000 {
1001 cout<<"Newform #"<<(i+1)<<": bplus = "<<bas<<endl;
1002 cout<<" fac = "<<nflist[i].fac<<endl;
1003 }
1004 }
1005 if(verbose) cout<<"jlist="<<jlist<<endl;
1006 }
1007 }
1008
1009 // Compute homspace::projcoord, so proj_coords can be used
1010 // Replaces coord_vecs of homspace with projections onto eigenspaces
1011 // NB if #newforms>1 this MUST be re-called after any sorting of newforms
make_projcoord()1012 void newforms::make_projcoord()
1013 {
1014 h1->projcoord.init(h1->coord_vecs.size()-1,n1ds);
1015 int j;
1016 if(sign==-1)
1017 for (j=1; j<=n1ds; j++)
1018 h1->projcoord.setcol(j, nflist[j-1].coordsminus);
1019 else
1020 for (j=1; j<=n1ds; j++)
1021 h1->projcoord.setcol(j, nflist[j-1].coordsplus);
1022 }
1023
dimoldpart(const vector<long> l)1024 long newforms::dimoldpart(const vector<long> l)
1025 {
1026 return of->dimoldpart(l);
1027 }
1028
1029 // if(!cuspidal) we really should check here that the basis vector b1
1030 // is in ker(delta), by checking that b1*h1->deltamat == 0
1031
use(const vec & b1,const vec & b2,const vector<long> aplist)1032 void newforms::use(const vec& b1, const vec& b2, const vector<long> aplist)
1033 {
1034 if(basisflag) // we already have all the data except the
1035 // basis vector, so not much needs doing
1036 {
1037 int n = nf_subset[j1ds++];
1038 newform& nf = nflist[n];
1039 if(verbose)
1040 cout<<"Filling in data for for newform #"<<(n+1)<<": bases..."<<flush;
1041 nf.sign=sign;
1042 if(sign==+1)
1043 nf.bplus=b1;
1044 if(sign==-1)
1045 nf.bminus=b1; // formfinder puts the basis vector in b1
1046 if(sign==0)
1047 {
1048 nf.bplus=b1;
1049 nf.bminus=b2;
1050 }
1051 if(verbose)
1052 cout<<"type and cuspidal factors..."<<flush;
1053 nf.find_cuspidal_factors();
1054 if(verbose)
1055 cout<<"coords..."<<flush;
1056 nf.find_coords_plus_minus();
1057 if(sign==0)
1058 {
1059 if(verbose)
1060 cout<<"twisting primes..."<<flush;
1061 nf.find_twisting_primes();
1062 if(verbose)
1063 cout<<"matrix..."<<flush;
1064 nf.find_matrix();
1065 }
1066 if(verbose)
1067 cout<<"done."<<endl;
1068 if(verbose)
1069 cout<<"Finished filling in data for newform #"<<(n+1)<<endl;
1070 return;
1071 }
1072
1073 // Code for initial newform construction
1074
1075 // We use the newform constructor to do all the work, given the basis vector(s) and aplist:
1076
1077 n1ds++;
1078 if(verbose)
1079 {
1080 cout<<"Constructing newform #"<<n1ds<<" with eigs ";
1081 vec_out(cout,aplist,10);
1082 cout<<endl;
1083 }
1084 if(sign==-1)
1085 nflist.push_back(newform(b1,b1,aplist,this)); // only 2nd vector used
1086 else
1087 nflist.push_back(newform(b1,b2,aplist,this));
1088 if(verbose)
1089 cout<<"Finished constructing newform #"<<n1ds<<" with sign = "<<sign<<endl;
1090 }
1091
1092 // Sort newforms
sort(int oldorder)1093 void newforms::sort(int oldorder)
1094 {
1095 if(oldorder)
1096 ::sort(nflist.begin(),nflist.end(),less_newform_old());
1097 else
1098 ::sort(nflist.begin(),nflist.end(),less_newform_new());
1099 }
1100
1101 // Before recovering eigenbases, we need to put back the aq into the
1102 // aplist (and resort, for efficiency).
unfix_eigs()1103 void newforms::unfix_eigs()
1104 {
1105 for(int i=0; i<n1ds; i++)
1106 nflist[i].unfix_eigs();
1107 }
1108
1109 // After recovering eigenbases, we need to refix the aplist
refix_eigs()1110 void newforms::refix_eigs()
1111 {
1112 for(int i=0; i<n1ds; i++)
1113 nflist[i].refix_eigs();
1114 }
1115
display(void) const1116 void newforms::display(void) const
1117 {
1118 if (n1ds==0) {cout<<"No newforms."<<endl; return;}
1119 cout << "\n"<<n1ds<<" newform(s) at level " << modulus << ":" << endl;
1120 cout<<"p0="<<p0<<endl;
1121 // if(dim(mvp)!=0) cout<<"mvp="<<mvp<<endl;
1122 cout<<"#ap=\t"<<nflist[0].aplist.size()<<endl;
1123 long i;
1124 for(i=0; i<n1ds; i++)
1125 {cout<<i+1<<":\t";
1126 nflist[i].display();
1127 }
1128 }
1129
1130 #define DEBUG_SCALING
1131
display_modular_symbol_map(int check) const1132 void newforms::display_modular_symbol_map(int check) const
1133 {
1134 long i,j,k;
1135 rational rplus, rminus;
1136 vector<bigfloat> x0s;
1137 vector<bigfloat> y0s;
1138 bigfloat x0,y0;
1139 if (check)
1140 for(k=0; k<n1ds; k++)
1141 {
1142 cout<<"getting period(s) for newform # "<<(k+1)<<endl;
1143 get_both_periods(k,x0,y0);
1144 cout<<"x0="<<x0<<", y0="<<y0<<endl;
1145 x0s.push_back(x0);
1146 y0s.push_back(y0);
1147 }
1148
1149 for(i=0; i<h1->nsymb; i++)
1150 {
1151 symb s = h1->symbol(i);
1152 modsym ms = modsym(s);
1153 cout<<s<<" = "<<ms<<" -> ";
1154 rational alpha = ms.alpha(), beta = ms.beta();
1155 long a, b, c, d, g=0;
1156 if (num(alpha)==0)
1157 {
1158 b = num(beta);
1159 d = den(beta);
1160 g = bezout(-modulus*b,d,c,a); // so g=a*d-b*N*c
1161 }
1162 j=h1->coordindex[i];
1163 long sg=::sign(j); j=abs(j);
1164 // cout<<"j="<<j<<"("<<sg<<")"<<endl;
1165 if(j==0)
1166 for(k=0; k<n1ds; k++)
1167 if(sign!=0)
1168 cout<<"0 ";
1169 else
1170 cout<<"(0,0) ";
1171 else
1172 for(k=0; k<n1ds; k++)
1173 {
1174 if (check && (g==1))
1175 {
1176 periods_direct integrator(this,&(nflist[k]));
1177 integrator.compute(a,b,c,d);
1178 x0 = integrator.rper();
1179 y0 = integrator.iper();
1180 }
1181 if(sign!=-1)
1182 {
1183 long nrplus = sg*nflist[k].coordsplus[j];
1184 long drplus = nflist[k].cuspidalfactorplus;
1185 rplus = rational(nrplus,drplus);
1186 if (check && (g==1))
1187 {
1188 long n = I2long(Iround(drplus *x0/x0s[k])); // should = nrplus
1189 if (n != nrplus)
1190 {
1191 cout << "plus check fails: rplus = "<<nrplus<<"/"<<drplus<<" = "<<rplus<<endl;
1192 cout << "real part = "<<x0<<endl;
1193 cout << "x0 = "<<x0s[k]<<endl;
1194 cout << "ratio = "<<x0/x0s[k]<<endl;
1195 cout << "scaled ratio = "<<n<<endl;
1196 }
1197 }
1198 rplus *= nflist[k].optimalityfactorplus;
1199 }
1200 if(sign!=+1)
1201 {
1202 long nrminus = sg*nflist[k].coordsminus[j];
1203 long drminus = nflist[k].cuspidalfactorminus;
1204 rminus = rational(nrminus,drminus);
1205 if (check && (g==1))
1206 {
1207 long n = I2long(Iround(drminus *y0/y0s[k])); // should = nrminus
1208 if (n != nrminus)
1209 {
1210 cout << "minus check fails: rminus = "<<rminus<<endl;
1211 cout << "imag part = "<<y0<<endl;
1212 cout << "y0 = "<<y0s[k]<<endl;
1213 cout << "ratio = "<<y0/y0s[k]<<endl;
1214 cout << "scaled ratio = "<<n<<endl;
1215 }
1216 }
1217 rminus *= nflist[k].optimalityfactorminus;
1218 }
1219 if(sign==+1)
1220 cout<<rplus<<" ";
1221 else if(sign==-1)
1222 cout<<rminus<<" ";
1223 else
1224 cout<<"("<<rplus<<","<<rminus<<") ";
1225 }
1226 cout<<endl;
1227 }
1228 }
1229
display(void) const1230 void newform::display(void) const
1231 {
1232 cout << "aplist = ";
1233 vec_out(cout,aplist,20); // outputs at most 20 eigs.
1234 cout<< endl;
1235 // cout << "basis = " << bplus << endl;
1236 cout << "aq = " << aqlist<<endl;
1237 cout << "ap0 = " << ap0
1238 <<", dp0 = " << dp0
1239 <<", np0 = " << np0;
1240 if(pdot!=0) cout <<", pdot = " << pdot;
1241 cout <<endl;
1242 cout << "SFE = " << sfe << ",\tL/P = " << loverp << endl;
1243 if(lplus>0) cout << "lplus = " << lplus << ", mplus = " << mplus << endl;
1244 if(lminus>0) cout << "lminus = " << lminus << ", mminus = " << mminus << endl;
1245 if(a!=0)
1246 {
1247 cout << "[(" <<a<<","<<b<<";"<<c
1248 <<","<<d<<"),"<<dotplus<<","<<dotminus
1249 <<";";
1250 if(type)
1251 cout<<type;
1252 else
1253 cout<<"?";
1254 cout<<"]"<<endl;
1255 }
1256
1257 if(index!=-1)cout << "Splitting index = " << index << endl;
1258 }
1259
putout(ofstream & of,short a,int binflag)1260 void putout(ofstream& of, short a, int binflag)
1261 {
1262 if(binflag)
1263 of.write((char*)&a,sizeof(short));
1264 else
1265 of<<setw(5)<<a;
1266 }
putout(ofstream & of,int a,int binflag)1267 void putout(ofstream& of, int a, int binflag)
1268 {
1269 if(binflag)
1270 of.write((char*)&a,sizeof(int));
1271 else
1272 of<<setw(10)<<a;
1273 }
putout(ofstream & of,long a,int binflag)1274 void putout(ofstream& of, long a, int binflag)
1275 {
1276 if(binflag)
1277 of.write((char*)&a,sizeof(long));
1278 else
1279 of<<setw(15)<<a;
1280 }
nl(ofstream & of,int binflag)1281 void nl(ofstream& of, int binflag)
1282 {if(!binflag) of<<"\n";}
1283
output_to_file(int binflag,int smallflag) const1284 void newforms::output_to_file(int binflag, int smallflag) const
1285 {
1286 long i,j;
1287 char prefix = 'e';
1288 if(binflag) prefix = 'x';
1289 string name = smallflag ? small_nf_filename(modulus,prefix)
1290 : nf_filename(modulus,prefix);
1291 ofstream out(name.c_str());
1292 if(!out)
1293 {
1294 cerr<<"Unable to open file "<<name<<" for newform output"<<endl;
1295 return;
1296 }
1297 // else
1298 // {
1299 // cout<<"--outputting newforms data to "<<name<<" (smallflag="<<smallflag<<")"<<endl;
1300 // display();
1301 // }
1302
1303 if(n1ds==0)
1304 {
1305 putout(out,(int)0,binflag);
1306 putout(out,(int)0,binflag);
1307 putout(out,(int)0,binflag);
1308 out.close();
1309 return;
1310 }
1311
1312 // Line 1: #newforms, #aq, #ap
1313 int nap = nflist[0].aplist.size();
1314 if (smallflag)
1315 {
1316 if (nap>=DEFAULT_SMALL_NAP) nap=DEFAULT_SMALL_NAP;
1317 else
1318 {
1319 cout<<"Warning: small newforms output will only have" << nap
1320 << "a_p (at least " << DEFAULT_SMALL_NAP <<"required" << endl;
1321 }
1322 }
1323 putout(out,(int)n1ds,binflag);
1324 putout(out,(int)nflist[0].aqlist.size(),binflag);
1325 putout(out, nap,binflag);
1326 nl(out,binflag);
1327 // Line 2: blank line
1328 nl(out,binflag);
1329 // Line 3: sign of f.e. for each newform
1330 for(i=0; i<n1ds; i++) putout(out,(int)nflist[i].sfe,binflag);
1331 nl(out,binflag);
1332 // Line 4: ap0 for each newform
1333 for(i=0; i<n1ds; i++) putout(out,(int)nflist[i].ap0,binflag);
1334 nl(out,binflag);
1335 // Line 5: np0 for each newform
1336 for(i=0; i<n1ds; i++) putout(out,(int)nflist[i].np0,binflag);
1337 nl(out,binflag);
1338 // Line 6: dp0 for each newform
1339 for(i=0; i<n1ds; i++) putout(out,(int)nflist[i].dp0,binflag);
1340 nl(out,binflag);
1341 // Line 7: lplus for each newform
1342 for(i=0; i<n1ds; i++) putout(out,(int)nflist[i].lplus,binflag);
1343 nl(out,binflag);
1344 // Line 8: mplus each newform
1345 for(i=0; i<n1ds; i++) putout(out,(int)nflist[i].mplus,binflag);
1346 nl(out,binflag);
1347 // Line 9: lminus each newform
1348 for(i=0; i<n1ds; i++) putout(out,(int)nflist[i].lminus,binflag);
1349 nl(out,binflag);
1350 // Line 10: mminus for each newform
1351 for(i=0; i<n1ds; i++) putout(out,(int)nflist[i].mminus,binflag);
1352 nl(out,binflag);
1353 // Line 11: matrix entry a for each newform
1354 for(i=0; i<n1ds; i++) putout(out,(int)nflist[i].a,binflag);
1355 nl(out,binflag);
1356 // Line 12: matrix entry b for each newform
1357 for(i=0; i<n1ds; i++) putout(out,(int)nflist[i].b,binflag);
1358 nl(out,binflag);
1359 // Line 13: matrix entry c for each newform
1360 for(i=0; i<n1ds; i++) putout(out,(int)nflist[i].c,binflag);
1361 nl(out,binflag);
1362 // Line 14: matrix entry d for each newform
1363 for(i=0; i<n1ds; i++) putout(out,(int)nflist[i].d,binflag);
1364 nl(out,binflag);
1365 // Line 15: dotplus for each newform
1366 for(i=0; i<n1ds; i++) putout(out,(int)nflist[i].dotplus,binflag);
1367 nl(out,binflag);
1368 // Line 16: dotminus for each newform
1369 for(i=0; i<n1ds; i++) putout(out,(int)nflist[i].dotminus,binflag);
1370 nl(out,binflag);
1371 // Line 17: lattice type for each newform
1372 for(i=0; i<n1ds; i++) putout(out,(int)nflist[i].type,binflag);
1373 nl(out,binflag);
1374 // Line 18: deg(phi) for each newform
1375 for(i=0; i<n1ds; i++) putout(out,(int)nflist[i].degphi,binflag);
1376 nl(out,binflag);
1377 // Line 19: blank line
1378 nl(out,binflag);
1379
1380 // Lines 20-(20+#aq): aq for each newform; then blank line
1381 for(j=0; j<int(nflist[0].aqlist.size()); j++)
1382 {
1383 for(i=0; i<n1ds; i++) putout(out,(short)nflist[i].aqlist[j],binflag);
1384 nl(out,binflag);
1385 }
1386 nl(out,binflag);
1387
1388 // Lines (21+#aq)-(20+#aq+#ap): ap for each newform
1389 for(j=0; j<nap; j++)
1390 {
1391 for(i=0; i<n1ds; i++) putout(out,(short)nflist[i].aplist[j],binflag);
1392 nl(out,binflag);
1393 }
1394 out.close();
1395 }
1396
1397 // Read in newform data from file NF_DIR/xN
1398
createfromdata(int s,long ntp,int create_from_scratch_if_absent,int small_data_ok)1399 void newforms::createfromdata(int s, long ntp,
1400 int create_from_scratch_if_absent,
1401 int small_data_ok)
1402 {
1403 sign = s;
1404 long i, j, n = modulus;
1405 if(verbose) cout << "Retrieving newform data for N = " << n << endl;
1406
1407 string name = nf_filename(modulus,'x');
1408 ifstream datafile(name.c_str());
1409 if(small_data_ok && !datafile.is_open())
1410 {
1411 if (verbose)
1412 cout << "No file "<<name<<" exists, trying ";
1413 name = small_nf_filename(modulus,'x');
1414 if (verbose)
1415 cout << name << "instead..."<<endl;
1416 datafile.open(name.c_str());
1417 }
1418 if(!datafile.is_open())
1419 {
1420 if(verbose) cout<<"Unable to open file "<<name<<" for oldform input"<<endl;
1421 if(create_from_scratch_if_absent)
1422 {
1423 if(verbose) cout<<"Creating from scratch instead"<<endl;
1424 createfromscratch(sign, ntp);
1425 output_to_file(1,0); // full newform data
1426 output_to_file(1,1); // short newform data (only DEFAULT_SMALL_NAP ap)
1427 if(verbose) cout << "Finished creating newform data for N = " << n << endl;
1428 if(verbose) display();
1429 }
1430 return;
1431 }
1432
1433 int temp_int;
1434 datafile.read((char*)&temp_int,sizeof(int)); // = number of newforms
1435 n1ds=temp_int;
1436 datafile.read((char*)&temp_int,sizeof(int)); // = number of bad primes
1437 datafile.read((char*)&temp_int,sizeof(int)); // = number of eigs
1438 nap=temp_int;
1439 if(n1ds==0)
1440 {
1441 if(verbose) cout << "No newforms at level " << n << endl;
1442 datafile.close();
1443 return;
1444 }
1445
1446 vector<int> * data = new vector<int>[n1ds];
1447 vector<long> * aq = new vector<long>[n1ds];
1448 vector<long> * ap = new vector<long>[n1ds];
1449
1450 // read extra data for each newform
1451 for(i=0; i<n1ds; i++) data[i].resize(16);
1452 long ntotal = 16*n1ds;
1453 int* batch_i = new int[ntotal];
1454 datafile.read((char*)batch_i,ntotal*sizeof(int));
1455 int *batch_i_ptr = batch_i;
1456 for(j=0; j<16; j++)
1457 for(i=0; i<n1ds; i++)
1458 data[i][j]=*batch_i_ptr++;
1459 delete[] batch_i;
1460 // cout<<"Raw data:\n"; for(i=0; i<n1ds; i++) cout<<data[i]<<endl;
1461
1462 // read aq for each newform
1463 for(i=0; i<n1ds; i++) aq[i].resize(npdivs);
1464 ntotal = npdivs*n1ds;
1465 short* batch = new short[ntotal];
1466 datafile.read((char*)batch,ntotal*sizeof(short));
1467 short *batchptr = batch;
1468 for(j=0; j<npdivs; j++)
1469 for(i=0; i<n1ds; i++)
1470 aq[i][j]=*batchptr++;
1471 // cout<<"Raw aq:\n"; for(i=0; i<n1ds; i++) cout<<aq[i]<<endl;
1472
1473 // read ap for each newform
1474 for(i=0; i<n1ds; i++) ap[i].resize(nap);
1475 ntotal = nap*n1ds;
1476 delete[] batch;
1477 batch = new short[ntotal];
1478 datafile.read((char*)batch,ntotal*sizeof(short));
1479 batchptr = batch;
1480 for(j=0; j<nap; j++)
1481 for(i=0; i<n1ds; i++)
1482 ap[i][j]=*batchptr++;
1483 // cout<<"Raw ap:\n"; for(i=0; i<n1ds; i++) cout<<ap[i]<<endl;
1484 delete[] batch;
1485 datafile.close();
1486
1487 // construct the newforms from this data
1488 nflist.reserve(n1ds);
1489 for(i=0; i<n1ds; i++)
1490 nflist.push_back(newform(data[i],aq[i],ap[i],this));
1491
1492 delete[] ap; delete[] aq; delete[] data;
1493
1494 if(verbose)
1495 {
1496 cout << "Finished reading newform data for N = " << n << endl;
1497 if(verbose>1) display();
1498 }
1499 }
1500
1501
1502 // Create from a list of Hecke eigenvalues from an elliptic curve
1503
eiglist(CurveRed & C,int nap)1504 vector<long> eiglist(CurveRed& C, int nap)
1505 {
1506 long N = I2long(getconductor(C));
1507 long p; bigint pp;
1508 vector<long> ans;
1509 for(primevar pr(nap); pr.ok(); pr++)
1510 {
1511 p=pr; pp=BIGINT(p);
1512 if(N%p==0)
1513 ans.push_back(LocalRootNumber(C,pp));
1514 else
1515 ans.push_back(I2long(Trace_Frob(C,pp)));
1516 }
1517 // cout<<"eiglist("<<(Curve)C<<") = "<<ans<<endl;
1518 return ans;
1519 }
1520
1521 // extract the eigenvalues for bad primes
aqlist(vector<long> aplist,long N)1522 vector<long> aqlist(vector<long> aplist, long N)
1523 {
1524 long iq=0, p, naq = pdivs(N).size();
1525 //cout << "Setting aq of size "<<naq<<endl;
1526 vector<long>::const_iterator api = aplist.begin();
1527 vector<long> aq(naq);
1528 for(primevar pr; (iq<naq)&&pr.ok(); pr++)
1529 {
1530 p=pr;
1531 if (N%p==0)
1532 {
1533 //cout << "Setting aq["<<p<<"] = "<<*api<<endl;
1534 aq[iq++] = *api;
1535 }
1536 api++;
1537 }
1538 return aq;
1539 }
1540
1541 // Create from a list of elliptic curves of the right conductor:
1542
createfromcurve(int s,const CurveRed & C,int nap)1543 void newforms::createfromcurve(int s, const CurveRed& C, int nap)
1544 {
1545 vector<CurveRed> Clist; Clist.push_back(C);
1546 return createfromcurves(s,Clist,nap);
1547 }
createfromcurves(int s,vector<CurveRed> Clist,int nap)1548 void newforms::createfromcurves(int s, vector<CurveRed> Clist, int nap)
1549 {
1550 if(verbose) cout << "In newforms::createfromcurves()..."<<endl;
1551 sign=s;
1552 int ncurves = Clist.size();
1553 if(ncurves==0) return;
1554 if(verbose) cout << "Making homspace..."<<flush;
1555 makeh1(sign);
1556 if(verbose) cout << "done." << endl;
1557 mvp=h1->maninvector(p0);
1558 long nap_default = long(100*sqrt(modulus));
1559 if (nap<nap_default)
1560 {
1561 if(verbose)
1562 {
1563 cout << "--increasing nap from "<< nap << " to " << nap_default << endl;
1564 }
1565 nap=nap_default;
1566 }
1567 if(verbose) cout << "Making form_finder (nap="<<nap<<")..."<<flush;
1568 form_finder splitspace(this, (sign!=0), nap, 0, 1, 0, verbose);
1569 if(verbose) cout << "Recovering eigenspace bases with form_finder..."<<endl;
1570 // j1ds counts through the newforms as they are found
1571 basisflag=0; j1ds=0;
1572 vector< vector<long> > eigs(ncurves);
1573 int i;
1574
1575 for(i=0; i<ncurves; i++)
1576 eigs[i]=eiglist(Clist[i],nap);
1577 n1ds=0; nflist.resize(0);
1578 splitspace.recover(eigs); // NB newforms::use() determines what is
1579 // done with each one as it is found;
1580 // this depends on basisflag and sign
1581
1582 // Get period lattice of each newform (and hence of optimal curves)
1583 for(i=0; i<ncurves; i++)
1584 {
1585 if(verbose)
1586 cout<<"Finding optimality scaling factors using approximate periods"<<endl;
1587 nflist[i].find_optimality_factors(Clist[i], i);
1588 if(verbose)
1589 cout<<"Fixing sign normalization using approximate periods"<<endl;
1590 nflist[i].sign_normalize();
1591 }
1592 if(verbose) cout << "...done."<<endl;
1593 }
1594
createfromcurves_mini(vector<CurveRed> Clist,int nap)1595 void newforms::createfromcurves_mini(vector<CurveRed> Clist, int nap)
1596 {
1597 if(verbose) cout << "In newforms::createfromcurves_mini()..."<<endl;
1598 int i; long N;
1599 n1ds = Clist.size();
1600 nflist.reserve(n1ds);
1601
1602 if (n1ds>0) // construct the ap and aq vectors from the curves
1603 {
1604 N = I2long(getconductor(Clist[0]));
1605 for(i=0; i<n1ds; i++)
1606 {
1607 vector<long> ap=eiglist(Clist[i],nap);
1608 vector<long> aq=aqlist(ap,N);
1609 // dummy data -- these fields are not set and will not be used
1610 vector<int> data(16,0);
1611 newform nf(data,aq,ap,this);
1612 if (verbose)
1613 {
1614 cout<<"adding this newform: "<<endl;
1615 nf.display();
1616 }
1617 nflist.push_back(newform(data,aq,ap,this));
1618 }
1619 }
1620 if(verbose) cout << "...done."<<endl;
1621 }
1622
1623 // Read in newform data from old-style data files eigs/xN and intdata/[ep]N
1624
createfromolddata()1625 void newforms::createfromolddata()
1626 {
1627 long i, j, n = modulus;
1628 if(verbose)
1629 cout << "Retrieving old-style newform data for N = " << n << endl;
1630
1631 stringstream eigsname;
1632 eigsname << "eigs/x" << n;
1633 ifstream eigsfile(eigsname.str().c_str());
1634
1635 if(!eigsfile.is_open())
1636 {
1637 cout<<"Unable to open file "<<eigsname.str()<<" for eigs input"<<endl;
1638 return;
1639 }
1640
1641 short temp;
1642 eigsfile.read((char*)&temp,sizeof(short)); // # newforms
1643 n1ds=temp;
1644 eigsfile.read((char*)&temp,sizeof(short)); // # irrational newforms
1645 eigsfile.read((char*)&temp,sizeof(short)); // # ap
1646 nap=temp;
1647 if(n1ds==0)
1648 {
1649 if(verbose) cout << "No newforms at level " << n << endl;
1650 eigsfile.close();
1651 return;
1652 }
1653
1654 // read ap for each newform
1655 vector<long> * ap = new vector<long>[n1ds];
1656 for(i=0; i<n1ds; i++) ap[i].resize(nap);
1657 long ntotal = nap*n1ds;
1658 short* batch = new short[ntotal];
1659 eigsfile.read((char*)batch,ntotal*sizeof(short));
1660 eigsfile.close();
1661 short* batchptr = batch;
1662 for(j=0; j<nap; j++)
1663 for(i=0; i<n1ds; i++)
1664 ap[i][j]=*batchptr++;
1665 // cout<<"Raw ap:\n"; for(i=0; i<n1ds; i++) cout<<ap[i]<<endl;
1666 delete[] batch;
1667
1668 // extract aq for each newform
1669 vector<long> * aq = new vector<long>[n1ds];
1670 for(i=0; i<n1ds; i++) aq[i].resize(npdivs);
1671 primevar pr; long q, k, a;
1672 for(k=0, j=0; j<int(plist.size()); j++)
1673 {
1674 q = plist[j];
1675 int q2divN = ::divides(q*q,modulus);
1676 while((long)pr!=q) {pr++; k++;}
1677 for(i=0; i<n1ds; i++)
1678 {
1679 a = ap[i][k];
1680 aq[i][j] = a;
1681 ap[i][k] = (q2divN? 0: -a);
1682 }
1683 }
1684 // cout<<"Raw aq:\n"; for(i=0; i<n1ds; i++) cout<<aq[i]<<endl;
1685
1686 // read extra data for each newform
1687 vector<int> * data = new vector<int>[n1ds];
1688 for(i=0; i<n1ds; i++) data[i].resize(16);
1689 stringstream intdataname;
1690 intdataname << "intdata/e" << n;
1691 ifstream intdatafile(intdataname.str().c_str());
1692 if(!intdatafile.is_open())
1693 {
1694 cout<<"Unable to open data file "<<intdataname.str()<<" for data input"<<endl;
1695 return;
1696 }
1697
1698 long nloverp, dloverp, dp0, np0;
1699 for(i=0; i<n1ds; i++)
1700 {
1701 // cout<<"Reading intdata for form #"<<(i+1)<<endl;
1702 intdatafile >> data[i][15]; // degphi
1703 // cout<<"degphi = "<<data[i][15]<<endl;
1704 intdatafile >> data[i][0]; // sfe
1705 // cout<<"sfe = "<<data[i][0]<<endl;
1706 intdatafile >> nloverp; // num(L/P)
1707 intdatafile >> dloverp; // den(L/P)
1708 intdatafile >> data[i][14]; // type
1709 // cout<<"type = "<<data[i][14]<<endl;
1710 intdatafile >> dp0; data[i][3]=dp0; // dp0
1711 // cout<<"dp0 = "<<data[i]=dp0[3]<<endl;
1712 intdatafile >> np0; data[i][2]=np0; // np0
1713 // cout<<"np0 = "<<data[i][2]<<endl;
1714 if(dp0==0) data[i][3]=(2*nloverp*np0)/dloverp;
1715 data[i][1]=1+p0-data[i][2]; // ap0 (not in intdata file)
1716 // cout<<"ap0 = "<<data[i][1]<<endl;
1717 intdatafile >> data[i][4]; // lplus
1718 intdatafile >> data[i][5]; // mplus
1719 intdatafile >> data[i][6]; // lminus
1720 intdatafile >> data[i][7]; // mminus
1721 intdatafile >> data[i][8]; // a
1722 intdatafile >> data[i][9]; // b
1723 intdatafile >> data[i][10]; // c
1724 intdatafile >> data[i][11]; // d
1725 intdatafile >> data[i][12]; // dotplus
1726 intdatafile >> data[i][13]; // dotminus
1727 }
1728 intdatafile.close();
1729 // cout<<"Raw data:\n"; for(i=0; i<n1ds; i++) cout<<data[i]<<endl;
1730
1731
1732 // construct the newforms from this data
1733 nflist.reserve(n1ds);
1734 for(i=0; i<n1ds; i++)
1735 nflist.push_back(newform(data[i],aq[i],ap[i],this));
1736 // delete ap; delete aq; delete data;
1737
1738 if(verbose)
1739 {
1740 cout << "Finished reading oldstyle newform data for N = " << n << endl;
1741 display();
1742 }
1743 }
1744
1745 // Construct bases (homology eigenvectors) from eigenvalue lists:
makebases(int flag,int all_nf)1746 void newforms::makebases(int flag, int all_nf)
1747 {
1748 if(n1ds==0) return;
1749
1750 if(((sign==-1)||(dim(nflist[0].bplus)>0)) &&
1751 ((sign==+1)||(dim(nflist[0].bminus)>0)))
1752 return;
1753 if(verbose) cout << "Making homspace..."<<flush;
1754 makeh1(sign);
1755 if(verbose) cout << "done." << endl;
1756 mvp=h1->maninvector(p0);
1757 if(verbose) cout << "Making form_finder (nap="<<nap<<")..."<<flush;
1758 form_finder splitspace(this, (sign!=0), nap, 0, 1, 0, verbose);
1759 if(verbose) cout << "Recovering eigenspace bases with form_finder..."<<endl;
1760 // basisflag controls what ::use() does with the nfs when found
1761 // j1ds counts through the newforms as they are found
1762 basisflag=flag;
1763 int i;
1764 j1ds = 0; // counts through newforms as they are recovered
1765 vector< vector<long> > eigs;
1766
1767 if (all_nf)
1768 {
1769 nf_subset.clear();
1770 for(i=0; i<n1ds; i++)
1771 nf_subset.push_back(i);
1772 }
1773
1774 unfix_eigs();
1775 sort();
1776 for(i=0; i<nf_subset.size(); i++)
1777 eigs.push_back(nflist[nf_subset[i]].aplist);
1778
1779 splitspace.recover(eigs); // NB newforms::use() determines what is
1780 // done with each one as it is found;
1781 // this depends on basisflag and sign
1782 if(verbose) cout << "...done."<<endl;
1783 refix_eigs();
1784 if(verbose>1) cout<<"Reordering newforms after recovery"<<endl;
1785 if(verbose>1) {cout<<"Before sorting:\n"; display();}
1786 sort(int(modulus<130000)); // old order for N<130000, else new order
1787 if(verbose>1) {cout<<"After sorting:\n"; display();}
1788 }
1789
merge(int all_nf)1790 void newforms::merge(int all_nf)
1791 {
1792 if(n1ds==0) return;
1793 if(verbose) cout << "Making homspace..."<<flush;
1794 makeh1(0);
1795 if(verbose) cout << "done." << endl;
1796 vec bplus, bminus;
1797 j1ds = 0;
1798 basisflag = 1;
1799 mvlplusvecs.clear();
1800 mvlminusvecs.clear();
1801 if (verbose>1)
1802 cout<<"merging newforms " << nf_subset << endl;
1803 int inf, jnf;
1804 unfix_eigs();
1805 sort();
1806 for(jnf=0; jnf<nf_subset.size(); jnf++)
1807 {
1808 inf = nf_subset[jnf];
1809 if(verbose) cout << "Newform #"<<(inf+1)<<":"<<endl;
1810 if(verbose) cout << "-about to extend bplus,bminus..."<<flush;
1811 bplus.init(h1->nsymb);
1812 bminus.init(h1->nsymb);
1813 int i,j;
1814 for(i=1; i<=h1->nsymb; i++)
1815 {
1816 j = h1plus->coordindex[i-1];
1817 if (j==0) bplus[i] = 0;
1818 else if (j>0) bplus[i] = nflist[inf].coordsplus[j];
1819 else if (j<0) bplus[i] = -nflist[inf].coordsplus[-j];
1820 j = h1minus->coordindex[i-1];
1821 if (j==0) bminus[i] = 0;
1822 else if (j>0) bminus[i] = nflist[inf].coordsminus[j];
1823 else if (j<0) bminus[i] = -nflist[inf].coordsminus[-j];
1824 }
1825 if(verbose) cout<< "done, about to contract bplus,bminus..."<<flush;
1826 bplus = h1->contract_coords(bplus);
1827 bplus /= vecgcd(bplus);
1828 bminus = h1->contract_coords(bminus);
1829 bminus /= vecgcd(bminus);
1830 if(verbose) cout<< "done."<<endl;
1831 if(verbose>1)
1832 {
1833 cout << " new bplus = "<<bplus <<":"<<endl;
1834 cout << " new bminus = "<<bminus<<":"<<endl;
1835 }
1836 // These new dual eigenvectors are used to compute all
1837 // additional data needed for curve and modular symbol
1838 // computation (scaling and cuspidal factors and type)
1839 use(bplus, bminus, nflist[inf].aplist);
1840 }
1841 refix_eigs();
1842 sort(int(modulus<130000)); // old order for N<130000, else new order
1843 }
1844
1845
update(const mat & pcd,vec & imagej,long ind)1846 void update(const mat& pcd, vec& imagej, long ind)
1847 {
1848 if(ind>0) imagej.add_row(pcd,ind);
1849 else
1850 if(ind<0) imagej.sub_row(pcd,-ind);
1851 }
1852
apvec(long p)1853 vector<long> newforms::apvec(long p) // computes a[p] for each newform
1854 {
1855 //cout<<"In apvec with p = "<<p<<endl;
1856 vector<long> apv(n1ds);
1857 vec v;
1858 long i,j,iq,ap;
1859 if(::divides(p,modulus)) // we already have all the aq
1860 {
1861 if(::divides(p*p,modulus))
1862 for (i=0; i<n1ds; i++) apv[i] = 0;
1863 else
1864 {
1865 iq = find(plist.begin(),plist.end(),p)-plist.begin();
1866 for (i=0; i<n1ds; i++) apv[i] = -nflist[i].aqlist[iq];
1867 }
1868 return apv;
1869 }
1870
1871 // now p is a good prime
1872
1873 long maxap=(long)(2*sqrt((double)p)); // for validity check
1874
1875 map<long,vec> images; // [j,v] stores image of j'th M-symbol in v
1876 // (so we don't compute any more than once)
1877 vec bas, imagej;
1878 long p2=(p-1)>>1; // (p-1)/2
1879 long sg, a, b, c, q, r;
1880 long u1,u2,u3;
1881 long ind;
1882
1883 // Compute the image of the necessary M-symbols (hopefully only one)
1884 //cout<<"Computing images of M-symbols"<<endl<<flush;
1885 //cout<<"jlist = "<<jlist<<endl;
1886
1887 for(std::set<long>::const_iterator jj=jlist.begin(); jj!=jlist.end(); jj++)
1888 {
1889 imagej=vec(n1ds); // initialised to 0
1890 j=*jj;
1891 symb s = h1->symbol(h1->freegens[j-1]);
1892 //cout<<"Computing image of "<<j<<"'th M-symbol "<<s<<endl;
1893 //cout<<" = "<<s<<"..."<<flush;
1894 long u=s.cee(),v=s.dee();
1895 mat& pcd = h1->projcoord;
1896 //cout<<"projcoord = "<<pcd;
1897 // Matrix [1,0;0,p]
1898 ind = h1->coordindex[h1->index2(u,p*v)];
1899 update(pcd,imagej,ind);
1900 //cout<<"(1) (u1,u2)=("<<u<<","<<p*v<<") partial image index is "<<ind<<", subtotal="<<imagej<<endl;
1901 // Matrix [p,0;0,1]
1902 ind = h1->coordindex[h1->index2(p*u,v)];
1903 update(pcd,imagej,ind);
1904 //cout<<"(2) (u1,u2)=("<<p*u<<","<<v<<") partial image index is "<<ind<<", subtotal="<<imagej<<endl;
1905 // Other matrices
1906 for(sg=0; sg<2; sg++) // signs
1907 for(r=1; r<=p2; r++)
1908 {
1909 a = -p;
1910 b = sg ? -r : r ;
1911 u1=u*p; u2=v-u*b;
1912 ind = h1->coordindex[h1->index2(u1,u2)];
1913 update(pcd,imagej,ind);
1914 //cout<<"(3) (u1,u2)=("<<u1<<","<<u2<<") partial image index is "<<ind<<", subtotal="<<imagej<<endl;
1915 while(b!=0)
1916 {
1917 c=mod(a,b); q=(a-c)/b;
1918 if(q==1) {u3= u2-u1;} else {u3=q*u2-u1;}
1919 a=-b; b=c; u1=u2; u2=u3;
1920 ind = h1->coordindex[h1->index2(u1,u2)];
1921 update(pcd,imagej,ind);
1922 //cout<<"(4) (u1,u2)=("<<u1<<","<<u2<<") partial image index is "<<ind<<", subtotal="<<imagej<<endl;
1923 }
1924 }
1925 images[j]=imagej;
1926 //cout<<" image is "<<imagej<<endl;
1927 }
1928
1929 for (i=0; i<n1ds; i++)
1930 {
1931 // recover eigenvalue:
1932 //cout<<"Numerator = "<< images[nflist[i].j0][i+1] <<endl;
1933 //cout<<"Denominator = "<< nflist[i].fac << endl;
1934 ap = images[nflist[i].j0][i+1]/nflist[i].fac;
1935 ap *= (sign==-1? nflist[i].contminus: nflist[i].contplus);
1936 ap /= h1->h1denom();
1937 apv[i]=ap;
1938 // check it is in range:
1939 if((ap>maxap)||(-ap>maxap))
1940 {
1941 cout<<"Error: eigenvalue "<<ap<<" for p="<<p
1942 <<" for form # "<<(i+1)<<" is outside valid range "
1943 <<-maxap<<"..."<<maxap<<endl;
1944 break; // no point in trying to compute any more.
1945 }
1946 }
1947 return apv;
1948 }
1949
addap(long last)1950 void newforms::addap(long last) // adds ap for primes up to the last'th prime
1951 {
1952 if(n1ds==0) return;
1953 long i, j, p;
1954 j=0;
1955 if(verbose>1) // output the ap already known...
1956 for(primevar pr(nflist[0].aplist.size()); pr.ok(); pr++, j++)
1957 {
1958 p=(long)pr;
1959 if(ndivides(p,modulus)) cout<<"p="; else cout<<"q=";
1960 cout<<p<<":\t";
1961 {
1962 for (i=0; i<n1ds; i++) cout<<nflist[i].aplist[j]<<"\t";
1963 cout<<endl;
1964 }
1965 }
1966 // Now compute and output the rest of the ap...
1967 for(primevar pr(last,1+nflist[0].aplist.size()); pr.ok(); pr++)
1968 {
1969 p=(long)pr;
1970 vector<long> apv=apvec(p);
1971 if(verbose>1)
1972 {
1973 if(ndivides(p,modulus)) cout<<"p="; else cout<<"q=";
1974 cout<<p<<":\t";
1975 for (i=0; i<n1ds; i++) cout<<apv[i]<<"\t";
1976 cout<<endl;
1977 }
1978 for (long i=0; i<n1ds; i++) nflist[i].aplist.push_back(apv[i]);
1979 }
1980 }
1981
output_to_file_no_newforms(long n,int binflag,int smallflag)1982 void output_to_file_no_newforms(long n, int binflag, int smallflag)
1983 {
1984 char prefix = 'e';
1985 if(binflag) prefix = 'x';
1986 string name = smallflag ? small_nf_filename(n,prefix)
1987 : nf_filename(n,prefix);
1988 ofstream out(name.c_str());
1989 if(binflag)
1990 {
1991 int a=0;
1992 out.write((char*)&a,sizeof(int));
1993 out.write((char*)&a,sizeof(int));
1994 out.write((char*)&a,sizeof(int));
1995 }
1996 else
1997 {
1998 out<<"0 0 0\n";
1999 }
2000 out.close();
2001
2002 }
2003
2004 // for the i'th newform return the value of the modular symbol {0,r} (default) or {oo,r}
plus_modular_symbol(const rational & r,long i,int base_at_infinity) const2005 rational newforms::plus_modular_symbol(const rational& r, long i, int base_at_infinity) const
2006 {
2007 rational a(h1->nfproj_coords(num(r),den(r),nflist[i].coordsplus),
2008 nflist[i].cuspidalfactorplus);
2009 // {oo,r} = {0,r}+{oo,0} and loverp={oo,0} (not {0,oo}!)
2010 if (base_at_infinity) a+=nflist[i].loverp;
2011 a *= nflist[i].optimalityfactorplus;
2012 return a;
2013 }
2014
minus_modular_symbol(const rational & r,long i,int base_at_infinity) const2015 rational newforms::minus_modular_symbol(const rational& r, long i, int base_at_infinity) const
2016 {
2017 // Ignore the value of base_at_infinity as it does not affect the minus symbol
2018 rational a(h1->nfproj_coords(num(r),den(r),nflist[i].coordsminus),
2019 nflist[i].cuspidalfactorminus);
2020 a *= nflist[i].optimalityfactorminus;
2021 return a;
2022 }
2023
full_modular_symbol(const rational & r,long i,int base_at_infinity) const2024 pair<rational,rational> newforms::full_modular_symbol(const rational& r, long i, int base_at_infinity) const
2025 {
2026 mat m(h1->coord_vecs.size()-1,2);
2027 m.setcol(1,nflist[i].coordsplus);
2028 m.setcol(2,nflist[i].coordsminus);
2029 vec a = h1->proj_coords(num(r),den(r),m);
2030 rational a1(a[1],nflist[i].cuspidalfactorplus);
2031 // {oo,r} = {0,r}+{oo,0} and loverp={oo,0} (not {0,oo}!)
2032 if (base_at_infinity) a1 += nflist[i].loverp;
2033 a1 *= nflist[i].optimalityfactorplus;
2034 rational a2(a[2],nflist[i].cuspidalfactorminus);
2035 a2 *= nflist[i].optimalityfactorminus;
2036 return pair<rational,rational> ( a1, a2 );
2037 }
2038
2039 // Attempt to compute and display the elliptic curve for each
2040 // newform; return a list of newform indices where this failed.
showcurves(vector<int> forms,int verbose,string filename)2041 vector<int> newforms::showcurves(vector<int> forms, int verbose, string filename)
2042 {
2043 if((verbose>1)&&(sqfac>1)) cout<<"c4 factor " << sqfac << endl;
2044
2045 ofstream curve_out;
2046 int output_curves = (filename!="no");
2047 if (output_curves) curve_out.open(filename.c_str());
2048 bigfloat rperiod;
2049 bigint a1,a2,a3,a4,a6, N;
2050 vector<int> badcurves; // will hold the indices of forms for which we fail to find a curve
2051 vector<int>::const_iterator inf; // will iterate through the forms to be used
2052
2053 for(inf=forms.begin(); inf!=forms.end(); inf++)
2054 {
2055 if(verbose)
2056 cout<<"\n"<<"Form number "<<*inf+1<<"\n";
2057 else cout<<(*inf+1)<<" ";
2058
2059 if (output_curves)
2060 curve_out << modulus << " "<< codeletter(*inf) << " 1 ";
2061
2062 Curve C = getcurve(*inf,-1,rperiod,verbose);
2063 Curvedata CD(C,1); // The 1 causes minimalization
2064 if(verbose) cout << "\nCurve = \t";
2065 cout << (Curve)CD << "\t";
2066 CurveRed CR(CD);
2067 N = getconductor(CR);
2068 cout << "N = " << N << endl;
2069 if(verbose) cout<<endl;
2070
2071 if(N!=modulus)
2072 {
2073 cout<<"No curve found"<<endl;
2074 badcurves.push_back(*inf);
2075 if (output_curves)
2076 curve_out<<endl;
2077 }
2078 else
2079 if (output_curves)
2080 {
2081 C.getai(a1,a2,a3,a4,a6);
2082 curve_out<<"["<<a1<<","<<a2<<","<<a3<<","<<a4<<","<<a6<<"]";
2083 int nt = CD.get_ntorsion();
2084 int r = nflist[*inf].rank(); // analytic rank
2085 curve_out<<" "<<r<<" "<<nt<<" 0"<<endl;
2086 }
2087 }
2088 if (output_curves)
2089 curve_out.close();
2090
2091 return badcurves;
2092 }
2093