declare verbose Models, 1;

intrinsic FindEquation(f::RngSerLaurElt, g::RngSerLaurElt) -> RngUPolElt, RngSerLaurElt
{}

   vprintf Models : "f = %o.\ng = %o.\n", f,g;
   q := Parent(f).1;
   X := f/g;
   Y := q*Derivative(X)/g;  

   R<x>:=PolynomialRing(Rationals());
 
   h := x^6;
   s := X^6;
   Y2 := Y^2;
   for n in [5,4,3,2,1,0] do
      a := Coefficient(Y2-s,-n);
      h +:= a*x^n;  
      s +:= a*X^n;
   end for;

   return h, Y2-s;
end intrinsic;

intrinsic FindModel(A::ModSym) -> RngUPolElt, RngSerLaurElt
{}
   prec := 37;

   if Sign(A) ne 0 and Dimension(A) eq 2 or Sign(A) eq 0 and Dimension(A) eq 4 then
      R<q>:=LaurentSeriesRing(Rationals());
      E := qExpansionBasis(A,prec+1);
      t,f := IsCoercible(R,E[1]);
      if not t then 
         return PolynomialRing(Rationals())!0, R!0;
      end if;
      g := R!E[2];
      return FindEquation(f,g);
   end if;

   F := qEigenform(A,prec+1);
   coef := [Eltseq(Coefficient(F,n)): n in  [1..prec]];
   n := #coef[1];
   nz := [];
   for c in coef do
      for j in [1..#c] do
         if c[j] ne 0 then
            Append(~nz,j);
         end if;
      end for;
      if #nz ge 2 then
         break;
      end if;
   end for;
   nz := Setseq(Set(nz));


   R<q>:=LaurentSeriesRing(Rationals());
   if #nz ne 2 or Type(coef[n][nz[1]]) ne FldRatElt then 
      return PolynomialRing(Rationals())!0, R!0;
   end if;

   f := &+[coef[n][nz[1]]*q^n : n in [1..prec]];
   g := &+[coef[n][nz[2]]*q^n : n in [1..prec]];
   if Coefficient(f,1) eq 0 then
      t := g;
      g := f;
      f := g;
   end if;
  
   if Coefficient(g,2) eq 0 then
      return PolynomialRing(Rationals())!0, R!0;
   end if;
 
   g /:= R!Coefficient(g,2);

   return FindEquation(f,g);
end intrinsic;

intrinsic Qchi(N::RngIntElt) -> FldCyc
{The field Q(zeta_n), where n = exponent of (Z/NZ)^*.}
   return CyclotomicField(Exponent(UnitGroup(IntegerRing(N))));
end intrinsic;

intrinsic FullDirichlet(N::RngIntElt) -> GrpDrch
{}
//   G := DirichletGroup(N,Qchi(N));
   G := DirichletGroup(N,CyclotomicField(12));
   names := ["e" cat IntegerToString(i) : i in [1..Ngens(G)]];
   AssignNames(~G,names);
   return G;
end intrinsic;

intrinsic Analyze(N::RngIntElt) -> SeqEnum, SeqEnum
{}
   G := FullDirichlet(N);
   E := Elements(G);
   
   S := [<DimensionCuspForms(e,2),e,e^(-1)> : e in Elements(G) | Order(e) in [1,2,3,4,6] and not IsTrivial(e)];

   M := &cat[Decomposition(NewSubspace(CuspidalSubspace(ModularSymbols(s[2],2,+1))),23) : s in S | s[1] gt 0];
   return S, [A : A in M | Dimension(A) eq 1 or Order(DirichletCharacter(A)) eq 2 and Dimension(A) eq 2];
end intrinsic;

intrinsic DescribeLevel(N::RngIntElt) -> SeqEnum
{}

 printf " \n\n>>>>>>>>>> X_1(%o) <<<<<<<<\n\n",N;
 S,M:=Analyze(N);    
 S; 
 for i in [1..#M] do
    F,dif := FindModel(M[i]); 
    if F ne 0 then
       printf "\n\n S_2[I]=%o\n", qExpansionBasis(M[i],8); 
       F;
       dif;
    end if;
 end for;
 return M;
end intrinsic;
