Sharedwww / Tables / hijikata.mOpen in CoCalc
Author: William A. Stein
1
/*****************************************************
2
tr.m
3
4
AUTHOR -- William A. Stein
5
CREATED: 14 January 2001
6
7
DESCRIPTION: Hijikata trace formula for tr(T_n)
8
on S_k(Gamma_0(N)), any n>=1, except if n|N,
9
then n must be prime.
10
11
*****************************************************/
12
13
//classgroup_algorithm := "ClassGroup";
14
//classgroup_algorithm := "Shanks"; // heuristic -- can be *wrong*!!
15
classgroup_algorithm := "ReducedForms";
16
17
function w(d)
18
case d:
19
when -4:
20
return 2;
21
when -3:
22
return 3;
23
end case;
24
return 1;
25
end function;
26
27
28
function tof(a)
29
t := &*[pn[1]^((pn[2]-(pn[2] mod 2)) div 2) :
30
pn in Factorization(a)];
31
32
if (Integers()!(a/t^2) mod 4) eq 1 then
33
return t;
34
else
35
return Integers()!(t/2);
36
end if;
37
end function;
38
39
40
function mu(N)
41
return N *
42
&*[Rationals()|1+1/p[1] : p in Factorization(N)];
43
end function;
44
45
46
function sig(n, N)
47
if N mod n eq 0 then
48
return n;
49
else
50
return &*[ (1 -p[1]^(p[2]+1))/(1-p[1]) : p in Factorization(n)];
51
end if;
52
end function;
53
54
55
function quadpoints(s,n,p,v)
56
// WARNING: This is a very very stupid algorithm!
57
pv := p^v;
58
return [x : x in [0..pv-1] | (x^2 - s*x + n) mod pv eq 0];
59
end function;
60
61
62
function A(s,f,n,p,v)
63
rho := Valuation(f,p);
64
sr := Set([x mod p^(v+rho) : x in quadpoints(s,n,p,v+2*rho)]);
65
return #[x : x in sr |
66
(2*x-s) mod p^rho eq 0
67
and ((n ne p) or ((n eq p) and x mod p ne 0))];
68
end function;
69
70
71
function B(s,f,n,p,v)
72
rho := Valuation(f,p);
73
return #Set([x mod p^(v+rho) : x in quadpoints(s,n,p,v+2*rho+1)]);
74
end function;
75
76
77
function cp(s,f,n,p,v)
78
ans := A(s,f,n,p,v);
79
if (Integers()!((s^2-4*n)/f^2)) mod p eq 0 then
80
ans +:= B(s,f,n,p,v);
81
end if;
82
return ans;
83
end function;
84
85
86
function c(s,f,n,N)
87
return &*[Integers()| cp(s,f,n,p[1],p[2]) : p in Factorization(N)];
88
end function;
89
90
91
function type_p(n,k,N)
92
if IsSquare(n) then
93
s := Floor(Sqrt(4*n));
94
ans := 1/4*(s/2)*n^((k div 2)-1)*(c(s,1,n,N) + (-1)^k*c(-s,1,n,N));
95
return ans;
96
else
97
return 0;
98
end if;
99
end function;
100
101
102
function absxy(s,n,k,N)
103
t := Floor(Sqrt(s^2-4*n));
104
x := (s-t)/2;
105
y := (s+t)/2;
106
ans := (Min([Abs(x), Abs(y)])^(k-1)/Abs(x-y)) * Sign(x)^k *
107
&+[1/2*EulerPhi(Integers()!(t/f))*c(s,f,n,N) : f in Divisors(t)];
108
return ans;
109
end function;
110
111
112
function type_h(n,k,N)
113
start := Ceiling(2*Sqrt(n));
114
if IsSquare(n) then
115
start +:= 1;
116
end if;
117
ans := &+[Rationals()|absxy(s,n,k,N) + absxy(-s,n,k,N) :
118
s in [start..4*n] | IsSquare(s^2-4*n)];
119
return ans;
120
end function;
121
122
function xy(s,n,k,N)
123
K<a> := QuadraticField(s^2-4*n);
124
D := Discriminant(K); // a = sqrt(D) or sqrt(D/4).
125
if D mod 4 eq 0 then
126
D := D div 4;
127
end if;
128
b := Round(Sqrt((s^2-4*n)/D));
129
x := (s+a*b)/2;
130
y := (s-a*b)/2;
131
ans := 1/2 * (x^(k-1)-y^(k-1))/(x-y)
132
* &+[ClassNumber(Integers()!((s^2-4*n)/f^2) : Al:=classgroup_algorithm)
133
/ w((s^2-4*n)/f^2) * c(s,f,n,N) :
134
f in Divisors(tof(s^2-4*n))];
135
return Rationals()!ans;
136
end function;
137
138
139
function type_e(n,k,N)
140
r := Floor(2*Sqrt(n));
141
if IsSquare(n) then
142
r -:= 1;
143
end if;
144
145
// return &+[xy(s,n,k,N) : s in [-r..r]];
146
// WARNING: I *conjecture* that the sum below is the same.
147
148
return xy(0,n,k,N) + 2* &+[xy(s,n,k,N) : s in [1..r]];
149
end function;
150
151
152
function sum_s(n,k,N)
153
return type_p(n,k,N)+type_h(n,k,N)+type_e(n,k,N);
154
end function;
155
156
157
intrinsic TraceHeckeOperator(n::RngIntElt,
158
k::RngIntElt) -> RngIntElt
159
{}
160
require n ge 1 : "Argument 1 must be at least 1.";
161
require k ge 2 : "Argument 2 must be at least 2.";
162
return TraceHeckeOperator(n,k,1);
163
end intrinsic;
164
165
intrinsic TraceHeckeOperator(n::RngIntElt,
166
k::RngIntElt,
167
N::RngIntElt) -> RngIntElt
168
{The trace of the Hecke operator T_n on S_k(Gamma_0(N))).
169
The only constraint on n is that if GCD(n,N) ne 1,
170
then n must be prime.}
171
require n ge 1 : "Argument 1 must be at least 1.";
172
require k ge 2 : "Argument 2 must be at least 2.";
173
require N ge 1 : "Argument 3 must be at least 1.";
174
175
if GCD(n,N) ne 1 then
176
require IsPrime(n) :
177
"Argument 1 must be prime when the GCD of arguments 1 and 3 is not 1.";
178
end if;
179
180
if k mod 2 ne 0 then
181
return 0;
182
end if;
183
if k eq 2 then
184
t := sig(n,N);
185
else
186
t := 0;
187
end if;
188
t +:= -sum_s(n,k,N);
189
if IsSquare(n) then
190
t +:= (k-1)*mu(N)/12*n^((k div 2)-1);
191
end if;
192
return t;
193
end intrinsic;
194
195
196
intrinsic TraceModularForm(k::RngIntElt, prec::RngIntElt)
197
-> RngSerPowElt
198
{The trace modular form sum Tr(T_n)*q^n to absolute precision prec,
199
where T_n is the nth Hecke operator on S_k(SL_2(Z)). The complexity
200
is almost entirely a function of prec, but *not* of k.}
201
202
R<q>:=PowerSeriesRing(Rationals());
203
return &+[TraceHeckeOperator(n,k,1)*q^n : n in [1..prec-1]] + O(q^(prec));
204
205
end intrinsic;
206
207
intrinsic TraceToFile(k::RngIntElt,
208
start::RngIntElt, stop::RngIntElt,
209
file::File)
210
{}
211
for n in [start..stop] do
212
t := Cputime();
213
tr := TraceHeckeOperator(n,k,1);
214
fprintf file, "+%o*q^%o // time = %o s\n", tr, n, Cputime(t);
215
Flush(file);
216
end for;
217
end intrinsic;
218
219
function Tp(p, r, k, f)
220
assert Type(p) eq RngIntElt;
221
assert Type(r) eq RngIntElt;
222
assert Type(f) eq RngSerPowElt;
223
224
if r gt 1 then
225
return Tp(p,1,k, Tp(p,r-1,k,f)) - p^(k-1)*Tp(p,r-2,k, f);
226
end if;
227
if r eq 1 then
228
R<q> := Parent(f);
229
prec := Floor(((AbsolutePrecision(f)-1)/p) + 1);
230
return &+[R|Coefficient(f,n*p)*q^n + p^(k-1)*Coefficient(f,n)*q^(n*p) :
231
n in [1..prec-1]] + O(q^prec);
232
end if;
233
if r eq 0 then
234
return f;
235
end if;
236
end function;
237
238
239
intrinsic HeckeOperator(n::RngIntElt, k::RngIntElt, f::RngSerPowElt)
240
-> RngSerPowElt
241
{The image T_n(f) of f under the Hecke operator T_n on level-1
242
weight-k modular forms.}
243
require n ge 1 : "Argument 1 must be at least 1.";
244
require k ge 2 : "Argument 2 must be at least 2.";
245
246
for p in Factorization(n) do
247
f := Tp(p[1],p[2],k,f);
248
end for;
249
return f;
250
251
end intrinsic;
252
253
intrinsic TraceFormulaBasis(k::RngIntElt, prec::RngIntElt) -> SeqEnum
254
{Basis for S_k(SL_2(Z)) computed using the trace formula.}
255
256
f := TraceModularForm(k, prec);
257
d := Integers()!Coefficient(f,1); // the dimension
258
B := [HeckeOperator(n,k,f) : n in [1..d]];
259
return B;
260
end intrinsic;
261
262
function BasisMatrix(B)
263
d := #B;
264
I := MatrixAlgebra(Rationals(),d)!0;
265
for r in [1..d] do
266
for c in [1..d] do
267
I[r,c] := Coefficient(B[r],c);
268
end for;
269
end for;
270
return I;
271
end function;
272
273
intrinsic HeckeOperatorMatrix(k::RngIntElt, n::RngIntElt) -> AlgMatElt
274
{Matrix representing the Hecke operator T_n (acting on the
275
right, as usual in MAGMA) on weight-k level-1 modular forms
276
with respect to the trace basis.}
277
d := DimensionCuspFormsGamma0(1,k);
278
B := TraceFormulaBasis(k, n*d^2+1);
279
return HeckeOperatorMatrix(k,n,B);
280
end intrinsic;
281
282
283
intrinsic HeckeOperatorMatrix(k::RngIntElt,
284
n::RngIntElt,
285
B::SeqEnum ) -> AlgMatElt
286
287
{Matrix representing the Hecke operator T_n (acting on the
288
right, as usual in MAGMA) on weight-k level-1 modular forms
289
with respect to the basis B.}
290
d := #B;
291
I := BasisMatrix(B)^(-1);
292
A := Parent(I)!&cat[ [Coefficient(g,i) : i in [1..d]]
293
where g is HeckeOperator(n,k,B[j]) : j in [1..d]];
294
return I*A;
295
end intrinsic;
296
297
298
/*
299
Here's a good test that the program is working:
300
301
for k in [i : i in [48..60] | IsEven(i)] do
302
B:=TraceFormBasis(k,k);
303
I := BasisMatrix(B);
304
D := DiscriminantOfHeckeAlgebra(CS(MS(1,k,+1)));
305
k, Determinant(I)-D;
306
end for;
307
308
*/
309
310