// This is an algorithm wrote by Haluk Sengun designed to compute abelianizations
// of congruence subgroups Gamma_0(I) of GL_2(O_K), where I is an ideal.
// In this specific example, we produce a function that computes abelianizations
//of congruence subgroups corresponding to ideals that are powers of (sqrt(-2))
// the only prime that lies above 2 in Q(sqrt(-2)).
// It can be easily adapted to work for any quadratic imaginary of class number 1;

Gamma0_abel_1:=function(power)
// We use generators and relations for GL_2(O_K) computed by Swan in
// Richard G Swan, Generators and relations for certain special linear groups,
// In Advances in Mathematics, Volume 6, Issue 1, 1971, Pages 1-77,
// ISSN 0001-8708, https://doi.org/10.1016/0001-8708(71)90027-2.
// (http://www.sciencedirect.com/science/article/pii/0001870871900272)
// This paper contains generators and relations for GL_2(O_K) of every quad.
// imag of trivial class class group.
      A<j,t,u,a,e>:=FreeGroup(5);
      BG<j,t,u,a,e>:=quo<A|j^2,t*u=u*t, a^2=j,(t*a)^3=j,(a*u^(-1)*a*u)^2=j, a*j=j*a, j*t=t*j,j*u=u*j,j,e*t*e=t^(-1), e*u*e=u^(-1), e*j=j*e, e*a*e = j*a>;

      K:=QuadraticField(-2);
      O<z>:=MaximalOrder(K);
// "I" is just a power of two here, but the code
// works of any ideal of O_K;     
      I := (Factorization(2*O)[1,1])^(power);

      J:=Matrix(O,2,2,[-1,0,0,-1]);
      T:=Matrix(O,2,2,[1,1,0,1]);
      U:=Matrix(O,2,2,[1,z,0,1]);
      A := Matrix(O,2,2,[0,-1,1,0]);
      E := Matrix(O,2,2,[-1,0,0,1]);

      PL,r:=ProjectiveLine(quo<O|I>);

      proj_action:=function(M,i)
        t,im :=r(PL[i]*M,true,false);
        return Index(PL,im);
      end function;

      SeqE:=[];
      for i in [1..#PL] do
        Append(~SeqE,proj_action(E,i));
      end for;

      SeqT:=[];
      for i in [1..#PL] do
        Append(~SeqT,proj_action(T,i));
      end for;

      SeqJ:=[];
      for i in [1..#PL] do
        Append(~SeqJ,proj_action(J,i));
      end for;

      SeqU:=[];
      for i in [1..#PL] do
        Append(~SeqU,proj_action(U,i));
      end for;

      SeqA:=[];
      for i in [1..#PL] do
        Append(~SeqA,proj_action(A,i));
      end for;

      Symm:=Sym(#PL);
      pJ:=Symm!SeqJ;
      pT:=Symm!SeqT;
      pU:=Symm!SeqU;
      pA := Symm!SeqA;
      pE := Symm!SeqE;

      f:=hom<BG -> Symm | pJ, pT, pU,pA,pE >;
      H:=sub<BG | f>;

      return AQInvariants(H);
end function;