CoCalc Shared Fileswww / mazur / katz / katz.m
Author: William A. Stein
1// \begin{verbatim}
2PI  := Pi(ComplexField());
3i := ComplexField().1;
4
5intrinsic 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]);
9end intrinsic;
10
11intrinsic Angle(p::RngIntElt, a::RngIntElt) -> FldPrElt
12{}
13   printf "[%o,%o]",p,a;
14   return Arccos(Kl(p,a)/(2*Sqrt(p)));
15end intrinsic;
16
17intrinsic 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;
25end intrinsic;
26
27intrinsic AnglePrim(p::RngIntElt) -> FldPrElt
28{}
29   a := SmallestPrimitiveRoot(p);
30   return Angle(p,a);
31end intrinsic;
32
33intrinsic RoundPrec(x::FldPrElt, prec::RngIntElt) -> FldRatElt
34{}
35   return Round(x*10^prec)/10^prec;
36end intrinsic;
37
38intrinsic Multiplicities(z::SeqEnum) -> SeqEnum
39{Return a sequence of pairs <x,multiplicity of x in z>,
40where the "multiplicities" are normalized so the sum
41of 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;
62end intrinsic;
63
64intrinsic 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;
82end intrinsic;
83
84
85intrinsic 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;
113end intrinsic;
114
115intrinsic 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;
124end intrinsic;
125
126intrinsic 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;
135end intrinsic;
136
137intrinsic 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;
145end intrinsic;
146
147intrinsic 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;
159end intrinsic;
160
161// \end{verbatim}
162
163
164////////////////////////////////////////////////////////////////
165
166// Up stuff
167
168/*
169
170P.S. While things seem to be on a roll, Here are is another completely
171different thing that I was vaguely hoping to display at this Katz lecture,
172but left to my own devices I probably wouldn't.   Consider j  the
173elliptical modular function.  Let  p=5. I'd love to compute (mod p^4 = 725,
174say) the power series
175
176         1/p U_p(j-744)
177
178as a polynomial divided by a power of j. Ditto for 1/p^2 U_p^2(j-744).
179Ditto 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
184exists 5-adically (it is the eigenfunction that Gouvea and I call f_1 in
185our paper "Searching for eigenfunctions") and has an expansion as a
186polynomial in j plus a power series in j^{-1}.
187
188*/
189
190p := 5;
192pow := 10;
193K := IntegerRing(p^pow);
194R<q> := LaurentSeriesRing(K);
195j := jInvariant(q+O(q^2000));
196
197intrinsic 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));
201end intrinsic;
202
203intrinsic Val(x::RngIntResElt) -> RngIntElt
204{}
205   return Max([i : i in [0..pow] | Integers(p^i)!x eq 0]);
206end intrinsic;
207
208intrinsic 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;
218print "h = ", h;
219print "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);
223print "F = ", F;
224         h := h - Coefficient(h,i)*j^(-i);
225      end for;
226      prec := AbsolutePrecision(h)-1;
227      print "prec = ", prec;
228print "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;
236end intrinsic;
237
238intrinsic 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;
256end intrinsic;
257