/* List of functions:
   getMaxRootOfUnityToAddToCharZeroDeg2PID(R)
   getI(R)
   getMaxPossibleCosArg(K)
   getListOfValuesToCheckForCosArg(K)
   getCyclicsAsMG(K)
   getDihedralsAsMG(K, isByConj)
   getMatrixA_forA4_A5_andS4Star(R)
   getA4starAsMG(R)
   getA5starAsMG(R) 
   getS4starAsMG(R) 
   getFiniteInsolubleMatrixGroupsOfDimTwoOverPIDOfCharZero(R, isGL, isByConj)
   getFiniteSolubleMatrixGroupsOfDimTwoOverPIDOfCharZero(R, isGL, isByConj)   
   getFiniteMatrixGroupsOfDimTwoOverPIDOfCharZero(R, isGL, isByConj)	
**************************/

getMaxRootOfUnityToAddToCharZeroDeg2PID:= function(R)
	local K, u;
	
	if (Category(R) eq RngUPol) then
		while (Category(R) eq RngUPol) do // find ground ring of univariate polynomial ring
			R:= BaseRing(R);  
		end while;		
	end if;

	if (Category(R) in [FldRe, FldCom, FldRat, RngInt]) then return 2; end if;

	if (Category(R) eq FldNum) then R:=AbsoluteField(R); end if;
		
	K:= Integers(R);
	u:= UnitGroup(K).1;
	return Order(u);

end function;

getI:= function(R)
	return getRoots("x^2 + 1", R);
end function;

getMaxPossibleCosArg:= function(K)
	// uses following facts:
	//  deg(cos(2pi/n)) = 1  (n=1 or 2), phi(n)/2 (n > 2)
	//  bound of phi(n) <= sqrt(n) (n neq 2 or 6)
	// cf. Mitrinovic & Sandor, 1995, p. 9
		
	if (Category(K) eq RngUPol) then
		while (Category(K) eq RngUPol) do // find ground ring of univariate polynomial ring
			K:= BaseRing(K);  
		end while;		
	end if;
	
	return Max(4*Degree(K)^2, 6);
end function;

getListOfValuesToCheckForCosArg:= function(K)
	return [x : x in [3..getMaxPossibleCosArg(K)] | not(x in Divisors(getMaxRootOfUnityToAddToCharZeroDeg2PID(K)))];
end function;

getCyclicsAsMG:= function(K)	
	/* Generators as per Theorem 7.2.2 (i) */

	// ***Declarations***************************************** //

	local f, P, x, rts;
	
	//Lists
	local lst_cyclics;
	
	// Matrices
	local mat;

	// Matrix Group
	local MG;

	// ******************************************************* //
	
	P<x> := PolynomialRing(K);
	lst_cyclics:= [];
		
	for s in getListOfValuesToCheckForCosArg(K) do

		twoCosTwoPiOver_s:= RootOfUnity(s)+RootOfUnity(s)^(-1);
		f:= P! Eltseq(MinimalPolynomial(twoCosTwoPiOver_s));
		rts:= Roots(f);

		if not(IsEmpty(rts)) then
			mat:= GL(2,K)! Matrix([[0,-1],[1,rts[1][1]]]); //any root will do, so pick 1st
			MG:= MatrixGroup<2, K | mat>; 
			MG`CentralQuotient:= "C_" cat IntegerToString(s);
			MG`Description:= MG`CentralQuotient;
			MG`Field:= K;
			MG`Properties:= "Monomial, Soluble, Abelian";		
			lst_cyclics:= lst_cyclics cat [MG];
		end if;
	end for;
	
	return lst_cyclics;
	
end function;

getDihedralsAsMG:= function(K, isByConj)
	/* Generators as per Theorem 7.2.2 (ii)  */

	// ***Declarations***************************************** //
	
	local f, P, x, two, rts;
	
	//Lists
	local lst_dihedrals;
	
	// Matrices
	local matA, matB;

	// Matrix Group
	local MG, MGnc;

	// ******************************************************* //
	
	P<x> := PolynomialRing(K);
	lst_dihedrals:= [];
	
	matA:= GL(2,K)! Matrix([[0,1],
							[1,0]]);					
							
	two:= K! 2;
	
	for s in getListOfValuesToCheckForCosArg(K) do
	
		twoCosTwoPiOver_s:= RootOfUnity(s)+RootOfUnity(s)^(-1);
		f:= P! Eltseq(MinimalPolynomial(twoCosTwoPiOver_s));
		rts:= Roots(f);
		
		if not(IsEmpty(rts)) then			
			matB:= GL(2,K)! Matrix([[0,-1],[1,rts[1][1]]]); //any root will do, so pick 1st	
			MG:= MatrixGroup<2, K | matA, matB>; 
			MG`CentralQuotient:= "D_" cat IntegerToString(2*s);
			MG`Description:= MG`CentralQuotient;
			MG`Field:= K;
			MG`Properties:= "Monomial, Soluble";		
			lst_dihedrals:= lst_dihedrals cat [MG];
			
			// applying condition as outlined in Corollory 7.2.17:
			if not(IsUnit(two)) and (isByConj) and (IsOdd(s)) then
				MGnc:= MatrixGroup<2, K | -matA, matB>; 
				cloneAttributesOfMG(~MGnc, MG);
				MGnc`Description:= MGnc`Description cat " (NonConj Copy)";
				lst_dihedrals:= lst_dihedrals cat [MGnc];
			end if;
			
		end if;
	end for;
	
	return lst_dihedrals;
	
