
Z:=Integers();
Zx<x>:=PolynomialRing(Z);

function PrimeIndex(n)
  i := 1;
  p := 2;
  while(p le n) do
    p := NextPrime(p);
    i := i + 1;
  end while;
  return i - 1;
end function;

function PNHeckePoly(p,n,k)
  if (n eq 0) then
    return 1;
  end if;
  if (n eq 1) then
    return Zx.1;
  else
    return Zx.1 * PNHeckePoly(p,n-1,k) - p^(k-1) * PNHeckePoly(p,n-2,k);
  end if;
end function; 

intrinsic HeckePoly(n::RngIntElt,T::RngMPol,k::RngIntElt) -> RngMPolElt
{Given n and a multivariate polynomial ring T, writes the nth Hecke
operator T_n as a polynomial in the generators of T, where T.1 = T_2,
T.2 = T_3, T.3 = T_5, etc.}
  P := 1;
  F:=Factorization(n);
  for term in F do
    P1 := PNHeckePoly(term[1],term[2],k);
    P := P * MultivariatePolynomial(T,P1,PrimeIndex(term[1])); 
  end for;
  return P;
end intrinsic;

function HeckeMap(T,Tprime,k)
  L := [];
  for i := 1 to Rank(T) do
    L[i] := HeckePoly(i+1,Tprime,k);
  end for;
  return hom<T->Tprime|L>;
end function;

intrinsic IsGorenstein(I::RngMPol) -> BoolElt
{Given a zero-dimensional, primary, ideal I of a polynomial ring R, 
tests whether the local ring R/I is Gorenstein.}
  require IsZeroDimensional(I) : "Argument must be zero-dimensional";
  require IsPrimary(I) : "Argument must be primary";
  m:=Radical(I);
  J:=ColonIdeal(I,m);
  Jin:=LeadingMonomialIdeal(J);
  Iin:=LeadingMonomialIdeal(I);
  count:=0;
  for term in Basis(Jin) do
    if not (term in Iin) then
      count:=count + 1;
    end if;
  end for;
  assert (count gt 0);
  if (count eq 1) then
    return true;
  end if;
  return false;
end intrinsic; 

function MapIdeal(I,g)
  Iprime := ideal<Codomain(g)|0>;
  for m in Basis(I) do
    Iprime := Iprime + ideal<Codomain(g)|g(m)>;
  end for;
  return Iprime;
end function;
 
function LinearForm(T, S);
  p := S[1];
  for i := 2 to #S do
    p := p + S[i] * T.(i-1);
  end for;
  return p;
end function;

intrinsic HeckeIdeal(M::ModSym, Bound::RngIntElt) -> RngMPol,RngMPol
{Given a modular symbols space M and an integer Bound, constructs the
Z-algebra generated by the Hecke operators T_n with p less than Bound,
and returns the kernel of the map from the free Z-algebra on the T_n
to this Z-algebra.}
  require IsTrivial(DirichletCharacter(M)) :
  "Space must have trivial character.";
  require IsCuspidal(M) : "Space must be cuspidal.";
  d:=Dimension(M);
  Matrices := [];
  for n := 1 to Bound do 
    Tn := IntegralHeckeOperator(M,n);
    if (n eq 1) then
      Mat := Matrix(Integers(),1,d*d,ElementToSequence(Tn));
    else 
      M2 := Matrix(Integers(),1,d*d,ElementToSequence(Tn));
      Mat := VerticalJoin(Mat,M2);
    end if;
    Matrices[n] := Tn;
  end for;
  T:=PolynomialRing(Integers(),Bound-1);
  Tprime:=PolynomialRing(Integers(),PrimeIndex(Bound));
  K:=NullspaceMatrix(Mat); 
  assert(Bound eq NumberOfColumns(K));
  I:=ideal<T|0>;
  // Add linear relations to ideal;
  for i := 1 to NumberOfRows(K) do
    f := LinearForm(T, ElementToSequence(K[i]));
    I := I + ideal<T|f>;
  end for;
  // Construct quadratic relations;
  for i := 2 to Bound do
    for j := i to Bound do
      P := T.(i-1) *T.(j-1);
      V := Vector(d*d,ElementToSequence(Matrices[i]*Matrices[j]));
      f := LinearForm(T, ElementToSequence(Solution(Mat,V)));
      P := P - f;
      I := I + ideal<T|P>;
    end for;
  end for;
  Iprime := MapIdeal(I,HeckeMap(T,Tprime,Weight(M)));
  return Tprime,Iprime;
end intrinsic;
 
intrinsic HeckeIdealDecomp(M::ModSym,Bound:RngIntElt) -> RngMPol,RngMPol,
List
{Given a modular symbols space M, and a Bound, decomposes M into
subspace corresponding to its newforms, and returns, in order, a ring
T = Z[T_2,T_3, ..., T_p] where p is the largest prime less than or
equal to Bound, an ideal I equal to the annhilator of M in T,
the ideal corresponding to the annhilator of the Eisenstein subspace of M,
and a list of ideals which are the annhilators of the spaces appearing
in NewformDecomposition(CuspidalSubspace(M)), in sorted order.}
  require IsTrivial(DirichletCharacter(M)) : 
  "Space must have trivial character.";
  require IsCuspidal(M) : "Space must be cuspidal.";
  Forms:=SortDecomposition(NewformDecomposition(CuspidalSubspace(M))); 
  L := [* *];
  T := PolynomialRing(Integers(),PrimeIndex(Bound));
  I := ideal<T|1>;;
  for f in Forms do
    Tf,If := HeckeIdeal(f,Bound);
    h := Coercion(Tf,T);
    If := MapIdeal(If,h); 
    L := Append(L,If);
    I := I meet If;
  end for;
  return T,I,L;
end intrinsic; 
