// \begin{verbatim}
PI  := Pi(ComplexField());
i := ComplexField().1;

intrinsic Kl(p::RngIntElt, a::RngIntElt) -> FldPrElt
{}
   return Real(&+[ Exp(2*PI*i*(Integers()!x + Integers()!(a/x))/p) : 
                x in GF(p) | x ne 0]);
end intrinsic;

intrinsic Angle(p::RngIntElt, a::RngIntElt) -> FldPrElt
{}
   printf "[%o,%o]",p,a;
   return Arccos(Kl(p,a)/(2*Sqrt(p)));
end intrinsic;

intrinsic SmallestPrimitiveRoot(p::RngIntElt) -> RngIntElt
{}
   F := GF(p);
   for a in [1..p-1] do 
      if Order(F!a) eq p-1 then
         return a;
      end if;
   end for;
end intrinsic;

intrinsic AnglePrim(p::RngIntElt) -> FldPrElt
{}
   a := SmallestPrimitiveRoot(p);
   return Angle(p,a);
end intrinsic;

intrinsic RoundPrec(x::FldPrElt, prec::RngIntElt) -> FldRatElt
{}
   return Round(x*10^prec)/10^prec;
end intrinsic;

intrinsic Multiplicities(z::SeqEnum) -> SeqEnum
{Return a sequence of pairs <x,multiplicity of x in z>, 
where the "multiplicities" are normalized so the sum
of multiplicities is 1.}
   if #z eq 0 then
      return [];
   end if;
   z := Sort(z);
   ans := [];
   cur := z[1];
   m := 0;
   i := 1;
   while i le #z do
      while i le #z and z[i] eq cur do
         i := i + 1;
         m := m + 1;
      end while;
      Append(~ans, <cur, m/#z>);
      m := 0;
      if i le #z then
         cur := z[i];
      end if;
   end while;
   return ans;
end intrinsic;

intrinsic RoundStr(x::RngElt, prec::RngIntElt) -> MonStgElt
{}
   y := Sprintf("%o",x*1.0);
   j := #y;
   for i in [1..#y] do 
      if y[i] eq "." then
         j := i+prec;
         break;
      end if;
   end for;
   a := "";
   if j ge #y then
      return y;
   end if;
   for i in [1..j] do 
      a := a*y[i];
   end for;
   return a;
end intrinsic;


intrinsic MakeFrequencyGraph(data::Tup,
                 bins::RngIntElt) -> MonStgElt
{}
   if #data eq 3 then
      a,pmax,x := Explode(data);
      vary_a := false;
   else
      p, x  := Explode(data);
      vary_a := true;
   end if;

   binsize := PI/bins;
   // sequence of bin numbers of things in x;
   bin_nums := [Round(y/binsize) : y in x];
   m := Multiplicities(bin_nums);
   if vary_a then
      ans := Sprintf("\\graphVaryA{%o}{%o}{",p,bins);   
   else
      ans := Sprintf("\\graph{%o}{%o}{%o}{",a,pmax,bins);
   end if;
   for w in m do 
      ans := ans*Sprintf("\\bar{%o}{%o}{%o}", 
              RoundStr(w[1]*binsize,3),
              RoundStr((1/binsize)*w[2],3),
              RoundStr(binsize,3));
   end for;
   ans := ans*"}\n";
   return ans;
end intrinsic;

intrinsic MakeFrequencyGraph(a::RngIntElt, 
     pmax::RngIntElt,
     bins::RngIntElt) -> MonStgElt, SeqEnum
{}
   print "Computing data.";
   time x := [Angle(p,a) : p in PrimeSeq(2, pmax)];
   print "Done";
   data := < a, pmax, x >;
   return MakeFrequencyGraph(data, bins), data;
end intrinsic;
              
intrinsic MakeFrequencyGraphAPrim(
     pmax::RngIntElt,
     bins::RngIntElt) -> MonStgElt, SeqEnum
{}
   print "Computing data for a equal to first primitive root.";
   time x := [AnglePrim(p) : p in PrimeSeq(2, pmax)];
   print "Done";
   data := < "\\text{First primitive root}", pmax, x >;
   return MakeFrequencyGraph(data, bins), data;
end intrinsic;

intrinsic MakeFrequencyGraphVaryA(p::RngIntElt,
     bins::RngIntElt) -> MonStgElt, Tup
{}
   print "Computing data.";
   time x := [Angle(p,a) : a in [1..p-1] | GCD(a,p) eq 1];
   print "Done";
   data := < p, x >;
   return MakeFrequencyGraph(data, bins), data;
end intrinsic;

intrinsic SatoTate(sample_points::RngIntElt) -> MonStgElt
{pstricks graph of Sato-Tate distribution 2/pi*Sin^2(theta).}
   step := PI/sample_points;
   x := 0.0001;
   ans := "\\pscurve[linecolor=red]";
   for i in [1..sample_points] do 
      ans := ans*Sprintf("(%o,%o)", RoundStr(x,3), 
                RoundStr((2/PI)*Sin(x)^2, 3) );
      x := x + step;
   end for;
   ans := ans * "\n";
   return ans;
end intrinsic;

// \end{verbatim}


////////////////////////////////////////////////////////////////

// Up stuff

/*

P.S. While things seem to be on a roll, Here are is another completely
different thing that I was vaguely hoping to display at this Katz lecture,
but left to my own devices I probably wouldn't.   Consider j  the
elliptical modular function.  Let  p=5. I'd love to compute (mod p^4 = 725,
say) the power series

         1/p U_p(j-744)

as a polynomial divided by a power of j. Ditto for 1/p^2 U_p^2(j-744).
Ditto for
 1/p^3 U_p^3(j-744).  (Until they settle down mod 725.)  In fact, the limit

 Lim_n 1/p^n U_p^n(j-744)

exists 5-adically (it is the eigenfunction that Gouvea and I call f_1 in
our paper "Searching for eigenfunctions") and has an expansion as a
polynomial in j plus a power series in j^{-1}.

*/

p := 5;
//K := pAdicField(p);
pow := 10;
K := IntegerRing(p^pow);
R<q> := LaurentSeriesRing(K);
j := jInvariant(q+O(q^2000));

intrinsic Up(f::RngSerLaurElt) -> RngSerLaurElt
{}
   d := Floor((AbsolutePrecision(f)-1)/p);
   return &+[Coefficient(f,n*p)*q^n : n in [1..d]] + O(q^(d+1));
end intrinsic;

intrinsic Val(x::RngIntResElt) -> RngIntElt
{}
   return Max([i : i in [0..pow] | Integers(p^i)!x eq 0]);
end intrinsic;

intrinsic InJ(g::RngSerLaurElt, precision::RngIntElt) -> RngIntElt, RngUPolElt
{}
   assert Valuation(g) gt 0;
   r := 1;
   while true do
      print r;
      h := j^r*g;
      v := [];
      S<X> := PolynomialRing(K);
      F := 0;
print "h = ", h;
print "Valuation(h) = ", Valuation(h);
      for i in [Valuation(h)..0] do 
         Append(~v,Coefficient(h,i));
         F := F + Coefficient(h,i)*X^(-i);
print "F = ", F;
         h := h - Coefficient(h,i)*j^(-i);
      end for;
      prec := AbsolutePrecision(h)-1;
      print "prec = ", prec;
print "h = ", h;
      z := Min([Val(Coefficient(h,n)) : n in [1..prec]]);
      print [Val(Coefficient(h,n)) : n in [1..prec]];
      if z ge precision then
         return F, r;
      end if;
      r := r + 1;
   end while;
end intrinsic;

intrinsic f1(n::RngIntElt) -> ., .
{}
   f := j - 744;
   for i in [1..n] do 
      f := Up(f);
   end for;
   r := p^n;
   prec := AbsolutePrecision(f);
   f := &+[(Integers()!(Coefficient(f,m)) div r)*q^m : m in 
             [Valuation(f)..prec-1]] + O(q^prec);
   printf "1/p^%o Up^%o(j-744) = %o", n, n, f;
   print "146560/j = ", 146560/j + O(q^100);
//h := 146560/j + O(q^100) - f;
//   print "dif = ", h;
//   z := [Val(Coefficient(h,n)) : n in [1..30]];
//print "z = ", z;
   F, r :=  InJ(f,4);
   return F, r;
end intrinsic;