end function;

getMatrixA_forA4_A5_andS4Star:= function(R)
	local lstRts, rt, rt_inv, matA;
	
	if not(IsUnit(R! 2)) then return false, false; end if;

	lstRts:= getRoots("x^2 - x + 2^(-1)", R);
	if IsEmpty(lstRts) then return false, false; end if;

	if IsEmpty(lstRts) then 
		return false, false;
	else
		rt:= lstRts[1];
		rt_inv:= lstRts[2];	
		matA:= GL(2,R)! ( Matrix([[rt_inv, -rt],
								  [rt_inv,  rt]]));	
		return true, matA;
	end if;
end function;

getA4starAsMG:= function(R)
	/* Generators as per Lemma 3.0.6. (iii) */

	// ***Declarations***************************************** //

	local lstI;
	local i, i_inv;
	local isMatA;
	
	// Matrices
	local matA, matB;

	// Matrix Group
	local MG_A4star;

	// ******************************************************* //
	
	lstI:= getI(R);
	if IsEmpty(lstI) then return false, false; end if;

	isMatA, matA:= getMatrixA_forA4_A5_andS4Star(R);
	if not(isMatA) then return false, false; end if;	
	
	i:=lstI[1];
	i_inv:=lstI[2];	

	matB:= GL(2,R)! Matrix([[i_inv,0],
							[    0,i]]);	

	MG_A4star:= MatrixGroup<2, R | matA, matB>; 
	MG_A4star`CentralQuotient := "Alt(4)";
	MG_A4star`Description := "A4*";
	MG_A4star`Field:= R;
	MG_A4star`Properties:= "Primitive, Soluble";

	return true, MG_A4star;
end function;

getA5starAsMG:= function(R)
	/* Generators as per Lemma 3.0.6. (v) */

	// ***Declarations***************************************** //

	local P, g, y, m;
	local i, i_inv;
	local s, s_inv;
	
	// Matrices
	local matC;

	// Matrix Group
	local MG_A4star, MG_A5star;

	// ******************************************************* //
	
	isA4Star, MG_A4star:= getA4starAsMG(R);
	if not(isA4Star) then return false, false; end if;
	if not(IsUnit(R! 4)) then return false, false; end if;
		
	lstGratio:= getRoots("x^2 - x - 1", R);
	if IsEmpty(lstGratio) then return false, false; end if;
	s:= lstGratio[1];
	s_inv:= lstGratio[2];
	
	P<y>:= PolynomialRing(R);
	g:= y^2 - s*y + 3*(4^(-1));
	if (Degree(P! g) lt 2) then return false, false; end if;
	m:= Roots(g);

	lstI:= getI(R);
	i:= lstI[1];
	i_inv:= lstI[2];
	
	matC:= GL(2,R)! (Matrix([[  i_inv/2, m[2][1]*i],
					         [m[1][1]*i, i/2]]));	
									  
	MG_A5star:= MatrixGroup<2, R | MG_A4star.1, MG_A4star.2, matC>;
	MG_A5star`CentralQuotient:= "Alt(5)";
	MG_A5star`Description:= "A5*";
	MG_A5star`Field:= R;
	MG_A5star`Properties:= "Primitive, Soluble";

	return true, MG_A5star;
