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*!!
15classgroup_algorithm := "ReducedForms";
16
17function w(d)
18   case d:
19      when -4:
20         return 2;
21      when -3:
22         return 3;
23   end case;
24   return 1;
25end function;
26
27
28function 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;
37end function;
38
39
40function mu(N)
41   return N *
42           &*[Rationals()|1+1/p[1] : p in Factorization(N)];
43end function;
44
45
46function 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;
52end function;
53
54
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];
59end function;
60
61
62function 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))];
68end function;
69
70
71function 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)]);
74end function;
75
76
77function 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;
83end function;
84
85
86function c(s,f,n,N)
87   return &*[Integers()| cp(s,f,n,p[1],p[2]) : p in Factorization(N)];
88end function;
89
90
91function 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;
99end function;
100
101
102function 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;
109end function;
110
111
112function 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;
120end function;
121
122function xy(s,n,k,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;
136end function;
137
138
139function 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]];
149end function;
150
151
152function sum_s(n,k,N)
153   return type_p(n,k,N)+type_h(n,k,N)+type_e(n,k,N);
154end function;
155
156
157intrinsic 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);
163end intrinsic;
164
165intrinsic 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))).
169The only constraint on n is that if GCD(n,N) ne 1,
170then 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;
193end intrinsic;
194
195
196intrinsic TraceModularForm(k::RngIntElt, prec::RngIntElt)
197      -> RngSerPowElt
198{The trace modular form sum Tr(T_n)*q^n to absolute precision prec,
199where T_n is the nth Hecke operator on S_k(SL_2(Z)).  The complexity
200is 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
205end intrinsic;
206
207intrinsic 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;
217end intrinsic;
218
219function 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;
236end function;
237
238
239intrinsic 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
242weight-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
251end intrinsic;
252
253intrinsic 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;
260end intrinsic;
261
262function 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;
271end function;
272
273intrinsic 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);
280end intrinsic;
281
282
283intrinsic HeckeOperatorMatrix(k::RngIntElt,
284                              n::RngIntElt,
285                              B::SeqEnum ) -> AlgMatElt
286
287{Matrix representing the Hecke operator T_n (acting on the
288right, as usual in MAGMA) on weight-k level-1 modular forms
289with 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;
295end intrinsic;
296
297
298/*
299Here'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