1 // FILE HOMSPACE.CC: Implemention of class homspace
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 #define USE_SMATS // Warning: no longer testing without this switched on!
25
26 #include <eclib/cusp.h>
27 #include <eclib/homspace.h>
28 #include <eclib/timer.h>
29
30 #ifdef USE_SMATS
31 #include <eclib/smatrix_elim.h>
32 #endif
33
34 const string W_opname("W");
35 const string T_opname("T");
36
matop(long a,long b,long c,long d)37 matop::matop(long a, long b, long c, long d)
38 {
39 mats.push_back(mat22(a,b,c,d));
40 }
41
matop(long p,long n)42 matop::matop(long p, long n)
43 {
44 if (p==n)
45 {
46 mats.push_back(mat22(0,-1,n,0));
47 return;
48 }
49 if ((n%p)==0) // W involution, 1 term
50 {
51 long u,v,a,b;
52 for (u=1, v=n; v%p==0; v/=p, u*=p) ;
53 bezout(u,v,a,b);
54 mats.push_back(mat22(u*a,-b,n,u));
55 return;
56 }
57 // else Hecke operator, p+1 terms
58 {
59 mats.resize(p+1);
60 long j, p2 = p>>1;
61 for (j=0; j<p; j++) mats[j] = mat22(1,j-p2,0,p);
62 mats[p] = mat22(p,0,0,1);
63 }
64 }
65
homspace(long n,int hp,int hcusp,int verbose)66 homspace::homspace(long n, int hp, int hcusp, int verbose) :symbdata(n)
67 {
68 init_time();
69 plusflag=hp;
70 cuspidal=hcusp;
71 long i,j,k;
72 coordindex = new int[nsymb];
73 if (!coordindex)
74 {
75 out_of_memory_error("coordindex");
76 return;
77 }
78 int* check = new int[nsymb];
79 if (!check)
80 {
81 out_of_memory_error("check");
82 return;
83 }
84 i=nsymb; while (i--) check[i]=0;
85 long ngens=0;
86 int* gens = new int[nsymb+1]; // N.B. Start of gens array is at 1 not 0
87 if (!gens)
88 {
89 out_of_memory_error("gens");
90 return;
91 }
92 long* rel = new long[4];
93 if (!rel)
94 {
95 out_of_memory_error("rel");
96 return;
97 }
98
99 // 2-term relations:
100
101 // if (plusflag==1)
102 // for (j=0; j<nsymb; j++)
103 // {if (check[j]==0)
104 // { rel[0]=j;
105 // rel[1]=rsof(j);
106 // rel[2]=sof(j);
107 // rel[3]=sof(rel[1]);
108 // if (verbose>1)
109 // cout << "Relation: " << rel[0]<<" "<<rel[1]<<" "<<rel[2]<<" "<<rel[3]<<endl;
110 // for (k=0; k<4; k++) check[rel[k]]=1;
111 // if ( (j==rel[2]) || (j==rel[3]) )
112 // for (k=0; k<4; k++) coordindex[rel[k]]=0;
113 // else
114 // { ngens++;
115 // gens[ngens] = j;
116 // if (verbose>1) cout << "gens["<<ngens<<"]="<<j<<endl;
117 // coordindex[rel[0]] = ngens;
118 // coordindex[rel[1]] = ngens;
119 // coordindex[rel[2]] = -ngens;
120 // coordindex[rel[3]] = -ngens;
121 // }
122 // }
123 // }
124 if (plusflag!=0)
125 for (j=0; j<nsymb; j++)
126 {if (check[j]==0)
127 { rel[0]=j;
128 if (plusflag==-1)
129 rel[1]=rof(j);
130 else
131 rel[1]=rsof(j);
132 rel[2]=sof(j);
133 rel[3]=sof(rel[1]);
134 if (verbose>1)
135 cout << "Relation: " << rel[0]<<" "<<rel[1]<<" "<<rel[2]<<" "<<rel[3]<<endl;
136 for (k=0; k<4; k++) check[rel[k]]=1;
137 if ( (j==rel[2]) || (j==rel[3]) )
138 for (k=0; k<4; k++) coordindex[rel[k]]=0;
139 else
140 { ngens++;
141 gens[ngens] = j;
142 if (verbose>1) cout << "gens["<<ngens<<"]="<<j<<endl;
143 coordindex[rel[0]] = ngens;
144 coordindex[rel[1]] = ngens;
145 coordindex[rel[2]] = -ngens;
146 coordindex[rel[3]] = -ngens;
147 }
148 }
149 }
150 if (plusflag==0)
151 {for (j=0; j<nsymb; j++)
152 {if (check[j]==0)
153 {rel[0]=j;
154 rel[1]=sof(j);
155 check[rel[0]] = check[rel[1]] = 1;
156 if (j==rel[1])
157 for (k=0; k<2; k++) coordindex[rel[k]]=0;
158 else
159 { ngens++;
160 gens[ngens] = j;
161 coordindex[rel[0]] = ngens;
162 coordindex[rel[1]] = -ngens;
163 }
164 }
165 }
166 }
167
168
169 // end of 2-term relations
170 if (verbose)
171 {cout << "After 2-term relations, ngens = "<<ngens<<"\n";
172 // Compare with predicted value:
173 /*
174 int nu2=(::divides((long)4,modulus)?0:1);
175 static int nu2table[4] = {0,2,1,0};
176 for(i=0; nu2&&(i<npdivs); i++) nu2 *= nu2table[plist[i]%4];
177 int ngens0=(nsymb-nu2)/2;
178 cout<<"predicted value of ngens = "<<ngens0;
179 if(!plusflag) if(ngens!=ngens0) cout<<" --WRONG!";
180 cout<<endl;
181 */
182 if (verbose>1)
183 {
184 cout << "gens = ";
185 for (i=1; i<=ngens; i++) cout << gens[i] << " ";
186 cout << "\n";
187 cout << "coordindex = ";
188 for (i=0; i<nsymb; i++) cout << coordindex[i] << " ";
189 cout << "\n";
190 }}
191 //
192 // 3-term relations
193
194 // long maxnumrel = 20+(2*ngens)/3;
195 long maxnumrel = ngens+10;
196
197 if (verbose)
198 {
199 cout << "ngens = "<<ngens<<endl;
200 cout << "maxnumrel = "<<maxnumrel<<endl;
201 cout << "relation matrix has = "<<(maxnumrel*ngens)<<" entries..."<<endl;
202 }
203 {
204 #ifdef USE_SMATS
205 smat relmat((int)maxnumrel,(int)ngens);
206 svec newrel(ngens);
207 #else
208 mat relmat(maxnumrel,ngens);
209 vec newrel(ngens);
210 if (verbose) cout << "successfully allocated "<<endl;
211 #endif
212 int numrel = 0;
213 long ij; int fix;
214
215 for (i=0; i<nsymb; i++) check[i]=0;
216 for (k=0; k<nsymb; k++)
217 {
218 if (check[k]==0)
219 {
220 #ifdef USE_SMATS
221 newrel.clear();
222 #else
223 for (j=1; j<=ngens; j++) newrel[j]=0;
224 #endif
225 rel[2]=tof(rel[1]=tof(rel[0]=k));
226 if (verbose>1) cout << "Relation: " << rel[0]<<" "<<rel[1]<<" "<<rel[2]<<" --> ";
227 for (j=0; j<3; j++)
228 {ij = rel[j];
229 check[ij] = 1;
230 if (plusflag) check[rof(ij)] = 1;
231 fix = coordindex[ij];
232 #ifdef USE_SMATS
233 if(fix) newrel.add(abs(fix),(fix>0?1:-1));
234 #else
235 if (verbose>1) cout << fix << " ";
236 if (fix!=0) newrel[abs(fix)] += sign(fix);
237 #endif
238 }
239 #ifdef USE_SMATS
240 if(verbose>1) cout<<newrel<<"\n";
241 #else
242 if(verbose>1) cout<<"\n";
243 #endif
244 #ifdef USE_SMATS
245 if(newrel.size()!=0)
246 {
247 numrel++;
248 make_primitive(newrel);
249 if(numrel<=maxnumrel)
250 relmat.setrow(numrel,newrel);
251 else
252 cout<<"Too many 3-term relations (numrel="<<numrel
253 <<", maxnumrel="<<maxnumrel<<")"<<endl;
254 }
255 #else
256 if (verbose>1) cout << newrel << "\n";
257 long h = vecgcd(newrel);
258 if (h!=0)
259 { if (h>1) newrel/=h;
260 if(numrel<maxnumrel)
261 relmat.setrow(++numrel,newrel);
262 else cout<<"Too many 3-term relations!"<<endl;
263 }
264 #endif
265 }
266 }
267 if (verbose)
268 {
269 cout << "Finished 3-term relations: numrel = "<<numrel<<" ( maxnumrel = "<<maxnumrel<<")"<<endl;
270 // Compare with predicted value (for full H_1 only):
271 /*
272 if(!plusflag)
273 {
274 int nu3=(::divides((long)9,modulus)?0:1);
275 static int nu3table[3] = {1,2,0};
276 for(i=0; nu3&&(i<npdivs); i++) nu3 *= nu3table[plist[i]%3];
277 int ntriangles=(nsymb+2*nu3)/3;
278 cout<<"predicted value of ntriangles = "<<ntriangles;
279 if(ntriangles!=numrel) cout<<" --WRONG!";
280 cout<<endl;
281 }
282 */
283 }
284
285
286 // end of 3-term relations
287
288 if(verbose)
289 {
290 #ifdef USE_SMATS
291 cout << "relmat has "<< get_population(relmat)<<" nonzero entries (density = "<<density(relmat)<<")"<<endl;
292 #endif
293
294 if(verbose>1)
295 cout << "relmat = " << relmat.as_mat().slice(numrel,ngens) << endl;
296 cout << "Computing kernel..."<<endl;
297 }
298 vec pivs, npivs;
299 #ifdef USE_SMATS
300 smat_elim sme(relmat);
301 //relmat=smat(0,0); // clear space
302 int d1;
303 start_time();
304 smat sp;
305 int ok = liftmat(sme.kernel(npivs,pivs),MODULUS,sp,d1);
306 stop_time();
307 if (!ok)
308 cout << "**!!!** failed to lift modular kernel of relation matrix\n" << endl;
309 else
310 if (verbose) {cout<<"time to compute kernel = "; show_time(); cout<<endl;}
311 denom1=d1;
312 if(verbose>1)
313 {
314 cout << "kernel of relmat = " << sp.as_mat() << endl;
315 cout << "pivots = "<<pivs << endl;
316 cout << "denom = "<<d1 << endl;
317 }
318 rk = sp.ncols();
319 coord_vecs.resize(ngens+1); // 0'th is unused
320 for(i=1; i<=ngens; i++)
321 coord_vecs[i]=sp.row(i);
322 //sp=smat(0,0); // clear space
323 #else
324 subspace sp = kernel(relmat);
325 rk = dim(sp);
326 coord = basis(sp);
327 pivs = pivots(sp);
328 denom1 = denom(sp);
329 for(i=1; i<=ngens; i++)
330 coord_vecs[i]=svec(coord.row(i));
331 sp.clear();
332 relmat.init(); newrel.init();
333 #endif
334
335 // cout<<"ngens = "<<ngens<<endl;
336 if (verbose)
337 {
338 cout << "rk = " << rk << endl;
339 if (verbose>1)
340 {
341 // cout << "coord:\n" ; coord.output_pretty();
342 cout << "coord_vecs:\n";
343 for(i=1; i<=ngens; i++)
344 cout<<i<<": "<<coord_vecs[i].as_vec()<<"\n";
345 cout << "pivots = " << pivs <<endl;
346 }
347 }
348 freegens = new int[rk];
349 if (!freegens)
350 {
351 out_of_memory_error("freegens");
352 return;
353 }
354 if (rk>0)
355 {
356 for (i=0; i<rk; i++) freegens[i] = gens[pivs[i+1]];
357 if (verbose>1)
358 { cout << "freegens: ";
359 for (i=0; i<rk; i++) cout << freegens[i] << " ";
360 cout << "\n";
361 }
362 }
363 delete[] check; delete[] gens; delete[] rel;
364 pivs.init(); npivs.init();
365 }
366
367 // Compute the number of cusps
368 long maxncusps =0, dd, pp, nc;
369 vector<long>::const_iterator d,p;
370 for(d=(dlist).begin();d!=(dlist).end();d++)
371 {
372 dd = *d;
373 nc = ::gcd(dd,modulus/dd);
374 for(p=plist.begin();p!=plist.end();p++) // computing phi(dd)
375 {
376 pp = *p;
377 if ((nc%pp)==0)
378 nc = nc*(pp-1)/pp;
379 }
380 maxncusps += nc;
381 }
382 if (verbose) cout << "Number of cusps is "<<maxncusps<<endl;
383
384 cusplist cusps(maxncusps, this);
385
386 for (i=0; i<rk; i++)
387 {
388 modsym m(symbol(freegens[i]));
389 for (j=1; j>-3; j-=2)
390 {
391 rational c = (j==1 ? m.beta() : m.alpha());
392 if (plusflag==-1)
393 k = cusps.index_1(c); //adds automatically if new, ignores if [c]=[-c]
394 else
395 k = cusps.index(c); //adds automatically if new
396 }
397 }
398 ncusps=cusps.count();
399
400 if(verbose)
401 {
402 cout << "ncusps = " << ncusps << endl;
403 if(verbose>1) {cusps.display(); cout<<endl;}
404 }
405
406 if (verbose) cout << "About to compute matrix of delta"<<endl;
407 mat deltamat=mat(ncusps,rk); // should make this sparse
408
409 for (i=0; i<rk; i++)
410 {
411 modsym m(symbol(freegens[i]));
412 for (j=1; j>-3; j-=2)
413 {
414 rational c = (j==1 ? m.beta() : m.alpha());
415 if (plusflag==-1)
416 k = cusps.index_1(c);
417 else
418 k = cusps.index(c);
419 if (k>0)
420 deltamat(k,i+1) += j;
421 if (k<0)
422 deltamat(-k,i+1) -= j;
423 }
424 }
425 if (verbose)
426 {
427 cout << "delta matrix done: size "<<deltamat.nrows()<<"x"<<deltamat.ncols()<<". "<<endl;
428 if(verbose>1)
429 cout<<"deltamat = "<<deltamat<<endl;
430 cout << "About to compute kernel of delta"<<endl;
431 }
432
433 kern=kernel(smat(deltamat));
434 vec pivs, npivs;
435 int d2;
436 smat sk;
437 int ok = liftmat(smat_elim(deltamat).kernel(npivs,pivs),MODULUS,sk,d2);
438 if (!ok)
439 cout << "**!!!** failed to lift modular kernel of delta matrix\n" << endl;
440 denom2=d2;
441 tkernbas = transpose(kern.bas()); // dim(kern) x rank
442 deltamat.init(0); // clear space.
443 if(verbose>1)
444 {
445 cout<<"tkernbas = "<<tkernbas.as_mat()<<endl;
446 }
447
448 if (verbose) cout << "done "<<endl;
449
450 if(cuspidal)
451 dimension = dim(kern);
452 else
453 dimension = rk;
454 // denom2 = 1;
455 denom3 = denom1*denom2;
456
457 freemods = new modsym[rk];
458 if (!freemods)
459 {
460 out_of_memory_error("freemods");
461 return;
462 }
463 needed = new int[rk];
464 if (!needed)
465 {
466 out_of_memory_error("needed");
467 return;
468 }
469
470 if (dimension>0)
471 {
472 if (verbose>1) cout << "Freemods:\n";
473 for (i=0; i<rk; i++)
474 {
475 freemods[i] = modsym(symbol(freegens[i])) ;
476 needed[i] = (cuspidal? ! trivial(kern.bas().row(i+1).as_vec())
477 : 1);
478 if (verbose>1)
479 {
480 cout << i << ": " << freemods[i];
481 if (!needed[i]) cout << " (not needed)";
482 cout << "\n";
483 }
484 }
485 if ((verbose>1)&&cuspidal)
486 {
487 cout << "Basis of ker(delta):\n";
488 cout << kern.bas().as_mat()<<endl;
489 cout << "pivots: " << pivots(kern) << endl;
490 }
491 }
492 if (verbose) cout << "Finished constructing homspace." << endl;
493 }
494
~homspace()495 homspace::~homspace()
496 {
497 if (coordindex) delete[] coordindex;
498 if (needed) delete[] needed;
499 if (freegens) delete[] freegens;
500 if (freemods) delete[] freemods;
501 }
502
503 // Extend a dual vector of length rk to one of length nsymb:
extend_coords(const vec & v)504 vec homspace::extend_coords(const vec& v)
505 {
506 // cout<<"Extending vector "<<v<<endl;
507 vec ans(nsymb);
508 int i,j;
509 for(i=1; i<=nsymb; i++)
510 {
511 j = coordindex[i-1];
512 if (j==0) ans[i] = 0;
513 else if (j>0) ans[i] = v*coord_vecs[j];
514 else if (j<0) ans[i] = -v*coord_vecs[-j];
515 }
516 // cout<<"returning "<<ans<<endl;
517 return ans;
518 }
519
520 // Contract a dual vector of length nsymb to one of length rk:
contract_coords(const vec & v)521 vec homspace::contract_coords(const vec& v)
522 {
523 // cout<<"Contracting vector "<<v<<endl;
524 vec ans(rk);
525 int i;
526 for(i=1; i<=rk; i++)
527 ans[i] = v[1+freegens[i-1]];
528 // cout<<"returning "<<ans<<endl;
529 return ans;
530 }
531
coords_from_index(int ind) const532 svec homspace::coords_from_index(int ind) const
533 {
534 long i= coordindex[ind];
535 if (i>0) return coord_vecs[i];
536 if (i<0) return -coord_vecs[-i];
537 return zero_coords();
538 }
539
proj_coords_from_index(int ind,const mat & m) const540 vec homspace::proj_coords_from_index(int ind, const mat& m) const
541 {
542 long i= coordindex[ind];
543 if (i>0) return m.row(i);
544 if (i<0) return -m.row(-i);
545 return vec(m.ncols());
546 }
547
nfproj_coords_from_index(int ind,const vec & bas) const548 long homspace::nfproj_coords_from_index(int ind, const vec& bas) const
549 {
550 long i= coordindex[ind];
551 if (i>0) return bas[i];
552 if (i<0) return -bas[-i];
553 return 0;
554 }
555
coords(const symb & s) const556 svec homspace::coords(const symb& s) const
557 {
558 return coords_from_index(index(s));
559 }
560
add_coords(svec & v,const symb & s) const561 void homspace::add_coords(svec& v, const symb& s) const
562 {
563 v += coords_from_index(index(s));
564 }
565
coords_cd(long c,long d) const566 svec homspace::coords_cd(long c, long d) const
567 {
568 return coords_from_index(index2(c,d));
569 }
570
add_coords_cd(svec & v,long c,long d) const571 void homspace::add_coords_cd(svec& v, long c, long d) const
572 {
573 v += coords_from_index(index2(c,d));
574 }
575
proj_coords_cd(long c,long d,const mat & m) const576 vec homspace::proj_coords_cd(long c, long d, const mat& m) const
577 {
578 return proj_coords_from_index(index2(c,d), m);
579 }
580
add_proj_coords_cd(vec & v,long c,long d,const mat & m) const581 void homspace::add_proj_coords_cd(vec& v, long c, long d, const mat& m) const
582 {
583 long n = coordindex[index2(c,d)];
584 if (n>0) v.add_row(m,n);
585 else if (n<0) v.sub_row(m,-n);
586 }
587
nfproj_coords_cd(long c,long d,const vec & bas) const588 long homspace::nfproj_coords_cd(long c, long d, const vec& bas) const
589 {
590 return nfproj_coords_from_index(index2(c,d), bas);
591 }
592
add_nfproj_coords_cd(long & a,long c,long d,const vec & bas) const593 void homspace::add_nfproj_coords_cd(long& a, long c, long d, const vec& bas) const
594 {
595 a += nfproj_coords_from_index(index2(c,d), bas);
596 }
597
coords(long nn,long dd) const598 svec homspace::coords(long nn, long dd) const
599 {
600 svec ans = zero_coords();
601 add_coords(ans, nn, dd);
602 return ans;
603 }
604
coords(const modsym & m) const605 svec homspace::coords(const modsym& m) const
606 {
607 svec ans = zero_coords();
608 add_coords(ans, m);
609 return ans;
610 }
611
add_coords(svec & vv,const modsym & m) const612 void homspace::add_coords(svec& vv, const modsym& m) const
613 {
614 rational al = m.alpha(), be=m.beta();
615 long a=num(be), b=num(al), c=den(be), d=den(al), u, v;
616 long de = a*d-b*c;
617 if (de<0) {de=-de; b=-b; d=-d;}
618 if (de==1)
619 {
620 add_coords_cd(vv, c, d);
621 return;
622 }
623 // now de>1
624 bezout(a,c, u,v); // =1
625 long nu = b*u+v*d;
626 //
627 // now m = M{0,infinity} = U.{nu/de,infinity} where U=[a,-v;c,u], so
628 // we find the CF expansion of nu/de.
629 //
630 long C=c, D=u, r=nu, s=de, q, t, e;
631 while (s)
632 {
633 t=s; s=mod(r,s); q=(r-s)/t; r=-t;
634 e=D; D=-C; C= q*C+e;
635 add_coords_cd(vv, -D, C);
636 }
637 }
638
639
proj_coords(long nn,long dd,const mat & m) const640 vec homspace::proj_coords(long nn, long dd, const mat& m) const
641 {
642 vec ans = vec(m.ncols());
643 add_proj_coords(ans, nn, dd, m);
644 return ans;
645 }
646
nfproj_coords(long nn,long dd,const vec & bas) const647 long homspace::nfproj_coords(long nn, long dd, const vec& bas) const
648 {
649 long ans = 0;
650 add_nfproj_coords(ans, nn, dd, bas);
651 return ans;
652 }
653
add_coords(svec & v,long nn,long dd) const654 void homspace::add_coords(svec& v, long nn, long dd) const
655 {
656 add_coords_cd(v,0,1);
657 long c=0, d=1, e, a=nn, b=dd, q, f;
658 while (b)
659 {
660 f=b; b=mod(a,b); q=(a-b)/f; a= -f;
661 e=d; d=-c; c=q*c+e;
662 add_coords_cd(v,c,d);
663 }
664 }
665
add_proj_coords(vec & v,long nn,long dd,const mat & m) const666 void homspace::add_proj_coords(vec& v, long nn, long dd, const mat& m) const
667 {
668 add_proj_coords_cd(v,0,1,m);
669 long c=0, d=1, e, a=nn, b=dd, q, f;
670 while (b)
671 {
672 f=b; b=mod(a,b); q=(a-b)/f; a= -f;
673 e=d; d=-c; c=q*c+e;
674 add_proj_coords_cd(v,c,d,m);
675 }
676 }
677
add_nfproj_coords(long & aa,long nn,long dd,const vec & bas) const678 void homspace::add_nfproj_coords(long& aa, long nn, long dd, const vec& bas) const
679 {
680 add_nfproj_coords_cd(aa,0,1,bas);
681 long c=0, d=1, e, a=nn, b=dd, q, f;
682 while (b)
683 {
684 f=b; b=mod(a,b); q=(a-b)/f; a= -f;
685 e=d; d=-c; c=q*c+e;
686 add_nfproj_coords_cd(aa,c,d,bas);
687 }
688 }
689
applyop(const matop & mlist,const rational & q) const690 svec homspace::applyop(const matop& mlist, const rational& q) const
691 { svec ans(rk);
692 long i=mlist.size();
693 while (i--) add_coords(ans,mlist[i](q));
694 return ans;
695 }
696
applyop(const matop & mlist,const modsym & m) const697 svec homspace::applyop(const matop& mlist, const modsym& m) const
698 { svec ans(rk);
699 long i=mlist.size();
700 while (i--) ans += coords((mlist[i])(m));
701 return ans;
702 }
703
704 // copy of routine from ../procs/xsplit.cc:
705 mat sparse_restrict(const mat& m, const subspace& s);
706 smat restrict_mat(const smat& m, const subspace& s);
707
calcop_restricted(string opname,long p,const matop & mlist,const subspace & s,int dual,int display) const708 mat homspace::calcop_restricted(string opname, long p, const matop& mlist,
709 const subspace& s, int dual, int display) const
710 {
711 long d=dim(s);
712 mat m(d,rk);
713 for (long j=0; j<d; j++)
714 {
715 long jj = pivots(s)[j+1]-1;
716 svec colj = applyop(mlist,freemods[jj]);
717 m.setrow(j+1,colj.as_vec());
718 }
719 m = (smat(m)*smat(basis(s))).as_mat();
720 if(!dual) m=transpose(m); // dual is default for restricted ops
721 if (display)
722 {
723 cout << "Matrix of " << opname << "(" << p << ") = ";
724 if (dimension>1) cout << "\n";
725 m.output_pretty();
726 }
727 return m;
728 }
729
s_calcop_restricted(string opname,long p,const matop & mlist,const ssubspace & s,int dual,int display) const730 smat homspace::s_calcop_restricted(string opname, long p, const matop& mlist,
731 const ssubspace& s, int dual, int display) const
732 {
733 long d=dim(s);
734 smat m(d,rk);
735 for (long j=1; j<=d; j++)
736 {
737 long jj = pivots(s)[j];
738 svec colj = applyop(mlist,freemods[jj-1]);
739 m.setrow(j,colj);
740 }
741 // m = m*basis(s);
742 m = mult_mod_p(m,basis(s),MODULUS);
743 if(!dual) m=transpose(m); // dual is default for restricted ops
744 if (display)
745 {
746 cout << "Matrix of " << opname << "(" << p << ") = ";
747 if (dimension>1) cout << "\n";
748 cout<<m<<endl;
749 }
750 return m;
751 }
752
calcop(string opname,long p,const matop & mlist,int dual,int display) const753 mat homspace::calcop(string opname, long p, const matop& mlist,
754 int dual, int display) const
755 {
756 mat m(rk,rk);
757 for (long j=0; j<rk; j++) if (needed[j])
758 { svec colj = applyop(mlist,freemods[j]);
759 m.setcol(j+1,colj.as_vec());
760 }
761 if(cuspidal) m = restrict_mat(smat(m),kern).as_mat();
762 if(dual) { m=transpose(m);}
763 if (display)
764 {
765 cout << "Matrix of " << opname << "(" << p << ") = ";
766 if (dimension>1) cout << "\n";
767 m.output_pretty();
768 }
769 return m;
770 }
771
calcop_col(string opname,long p,int j,const matop & mlist,int display) const772 vec homspace::calcop_col(string opname, long p, int j, const matop& mlist,
773 int display) const
774 {
775 // j counts from 1
776 vec colj = applyop(mlist,freemods[j-1]).as_vec();
777 if (display)
778 {
779 cout << "Image of "<<j<<"-th generator under " << opname << "(" << p << ") = "
780 << colj << endl;
781 }
782 return colj;
783 }
784
calcop_cols(string opname,long p,const vec & jlist,const matop & mlist,int display) const785 mat homspace::calcop_cols(string opname, long p, const vec& jlist, const matop& mlist,
786 int display) const
787 {
788 int i, d = dim(jlist);
789 mat m(d,rk);
790 for (i=1; i<=d; i++)
791 {
792 m.setrow(i, applyop(mlist,freemods[jlist[i]-1]).as_vec());
793 }
794 return m;
795 }
796
s_calcop_col(string opname,long p,int j,const matop & mlist,int display) const797 svec homspace::s_calcop_col(string opname, long p, int j, const matop& mlist,
798 int display) const
799 {
800 // j counts from 1
801 svec colj = applyop(mlist,freemods[j-1]);
802 if (display)
803 {
804 cout << "Image of "<<j<<"-th generator under " << opname << "(" << p << ") = "
805 << colj.as_vec() << endl;
806 }
807 return colj;
808 }
809
s_calcop_cols(string opname,long p,const vec & jlist,const matop & mlist,int display) const810 smat homspace::s_calcop_cols(string opname, long p, const vec& jlist, const matop& mlist,
811 int display) const
812 {
813 int i, d = dim(jlist);
814 smat m(d,rk);
815 for (i=1; i<=d; i++)
816 {
817 m.setrow(i, applyop(mlist,freemods[jlist[i]-1]));
818 }
819 return m;
820 }
821
s_calcop(string opname,long p,const matop & mlist,int dual,int display) const822 smat homspace::s_calcop(string opname, long p, const matop& mlist,
823 int dual, int display) const
824 {
825 smat m(rk,rk);
826 for (long j=0; j<rk; j++) if (needed[j])
827 { svec colj = applyop(mlist,freemods[j]);
828 m.setrow(j+1,colj);
829 }
830 if(cuspidal)
831 {
832 m = restrict_mat(transpose(m),kern);
833 if(dual) m = transpose(m);
834 }
835 else if(!dual) {m=transpose(m);}
836 if (display)
837 {
838 cout << "Matrix of " << opname << "(" << p << ") = ";
839 if (dimension>1) cout << "\n";
840 cout<<m.as_mat();
841 }
842 return m;
843 }
844
845 /* NOTE: heckeop actually only returns the Hecke operator for p not dividing
846 the level. For p that divide the level it actually computes the Atkin-Lehner
847 operator.
848 */
849
heckeop(long p,int dual,int display) const850 mat homspace::heckeop(long p, int dual, int display) const
851 {
852 matop matlist(p,modulus);
853 string name = ((modulus%p) ? T_opname : W_opname);
854 return calcop(name,p,matlist,dual,display);
855 }
856
s_heckeop(long p,int dual,int display) const857 smat homspace::s_heckeop(long p, int dual, int display) const
858 {
859 matop matlist(p,modulus);
860 string name = ((modulus%p) ? T_opname : W_opname);
861 return s_calcop(name,p,matlist,dual,display);
862 }
863
heckeop_col(long p,int j,int display) const864 vec homspace::heckeop_col(long p, int j, int display) const
865 {
866 matop matlist(p,modulus);
867 string name = ((modulus%p) ? T_opname : W_opname);
868 return calcop_col(name,p,j,matlist,display);
869 }
870
heckeop_cols(long p,const vec & jlist,int display) const871 mat homspace::heckeop_cols(long p, const vec& jlist, int display) const
872 {
873 matop matlist(p,modulus);
874 string name = ((modulus%p) ? T_opname : W_opname);
875 return calcop_cols(name,p,jlist,matlist,display);
876 }
877
s_heckeop_col(long p,int j,int display) const878 svec homspace::s_heckeop_col(long p, int j, int display) const
879 {
880 matop matlist(p,modulus);
881 string name = ((modulus%p) ? T_opname : W_opname);
882 return s_calcop_col(name,p,j,matlist,display);
883 }
884
s_heckeop_cols(long p,const vec & jlist,int display) const885 smat homspace::s_heckeop_cols(long p, const vec& jlist, int display) const
886 {
887 matop matlist(p,modulus);
888 string name = ((modulus%p) ? T_opname : W_opname);
889 return s_calcop_cols(name,p,jlist,matlist,display);
890 }
891
heckeop_restricted(long p,const subspace & s,int dual,int display) const892 mat homspace::heckeop_restricted(long p,const subspace& s,
893 int dual, int display) const
894 {
895 matop matlist(p,modulus);
896 string name = ((modulus%p) ? T_opname : W_opname);
897 return calcop_restricted(name,p,matlist,s,dual,display);
898 }
899
s_heckeop_restricted(long p,const ssubspace & s,int dual,int display) const900 smat homspace::s_heckeop_restricted(long p,const ssubspace& s,
901 int dual, int display) const
902 {
903 matop matlist(p,modulus);
904 string name = ((modulus%p) ? T_opname : W_opname);
905 return s_calcop_restricted(name,p,matlist,s,dual,display);
906 }
907
newheckeop(long p,int dual,int display) const908 mat homspace::newheckeop(long p, int dual, int display) const
909 {
910 if(::divides(p,modulus)) return wop(p,display);
911 matop hmats(p); // constructs H-matrices
912 long j, nmats=hmats.size();
913 svec colj(rk); mat m(rk,rk);
914 for (j=0; j<rk; j++) if (needed[j])
915 { symb s = symbol(freegens[j]);
916 colj = hmats[0](s,this);
917 colj+= hmats[1](s,this);
918 for(long i=2; i<nmats; i++) colj+= (hmats[i](s,this));
919 m.setcol(j+1,colj.as_vec());
920 }
921 if(cuspidal) m = restrict_mat(smat(m),kern).as_mat();
922 if(dual) m=transpose(m);
923 if (display)
924 {
925 cout << "Matrix of T(" << p << ") = ";
926 if (dimension>1) cout << "\n";
927 m.output_pretty();
928 }
929 return m;
930 }
931
wop(long q,int dual,int display) const932 mat homspace::wop(long q, int dual, int display) const
933 {
934 matop matlist(q,modulus);
935 return calcop(W_opname,q,matlist,dual,display);
936 }
937
s_wop(long q,int dual,int display) const938 smat homspace::s_wop(long q, int dual, int display) const
939 {
940 matop matlist(q,modulus);
941 return s_calcop(W_opname,q,matlist,dual,display);
942 }
943
conj(int dual,int display) const944 mat homspace::conj(int dual, int display) const
945 {
946 mat m(rk,rk);
947 for (long j=1; j<=rk; j++) if (needed[j-1])
948 { symb s = symbol(freegens[j-1]);
949 svec colj = coords_cd(-s.cee(),s.dee());
950 m.setcol(j,colj.as_vec());
951 }
952 if(cuspidal) m = restrict_mat(smat(m),kern).as_mat();
953 if(dual) m=transpose(m);
954 if (display) cout << "Matrix of conjugation = " << m;
955 return m;
956 }
957
conj_col(int j,int display) const958 vec homspace::conj_col(int j, int display) const
959 {
960 // j counts from 1
961 symb s = symbol(freegens[j-1]);
962 vec colj = coords_cd(-s.cee(),s.dee()).as_vec();
963 if (display) cout << "Column "<<j<<" of matrix of conjugation = " << colj << endl;
964 return colj;
965 }
966
conj_cols(const vec & jlist,int display) const967 mat homspace::conj_cols(const vec& jlist, int display) const
968 {
969 int d = dim(jlist);
970 mat m(d,rk);
971 int i;
972 for (i=1; i<=d; i++)
973 {
974 symb s = symbol(freegens[jlist[i]-1]);
975 m.setrow(i, coords_cd(-s.cee(),s.dee()).as_vec());
976 }
977 return m;
978 }
979
s_conj_col(int j,int display) const980 svec homspace::s_conj_col(int j, int display) const
981 {
982 // j counts from 1
983 symb s = symbol(freegens[j-1]);
984 svec colj = coords_cd(-s.cee(),s.dee());
985 if (display) cout << "Column "<<j<<" of matrix of conjugation = " << colj.as_vec() << endl;
986 return colj;
987 }
988
s_conj_cols(const vec & jlist,int display) const989 smat homspace::s_conj_cols(const vec& jlist, int display) const
990 {
991 int d = dim(jlist);
992 smat m(d,rk);
993 int i;
994 for (i=1; i<=d; i++)
995 {
996 symb s = symbol(freegens[jlist[i]-1]);
997 m.setrow(i, coords_cd(-s.cee(),s.dee()));
998 }
999 return m;
1000 }
1001
s_conj(int dual,int display) const1002 smat homspace::s_conj(int dual, int display) const
1003 {
1004 smat m(rk,rk);
1005 for (long j=1; j<=rk; j++) if (needed[j-1])
1006 { symb s = symbol(freegens[j-1]);
1007 svec colj = coords_cd(-s.cee(),s.dee());
1008 m.setrow(j,colj);
1009 }
1010 if(cuspidal)
1011 {
1012 m = restrict_mat(transpose(m),kern);
1013 if(dual) m = transpose(m);
1014 }
1015 else if(!dual) {m=transpose(m);}
1016 if (display) cout << "Matrix of conjugation = " << m;
1017 return m;
1018 }
1019
1020 // Computes matrix of conjugation restricted to a DUAL subspace; here
1021 // the ambient space of s must be H_1(-;cusps) of dimension rk
conj_restricted(const subspace & s,int dual,int display) const1022 mat homspace::conj_restricted(const subspace& s,
1023 int dual, int display) const
1024 {
1025 long d = dim(s);
1026 mat m(d,rk);
1027 for (long j=1; j<=d; j++)
1028 {
1029 long jj=pivots(s)[j];
1030 symb sy = symbol(freegens[jj-1]);
1031 svec colj = coords_cd(-sy.cee(),sy.dee());
1032 m.setrow(j,colj.as_vec());
1033 }
1034 m = matmulmodp(m,basis(s),MODULUS);
1035 if(!dual) m=transpose(m); // dual is default for restricted ops
1036 if (display) cout << "Matrix of conjugation = " << m;
1037 return m;
1038 }
1039
s_conj_restricted(const ssubspace & s,int dual,int display) const1040 smat homspace::s_conj_restricted(const ssubspace& s,
1041 int dual, int display) const
1042 {
1043 long d = dim(s);
1044 smat m(d,rk);
1045 for (long j=1; j<=d; j++)
1046 {
1047 long jj=pivots(s)[j];
1048 symb sy = symbol(freegens[jj-1]);
1049 svec colj = coords_cd(-sy.cee(),sy.dee());
1050 m.setrow(j,colj);
1051 }
1052 // cout<<"m = "<<m<<" = "<<m.as_mat()<<endl;
1053 m = mult_mod_p(m,basis(s),MODULUS);
1054 if(!dual) m=transpose(m); // dual is default for restricted ops
1055 if (display) cout << "Matrix of conjugation = " << m.as_mat();
1056 return m;
1057 }
1058
fricke(int dual,int display) const1059 mat homspace::fricke(int dual, int display) const
1060 {
1061 matop frickelist(modulus,modulus);
1062 return calcop(W_opname,modulus,frickelist,dual,display);
1063 }
1064
op_prime(int i)1065 long homspace::op_prime(int i) // the i'th operator prime for Tp or Wq
1066 {
1067 #ifdef NEW_OP_ORDER
1068 long p = prime_number(i+1);
1069 // cout<<"opmat("<<i<<") using p="<<p<<endl;
1070 #else
1071 long p = primelist[i];
1072 #endif
1073 return p;
1074 }
1075
opmat(int i,int dual,int v)1076 mat homspace::opmat(int i, int dual, int v)
1077 {
1078 if(i==-1) return conj(dual,v);
1079 if((i<0)||(i>=nap))
1080 {
1081 cerr<<"Error in homspace::opmat(): called with i = " << i << endl;
1082 return mat(dimension); // shouldn't happen
1083 }
1084 long p = op_prime(i);
1085 if(v)
1086 {
1087 cout<<"Computing " << ((::divides(p,modulus)) ? W_opname : T_opname) <<"("<<p<<")..."<<flush;
1088 mat ans = heckeop(p,dual,0); // Automatically chooses W or T
1089 cout<<"done."<<endl;
1090 return ans;
1091 }
1092 else return heckeop(p,dual,0); // Automatically chooses W or T
1093 }
1094
opmat_col(int i,int j,int v)1095 vec homspace::opmat_col(int i, int j, int v)
1096 {
1097 if(i==-1) return conj_col(j,v);
1098 if((i<0)||(i>=nap))
1099 {
1100 cerr<<"Error in homspace::opmat_col(): called with i = " << i << endl;
1101 return vec(dimension); // shouldn't happen
1102 }
1103 long p = op_prime(i);
1104 if(v)
1105 {
1106 cout<<"Computing col "<<j<<" of " << ((::divides(p,modulus)) ? W_opname : T_opname)
1107 <<"("<<p<<")..."<<flush;
1108 }
1109 vec ans = heckeop_col(p,j,0); // Automatically chooses W or T
1110 if(v)
1111 {
1112 cout<<"done."<<endl;
1113 }
1114 return ans;
1115 }
1116
opmat_cols(int i,const vec & jlist,int v)1117 mat homspace::opmat_cols(int i, const vec& jlist, int v)
1118 {
1119 if(i==-1) return conj_cols(jlist,v);
1120 int d = dim(jlist);
1121 if((i<0)||(i>=nap))
1122 {
1123 cerr<<"Error in homspace::opmat_cols(): called with i = " << i << endl;
1124 return mat(d,rk); // shouldn't happen
1125 }
1126 long p = op_prime(i);
1127 if(v)
1128 {
1129 cout<<"Computing "<<d<<" cols of " << ((::divides(p,modulus)) ? W_opname : T_opname)
1130 <<"("<<p<<")..."<<flush;
1131 }
1132 mat ans = heckeop_cols(p,jlist,0); // Automatically chooses W or T
1133 if(v)
1134 {
1135 cout<<"done."<<endl;
1136 }
1137 return ans;
1138 }
1139
s_opmat_col(int i,int j,int v)1140 svec homspace::s_opmat_col(int i, int j, int v)
1141 {
1142 if(i==-1) return s_conj_col(j,v);
1143 if((i<0)||(i>=nap))
1144 {
1145 cerr<<"Error in homspace::opmat_col(): called with i = " << i << endl;
1146 return svec(dimension); // shouldn't happen
1147 }
1148 long p = op_prime(i);
1149 if(v)
1150 {
1151 cout<<"Computing col "<<j<<" of " << ((::divides(p,modulus)) ? W_opname : T_opname)
1152 <<"("<<p<<")..."<<flush;
1153 }
1154 svec ans = s_heckeop_col(p,j,0); // Automatically chooses W or T
1155 if(v)
1156 {
1157 cout<<"done."<<endl;
1158 }
1159 return ans;
1160 }
1161
s_opmat_cols(int i,const vec & jlist,int v)1162 smat homspace::s_opmat_cols(int i, const vec& jlist, int v)
1163 {
1164 if(i==-1) return s_conj_cols(jlist,v);
1165 int d = dim(jlist);
1166 if((i<0)||(i>=nap))
1167 {
1168 cerr<<"Error in homspace::opmat_col(): called with i = " << i << endl;
1169 return smat(d,rk); // shouldn't happen
1170 }
1171 long p = op_prime(i);
1172 if(v)
1173 {
1174 cout<<"Computing "<<d<<" cols of " << ((::divides(p,modulus)) ? W_opname : T_opname)
1175 <<"("<<p<<")..."<<flush;
1176 }
1177 smat ans = s_heckeop_cols(p,jlist,0); // Automatically chooses W or T
1178 if(v)
1179 {
1180 cout<<"done."<<endl;
1181 }
1182 return ans;
1183 }
1184
s_opmat(int i,int dual,int v)1185 smat homspace::s_opmat(int i, int dual, int v)
1186 {
1187 if(i==-1) return s_conj(dual,v);
1188 if((i<0)||(i>=nap))
1189 {
1190 cerr<<"Error in homspace::s_opmat(): called with i = " << i << endl;
1191 return smat(dimension); // shouldn't happen
1192 }
1193 long p = op_prime(i);
1194 if(v)
1195 {
1196 cout<<"Computing " << ((::divides(p,modulus)) ? W_opname : T_opname) <<"("<<p<<")..."<<flush;
1197 smat ans = s_heckeop(p,dual,0); // Automatically chooses W or T
1198 cout<<"done."<<endl;
1199 return ans;
1200 }
1201 else return s_heckeop(p,dual,0); // Automatically chooses W or T
1202 }
1203
opmat_restricted(int i,const subspace & s,int dual,int v)1204 mat homspace::opmat_restricted(int i, const subspace& s, int dual, int v)
1205 {
1206 if(i==-1) return conj_restricted(s,dual,v);
1207 if((i<0)||(i>=nap))
1208 {
1209 cerr<<"Error in homspace::opmat_restricted(): called with i = "
1210 << i << endl;
1211 return mat(dim(s)); // shouldn't happen
1212 }
1213 long p = op_prime(i);
1214 if(v)
1215 {
1216 cout<<"Computing " << ((::divides(p,modulus)) ? W_opname : T_opname) <<"("<<p
1217 <<") restricted to subspace of dimension "<<dim(s)<<" ..."<<flush;
1218 mat ans = heckeop_restricted(p,s,dual,0); // Automatically chooses W or T
1219 cout<<"done."<<endl;
1220 return ans;
1221 }
1222 else return heckeop_restricted(p,s,dual,0); // Automatically chooses W or T
1223 }
1224
s_opmat_restricted(int i,const ssubspace & s,int dual,int v)1225 smat homspace::s_opmat_restricted(int i, const ssubspace& s, int dual, int v)
1226 {
1227 if(i==-1) return s_conj_restricted(s,dual,v);
1228 if((i<0)||(i>=nap))
1229 {
1230 cerr<<"Error in homspace::s_opmat_restricted(): called with i = "
1231 << i << endl;
1232 return smat(dim(s)); // shouldn't happen
1233 }
1234 long p = op_prime(i);
1235 if(v)
1236 {
1237 cout<<"Computing " << ((::divides(p,modulus)) ? W_opname : T_opname) <<"("<<p
1238 <<") restricted to subspace of dimension "<<dim(s)<<" ..."<<flush;
1239 smat ans = s_heckeop_restricted(p,s,dual,(v>2)); // Automatically chooses W or T
1240 cout<<"done."<<endl;
1241 return ans;
1242 }
1243 else return s_heckeop_restricted(p,s,dual,0); // Automatically chooses W or T
1244 }
1245
1246 #ifdef OLD_EIG_ORDER
1247
1248 static long*pm1={1,-1};
1249
T_eigrange(long p)1250 vector<long> T_eigrange(long p)
1251 {
1252 vector<long> ans;
1253 ans.push_back(0);
1254 long aplim=1;
1255 while (aplim*aplim<=4*p) aplim++; aplim--;
1256 for(long ap=1; ap<=aplim; ap++)
1257 {
1258 ans.push_back(ap);
1259 ans.push_back(-ap);
1260 }
1261 return ans;
1262 }
1263
1264 #else // new eig order, in strict numerical order
1265
1266 static long pm1[]={-1,1};
1267
T_eigrange(long p)1268 vector<long> T_eigrange(long p)
1269 {
1270 long aplim=3, four_p=p<<2;
1271 while (aplim*aplim<=four_p) aplim++;
1272 aplim--;
1273 vector<long> ans(1+2*aplim);
1274 iota(ans.begin(),ans.end(),-aplim);
1275 return ans;
1276 }
1277 #endif
1278
eigrange(long i)1279 vector<long> homspace::eigrange(long i)
1280 {
1281 if((i<0)||(i>=nap)) return vector<long>(0); // shouldn't happen
1282 long p = op_prime(i);
1283 if(::divides(p,modulus)) return vector<long>(pm1,pm1+2);
1284 return T_eigrange(p);
1285 }
1286
maninvector(long p) const1287 vec homspace::maninvector(long p) const
1288 {
1289 long i,p2;
1290 svec tvec = coords(0,p); // =0, but sets the right length.
1291 if (plusflag!=-1)
1292 {
1293 if (p==2)
1294 add_coords(tvec,1,2);
1295 else
1296 {
1297 p2=(p-1)>>1;
1298 for (i=1; i<=p2; i++) { add_coords(tvec,i,p); }
1299 if(plusflag)
1300 tvec *=2;
1301 else
1302 for (i=1; i<=p2; i++) { add_coords(tvec,-i,p); }
1303 }
1304 }
1305 if(cuspidal)
1306 return cuspidalpart(tvec.as_vec());
1307 else
1308 return tvec.as_vec();
1309 }
1310
manintwist(long p) const1311 vec homspace::manintwist(long p) const
1312 {
1313 svec sum = coords(0,p); // =0, but sets the right length.
1314 for (long i=1; i<p; i++) sum += legendre(i,p)*coords(i,p);
1315 if(cuspidal) return cuspidalpart(sum.as_vec());
1316 else return sum.as_vec();
1317 }
1318
1319 #if(0) // no longer used
projmaninvector(long p) const1320 vec homspace::projmaninvector(long p) const // Will only work after "proj"
1321 {
1322 long i,p2;
1323 vec tvec = proj_coords(0,p); // =0, but sets the right length.
1324 if (plusflag==-1) return tvec;
1325 if (p==2)
1326 add_proj_coords(tvec,1,2);
1327 else
1328 {
1329 p2=(p-1)>>1;
1330 for (i=1; i<=p2; i++) { add_proj_coords(tvec,i,p); }
1331 if(plusflag)
1332 tvec *=2;
1333 else
1334 for (i=1; i<=p2; i++) { add_proj_coords(tvec,-i,p); }
1335 }
1336 return tvec;
1337 }
1338
projmaninvector(long p,const mat & m) const1339 vec homspace::projmaninvector(long p, const mat& m) const
1340 {
1341 long i,p2;
1342 vec tvec = proj_coords(0,p,m); // =0, but sets the right length.
1343 if (plusflag==-1) return tvec;
1344 if (p==2)
1345 add_proj_coords(tvec,1,2,m);
1346 else
1347 {
1348 p2=(p-1)>>1;
1349 for (i=1; i<=p2; i++) { add_proj_coords(tvec,i,p,m); }
1350 if(plusflag)
1351 tvec *=2;
1352 else
1353 for (i=1; i<=p2; i++) { add_proj_coords(tvec,-i,p,m); }
1354 }
1355 return tvec;
1356 }
1357
newhecke(long p,long n,long d) const1358 vec homspace::newhecke(long p, long n, long d) const
1359 // Will only work after "proj"
1360 { vec tvec = proj_coords(p*n,d);
1361 // cout<<"newhecke starts: tvec = "<<tvec<<endl;
1362 long p2 = p>>1, dp = d*p, k;
1363 for (k=1+p2-p; k<=p2; k++)
1364 {
1365 add_proj_coords(tvec,n+d*k,dp);
1366 // cout<<"tvec = "<<tvec<<endl;
1367 }
1368 // cout<<"newhecke returns: tvec = "<<tvec<<endl;
1369 return tvec;
1370 }
1371 #endif // no longer used
1372
matop(long p)1373 matop::matop(long p)
1374 {
1375
1376 // case 31: nmats = 106; break;
1377 // case 37: nmats = 128; break;
1378 // case 41: nmats = 146; break;
1379 // case 43: nmats = 154; break;
1380 // case 47: nmats = 170; break;
1381
1382 switch (p) {
1383 case 2:
1384 mats.resize(4);
1385 mats[0]=mat22(2,0,0,1); mats[1]=mat22(2,1,0,1);
1386 mats[2]=mat22(1,0,1,2); mats[3]=mat22(1,0,0,2);
1387 return;
1388 case 3:
1389 mats.resize(6);
1390 mats[0]=mat22(1,0,0,3);
1391 mats[1]=mat22(3,1,0,1);
1392 mats[2]=mat22(1,0,1,3);
1393 mats[3]=mat22(3,0,0,1);
1394 mats[4]=mat22(3,-1,0,1);
1395 mats[5]=mat22(-1,0,1,-3);
1396 return;
1397 case 5:
1398 mats.resize(12);
1399 mats[0]=mat22(1,0,0,5);
1400 mats[1]=mat22(5,2,0,1);
1401 mats[2]=mat22(2,1,1,3);
1402 mats[3]=mat22(1,0,3,5);
1403 mats[4]=mat22(5,1,0,1);
1404 mats[5]=mat22(1,0,1,5);
1405 mats[6]=mat22(5,0,0,1);
1406 mats[7]=mat22(5,-1,0,1);
1407 mats[8]=mat22(-1,0,1,-5);
1408 mats[9]=mat22(5,-2,0,1);
1409 mats[10]=mat22(-2,1,1,-3);
1410 mats[11]=mat22(1,0,-3,5);
1411 return;
1412 case 7:
1413 mats.resize(18);
1414 mats[0]=mat22(1,0,0,7);
1415 mats[1]=mat22(7,3,0,1);
1416 mats[2]=mat22(3,-1,1,2);
1417 mats[3]=mat22(-1,0,2,-7);
1418 mats[4]=mat22(7,2,0,1);
1419 mats[5]=mat22(2,1,1,4);
1420 mats[6]=mat22(1,0,4,7);
1421 mats[7]=mat22(7,1,0,1);
1422 mats[8]=mat22(1,0,1,7);
1423 mats[9]=mat22(7,0,0,1);
1424 mats[10]=mat22(7,-1,0,1);
1425 mats[11]=mat22(-1,0,1,-7);
1426 mats[12]=mat22(7,-2,0,1);
1427 mats[13]=mat22(-2,1,1,-4);
1428 mats[14]=mat22(1,0,-4,7);
1429 mats[15]=mat22(7,-3,0,1);
1430 mats[16]=mat22(-3,-1,1,-2);
1431 mats[17]=mat22(-1,0,-2,-7);
1432 return;
1433 case 11:
1434 mats.resize(30);
1435 mats[0]=mat22(1,0,0,11);
1436 mats[1]=mat22(11,5,0,1);
1437 mats[2]=mat22(5,-1,1,2);
1438 mats[3]=mat22(-1,0,2,-11);
1439 mats[4]=mat22(11,4,0,1);
1440 mats[5]=mat22(4,1,1,3);
1441 mats[6]=mat22(1,0,3,11);
1442 mats[7]=mat22(11,3,0,1);
1443 mats[8]=mat22(3,1,1,4);
1444 mats[9]=mat22(1,0,4,11);
1445 mats[10]=mat22(11,2,0,1);
1446 mats[11]=mat22(2,1,1,6);
1447 mats[12]=mat22(1,0,6,11);
1448 mats[13]=mat22(11,1,0,1);
1449 mats[14]=mat22(1,0,1,11);
1450 mats[15]=mat22(11,0,0,1);
1451 mats[16]=mat22(11,-1,0,1);
1452 mats[17]=mat22(-1,0,1,-11);
1453 mats[18]=mat22(11,-2,0,1);
1454 mats[19]=mat22(-2,1,1,-6);
1455 mats[20]=mat22(1,0,-6,11);
1456 mats[21]=mat22(11,-3,0,1);
1457 mats[22]=mat22(-3,1,1,-4);
1458 mats[23]=mat22(1,0,-4,11);
1459 mats[24]=mat22(11,-4,0,1);
1460 mats[25]=mat22(-4,1,1,-3);
1461 mats[26]=mat22(1,0,-3,11);
1462 mats[27]=mat22(11,-5,0,1);
1463 mats[28]=mat22(-5,-1,1,-2);
1464 mats[29]=mat22(-1,0,-2,-11);
1465 return;
1466 case 13:
1467 mats.resize(38);
1468 mats[0]=mat22(1,0,0,13);
1469 mats[1]=mat22(13,6,0,1);
1470 mats[2]=mat22(6,-1,1,2);
1471 mats[3]=mat22(-1,0,2,-13);
1472 mats[4]=mat22(13,5,0,1);
1473 mats[5]=mat22(5,2,1,3);
1474 mats[6]=mat22(2,-1,3,5);
1475 mats[7]=mat22(-1,0,5,-13);
1476 mats[8]=mat22(13,4,0,1);
1477 mats[9]=mat22(4,-1,1,3);
1478 mats[10]=mat22(-1,0,3,-13);
1479 mats[11]=mat22(13,3,0,1);
1480 mats[12]=mat22(3,-1,1,4);
1481 mats[13]=mat22(-1,0,4,-13);
1482 mats[14]=mat22(13,2,0,1);
1483 mats[15]=mat22(2,1,1,7);
1484 mats[16]=mat22(1,0,7,13);
1485 mats[17]=mat22(13,1,0,1);
1486 mats[18]=mat22(1,0,1,13);
1487 mats[19]=mat22(13,0,0,1);
1488 mats[20]=mat22(13,-1,0,1);
1489 mats[21]=mat22(-1,0,1,-13);
1490 mats[22]=mat22(13,-2,0,1);
1491 mats[23]=mat22(-2,1,1,-7);
1492 mats[24]=mat22(1,0,-7,13);
1493 mats[25]=mat22(13,-3,0,1);
1494 mats[26]=mat22(-3,-1,1,-4);
1495 mats[27]=mat22(-1,0,-4,-13);
1496 mats[28]=mat22(13,-4,0,1);
1497 mats[29]=mat22(-4,-1,1,-3);
1498 mats[30]=mat22(-1,0,-3,-13);
1499 mats[31]=mat22(13,-5,0,1);
1500 mats[32]=mat22(-5,2,1,-3);
1501 mats[33]=mat22(2,-1,-3,8);
1502 mats[34]=mat22(-1,0,8,-13);
1503 mats[35]=mat22(13,-6,0,1);
1504 mats[36]=mat22(-6,-1,1,-2);
1505 mats[37]=mat22(-1,0,-2,-13);
1506 return;
1507 case 17:
1508 mats.resize(52);
1509 mats[0]=mat22(1,0,0,17);
1510 mats[1]=mat22(17,8,0,1);
1511 mats[2]=mat22(8,-1,1,2);
1512 mats[3]=mat22(-1,0,2,-17);
1513 mats[4]=mat22(17,7,0,1);
1514 mats[5]=mat22(7,-3,1,2);
1515 mats[6]=mat22(-3,-1,2,-5);
1516 mats[7]=mat22(-1,0,-5,-17);
1517 mats[8]=mat22(17,6,0,1);
1518 mats[9]=mat22(6,1,1,3);
1519 mats[10]=mat22(1,0,3,17);
1520 mats[11]=mat22(17,5,0,1);
1521 mats[12]=mat22(5,-2,1,3);
1522 mats[13]=mat22(-2,-1,3,-7);
1523 mats[14]=mat22(-1,0,-7,-17);
1524 mats[15]=mat22(17,4,0,1);
1525 mats[16]=mat22(4,-1,1,4);
1526 mats[17]=mat22(-1,0,4,-17);
1527 mats[18]=mat22(17,3,0,1);
1528 mats[19]=mat22(3,1,1,6);
1529 mats[20]=mat22(1,0,6,17);
1530 mats[21]=mat22(17,2,0,1);
1531 mats[22]=mat22(2,1,1,9);
1532 mats[23]=mat22(1,0,9,17);
1533 mats[24]=mat22(17,1,0,1);
1534 mats[25]=mat22(1,0,1,17);
1535 mats[26]=mat22(17,0,0,1);
1536 mats[27]=mat22(17,-1,0,1);
1537 mats[28]=mat22(-1,0,1,-17);
1538 mats[29]=mat22(17,-2,0,1);
1539 mats[30]=mat22(-2,1,1,-9);
1540 mats[31]=mat22(1,0,-9,17);
1541 mats[32]=mat22(17,-3,0,1);
1542 mats[33]=mat22(-3,1,1,-6);
1543 mats[34]=mat22(1,0,-6,17);
1544 mats[35]=mat22(17,-4,0,1);
1545 mats[36]=mat22(-4,-1,1,-4);
1546 mats[37]=mat22(-1,0,-4,-17);
1547 mats[38]=mat22(17,-5,0,1);
1548 mats[39]=mat22(-5,-2,1,-3);
1549 mats[40]=mat22(-2,-1,-3,-10);
1550 mats[41]=mat22(-1,0,-10,-17);
1551 mats[42]=mat22(17,-6,0,1);
1552 mats[43]=mat22(-6,1,1,-3);
1553 mats[44]=mat22(1,0,-3,17);
1554 mats[45]=mat22(17,-7,0,1);
1555 mats[46]=mat22(-7,-3,1,-2);
1556 mats[47]=mat22(-3,1,-2,-5);
1557 mats[48]=mat22(1,0,-5,17);
1558 mats[49]=mat22(17,-8,0,1);
1559 mats[50]=mat22(-8,-1,1,-2);
1560 mats[51]=mat22(-1,0,-2,-17);
1561 return;
1562 case 19:
1563 mats.resize(58);
1564 mats[0]=mat22(1,0,0,19);
1565 mats[1]=mat22(19,9,0,1);
1566 mats[2]=mat22(9,-1,1,2);
1567 mats[3]=mat22(-1,0,2,-19);
1568 mats[4]=mat22(19,8,0,1);
1569 mats[5]=mat22(8,-3,1,2);
1570 mats[6]=mat22(-3,1,2,-7);
1571 mats[7]=mat22(1,0,-7,19);
1572 mats[8]=mat22(19,7,0,1);
1573 mats[9]=mat22(7,2,1,3);
1574 mats[10]=mat22(2,-1,3,8);
1575 mats[11]=mat22(-1,0,8,-19);
1576 mats[12]=mat22(19,6,0,1);
1577 mats[13]=mat22(6,-1,1,3);
1578 mats[14]=mat22(-1,0,3,-19);
1579 mats[15]=mat22(19,5,0,1);
1580 mats[16]=mat22(5,1,1,4);
1581 mats[17]=mat22(1,0,4,19);
1582 mats[18]=mat22(19,4,0,1);
1583 mats[19]=mat22(4,1,1,5);
1584 mats[20]=mat22(1,0,5,19);
1585 mats[21]=mat22(19,3,0,1);
1586 mats[22]=mat22(3,-1,1,6);
1587 mats[23]=mat22(-1,0,6,-19);
1588 mats[24]=mat22(19,2,0,1);
1589 mats[25]=mat22(2,1,1,10);
1590 mats[26]=mat22(1,0,10,19);
1591 mats[27]=mat22(19,1,0,1);
1592 mats[28]=mat22(1,0,1,19);
1593 mats[29]=mat22(19,0,0,1);
1594 mats[30]=mat22(19,-1,0,1);
1595 mats[31]=mat22(-1,0,1,-19);
1596 mats[32]=mat22(19,-2,0,1);
1597 mats[33]=mat22(-2,1,1,-10);
1598 mats[34]=mat22(1,0,-10,19);
1599 mats[35]=mat22(19,-3,0,1);
1600 mats[36]=mat22(-3,-1,1,-6);
1601 mats[37]=mat22(-1,0,-6,-19);
1602 mats[38]=mat22(19,-4,0,1);
1603 mats[39]=mat22(-4,1,1,-5);
1604 mats[40]=mat22(1,0,-5,19);
1605 mats[41]=mat22(19,-5,0,1);
1606 mats[42]=mat22(-5,1,1,-4);
1607 mats[43]=mat22(1,0,-4,19);
1608 mats[44]=mat22(19,-6,0,1);
1609 mats[45]=mat22(-6,-1,1,-3);
1610 mats[46]=mat22(-1,0,-3,-19);
1611 mats[47]=mat22(19,-7,0,1);
1612 mats[48]=mat22(-7,2,1,-3);
1613 mats[49]=mat22(2,-1,-3,11);
1614 mats[50]=mat22(-1,0,11,-19);
1615 mats[51]=mat22(19,-8,0,1);
1616 mats[52]=mat22(-8,-3,1,-2);
1617 mats[53]=mat22(-3,-1,-2,-7);
1618 mats[54]=mat22(-1,0,-7,-19);
1619 mats[55]=mat22(19,-9,0,1);
1620 mats[56]=mat22(-9,-1,1,-2);
1621 mats[57]=mat22(-1,0,-2,-19);
1622 return;
1623 case 23:
1624 mats.resize(74);
1625 mats[0]=mat22(1,0,0,23);
1626 mats[1]=mat22(23,11,0,1);
1627 mats[2]=mat22(11,-1,1,2);
1628 mats[3]=mat22(-1,0,2,-23);
1629 mats[4]=mat22(23,10,0,1);
1630 mats[5]=mat22(10,-3,1,2);
1631 mats[6]=mat22(-3,-1,2,-7);
1632 mats[7]=mat22(-1,0,-7,-23);
1633 mats[8]=mat22(23,9,0,1);
1634 mats[9]=mat22(9,4,1,3);
1635 mats[10]=mat22(4,-1,3,5);
1636 mats[11]=mat22(-1,0,5,-23);
1637 mats[12]=mat22(23,8,0,1);
1638 mats[13]=mat22(8,1,1,3);
1639 mats[14]=mat22(1,0,3,23);
1640 mats[15]=mat22(23,7,0,1);
1641 mats[16]=mat22(7,-2,1,3);
1642 mats[17]=mat22(-2,-1,3,-10);
1643 mats[18]=mat22(-1,0,-10,-23);
1644 mats[19]=mat22(23,6,0,1);
1645 mats[20]=mat22(6,1,1,4);
1646 mats[21]=mat22(1,0,4,23);
1647 mats[22]=mat22(23,5,0,1);
1648 mats[23]=mat22(5,2,1,5);
1649 mats[24]=mat22(2,-1,5,9);
1650 mats[25]=mat22(-1,0,9,-23);
1651 mats[26]=mat22(23,4,0,1);
1652 mats[27]=mat22(4,1,1,6);
1653 mats[28]=mat22(1,0,6,23);
1654 mats[29]=mat22(23,3,0,1);
1655 mats[30]=mat22(3,1,1,8);
1656 mats[31]=mat22(1,0,8,23);
1657 mats[32]=mat22(23,2,0,1);
1658 mats[33]=mat22(2,1,1,12);
1659 mats[34]=mat22(1,0,12,23);
1660 mats[35]=mat22(23,1,0,1);
1661 mats[36]=mat22(1,0,1,23);
1662 mats[37]=mat22(23,0,0,1);
1663 mats[38]=mat22(23,-1,0,1);
1664 mats[39]=mat22(-1,0,1,-23);
1665 mats[40]=mat22(23,-2,0,1);
1666 mats[41]=mat22(-2,1,1,-12);
1667 mats[42]=mat22(1,0,-12,23);
1668 mats[43]=mat22(23,-3,0,1);
1669 mats[44]=mat22(-3,1,1,-8);
1670 mats[45]=mat22(1,0,-8,23);
1671 mats[46]=mat22(23,-4,0,1);
1672 mats[47]=mat22(-4,1,1,-6);
1673 mats[48]=mat22(1,0,-6,23);
1674 mats[49]=mat22(23,-5,0,1);
1675 mats[50]=mat22(-5,2,1,-5);
1676 mats[51]=mat22(2,-1,-5,14);
1677 mats[52]=mat22(-1,0,14,-23);
1678 mats[53]=mat22(23,-6,0,1);
1679 mats[54]=mat22(-6,1,1,-4);
1680 mats[55]=mat22(1,0,-4,23);
1681 mats[56]=mat22(23,-7,0,1);
1682 mats[57]=mat22(-7,-2,1,-3);
1683 mats[58]=mat22(-2,-1,-3,-13);
1684 mats[59]=mat22(-1,0,-13,-23);
1685 mats[60]=mat22(23,-8,0,1);
1686 mats[61]=mat22(-8,1,1,-3);
1687 mats[62]=mat22(1,0,-3,23);
1688 mats[63]=mat22(23,-9,0,1);
1689 mats[64]=mat22(-9,4,1,-3);
1690 mats[65]=mat22(4,1,-3,5);
1691 mats[66]=mat22(1,0,5,23);
1692 mats[67]=mat22(23,-10,0,1);
1693 mats[68]=mat22(-10,-3,1,-2);
1694 mats[69]=mat22(-3,1,-2,-7);
1695 mats[70]=mat22(1,0,-7,23);
1696 mats[71]=mat22(23,-11,0,1);
1697 mats[72]=mat22(-11,-1,1,-2);
1698 mats[73]=mat22(-1,0,-2,-23);
1699 return;
1700 case 29:
1701 mats.resize(96);
1702 mats[0]=mat22(1,0,0,29);
1703 mats[1]=mat22(29,14,0,1);
1704 mats[2]=mat22(14,-1,1,2);
1705 mats[3]=mat22(-1,0,2,-29);
1706 mats[4]=mat22(29,13,0,1);
1707 mats[5]=mat22(13,-3,1,2);
1708 mats[6]=mat22(-3,-1,2,-9);
1709 mats[7]=mat22(-1,0,-9,-29);
1710 mats[8]=mat22(29,12,0,1);
1711 mats[9]=mat22(12,-5,1,2);
1712 mats[10]=mat22(-5,-2,2,-5);
1713 mats[11]=mat22(-2,1,-5,-12);
1714 mats[12]=mat22(1,0,-12,29);
1715 mats[13]=mat22(29,11,0,1);
1716 mats[14]=mat22(11,4,1,3);
1717 mats[15]=mat22(4,1,3,8);
1718 mats[16]=mat22(1,0,8,29);
1719 mats[17]=mat22(29,10,0,1);
1720 mats[18]=mat22(10,1,1,3);
1721 mats[19]=mat22(1,0,3,29);
1722 mats[20]=mat22(29,9,0,1);
1723 mats[21]=mat22(9,-2,1,3);
1724 mats[22]=mat22(-2,-1,3,-13);
1725 mats[23]=mat22(-1,0,-13,-29);
1726 mats[24]=mat22(29,8,0,1);
1727 mats[25]=mat22(8,3,1,4);
1728 mats[26]=mat22(3,1,4,11);
1729 mats[27]=mat22(1,0,11,29);
1730 mats[28]=mat22(29,7,0,1);
1731 mats[29]=mat22(7,-1,1,4);
1732 mats[30]=mat22(-1,0,4,-29);
1733 mats[31]=mat22(29,6,0,1);
1734 mats[32]=mat22(6,1,1,5);
1735 mats[33]=mat22(1,0,5,29);
1736 mats[34]=mat22(29,5,0,1);
1737 mats[35]=mat22(5,1,1,6);
1738 mats[36]=mat22(1,0,6,29);
1739 mats[37]=mat22(29,4,0,1);
1740 mats[38]=mat22(4,-1,1,7);
1741 mats[39]=mat22(-1,0,7,-29);
1742 mats[40]=mat22(29,3,0,1);
1743 mats[41]=mat22(3,1,1,10);
1744 mats[42]=mat22(1,0,10,29);
1745 mats[43]=mat22(29,2,0,1);
1746 mats[44]=mat22(2,1,1,15);
1747 mats[45]=mat22(1,0,15,29);
1748 mats[46]=mat22(29,1,0,1);
1749 mats[47]=mat22(1,0,1,29);
1750 mats[48]=mat22(29,0,0,1);
1751 mats[49]=mat22(29,-1,0,1);
1752 mats[50]=mat22(-1,0,1,-29);
1753 mats[51]=mat22(29,-2,0,1);
1754 mats[52]=mat22(-2,1,1,-15);
1755 mats[53]=mat22(1,0,-15,29);
1756 mats[54]=mat22(29,-3,0,1);
1757 mats[55]=mat22(-3,1,1,-10);
1758 mats[56]=mat22(1,0,-10,29);
1759 mats[57]=mat22(29,-4,0,1);
1760 mats[58]=mat22(-4,-1,1,-7);
1761 mats[59]=mat22(-1,0,-7,-29);
1762 mats[60]=mat22(29,-5,0,1);
1763 mats[61]=mat22(-5,1,1,-6);
1764 mats[62]=mat22(1,0,-6,29);
1765 mats[63]=mat22(29,-6,0,1);
1766 mats[64]=mat22(-6,1,1,-5);
1767 mats[65]=mat22(1,0,-5,29);
1768 mats[66]=mat22(29,-7,0,1);
1769 mats[67]=mat22(-7,-1,1,-4);
1770 mats[68]=mat22(-1,0,-4,-29);
1771 mats[69]=mat22(29,-8,0,1);
1772 mats[70]=mat22(-8,3,1,-4);
1773 mats[71]=mat22(3,-1,-4,11);
1774 mats[72]=mat22(-1,0,11,-29);
1775 mats[73]=mat22(29,-9,0,1);
1776 mats[74]=mat22(-9,-2,1,-3);
1777 mats[75]=mat22(-2,-1,-3,-16);
1778 mats[76]=mat22(-1,0,-16,-29);
1779 mats[77]=mat22(29,-10,0,1);
1780 mats[78]=mat22(-10,1,1,-3);
1781 mats[79]=mat22(1,0,-3,29);
1782 mats[80]=mat22(29,-11,0,1);
1783 mats[81]=mat22(-11,4,1,-3);
1784 mats[82]=mat22(4,-1,-3,8);
1785 mats[83]=mat22(-1,0,8,-29);
1786 mats[84]=mat22(29,-12,0,1);
1787 mats[85]=mat22(-12,-5,1,-2);
1788 mats[86]=mat22(-5,2,-2,-5);
1789 mats[87]=mat22(2,1,-5,12);
1790 mats[88]=mat22(1,0,12,29);
1791 mats[89]=mat22(29,-13,0,1);
1792 mats[90]=mat22(-13,-3,1,-2);
1793 mats[91]=mat22(-3,1,-2,-9);
1794 mats[92]=mat22(1,0,-9,29);
1795 mats[93]=mat22(29,-14,0,1);
1796 mats[94]=mat22(-14,-1,1,-2);
1797 mats[95]=mat22(-1,0,-2,-29);
1798 return;
1799 default:
1800 long p2=(p-1)>>1;
1801 long sl, sg, x1, x2, x3, y1, y2, y3, a, b, c, q, r;
1802 mats.push_back(mat22(1,0,0,p));
1803 mats.push_back(mat22(p,0,0,1));
1804 sl = -2;
1805 for(sg=1; sg>sl; sg-=2) // i.e. sg = +1 then -1
1806 for(r=1; r<=p2; r++)
1807 {
1808 x1=p; x2=-sg*r; y1=0; y2=1; a=-p; b=sg*r;
1809 mats.push_back(mat22(x1,x2,y1,y2));
1810 while(b!=0)
1811 {
1812 c=mod(a,b); q=(a-c)/b;
1813 x3=q*x2-x1; y3=q*y2-y1;
1814 a=-b; b=c; x1=x2; x2=x3; y1=y2; y2=y3;
1815 mats.push_back(mat22(x1,x2,y1,y2));
1816 }
1817 }
1818 }
1819 }
1820