/* kani.m */

declare attributes ModTupFld :
   dim,
   p,     
   prec,
   k,
   g_oldforms,
   fundamental_newforms, 
   cm_forms,

   n0, n1bar, n2bar, n3,
   vg0, vg1, vg2, vg3,
   
   dirichlet
   ;

intrinsic Print(S::ModTupFld) 
{Print argument 1.}

   printf "S_%o(Gamma(%o)), prec=%o.\n", S`k, S`p, S`prec;

end intrinsic;


intrinsic SkGamma(p::RngIntElt, k::RngIntElt, prec::RngIntElt) -> ModTupFld
{Create the space S_k(Gamma(p)).}

   requirege k,2 ;
   requirege prec, 1 ;
   require IsPrime(p) : "Argument 1 must be a prime.";

   S := VectorSpaceWithBasis([VectorSpace(Rationals(),1)![1]]);
   S`k := k;
   S`p := p;
   S`prec := prec;

   G<eps> :=  DirichletGroup(p,CyclotomicField(EulerPhi(p)));
   S`dirichlet := G;

   return S;

end intrinsic;


intrinsic G_Oldforms(S::ModTupFld) -> List
{ [* 
     [* list of q-expansions of Galois-conjugacy classes of newforms on Gamma_1(p) *],
     [* corresponding characters *],
  *]
}

   if not assigned S`g_oldforms then

      k := S`k;
      p := S`p;
      prec := S`prec;      
 
      Dir<eps> := DirichletGroup(p, CyclotomicField(EulerPhi(p)));
      E := [ eps^i : i in [0..#Dir div 2] | IsEven(eps^i) ];

      Q_p := [* *];
      eps_p := [* *];
      for e in E do
         M := CuspidalSubspace(ModularSymbols(e,k,+1));
         if Dimension(M) eq 0 then
            continue;
         end if;
         D := NewformDecomposition(M);
         for f in D do
            Append(~Q_p, qEigenform(f,prec));
            Append(~eps_p, e);
         end for;
      end for;
      S`g_oldforms := [* Q_p, eps_p *];

   end if;
   return S`g_oldforms;

end intrinsic;


intrinsic FundamentalNewforms(S::ModTupFld) -> List
{ 
  [* list of q-expansions of the newforms on Gamma_0(p^2) that 
     are not twists of newforms on Gamma_1(p). *]
}

   if not assigned S`fundamental_newforms then

      k := S`k;
      p := S`p;
      prec := S`prec;      
 
      M := NewSubspace(CuspidalSubspace(ModularSymbols(p^2,k,CyclotomicField(EulerPhi(p)),+1)));
      Q_psquared := [* *];
      if Dimension(M) gt 0 then
         D := NewformDecomposition(M);
         for f in D do
            Append(~Q_psquared, qEigenform(f,prec));
         end for;
      end if;   

      oldforms, eps := Explode(G_Oldforms(S));
      if #eps eq 0 then
         S`fundamental_newforms := Q_psquared;
         return  Q_psquared;
      end if;

      G<g> := Parent(eps[1]);
      for i in [1..#G] do
         e := g^i;
         ee := (e^2)^(-1);
         if not IsTrivial(e) then
            for j in [1..#eps] do
               if ee eq eps[j] then

                  f := oldforms[j];
                  traceftwist := Trace(Twist(f,e));
                  found_one := false;
                  for k in [1..#Q_psquared] do
                     if k le #Q_psquared and traceftwist eq Trace(Q_psquared[k]) then

                        if found_one eq true then
                           error "Found two!  There is a bug in FundamentalNewforms, or your precision is too low.";
                        end if;
                        found_one := true;
                        Remove(~Q_psquared,k);

                     end if;
                  end for;

                  if found_one eq false then
                     error "Didn't find one!  There is a bug in FundamentalNewforms, or your precision is too low.";
                  end if;

               end if;
            end for;
         end if;
      end for;      

      S`fundamental_newforms := Q_psquared;
   end if;

   return S`fundamental_newforms;
end intrinsic;


intrinsic CM_Newforms(S::ModTupFld) -> SeqEnum
{The CM newforms.}

   if not assigned S`cm_forms then
      S`cm_forms := [];
      G<eps> := DirichletGroup(Level(S), Rationals());
      D := FundamentalNewforms(S);

      for f in D do
         is_cm := true;
         for ell in [l : l in [2..Precision(S)-1] | l ne Level(S) and IsPrime(l)] do
            if Coefficient(f,ell) ne Evaluate(eps,ell)*Coefficient(f,ell) then
               is_cm := false;
               continue;
            end if;           
         end for;
         Append(~S`cm_forms, is_cm);
      end for;

      assert #S`cm_forms eq #D;
   end if;

   D := FundamentalNewforms(S);

   return S`cm_forms, [D[i] : i in [1..#D] | S`cm_forms[i]];

end intrinsic;

procedure PrintForms(X, label)
   R := PolynomialRing(RationalField());
   i := 1;
   for f in X do
      printf "%o%o := %o;\n", label, i, f;
      K := BaseRing(Parent(f));
      if Type(K) ne FldRat and Type(K) ne FldCyc then
         printf "// Where a satisfies %o\n", Modulus(K);
      end if;
      i := i + 1;
   end for;
end procedure;

intrinsic Everything(S::ModTupFld) 
{}
   X := G_Oldforms(S)[1];
   Y := FundamentalNewforms(S);
   print "OLD FORMS:";
   PrintForms(X, "f");
   print "\nFUNDAMENTAL NEWFORMS:";
   PrintForms(Y, "g");
end intrinsic;

intrinsic Dim(S::ModTupFld) -> RngIntElt
{Dimension of S.}
   if not assigned S`dim then
      if Weight(S) gt 2 then
         error "Ernst doesn't have the formula on him, so it's not implemented.";
      end if;
      p := S`p;
      S`dim := 1/24*(p+2)*(p-3)*(p-5);
   end if;
   return S`dim;
end intrinsic;


intrinsic Weight(S::ModTupFld) -> RngIntElt
{Weight of S.}
   return S`k;
end intrinsic;


intrinsic Precision(S::ModTupFld) -> RngIntElt
{Precision of S.}
   return S`prec;
end intrinsic;


intrinsic Level(S::ModTupFld) -> RngIntElt
{Level of S.}
   return S`p;
end intrinsic;



intrinsic N_0(S::ModTupFld) -> SeqEnum
{}
   if not assigned S`n0 then
      O, eps := Explode(G_Oldforms(S));
      S`n0 := [* *];
      for i in [1..#eps] do
         if IsTrivial(eps[i]) then
            Append(~S`n0,O[i]);
         end if;
      end for;
   end if;

   return S`n0;
end intrinsic;

intrinsic N_1bar(S::ModTupFld) -> SeqEnum
{}
   if not assigned S`n1bar then
      O, eps := Explode(G_Oldforms(S));
      S`n1bar := [* *];
      for i in [1..#eps] do
         if not IsTrivial(eps[i]) then
            Append(~S`n1bar,O[i]);
         end if;
      end for;
   end if;

   return S`n1bar;
end intrinsic;

intrinsic N_2bar(S::ModTupFld) -> SeqEnum
{}
   if not assigned S`n2bar then
      iscm := CM_Newforms(S);
      qexp := FundamentalNewforms(S);
      S`n2bar := [* *];
      for i in [1..#iscm] do
         if not iscm[i] then
            Append(~S`n2bar,qexp[i]);
         end if;
      end for;

      G<eps>:=DirichletGroup(Level(S),Rationals());
      i := 1;
      while true do
         if i ge #S`n2bar then
            break;
         end if;
         f := Twist(S`n2bar[i],eps);
         for j in [i+1..#S`n2bar] do
            if Trace(S`n2bar[j]) eq Trace(f) then
               Remove(~S`n2bar,j);
               break;
            end if;  
         end for;             
         i +:= 1;
      end while;
   end if;

   return S`n2bar;
end intrinsic;

intrinsic N_3(S::ModTupFld) -> SeqEnum
{}
   if not assigned S`n3 then
      a, S`n3 := CM_Newforms(S);
   end if;

   return S`n3;
