// David Kohel's code for decomposing a space as a direct
// sum of Hecke stable subspaces. 

intrinsic HeckeDecomposition(TM::AlgMatElt) -> SeqEnum
   {Returns the sequence of vector subspaces stabilized by 
   the Hecke operators.}
   T := BrandtHecke(TM,2);
   return HeckeDecomposition(TM,BaseModule(Parent(T)));
end intrinsic;

intrinsic HeckeDecomposition(TM::AlgMatElt,V::ModTupFld) -> SeqEnum
   {Returns the sequence of vector subspaces of V stabilized 
   by the Hecke operators.}
   p := 2;
   Decomp := [ V ];
   if Dimension(Decomp[1]) eq 1 then 
      return Decomp;
   end if;
   Irred := [ false ];
   while not &and Irred do
      T := BrandtHecke(TM,p);
      for V in [ V : V in Decomp | not Irred[Index(Decomp,V)] ] do
         fact := Factorization(CharacteristicPolynomial(T,V));
         if #fact gt 1 then
            Remove(~Irred,Index(Decomp,V));
            Exclude(~Decomp,V);
            for f in fact do
               Append(~Decomp,V meet Kernel(Evaluate(f[1],T)));
               if f[2] eq 1 then
                  Append(~Irred,true);
               else
                  Append(~Irred,false);
               end if;
            end for;
         elif fact[1][2] eq 1 then
            Irred[Index(Decomp,V)] := true;
         end if;
      end for;
      p := NextPrime(p);
   end while;
   return Decomp;
end intrinsic;

intrinsic HeckeLinearDecomposition(TM::AlgMatElt) -> SeqEnum
   {Returns the sequence of vector subspaces stabilized by 
   the Hecke operators.}
   T := BrandtHecke(TM,2);
   return HeckeLinearDecomposition(TM,BaseModule(Parent(T)));
end intrinsic;

