/* Code to compute the operator U_p acting on overconvergent $p$-adic modular forms of level 1 and weight 0.

Version 0.1, 14th February 2008.

Written by L. J. P. Kilford, 2007-2008, based on ideas of Kevin Buzzard.

Also includes William Stein's port of newtonslopes from PARI

*************

Examples of use:

> Attach("Up_acting_on_overconvergent_forms_of_genus_1.m");
> PR<x>:=PolynomialRing(Integers());
> time U17_dim6:=UpMatrix(17,6,17^13,13);
Time: 1.790
> NewtonSlopes(PR!CharacteristicPolynomial(U17_dim6),17);
[* 0, 0, 1, 2, 2, 3, 4, Infinity, Infinity, Infinity, Infinity, Infinity, Infinity*]
> time U19_dim5:=UpMatrix(19,5,19^11,11);
Time: 1.390
> NewtonSlopes(PR!CharacteristicPolynomial(U19_dim5),19);
[* 0, 0, 1, 1, 2, 3, 3, Infinity, Infinity, Infinity, Infinity*]

*/

intrinsic 
MakeTestMatrix(p::RngIntElt, r::RngIntElt,N::RngIntElt,acc::RngIntElt,z::RngSerPuisElt,Z::RngSerPuisElt : verbose:=false) -> Any
{This computes a matrix version of the basis [z,Z,z^2,Z^2,...]}
if verbose then
print "now running MakeTestMatrix";
end if;
psr<qq>:=Parent(z);
ZmodN:=quo<Integers()|ideal<Integers()|N>>;
PSRN<qN>:=PuiseuxSeriesRing(ZmodN,p*5*acc);
elts:=[];
elts:=elts cat [1] cat Eltseq(RMatrixSpace(ZmodN,acc-1,1)!0);
for i:=1 to r do
elts:=elts cat Eltseq(PSRN!(z^i+O(qq^acc)));
elts:=elts cat Eltseq(PSRN!((Z/p)^i+O(qq^acc)));
end for;
return elts;
end intrinsic;

intrinsic U_p(powser::RngSerPuisElt,acc::RngIntElt,p::RngIntElt)->RngSerPuisElt
{This computes the action of the U_p operator on power series}
psr:=Parent(powser);
temp:=0;
for i:=Valuation(powser) to acc do
if i mod p eq 0 then
temp:=temp+Coefficient(powser,i)*(psr.1)^(i/p);
end if;
end for;
return temp;
end intrinsic; 

intrinsic 
GetRepresentation_z(i::RngIntElt,r::RngIntElt,z::Any,Z::Any,mTM::Any,acc::RngIntElt : p) -> Any
{This computes z^i in terms of z^i and (p/z)^j}
rms_acc:=RMatrixSpace(Parent(mTM[1]),2*r+2,acc);
full_mat:=Eltseq(PuiseuxSeriesRing(Parent(mTM[1]),p*5*acc)!U_p(z^i,p*(acc-1),p)) cat mTM;
ns:=Basis(NullSpace(rms_acc!full_mat));
if #ns ne 1 then
error "Error: Not enough accuracy to compute the nullspace";
end if;
if ns[1][1] ne 1 then
print "Warning: not enough p-adic accuracy in GetRepresentation_z";
end if;
return ns[1];
end intrinsic;

intrinsic
GetRepresentation_Z(i::RngIntElt,r::RngIntElt,z::Any,Z::Any,mTM::Any,acc::RngIntElt : p) -> Any
{This computes (p/z)^i in terms of z^i and (p/z)^j}
rms_acc:=RMatrixSpace(Parent(mTM[1]),2*r+2,acc);
full_mat:=Eltseq(PuiseuxSeriesRing(Parent(mTM[1]),11*5*acc)!U_p((Z/p)^i,p*(acc-1),p)) cat mTM;
ns:=Basis(NullSpace(rms_acc!full_mat));
if #ns ne 1 then
error "Error: Not enough accuracy to compute the nullspace";
end if;
if ns[1][1] ne 1 then
print "Warning: not enough p-adic accuracy in GetRepresentation_Z";
end if;
return ns[1];
end intrinsic;

intrinsic UpMatrix(p::RngIntElt,dim::RngIntElt,N::RngIntElt,acc::RngIntElt : verbose:=false, 
degree:=500)->Any
{This computes the action of the U_p operator on overconvergent modular forms of level 1 and 
weight 0, for p in [11,17,19]}
assert p in [11,17,19];
ZmodN:=quo<Integers()|ideal<Integers()|N>>;
MA:=MatrixAlgebra(ZmodN,2*dim+1)!0;
MA[1,1]:=1;
PSR<qq>:=PuiseuxSeriesRing(Rationals(),degree);
j:=jInvariant(qq);
if p eq 11 then
a:=1;b:=0;
elif p eq 17 then
a:=8;b:=0;
else
a:=7;b:=18;
end if;
z:=(j-a)/(j-b);
Z:=p/z;
mTM:=MakeTestMatrix(p, dim,N,acc,z,Z: verbose:=verbose);
for i:=2 to dim+1 do
col:=GetRepresentation_z(i-1,dim,z,Z,mTM,acc: p:=p);
col_Z:=GetRepresentation_Z(i-1,dim,z,Z,mTM,acc: p:=p);
if verbose then
print col;
print col_Z;
end if;
for j:=1 to 2*dim do
MA[j,2*i-2]:=col[j+1];
MA[j,2*i+1-2]:=col_Z[j+1];
end for;
end for;
return MA;
end intrinsic;

/* 

********************************

Code to compute Newton slopes of polynomials follows, with attribution

********************************

*/

///////////////////////////////////////////////////////////////////
// Newton Polygons                                               //
// Ported from PARI by William Stein  Feb 12, 2000               //
///////////////////////////////////////////////////////////////////

intrinsic NewtonSlopes(f::RngUPolElt, p::RngIntElt) -> List
{The slopes of the Newton polygon of f with at the prime p, where the slopes are the valuations of the p-adic roots 
of f.}
   assert IsPrime(p);
   n := Degree(f);
   if n le 0 then
      return [];
   end if;

   y := [* *];
   vval := [* *];
   for i in [0..n] do 
      Append(~vval,Type(Parent(Coefficient(f,0))) in [RngInt,FldRat] select
                       Valuation(Coefficient(f,i),p)
                   else
   		       Valuation(Coefficient(f,i)));
      Append(~y,0);
   end for;

   a := 1; ind := 2; 
   while a le n do
      if vval[a] ne Infinity() then
         break;
      end if;
      y[ind] := Infinity();   
      ind +:= 1;
      a +:= 1;
   end while;

   b := a+1;
   while b le n+1 do
      while vval[b] eq Infinity() do
         b +:= 1;
      end while;
      u1 := vval[a] - vval[b];
      u2 := b - a;
      c := b+1;
      while c le n+1 do
         if vval[c] eq Infinity() then
            c +:= 1;
            continue;
         end if;
         r1 := vval[a] - vval[c];
         r2 := c - a;
         if u1*r2 le u2*r1 then
            u1 := r1;
            u2 := r2;
             b := c;
         end if;
         c +:= 1;
      end while;

      while ind le b do 
         if u1 eq 0 then
            y[ind] := 0;
         else
            y[ind] := u1 / u2; 
         end if;
         ind +:= 1;
      end while;
      a := b;
      b := a+1;
   end while;
   Remove(~y,1);
   z := [* *];
   for i in [1..#y] do
      Append(~z,y[#y-i+1]);
   end for;
   return z;
end intrinsic;