end intrinsic;


function AllTwists(qexp, G)

   if #qexp eq 0 then
      return [* *];
   end if;
   R<qn>  := Parent(qexp[1]);

   ans := [* *];
   traces := [];
   for f in qexp do
      for i in [1..#G] do   
         eps := G.1^i;
         ftwist := Twist(f,eps);
         tracetwist := Trace(ftwist);
         if Index(traces, tracetwist) eq 0 then
            Append(~ans, ftwist);
            Append(~traces, tracetwist);
         end if;
      end for;
   end for;
   return ans;

end function;

intrinsic VG_0(S::ModTupFld) -> List
{}
   if not assigned S`vg0 then

      qexp := N_0(S);
      if #qexp gt 0 then
         S`vg0 := AllTwists(qexp, S`dirichlet);

         R<qp> := Parent(S`vg0[1]);
         for f in qexp do
            Append(~S`vg0,Evaluate(f,qp^Level(S))+O(qp^(S`prec)));
         end for;
      else
         S`vg0 := [* *];
      end if;

   end if;

   return S`vg0;
end intrinsic;

intrinsic VG_1(S::ModTupFld) -> List
{}
   if not assigned S`vg1 then

      qexp := N_1bar(S);
      if #qexp gt 0 then

         S`vg1 := AllTwists(qexp, S`dirichlet);

         p := Level(S);
         R<qp> := Parent(S`vg1[1]);
         for f in qexp do
            fqp := Evaluate(f,qp^p) + O(qp^S`prec);
            Append(~S`vg1,Evaluate(f,qp^p));
            Append(~S`vg1,ComplexConjugate(fqp));
         end for;
      else
         S`vg1 := [* *];
      end if;
      
   end if;

   return S`vg1;
end intrinsic;


intrinsic VG_2(S::ModTupFld) -> List
{}
   if not assigned S`vg2 then
      qexp := N_2bar(S);
      S`vg2 := AllTwists(qexp, S`dirichlet);
   end if;
 
   return S`vg2;      
end intrinsic;

intrinsic VG_3(S::ModTupFld) -> List
{}
   if not assigned S`vg3 then
      qexp := N_3(S);
      S`vg3 := AllTwists(qexp, S`dirichlet);
   end if;
 
   return S`vg3;      
end intrinsic;



intrinsic OldEigenformBasis(S::ModTupFld) -> SeqEnum
{Collection of q-expansions, which are eigenforms for the anemic Hecke algebra, 
 and whose Gal(Qbar/Q(zeta_(p-1)))-conjugates generate the old subspace.}

end intrinsic;

intrinsic OldspaceBasis(S::ModTupFld) -> SeqEnum
{Basis for the old space as a vector space over Q(zeta_(p-1)).}

end intrinsic;


intrinsic NewformBasis(S::ModTupFld) -> SeqEnum
{Collection of q-expansions, which are eigenforms for the anemic Hecke algebra, 
 and whose Gal(Qbar/Q(zeta_(p-1)))-conjugates generate the new subspace.}

end intrinsic;

intrinsic NewspaceBasis(S::ModTupFld) -> SeqEnum
{Basis for the new space as a vector space over Q(zeta_(p-1)).}

end intrinsic;


///////////////////////////////////////// Support functions ////////////////////////////

intrinsic Twist(f::RngSerPowElt, eps:GrpDrchElt) -> RngSerPowElt
{Sum a_n * eps(n) * q^n.}
   R<q> := Parent(f);
   prec := AbsolutePrecision(f);
   return &+[Coefficient(f,n)*q^n*Evaluate(eps,n) : n in [0..AbsolutePrecision(f)-1]] + O(q^prec);
end intrinsic;


intrinsic Trace(f::RngSerPowElt) -> RngSerPowElt
{Sum trace(a_n) * q^n.}
   R<q> := PowerSeriesRing(Parent(tracemod(Coefficient(f,1))));
   prec := AbsolutePrecision(f);
   return &+[tracemod(Coefficient(f,n))*q^n : n in [0..AbsolutePrecision(f)-1]] + O(q^prec);
end intrinsic;

intrinsic ComplexConjugate(f::RngSerPowElt) -> RngSerPowElt
{Sum trace(bar(a_n))*q^n. }
   R<q> := Parent(f);
   F    := BaseRing(R);
   prec := AbsolutePrecision(f);
 
   if Type(F) eq FldCyc then
      return &+[ComplexConjugate(Coefficient(f,n))*q^n : n in [0..AbsolutePrecision(f)-1]] + O(q^prec);
   end if;

   f    := Modulus(F);
   S<x> := Parent(f);
   fbar := &+[x^i*ComplexConjugate(Coefficient(f,i)) : i in [0..Degree(f)]];
   Fbar := quo<S|fbar>;
   Rbar<q> := PowerSeriesRing(Fbar);

   return &+[(Fbar!(S!Coefficient(f,n)))*q^n : n in [0..AbsolutePrecision(f)-1]] + O(q^prec);

end intrinsic;

///////////////////////////////////////// stuff that shouldn't be here /////////////////

 
intrinsic charpolymod(f::RngUPolElt, g::RngUPolElt) -> RngUPolElt
{Uses resultants to compute the characteristic polynomial of f(x)
 in K[x]/(g(x)), where f and g are polynomials in K[x], and
 g is assumed irreducible.}
/* This algorithm is from page 162 of Cohen's "A course in computational
   algebraic number theory". */
   if Parent(f) ne Parent(g) then
      error "charpolymod: f and g must have the same parent.";
   end if;
   R<x> := Parent(f);
   K := Parent(Coefficient(f,0));
   S<z,y> := PolynomialRing(K,2);
   ff := &+[S|y^i * Coefficient(f,i) : i in [0..Degree(f)]];
   g /:= Coefficient(g,Degree(g));       // make g monic
   gg := &+[S|y^i * Coefficient(g,i) : i in [0..Degree(g)]];
   c  := Resultant(gg,z-ff,2);
   return &+[R|x^i * (K!LeadingCoefficient(Coefficient(c,1,i)))
                    : i in [0..Degree(c,1)]];
end intrinsic;
 
intrinsic charpolymod(f::RngUPolResElt) -> RngUPolElt
{Compute the characteristic polynomial of f in K[x]/(g(x)).  The built in
Magma algorithm for doing this doesn't work.}
   R := Parent(f);
   g := Modulus(R);
   K := Parent(Coefficient(g,0));
   liftf := PolynomialRing(K)!f;
   return charpolymod(liftf, g);
end intrinsic;
 
intrinsic charpolymod(f::FldPadElt) -> RngUPolElt
{Compute the characteristic polynomial of f in K[x]/(g(x)).  The built in
Magma algorithm for doing this doesn't work.}
   R<x> := PolynomialRing(Parent(f));
   return x-f;
end intrinsic;      

intrinsic tracemod(f::RngUPolResElt) -> FldCycElt
{Trace down to K of f in K[x]/(g(x)).}
   h := charpolymod(f);
   return -Coefficient(h,Degree(h)-1);
end intrinsic;

intrinsic tracemod(x::FldCycElt) -> FldCycElt
{Just x again.}
   return x;
end intrinsic;