CoCalc Shared Fileswww / mazur / katz / katz.mOpen in CoCalc with one click!
Author: William A. Stein
1
// \begin{verbatim}
2
PI := Pi(ComplexField());
3
i := ComplexField().1;
4
5
intrinsic Kl(p::RngIntElt, a::RngIntElt) -> FldPrElt
6
{}
7
return Real(&+[ Exp(2*PI*i*(Integers()!x + Integers()!(a/x))/p) :
8
x in GF(p) | x ne 0]);
9
end intrinsic;
10
11
intrinsic Angle(p::RngIntElt, a::RngIntElt) -> FldPrElt
12
{}
13
printf "[%o,%o]",p,a;
14
return Arccos(Kl(p,a)/(2*Sqrt(p)));
15
end intrinsic;
16
17
intrinsic SmallestPrimitiveRoot(p::RngIntElt) -> RngIntElt
18
{}
19
F := GF(p);
20
for a in [1..p-1] do
21
if Order(F!a) eq p-1 then
22
return a;
23
end if;
24
end for;
25
end intrinsic;
26
27
intrinsic AnglePrim(p::RngIntElt) -> FldPrElt
28
{}
29
a := SmallestPrimitiveRoot(p);
30
return Angle(p,a);
31
end intrinsic;
32
33
intrinsic RoundPrec(x::FldPrElt, prec::RngIntElt) -> FldRatElt
34
{}
35
return Round(x*10^prec)/10^prec;
36
end intrinsic;
37
38
intrinsic Multiplicities(z::SeqEnum) -> SeqEnum
39
{Return a sequence of pairs <x,multiplicity of x in z>,
40
where the "multiplicities" are normalized so the sum
41
of multiplicities is 1.}
42
if #z eq 0 then
43
return [];
44
end if;
45
z := Sort(z);
46
ans := [];
47
cur := z[1];
48
m := 0;
49
i := 1;
50
while i le #z do
51
while i le #z and z[i] eq cur do
52
i := i + 1;
53
m := m + 1;
54
end while;
55
Append(~ans, <cur, m/#z>);
56
m := 0;
57
if i le #z then
58
cur := z[i];
59
end if;
60
end while;
61
return ans;
62
end intrinsic;
63
64
intrinsic RoundStr(x::RngElt, prec::RngIntElt) -> MonStgElt
65
{}
66
y := Sprintf("%o",x*1.0);
67
j := #y;
68
for i in [1..#y] do
69
if y[i] eq "." then
70
j := i+prec;
71
break;
72
end if;
73
end for;
74
a := "";
75
if j ge #y then
76
return y;
77
end if;
78
for i in [1..j] do
79
a := a*y[i];
80
end for;
81
return a;
82
end intrinsic;
83
84
85
intrinsic MakeFrequencyGraph(data::Tup,
86
bins::RngIntElt) -> MonStgElt
87
{}
88
if #data eq 3 then
89
a,pmax,x := Explode(data);
90
vary_a := false;
91
else
92
p, x := Explode(data);
93
vary_a := true;
94
end if;
95
96
binsize := PI/bins;
97
// sequence of bin numbers of things in x;
98
bin_nums := [Round(y/binsize) : y in x];
99
m := Multiplicities(bin_nums);
100
if vary_a then
101
ans := Sprintf("\\graphVaryA{%o}{%o}{",p,bins);
102
else
103
ans := Sprintf("\\graph{%o}{%o}{%o}{",a,pmax,bins);
104
end if;
105
for w in m do
106
ans := ans*Sprintf("\\bar{%o}{%o}{%o}",
107
RoundStr(w[1]*binsize,3),
108
RoundStr((1/binsize)*w[2],3),
109
RoundStr(binsize,3));
110
end for;
111
ans := ans*"}\n";
112
return ans;
113
end intrinsic;
114
115
intrinsic MakeFrequencyGraph(a::RngIntElt,
116
pmax::RngIntElt,
117
bins::RngIntElt) -> MonStgElt, SeqEnum
118
{}
119
print "Computing data.";
120
time x := [Angle(p,a) : p in PrimeSeq(2, pmax)];
121
print "Done";
122
data := < a, pmax, x >;
123
return MakeFrequencyGraph(data, bins), data;
124
end intrinsic;
125
126
intrinsic MakeFrequencyGraphAPrim(
127
pmax::RngIntElt,
128
bins::RngIntElt) -> MonStgElt, SeqEnum
129
{}
130
print "Computing data for a equal to first primitive root.";
131
time x := [AnglePrim(p) : p in PrimeSeq(2, pmax)];
132
print "Done";
133
data := < "\\text{First primitive root}", pmax, x >;
134
return MakeFrequencyGraph(data, bins), data;
135
end intrinsic;
136
137
intrinsic MakeFrequencyGraphVaryA(p::RngIntElt,
138
bins::RngIntElt) -> MonStgElt, Tup
139
{}
140
print "Computing data.";
141
time x := [Angle(p,a) : a in [1..p-1] | GCD(a,p) eq 1];
142
print "Done";
143
data := < p, x >;
144
return MakeFrequencyGraph(data, bins), data;
145
end intrinsic;
146
147
intrinsic SatoTate(sample_points::RngIntElt) -> MonStgElt
148
{pstricks graph of Sato-Tate distribution 2/pi*Sin^2(theta).}
149
step := PI/sample_points;
150
x := 0.0001;
151
ans := "\\pscurve[linecolor=red]";
152
for i in [1..sample_points] do
153
ans := ans*Sprintf("(%o,%o)", RoundStr(x,3),
154
RoundStr((2/PI)*Sin(x)^2, 3) );
155
x := x + step;
156
end for;
157
ans := ans * "\n";
158
return ans;
159
end intrinsic;
160
161
// \end{verbatim}
162
163
164
////////////////////////////////////////////////////////////////
165
166
// Up stuff
167
168
/*
169
170
P.S. While things seem to be on a roll, Here are is another completely
171
different thing that I was vaguely hoping to display at this Katz lecture,
172
but left to my own devices I probably wouldn't. Consider j the
173
elliptical modular function. Let p=5. I'd love to compute (mod p^4 = 725,
174
say) the power series
175
176
1/p U_p(j-744)
177
178
as a polynomial divided by a power of j. Ditto for 1/p^2 U_p^2(j-744).
179
Ditto for
180
1/p^3 U_p^3(j-744). (Until they settle down mod 725.) In fact, the limit
181
182
Lim_n 1/p^n U_p^n(j-744)
183
184
exists 5-adically (it is the eigenfunction that Gouvea and I call f_1 in
185
our paper "Searching for eigenfunctions") and has an expansion as a
186
polynomial in j plus a power series in j^{-1}.
187
188
*/
189
190
p := 5;
191
//K := pAdicField(p);
192
pow := 10;
193
K := IntegerRing(p^pow);
194
R<q> := LaurentSeriesRing(K);
195
j := jInvariant(q+O(q^2000));
196
197
intrinsic Up(f::RngSerLaurElt) -> RngSerLaurElt
198
{}
199
d := Floor((AbsolutePrecision(f)-1)/p);
200
return &+[Coefficient(f,n*p)*q^n : n in [1..d]] + O(q^(d+1));
201
end intrinsic;
202
203
intrinsic Val(x::RngIntResElt) -> RngIntElt
204
{}
205
return Max([i : i in [0..pow] | Integers(p^i)!x eq 0]);
206
end intrinsic;
207
208
intrinsic InJ(g::RngSerLaurElt, precision::RngIntElt) -> RngIntElt, RngUPolElt
209
{}
210
assert Valuation(g) gt 0;
211
r := 1;
212
while true do
213
print r;
214
h := j^r*g;
215
v := [];
216
S<X> := PolynomialRing(K);
217
F := 0;
218
print "h = ", h;
219
print "Valuation(h) = ", Valuation(h);
220
for i in [Valuation(h)..0] do
221
Append(~v,Coefficient(h,i));
222
F := F + Coefficient(h,i)*X^(-i);
223
print "F = ", F;
224
h := h - Coefficient(h,i)*j^(-i);
225
end for;
226
prec := AbsolutePrecision(h)-1;
227
print "prec = ", prec;
228
print "h = ", h;
229
z := Min([Val(Coefficient(h,n)) : n in [1..prec]]);
230
print [Val(Coefficient(h,n)) : n in [1..prec]];
231
if z ge precision then
232
return F, r;
233
end if;
234
r := r + 1;
235
end while;
236
end intrinsic;
237
238
intrinsic f1(n::RngIntElt) -> ., .
239
{}
240
f := j - 744;
241
for i in [1..n] do
242
f := Up(f);
243
end for;
244
r := p^n;
245
prec := AbsolutePrecision(f);
246
f := &+[(Integers()!(Coefficient(f,m)) div r)*q^m : m in
247
[Valuation(f)..prec-1]] + O(q^prec);
248
printf "1/p^%o Up^%o(j-744) = %o", n, n, f;
249
print "146560/j = ", 146560/j + O(q^100);
250
//h := 146560/j + O(q^100) - f;
251
// print "dif = ", h;
252
// z := [Val(Coefficient(h,n)) : n in [1..30]];
253
//print "z = ", z;
254
F, r := InJ(f,4);
255
return F, r;
256
end intrinsic;
257