/* eigenrep.m -- Decompose a modular form in terms of eigenforms. */

forward ApproxCoordinates;
complex_prec := 500;

intrinsic EigenDecomp(M::ModFrm, f::ModFrmElt, prec::RngIntElt) -> SeqEnum, SeqEnum
{Decompose a modular form in terms of eigenforms.} 
   if prec lt Dimension(M) + 10 then
      print "WARNING: Precision is probably too low.";
   end if;
   return EigenDecomp(M, PowerSeries(f,prec), prec);
end intrinsic;

intrinsic EigenDecomp(M::ModFrm, f::RngSerElt, prec::RngIntElt) -> SeqEnum, SeqEnum
{}
   // Compute a basis of power series q-expansions of eigenforms for M. 
   E := Eigenbasis(M,prec);
   // Write f in terms of this basis
   V := VectorSpace(ComplexField(complex_prec),prec);
   w := ApproxCoordinates([V | [Coefficient(g,n) : n in [0..prec-1]] : g in E], 
          V![Coefficient(f,n) : n in [0..prec-1]], M);
   return w, E;
end intrinsic;

function Frob(f, n)
   // return f(q^n).
   R := Parent(f);
   return Evaluate(f,R.1^n) + O(R.1^AbsolutePrecision(f));
end function;

function NewformsImage(M, prec, t) 
   // Compute basis for newforms of level M*d coming from the newforms in M.
   // that is, compute a list of the newforms in M and take their
   // images under f(q) |-> f(q^d) for each divisor d of t.
   ans := [];
   for i in [1..NumberOfNewformClasses(M)] do
      f := Newform(M,i);
      for g in ComplexEmbeddings(f)[1] do
         h := PowerSeries(g,prec);
         for d in Divisors(t) do
            Append(~ans, Frob(h,d));
         end for;
      end for;
   end for;
   return ans;
end function;

intrinsic Eigenbasis(M::ModFrm, prec::RngIntElt) -> List
{}
   ans := [];
   char := DirichletCharacters(M);
   require #char eq 1 : "DirichletCharacters(argument 1) must have size 1.";
   cond := Conductor(char[1]);
   for N in [x : x in Divisors(Level(M))| x mod cond eq 0] do
      M_N := ModularForms(M, N);
      if IsCuspidal(M) then 
         M_N := CuspidalSubspace(M_N);
      end if;
      ans := ans cat NewformsImage(M_N, prec, Level(M) div N);
   end for;
   return ans;
end intrinsic;


function Pivots(B)   // find pivots of reduced basis.
   return [Min(Support(b)): b in B];
end function;

function Valuation(f)
   // temp := PowerSeries(f,Dimension(Parent(f))+10);
   n := 0;
   while Coefficient(f,n) eq 0 do
      n := n + 1;
   end while;
   return n;
end function;

function ApproxCoordinates(B, x, M)
   // B is a sequence of elements of a vector space
   // x is an element that is approximately in that vector space.

   V := Parent(B[1]);

   // 1. find pivot entries of B,

   // R := Basis(sub<V|B>);
   mat := RMatrixSpace(BaseRing(V),#B,Degree(V))!(&cat[Eltseq(v) : v in B]);
   ech := EchelonForm(mat);
   P := Pivots([ech[i] : i in [1..#B]]);

   // P := [Valuation(f)+1 : f in Basis(M)];
   // print "Pivots =", P;

   // 2. make matrix whose columns are just those entries of B,

   n := #B;
   A := MatrixAlgebra(ComplexField(complex_prec),n)!(&cat[[b[i] : i in P] : b in B]);
   A := Transpose(A);

   // 3. invert that matrix,
   print "A = ", A;
//   print "Determinant(A) = ", Determinant(A);
"now compute inverse!";
   C := A^(-1);
      
   // 4. apply inverted matrix to pivot subvector of x (viewed as a column vector).
   y := RMatrixSpace(ComplexField(complex_prec),n,1)![x[i] : i in P];

   return Eltseq(C*y);
end function;



/* 
Example:

> M := ModularForms(DirichletGroup(37).1,2);
> w, E := EigenDecomp(M, M.1+M.2+M.3, 10);
> &+[w[i]*E[i] : i in [1..4]];
1 + (0.99999999999999999999999999999999999999 +  0.E-39*i)*q + (1.0000000000000000000000000000000000000 + 8.8162076286971569060000000000000000000 E-39*i)*q^2 + (2.35098870098590850800000000000000000000000000000 E-39 +  0.E-39*i)*q^3 + (-1.0000000000000000000000000000000000000 +  0.E-39*i)*q^4 + (1.9999999999999999999999999999999999999 - 8.8162076286971569060000000000000000000 E-39*i)*q^5 + (1.9999999999999999999999999999999999999 - 8.8162076286971569060000000000000000000 E-39*i)*q^6 + (3.99999999999999999999999999999999999998883280366 +  0.E-39*i)*q^7 + (3.00000000000000000000000000000000000000000000000000000000 +  0.E-48*i)*q^8 + (1.00000000000000000000000000000000000000470197740 +  0.E-39*i)*q^9 + O(q^10)
> PowerSeries(M.1+M.2+M.3,10)-$1;
(5.8774717541114375398032809710000000000 E-39 +  0.E-39*i)*q + ( 0.E-38 - 8.8162076286971569060000000000000000000 E-39*i)*q^2 + (-2.35098870098590850800000000000000000000000000000 E-39 +  0.E-39*i)*q^3 + (4.7019774034944183618649214802461880450 E-39 + 8.8162076286971569060000000000000000000 E-39*i)*q^5 + (4.7019774034944183618649214802461880450 E-39 + 8.8162076286971569060000000000000000000 E-39*i)*q^6 + (1.1167196335001260175737169960000000000 E-38 +  0.E-39*i)*q^7 + (3.45446742139756679500000000000000000000000000000000000000 E-77 +  0.E-48*i)*q^8 + (-4.7019774019718170160000000000000000000 E-39 +  0.E-39*i)*q^9 + O(q^10)

*/
