forward VectorSpaceElement,
        ToLowerCaseLetter, 
        PrecisionAtWeight,
        qRelationsOfDegree,
        idxG0;

intrinsic HasCanonicalEmbeddingX0(N::RngIntElt) -> BoolElt
{True if and only if X_0(N) has genus > 2 and the canonical map associated to
X_0(N) is an embedding.}

   return N ge 11 and DimensionCuspFormsGamma0(N,2) gt 2 and  
               not (N in {30,33,35,37,39,40,41,46,47,48,59,71});

end intrinsic;


intrinsic CanonicalEmbeddingX0(N::RngIntElt) -> SchProj
{The canonical embedding of X_0(N).}

   requirege N,1;
   require HasCanonicalEmbeddingX0(N) : "The canonical map is not an embedding.";

   n := DimensionCuspFormsGamma0(N,2);
   K := RationalField();
   M := CuspidalSubspace(ModularSymbols(N,2,K,+1));
   R := PolynomialRing(K,n);   
   AssignNames(~R,[ToLowerCaseLetter(i) : i in [1..n]]);
   Pn := ProjectiveSpace(R);
   deg := 0;
   C := Scheme(Pn,R!0);
  
   while true do
      deg +:= 1;
      rels := qRelationsOfDegree(M,PrecisionAtWeight(N,2*deg),deg);
      if #rels gt 0 then 
         rels := [R!r : r in rels];
         C := C meet Scheme(Pn, rels);
         if Dimension(C) eq 1 and IsIrreducible(C) then 
//         if Dimension(C) eq 1 then 
            DisownChildren(AmbientSpace(M));
            return C;
         end if;
      end if;
   end while;

end intrinsic;



///////////////////////////////////////////////////////////
//  The functions used in defining the above intrinsic.  //
///////////////////////////////////////////////////////////

function qRelationsOfDegree(M, prec, d) 
// The algebraic relations of degree d satisfied by the sequence of
// q-expansions qExpansions(M,prec) of a basis for the space of
// modular forms corresponding to M.

   assert Type(d) eq RngIntElt;
   assert Type(prec) eq RngIntElt;
   assert Type(M) eq ModSym; 

   assert IsCuspidal(M) eq true; 
   assert Dimension(M) gt 0 ;
   assert d ge 1;

   K := BaseField(M);
   B := qExpansionBasis(M,prec);
   S := Parent(B[1]);
   R := PolynomialRing(K,#B);
   AssignNames(~R,[ToLowerCaseLetter(n) : n in [1..#B]]);
   monoms := MonomialsOfDegree(R,d);
   forms  := [&*[B[i]^exp[i] : i in [1..#B]] where exp := Exponents(m) : m in monoms];
   V    := VectorSpace(K,prec);
   vecs := [VectorSpaceElement(V,f) : f in forms];
   A := RMatrixSpace(K,#forms,prec)!vecs;
   relvecs := Basis(Kernel(A));
   rels := [&+[v[i]*monoms[i] : i in [1..#monoms]] : v in relvecs];

   return rels;
end function;

alphabet := {@ "a", "b", "c", "d", "e", "f", "g", "h", "i", "j", "k", "l",
     "m", "n", "o", "p", "q", "r", "s", "t", "u", "v", "w", "x", "y", "z" @};

function ToLowerCaseLetter(n)
// Returns the n-th letter coding.  The coding goes a,b,c,...,z,aa,bb,...
   assert Type(n) eq RngIntElt and n gt 0;
   return alphabet[((n-1)mod 26)+1]^(((n-1) div 26) + 1);
end function;

function idxG0(n)
   assert Type(n) eq RngIntElt;
   return 
      &*[Integers()| t[1]^t[2] + t[1]^(t[2]-1) : t in Factorization(n)];
end function;

function PrecisionAtWeight(N,k)
   assert Type(N) eq RngIntElt and Type(k) eq RngIntElt;
   return Floor(idxG0(N) * k/12 + 1) + 5;  // the "+5" should be unnecessary.
end function;

function VectorSpaceElement(V, f)
   R<q>:=Parent(f);
   S := [0 : i in [1..Valuation(f)]] cat Eltseq(f+O(q^Degree(V)));
   S := S cat [0 : i in [1..Degree(V)-#S]];
   return V!S;
end function;

