forked from JohnCremona/eclib
-
Notifications
You must be signed in to change notification settings - Fork 0
/
hecketest.cc
389 lines (362 loc) · 10.9 KB
/
hecketest.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
// FILE HECKETEST.CC -- Test program for Hecke operators
//////////////////////////////////////////////////////////////////////////
//
// Copyright 1990-2012 John Cremona
//
// This file is part of the eclib package.
//
// eclib is free software; you can redistribute it and/or modify it
// under the terms of the GNU General Public License as published by the
// Free Software Foundation; either version 2 of the License, or (at your
// option) any later version.
//
// eclib is distributed in the hope that it will be useful, but WITHOUT
// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
// for more details.
//
// You should have received a copy of the GNU General Public License
// along with eclib; if not, write to the Free Software Foundation,
// Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA
//
//////////////////////////////////////////////////////////////////////////
//
#include <eclib/interface.h>
#include <eclib/timer.h>
#include <NTL/mat_ZZ.h>
#include <NTL/mat_poly_ZZ.h>
#include <NTL/ZZXFactoring.h>
#include <NTL/LLL.h>
#include <eclib/moddata.h>
#include <eclib/symb.h>
#include <eclib/cusp.h>
#include <eclib/homspace.h>
#include <eclib/smatrix_elim.h>
#include <eclib/mmatrix.h>
#include <eclib/msubspace.h>
//#define AUTOLOOP
//#define COMPARE_OLD
//#define CHECK_COMMUTE
//#define TEST_EIGS
double sparsity(const mat_m& m);
double sparsity(const mat& m);
vector<long> eigrange(long p)
{
long aplim=3, four_p=p<<2;
while (aplim*aplim<=four_p) aplim++;
aplim--;
vector<long> ans(1+2*aplim);
iota(ans.begin(),ans.end(),-aplim);
return ans;
}
int main(void)
{
cout << "Program hecketest. Using METHOD = " << METHOD << " to find newforms" << endl;
#ifdef MODULAR
cout << "MODULUS for linear algebra = " << MODULUS << endl;
#endif
init_time();
start_time();
int n=1;
int plus=1;
int verbose=0;
cout << "See the hecke matrices (0/1)? "; cin >> verbose;
cout << "Plus space (0/1)? "; cin >> plus;
int limit;
#ifdef AUTOLOOP
cout<<"Enter limit on level: ";cin>>limit;
while (n<limit) { n++;
#else
while (n>0) { cout<<"Enter level: "; cin>>n;
#endif
if (n>0)
{
cout << ">>>Level " << n << "\t";
homspace hplus(n,plus,0,0);
int genus = hplus.h1dim();
long den = hplus.h1denom();
bigint den2; den2 = den*den;
cout << "Dimension = " << genus << "\n";
cout << "denominator = " << den << "\n";
vector<long> badprimes = hplus.plist;
int nq = badprimes.size(); int firstq=0; // =0 for all W's
if (genus>0)
{
mat_m id = idmat(genus);
mat_m id2 = den2*id;
mat_m* wqlist = new mat_m[nq];
cout << "Computing conjmat... " << flush;
smat conjmat = hplus.s_conj(1);
cout<<" done."<<endl;
cout << "Computing +1 eigenspace... " << flush;
ssubspace h1plus = eigenspace(conjmat,den);
cout<<" done, dimension = "<<dim(h1plus)<<endl;
cout << "Computing -1 eigenspace... " << flush;
ssubspace h1minus = eigenspace(conjmat,-den);
cout<<" done, dimension = "<<dim(h1minus)<<endl;
int w_eigs=0;
cout<<"Compute W-eigenspaces? "; cin>>w_eigs;
for (int i=0; i<nq; i++)
{long q=badprimes[i]; if(i<firstq) continue;
cout << "Computing W("<<q<<")... " << flush;
mat wq = hplus.heckeop(q,verbose);
cout << "done, sparsity = "<<sparsity(wq)<<". " << endl;
if(verbose)
{
cout<<"Computed matrix ";
if(den>1) cout<<" (scaled by "<<den<<") ";
cout<<" = "<<wq<<endl;
// wq.output_pretty();
}
if(w_eigs) {
smat swq(wq);
int e; long mult;
for(e=1; e>-2; e-=2)
{
/*
cout<<"Using modular matrix code..."<<flush;
start_time();
mult=dim(peigenspace(wq,e*den,MODULUS));
stop_time();
show_time(cerr); cerr<<endl;
cout<<"Dimension of "<<e<<"-eigenspace="<<mult<<endl;
*/
cout<<"Using sparse matrix code..."<<endl;
start_time();
mult=dim(eigenspace(swq,e*den));
stop_time();
show_time(cerr); cerr<<endl;
cout<<"Dimension of "<<e<<"-eigenspace="<<mult<<endl;
}
}
wqlist[i]=wq;
#ifdef CHECK_COMMUTE
if (mult_mod_p(swq,swq,MODULUS)==den*den*sidmat(genus))
cout << "Involution!" << "\n";
else
cout << "NOT an involution...." << "\n";
#else
cout<<endl;
#endif
}
int np=5,ip=0;
cout<<"How many T_p? "; cin>>np;
mat_m* tplist = new mat_m[np];
for (primevar pr(np+nq); pr.ok()&&ip<np; pr++, ip++)
{while (n%pr==0) pr++;
int p=pr;
cout << "\nComputing T_p for p = " << p << "..." << flush;
#ifdef COMPARE_OLD
cout<<endl;
start_time();
mat_m temp = hplus.heckeop(p,verbose);
stop_time();
cout<<"Time for old method: "; show_time(cerr); cerr<<endl;
#endif // COMPARE_OLD
start_time();
mat tp = hplus.newheckeop(p,verbose);
if(verbose)
{
cout<<"Computed matrix ";
if(den>1) cout<<" (scaled by "<<den<<") ";
cout<<" = "<<tp<<endl;
// tp.output_pretty();
}
tplist[ip] = tp;
// cout<<"Copied mat_m = "<<tplist[ip]<<endl;
stop_time();
#ifdef COMPARE_OLD
cout<<"Time for new method: "; show_time(cerr); cerr<<endl;
if(temp!=tplist[ip]) cout<<"Matrices differ!\n";
#else
cout << "done, sparsity = "<<sparsity(tplist[ip])<<". " << endl;
#endif // COMPARE_OLD
#ifdef TEST_EIGS
vector<long> eigs = eigrange(p); // hplus.eigrange(nq+ip);
cout<<"\nChecking for eigenvalues from "<<eigs<<endl;
long i,j,k,n=genus,r;
long nulty, nulty1, totalmult=0;
SCALAR dummy;
mat m = tplist[ip].shorten(dummy);
// if(verbose) cout<<"shortened matrix: \n"<<m<<endl;
smat sm=smat(m);
for(k=0; k<eigs.size(); k++)
{
long e = eigs[k];
cout<<"\nTrying eigenvalue e = "<<e<<" ("<<e<<")"<<endl;
e *= den;
/*
cout<<"Computing nullity, using my (modular) matrix code..."<<flush;
start_time();
nulty=dim(peigenspace(m,e,MODULUS));
stop_time();
show_time(cerr); cerr<<endl;
cout<<" nullity="<<nulty<<endl;
*/
cout<<"Computing nullity, using my sparse matrix code..."<<flush;
start_time();
nulty=dim(eigenspace(sm,e));
stop_time();
show_time(cerr); cerr<<endl;
cout<<" nullity="<<nulty<<endl;
cout<<"\n"<<e<<" (scaled) ";
if(nulty>0) cout<<" IS "; else cout<<" is NOT ";
cout<<"an eigenvalue";
if(nulty>0) cout<<", of multiplicity "<<nulty;
cout<<endl;
}
mat_ZZ M;
M.SetDims(n,n);
for(i=1; i<=n; i++)
for(j=1; j<=n; j++)
M(i,j)=m(i,j);
// cout<<"NTL matrix = "<<M<<endl;
// Evaluate "rational" charpoly of T_p:
cout<<"Computing product of (T-a*I) over possible eigs a..."<<flush;
start_time();
mat mm=m;
mat MT=mm, m2=m*m;
for(k=1; k<(eigs.size()); k++)
{
long e = eigs[k]*den;
if(e>0)
MT = MT * addscalar(m2,(-e*e));
}
stop_time();
cout<<"...done, sparsity = "<<sparsity(MT); show_time(cerr); cerr<<endl;
cout<<"Computing kernel, using my (modular) matrix code..."<<flush;
start_time();
subspace ker = pkernel(MT,MODULUS);
nulty=dim(ker);
int denker=denom(ker); // =1 for modular method, but set below
stop_time();
cout<<"done, nulty = "<<nulty; show_time(cerr); cerr<<endl;
cout<<endl;
vector<long> eigs1;
long n1f = 0;
if(nulty>0)
{
cout<<"lifting kernel..."<<flush;
start_time();
mat MTR;
int ok = liftmat(prestrict(m,ker,MODULUS),MODULUS,MTR,denker);
if (!ok)
cout << "**!!!** failed to lift modular kernel\n" << endl;
stop_time();
cout<<"done, denom(ker)="<<denker; show_time(cerr); cout<<endl;
if(nulty<21)
cout<<"Restriction of Tp to relevant subspace = \n" << MTR << endl;
mat_ZZ Msub;
Msub.SetDims(nulty,nulty);
for(i=1; i<=nulty; i++)
for(j=1; j<=nulty; j++)
Msub(i,j)=MTR(i,j);
ZZX cptp; ZZ cont;
cout<<"computing char poly..."<<flush;
start_time();
CharPoly(cptp, Msub);
stop_time();
cout<<"done "; show_time(cerr); cerr<<endl;
vec_pair_ZZX_long factors;
cout<<"\nfactorizing char poly..."<<flush;
start_time();
factor(cont,factors,cptp);
stop_time();
cout<<"done "; show_time(cerr); cerr<<endl;
cout<<"\nFactors are:"<<endl;
long nf = factors.length();
for(i=0; i<nf; i++)
{
cout<<(i+1)<<":\t"<<factors[i].a
<<"\t(degree "<<deg(factors[i].a)<<")"
<<"\t to power "<<factors[i].b;
cout<<endl;
if(deg(factors[i].a)==1)
{
long ap = -I2long(coeff(factors[i].a,0))/denker;
cout<<"Adding eigenvalue "<<ap<<endl;
eigs1.push_back(ap); n1f++;
}
}
cout<<"Rational eigenvalues (scaled by "<<den<<") are "<<eigs1<<endl;
for(k=0; k<n1f; k++)
{
long e = eigs1[k];
cout<<"\nTrying eigenvalue e = "<<e<<" ("<<e<<")"<<endl;
e *= denker;
ZZ ee; ee = e;
for(i=1;i<=nulty;i++) Msub(i,i)=MTR(i,i)-ee;
if(nulty<21) cout<<"Adjusted NTL matrix = "<<Msub<<endl;
ZZ detM;
/*
start_time();
cout<<"Computing determinant, using randomized strategy..."<<flush;
detM = determinant(Msub);
stop_time();
show_time(cerr); cerr<<endl;
cout<<"det="<<detM<<endl;
*/
cout<<"Computing determinant, using deterministic strategy..."<<flush;
start_time();
detM = determinant(Msub,1);
stop_time();
show_time(cerr); cerr<<endl;
cout<<"det="<<detM<<endl;
/*
cout<<"Computing rank, using NTL's LLL..."<<flush;
start_time();
mat_ZZ M2; M2=Msub;
r=image(detM,M2); // NB image() changes its 2nd arg!
stop_time();
show_time(cerr); cerr<<endl;
nulty1 = nulty-r;
cout<<" nullity="<<nulty1<<endl;
*/
cout<<"Computing nullity, using my (modular) matrix code..."<<flush;
start_time();
nulty1=dim(peigenspace(MTR,e,MODULUS));
stop_time();
show_time(cerr); cerr<<endl;
cout<<" nullity="<<nulty1<<endl;
cout<<"\n"<<e<<" (scaled) ";
if(nulty>0) cout<<" IS "; else cout<<" is NOT ";
cout<<"an eigenvalue";
if(nulty>0) cout<<", of multiplicity "<<nulty1;
cout<<endl;
totalmult+=nulty1;
}
}
cout<<"Total multiplicity of rational eigenvalues = "<<totalmult<<endl;
#endif // TEST_EIGS
#ifdef CHECK_COMMUTE
bigint P = BIGINT(MODULUS);
for (int kp=firstq; kp<nq; kp++)
{if (matmulmodp(wqlist[kp],tplist[ip],P)!=matmulmodp(tplist[ip],wqlist[kp],P))
{
cout << "Problem: T_p matrix "<<ip<<" and W_q matrix "<<kp<<" do not commute!" << "\n";
}
}
for (int jp=0; jp<ip; jp++)
{if (matmulmodp(tplist[ip],tplist[jp],P)!=matmulmodp(tplist[jp],tplist[ip],P))
{
cout << "Problem: T_p matrices "<<ip<<" and "<<jp<<" do not commute!" << "\n";
}
}
#endif //CHECK_COMMUTE
} // loop on p
delete[] wqlist; delete[] tplist;
} // end of if(genus>0)
} // end of if(n)
} // end of while(n>0) or while(n<limit)
cout<<endl;
exit(0);
} // end of main()
double sparsity(const mat_m& m)
{
double count=0;
long i,j,nr=m.nrows(), nc=m.ncols();
for(i=0; i<nr; i++)
for(j=0; j<nc; j++)
if(!is_zero(m(i+1,j+1))) count=count+1;
return count/(nr*nc);
}