AttachSpec("klngpspec"); /* Attaches the program written by A. Page */



import "bianchi.m" : QuatToMatrix; /* Returns matrix generators */


//precomputation, common to every group (with the same precision) : useful for intensive computations
import "geometry/volumes.m" : ComputeZetas;
print "precomputing coefficients...";
zetas := ComputeZetas(100);
print "...done.\n\n";

// The function below is based on the algorithm of Haluk Sengun and computes
// abelianizations for congruence subgroups of certain Kleinian subgroups
// Input: - prime : a rational prime number;
// 				- order : a positive integer which, together with "prime", will determine
// 									pricesely the prime ideal in the factorization of (prime*O_K)
// 				- power : the power of the prime ideal determined above by the paird (prime, order);
// 									the triple (prime,order,power) determines a prime power ideal I of O_K
//					- G, Generators : a Kleinian subgroup (finitely presented) together with its Generators;
// 													 these are obtained via the algorithm and implementation of A. Page;
// Output : The abelianization of the \Gamma_0 congruence subgroup of G determined by the ideal I;
// 					the output is presented as (d_1,d_2, ..., d_n, 0,...,0) where
// 					d_1 | d_2 | ... | d_n are natural numbers and (0,..,0) are r copies of 0
//					correspond to the presentation of the abelianization a finitely generated abelian group
//					Z_{d_1}+...+Z_{d_2}+(Z)^r
Gamma0_abel:=function(prime, ord, power,G, Generators)
	Mat0 := [QuatToMatrix(g) : g in Generators];
	K := Parent(Mat0[1][1,1]);
	OK := MaximalOrder(K);
	R := MatrixAlgebra(OK,2);
	MatGen := [R ! m : m in Mat0];
	/* Append(~MatGen, (R ! Matrix(OK,2,2,[-1,0,0,1])));*/
	I := (Factorization(prime*OK)[ord,1])^(power);

	PL,r := ProjectiveLine(quo<OK|I>);
	proj_action := function(M,i)
		t,im :=r(PL[i]*M,true,false);
		return Index(PL,im);
	end function;

	Seq := [[] : m in MatGen];
	for j in [1..#Seq] do
		for i in [1..#PL] do
			Append(~Seq[j], proj_action(MatGen[j],i));
		end for;
	end for;

	Symm := Sym(#PL);
	p := [Symm!el : el in Seq];

	f := hom<G -> Symm | p>;
	H := sub<G | f>;

	return AQInvariants(H);
end function;

// function that uses Aurel Page's program to return generators of PGL2(O_K)
PresPGL2 := function(d)
	O := BianchiOrder(d);
	_, Faces, Edges := NormalizedBasis(O : GroupType :="Units", zetas := zetas);
	PG, PGenerators := Presentation(Faces, Edges, O);
	return PG, PGenerators;
end function;

Abelianization := function(d, prime, ord, power)
	PG, PGenerators := PresPGL2(d);
	return Gamma0_abel(prime,ord, power, PG, PGenerators);
end function;