// prop20.m in Merel's paper


intrinsic MyHeilbronn(n::RngIntElt) -> SeqEnum
{}

   ans0 := []; 
   A := MatrixAlgebra(Integers(),2);
   for a in [0..n] do
      for b in [0..n] do
         for c in [0..n] do
            for d in [0..n] do   
               if a gt b and d gt c and a*d-b*c eq n then
                  Append(~ans0, A![a,b,c,d]);
               end if;
            end for;
         end for;
      end for;
   end for;
   return ans0;
end intrinsic;

intrinsic MyHeilbronn(n::RngIntElt, ell::RngIntElt) -> SeqEnum
{}

   ans := []; 
   A := MatrixAlgebra(GF(ell),2);
   for a in [0..ell-1] do
      for b in [0..ell-1] do
         for c in [0..ell-1] do
            for d in [0..ell-1] do   
               if a gt b and d gt c and (a*d-b*c - n) mod ell eq 0 then
                  Append(~ans, A![a,b,c,d]);
               end if;
            end for;
         end for;
      end for;
   end for;
   return ans;
end intrinsic;

intrinsic Stat(HB::SeqEnum, ell::RngIntElt) -> SeqEnum
{}
   A := MatrixAlgebra(GF(ell),2);
   ans_ell := [];
   for a in A do
      cnt := #[i : i in [1..#HB] | HB[i] eq a];
      if cnt gt 0 then
         Append(~ans_ell, <Eltseq(a),cnt>);
      end if;
   end for;
   return Sort(ans_ell);   
end intrinsic;

intrinsic PrimeSeq(p1::RngIntElt,p2::RngIntElt)  -> SeqEnum
{The sequence of primes [p1 ... p2].}
   return [p : p in [p1..p2] | IsPrime(p)];
end intrinsic;

intrinsic MyCremona(p::RngIntElt) -> SeqEnum
{Return the Heilbronn matrices of determinant p, as defined by Cremona.}
   assert IsPrime(p); //: "Argument 1 must be prime.";
   A := MatrixAlgebra(Integers(),2);
   if p eq 2 then
      return [A|[1,0, 0,2], [2,0, 0,1], [2,1, 0,1], [1,0, 1,2]];
   end if;
   ans := [[1,0, 0,p]];
   for r in [Ceiling(-p/2)..Floor(p/2)] do
      x1:=p; x2:=-r; y1:=0; y2:=1; a:=-p; b:=r; c:=0; x3:=0; y3:=0; q:=0;
      Append(~ans, [x1,x2, y1,y2]);
      while b ne 0 do
         q := Round(a/b);
         c := a-b*q;
         a := -b;
         b := c;
         x3 := q*x2-x1; 
         x1 := x2; x2 := x3; y3 := q*y2-y1; y1 := y2; y2 := y3;
         Append(~ans, [x1,x2, y1, y2]);
      end while;
   end for;
   return [A!a : a in ans];
end intrinsic;




////////////////////////////////////

function MatrixAction(A, r)
   if r[2] ne 0 then
      z := r[1]/r[2];
      if A[2,1]*z + A[2,2] eq 0 then
         return [1,0];
      end if;
      s := (A[1,1]*z+A[1,2])/(A[2,1]*z + A[2,2]);
      return [Numerator(s), Denominator(s)];
   end if;
   return A[1,2] / A[2,2];
end function;

intrinsic Act(H::AlgMatElt, x::ModSymElt) -> ModSymElt
{}
   return Act([<1,H>], x);
end intrinsic;

intrinsic Act(H::[Tup], x::ModSymElt) -> ModSymElt
{H is is a sequence of pairs <n, A>, where n is an integer
and A is a 2x2 rational matrix.}
   M := Parent(x);
   require Weight(M) eq 2 : "Argument 2 must have weight 2.";
   ans := M!0;
   m := ModularSymbol(x);
   for s in H do
      A := s[2]; 
      for z in m do
         if s[1]*z[1] ne 0 then
            ans := ans + M!<s[1]*z[1], 
                  [MatrixAction(A,a): a in z[2]]>;
         end if;
      end for;
   end for;
   return ans;
end intrinsic;

intrinsic Hecke1(p::RngIntElt, x::ModSymElt) -> ModSymElt
{Compute Hecke using the usual matrices.}
   A := MatrixAlgebra(Rationals(),2);
   H := [<1,A![1,r,0,p]> : r in [0..p-1]] cat
        [<1,A![p,0,0,1]>];
return H;
end intrinsic;

intrinsic Hecke2(p::RngIntElt, x::ModSymElt) -> ModSymElt
{Compute Hecke using the usual matrices.}
   ell := Level(Parent(x));
   require IsPrime(ell) : "The level of argument 2 must be prime.";
   A := MatrixAlgebra(Rationals(),2);
   H := [<Floor((p-1)/ell)+1, A![1,0,0,p]>] cat
        [<Floor((p-1-a)/ell)+1, A![1,a,0,p]> : a in [1..Min(ell-1,p-1)]] cat
        [<1,A![p,0,0,1]>];
   return Act(H,x);
end intrinsic;

intrinsic Hecke3(p::RngIntElt, x::ModSymElt) -> ModSymElt
{}
   ell := Level(Parent(x));
   require IsPrime(ell) : "The level of argument 2 must be prime.";

   A := MatrixAlgebra(Rationals(),2);
   H := [<1,A![1,p*Integers()!(Integers(ell)!(r/p)),0,p]> : r in [0..p-1]] cat
        [<1,A![p,0,0,1]>];
   // Now massage:
   /* Replace [1,r,0,p] by [1,p*alpha,0,p], where
      alpha = r/p mod ell. 
   */
   return Act(H,x);
end intrinsic;