end function;

getS4starAsMG:= function(R)
	/* Generators as per Lemma 3.0.6. (iv) */

	// ***Declarations***************************************** //

	local lstRts;
	local i, ru8, ru8_inv;
	local isMatA;
	
	// Matrices
	local matA, matS;

	// Matrix Group
	local MG_S4star;

	// ******************************************************* //
	
	isMatA, matA:= getMatrixA_forA4_A5_andS4Star(R);
	if not(isMatA) then return false, false; end if;

	lstRts:= getRoots("x^4 + 1", R);
	if IsEmpty(lstRts) then return false, false; end if;
	ru8:= lstRts[1];
	ru8_inv:= ru8^(-1);
	
	matS:= GL(2,R)! Matrix([[ru8_inv,  0],
						    [      0,ru8]]);	

	MG_S4star:= MatrixGroup<2, R | matA, matS>; 
	MG_S4star`CentralQuotient := "Sym(4)";
	MG_S4star`Description := "S4*";
	MG_S4star`Field:= R;
	MG_S4star`Properties:= "Primitive, Soluble";

	return true, MG_S4star;
end function;

getFiniteInsolubleMatrixGroupsOfDimTwoOverPIDOfCharZero:= function(R, isGL, isByConj)
	/* cf. Theorem 7.3.6. (vi) */

	// ***Declarations***************************************** //

	local intMaxRu;
	local strSLExtraCond;
	local isA5star;

	// Lists
	local lst_a5star, lst_a5starxCr;
		
	// Matrix Groups
	local MG_A5star;

	// maps
	local ru, DM, DMoe;
	
	// ******************************************************* //
	
	P<x>:= PolynomialRing(R);
	two:= R! 2;
	
	if not(IsUnit(two)) then return []; end if;
	
	//A5*
	isA5star, MG_A5star:= getA5starAsMG(R);
	
	if (isA5star) then
	
		if (Category(R) eq FldCom) and (isGL) then
			return "Error: Cannot compute all insoluble irreducible finite matrix subgroups of GL(2,C) as list is infinite. \nYou can compute for SL(2,C) if you like.\n";
		end if;
		
		intMaxRu:= getMaxRootOfUnityToAddToCharZeroDeg2PID(R);
		ru:= getRootOfUnityMap(R, intMaxRu);

		DM:= map<Integers() -> GL(2,R) | r:-> GL(2,R)! DiagonalMatrix([ru(r),ru(r)])>;
		DMoe:= map<Integers() -> GL(2,R) | r:-> DM((2 - (r mod 2))*r)>;
		
		if (isGL) then
			strSLExtraCond := "true";
		else
			strSLExtraCond := "ru(r)^2 eq 1";
		end if;		
				
		lst_a5starxCr := [makeMG(R, [MG_A5star.1, MG_A5star.2, MG_A5star.3, DMoe(r)],
				<"A5*", "A5*" cat getCyclicDescCent(r,2), "(iii)", MG_A5star`Properties>) : r in Divisors(intMaxRu div 2) | r ne 1 and eval(strSLExtraCond)];
		lst_a5star := [MG_A5star] cat lst_a5starxCr;
		
		return [lst_a5star];
	else
		return [];
	end if;
end function;

getFiniteSolubleMatrixGroupsOfDimTwoOverPIDOfCharZero:= function(R, isGL, isByConj)
	/* cf. Theorem 7.3.6. (i)--(v) */

	// ***Declarations***************************************** //

	local intMaxRu;
	local strSLExtraCond;

	// Lists
	local lst_groups;
	local lst_cyclics, lst_dihedrals;
	local lst_dicxCr;
	local lst_a4star, lst_a4starxCr;
	local lst_s4star, lst_s4starxCr;
	
	// boolean
	local isA4star, isS4star;
	
	// Matrices
	local matDicB;
	
	// Matrix Groups
	local MG_A4star, MG_S4star;

	// maps
	local ru, DM, DMoe, SM;
	
	// ******************************************************* //
	
	if (Category(R) eq FldCom) and (isGL) then
		return "Error: Cannot compute all soluble irreducible finite matrix subgroups of GL(2,C) as list is infinite. \nYou can compute for SL(2,C) if you like.\n";
	end if;
	
	lst_groups := [];
	
	if (Category(R) eq FldRe) then
		printf "For each integer s>=3, there is an irreducible subgroup in ";
		if (isGL) then
			printf "GL";
		else
			printf "SL";
		end if;
		// C_s: (cf. Theorem 7.2.2. (i))
		printf "(2,R) isomorphic to C_s, each with a single generating matrix:\n\n";
		printf "   [ 0      -1    ]\n";
		printf "   [ 1 2cos(2pi/s)]";
		if (isGL) then
			// D_2s: (cf. Theorem 7.2.2. (ii))
			printf "\n\nand irreducible subgroup isomorphic to D_2s generated by matrices\n\n";
			printf "   [ 0  1]  [ 0      -1    ]\n";
			printf "   [ 1  0], [ 1 2cos(2pi/s)]";
		end if;
		printf ".\n\n";
	else 
		if not (Category(R) eq FldCom) then
			lst_groups := lst_groups cat [getCyclicsAsMG(R)];
			if (isGL) then
				lst_groups := lst_groups cat [getDihedralsAsMG(R, isByConj)];
			end if;
		end if;
	end if;
	
	intMaxRu := getMaxRootOfUnityToAddToCharZeroDeg2PID(R);
	ru:= getRootOfUnityMap(R, intMaxRu);

	DM:= map<Integers() -> GL(2,R) | r:-> GL(2,R)! DiagonalMatrix([ru(r),ru(r)])>;
	DMoe:= map<Integers() -> GL(2,R) | r:-> DM((2 - (r mod 2))*r)>;
	SM:= map<Integers() -> GL(2,R) | r:-> GL(2,R)! DiagonalMatrix([ru(r),ru(r)^(-1)])>;

	if (isGL) then
		strSLExtraCond := "true";
	else
		strSLExtraCond := "ru(r)^2 eq 1";
	end if;

	// Dicyclics (cf. Theorem 3.0.6. (ii))	
	if ((intMaxRu mod 2 eq 0) and (intMaxRu gt 2)) then
		matDicB:= GL(2,R)! Matrix([[0,-1],
								   [1, 0]]);
		lst_dicxCr := [makeMG(R, [(GL(2,R)! SM(2*n)), GL(2,R)! matDicB, DMoe(r)],
					<"Dic(" cat IntegerToString(n) cat ")", "Dic(" cat IntegerToString(n) cat ")" cat getCyclicDescCent(r,2), "(ii)", "Soluble">) : 
					  n in Divisors(intMaxRu div 2), r in Divisors(intMaxRu div 2) | n ge 2 and eval(strSLExtraCond)];
		lst_groups := lst_groups cat [lst_dicxCr];
	end if;
	
	//A4*
	isA4star, MG_A4star:= getA4starAsMG(R);	
	if (isA4star) then
		lst_a4starxCr := [makeMG(R, [MG_A4star.1, MG_A4star.2, DMoe(r)],
				<"A4*", "A4*" cat getCyclicDescCent(r,2), "(iii)", MG_A4star`Properties>) : r in Divisors(intMaxRu div 2) | r ne 1 and eval(strSLExtraCond)];
		lst_a4star := [MG_A4star] cat lst_a4starxCr;
		lst_groups := lst_groups cat [lst_a4star];
	end if;

	//S4*
	isS4star, MG_S4star:= getS4starAsMG(R);
	if (isS4star) then
		lst_s4starxCr := [makeMG(R, [MG_S4star.1, MG_S4star.2, DMoe(r)],
				<"S4*", "S4*" cat getCyclicDescCent(r,2), "(iv)", MG_S4star`Properties>) : r in Divisors(intMaxRu div 2) | r ne 1 and eval(strSLExtraCond)];
		lst_s4star := [MG_S4star] cat lst_s4starxCr;
		lst_groups := lst_groups cat [lst_s4star];
	end if;
		
	return lst_groups;
end function;

getFiniteMatrixGroupsOfDimTwoOverPIDOfCharZero:= function(R, isGL, isByConj)
	return getFiniteSolubleMatrixGroupsOfDimTwoOverPIDOfCharZero(R, isGL, isByConj) cat 
		   getFiniteInsolubleMatrixGroupsOfDimTwoOverPIDOfCharZero(R, isGL, isByConj);
end function;