intrinsic HeckeLinearDecomposition(TM::AlgMatElt,V::ModTupFld,
   N::RngIntElt) -> SeqEnum
   {The sequence of one dimensional vector subspaces stabilized by 
   the Hecke operators, priority given to the divisors of the level N.}
   p := 2;
   Decomp := [ V ];
   for q in PrimeDivisors(N) do
      T := BrandtHecke(TM,q);
      for V in Decomp do
         if Dimension(V) gt 1 then
            eigen := Eigenvalues(T,V);
            Exclude(~Decomp,V);
            for a in eigen do
               Append(~Decomp,V meet Kernel(T-a[1]));
            end for;
         end if;
      end for;
   end for;
   while &or[ Dimension(V) ne 1 : V in Decomp ] do
      T := BrandtHecke(TM,p);
      for V in Decomp do
         if Dimension(V) gt 1 then
            eigen := Eigenvalues(T,V);
            Exclude(~Decomp,V);
            for a in eigen do
               Append(~Decomp,V meet Kernel(T-a[1]));
            end for;
         end if;
      end for;
      p := NextPrime(p);
   end while;
   for i in [2..#Decomp] do
      s := &+[ Basis(Decomp[i])[1][j] : j in [1..Degree(Parent(T))] ];
      if s ne 0 then
         v := Decomp[1];
         Decomp[1] := Decomp[i];
         Decomp[i] := v;
      end if;
   end for;
   return Decomp;
end intrinsic;

intrinsic HeckeLinearDecomposition(TM::AlgMatElt,V::ModTupFld) -> SeqEnum
   {The sequence of linear Hecke subspaces.}
   p := 2;
   Decomp := [ V ];
   while &*[ Integers() | Dimension(V) : V in Decomp ] gt 1 do
      T := BrandtHecke(TM,p);
      for V in Decomp do
         if Dimension(V) gt 1 then
            eigen := Eigenvalues(T,V);
            Exclude(~Decomp,V);
            for a in eigen do
                // Should be:
                //   Append(~Decomp,V meet Kernel(T-a[1]));
                Decomp := Decomp cat [ V meet Kernel(T-a[1]) ];
            end for;
         end if;
      end for;
      p := NextPrime(p);
   end while;
   for i in [2..#Decomp] do
      s := &+[ Basis(Decomp[i])[1][j] : j in [1..Degree(Parent(T))] ];
      if s ne 0 then
         v := Decomp[1];
         Decomp[1] := Decomp[i];
         Decomp[i] := v;
      end if;
   end for;
   return Decomp;
end intrinsic;

intrinsic HeckeThetaDecomposition(TM::AlgMatElt) -> SeqEnum
   {}
   Decomp := HeckeDecomposition(TM);
   VM := BaseModule(Parent(TM));
   ThetaDecomp := [ ];
   for V in Decomp do
      BM := [ VM!v : v in Basis(V) ];
      Append(~ThetaDecomp,EchelonSeries(
         [ InnerProduct(v,v*TM) : v in BM ] ));
   end for;
   return ThetaDecomp;
end intrinsic;

intrinsic HeckeEigenvectors(TM::AlgMatElt) -> SeqEnum
   {Returns the sequence of generators for the linear Hecke-invariant
   subspaces of the Hecke module.}
   n := Degree(Parent(TM));
   U := BaseModule(Parent(BrandtHecke(TM,2)));
   BQ := [ U!(Basis(V)[1]) : V in HeckeLinearDecomposition(TM) ];
   return [ v*LCM([ Denominator(v[i]) : i in [1..n] ]) : v in BQ ];
end intrinsic;

intrinsic HeckeEigenvalues(TM::AlgMatElt,v::ModTupFldElt,n:RngIntElt) -> SeqEnum
   {Returns the eigenvalues of Hecke for the first n primes.}
   S := [ RationalField() | ];
   p := 1;
   for i in [1..n] do
      p := NextPrime(p);
      T := BrandtHecke(TM,p);
      Append(~S,InnerProduct(v,v*T)/InnerProduct(v,v));
   end for;
   return S;
end intrinsic;

intrinsic CremonaReference(E::CurveEll) -> MonStgElt
   {The Cremona database reference for E.}
   N := Conductor(E);  
   if N lt 5000 then 
      SQ := ["A","B","C","D","E","F","G","H","I","J","K","L","M","N"];
      DB := EllipticCurveDatabase(); 
      E := MinimalModel(E);
      t := NumberOfCurves(DB,N); i := 1;
      while t gt 0 do
         if E in EllipticCurves(DB,N,i) then
            return IntegerToString(N) cat SQ[i];
         end if;
         t -:= NumberOfCurves(DB,N,i); i +:= 1;
      end while;
   end if;
   return "Not in tables."; 
end intrinsic;

intrinsic CremonaReferences(TM::AlgMatElt,N::RngIntElt) -> SeqEnum
   {The sequence of elliptic curves corresponding to linear 
   Hecke subspaces.}
   DB := EllipticCurveDatabase();
   SQ := HeckeLinearDecomposition(TM);
   p := 2; while GCD(p,N) eq p do p := NextPrime(p); end while;
   Tp := BrandtHecke(TM,p);
   return [ CremonaReference(TM,V,N) : V in SQ | Trace(Tp,V) ne p+1 ];
end intrinsic;

intrinsic CremonaReference(TM::AlgMatElt,V::ModTupFld,N::RngIntElt) 
   -> MonStgElt
   {The Cremona database reference for linear Hecke eigenspace V of level N.}
   DB := EllipticCurveDatabase(); 
   SQ := ["A","B","C","D","E","F","G","H","I","J","K","L","M","N"];
   TstPrms := [ p : p in InitialPrimes(25) | p notin PrimeDivisors(6*N) ];
   for M in [ M : M in Divisors(N) | (N div M) gt 10 ] do
      N0 := N div M;
      CremRef := IntegerToString(N0);
      t := NumberOfCurves(DB,N0);
      s := 0; // The number of isogeny classes.
      while t gt 0 do
        s +:= 1;
        t -:= NumberOfCurves(DB,N0,s); 
      end while;
      EllCurvs := [ EllipticCurve(DB,N0,i,1) : i in [1..s] ];
      TstInds := [1..#EllCurvs];
      j := 1;
      while (#TstInds gt 0) and (#TstPrms ge j) do
         p := TstPrms[j]; 
         Tp := BrandtHecke(TM,p);
         for k in TstInds do
            if Trace(Tp,V) ne Trace(EllCurvs[k],p) then
               Exclude(~TstInds,k);
            end if;
         end for;
         j +:= 1;
      end while;
      if #TstInds eq 1 then 
         if N0 eq N then 
            age := "New";
         else age := "Old";
         end if;
         return IntegerToString(N0) cat SQ[TstInds[1]], age;
      elif #TstInds gt 1 then 
         return "Inconclusive.";
      end if;
   end for;
   return "???";
end intrinsic;
