/*
This Magma code example is from Davide Lombardo's page:

https://people.dm.unipi.it/lombardo/files/Lab.m
*/ 

/******
Elementary approach to enumerating groups of order 2^n, n \leq 6
*******/

/*
Compute the Cayley embedding of a finite group G
*/
function CayleyEmbedding(G)
	// Get the order of the group
	n := #G;

	// Define the symmetric group on n elements
	S := SymmetricGroup(n);
	
	// Get the elements of G
	elements := [g : g in G];

	// Create a map from G to S_n
	CayleyMap := hom<G -> S | g :-> S![Index(elements, g * elements[i]) : i in [1..n]]>;

	return Image(CayleyMap), CayleyMap;
end function;

/*
Construct the "double" of a permutation sigma
*/
function DoublePermutation(sigma)
	n := Degree(Parent(sigma));
	S2n := Sym(2*n);
	L := [i^sigma : i in [1..n]] cat [i^sigma + n : i in [1..n]];
	
    return S2n!L;
end function;

/*
Construct the "double" of a group G
*/
function DoubleGroup(G)
	H := CayleyEmbedding(G);
	n := #H;
	HDouble := sub< Sym(2*n) | [ DoublePermutation(h) : h in Generators(H) ] >;
	
    return HDouble;
end function;


/*
Given a group G, construct a list of groups H
such that G < H with index 2, and up to isomorphism,
every group H with this property appears in the list.

The first version is elementary; the second uses the
correspondence theorem for subgroups; and the third
one optimises by not repeating subgroups that are
obviously conjugate to one another and only considering
a 2-Sylow subgroup of the normaliser.
*/
function ConstructDoubleCovers(G)
	GDouble := DoubleGroup(G);
	N := Normaliser( Sym(2*#G), GDouble );
	subs := Subgroups(N : OrderEqual := 2*#G );
	subs := [sub`subgroup : sub in subs];
	subs := [H : H in subs | GDouble subset H ];
	
    return subs;
end function;

function ConstructDoubleCovers2(G)
	GDouble := DoubleGroup(G);
	N := Normaliser( Sym(2*#G), GDouble );
	Quoziente, Proiezione := N/GDouble;
	Sezione := Proiezione^(-1);
	Els2 := [q : q in Quoziente | Order(q) eq 2];
	subs := [ sub< N | GDouble, Sezione(q) > : q in Els2 ];
	
    return subs;
end function;

function ConstructDoubleCovers3(G)
	"Computing double";
	time GDouble := DoubleGroup(G);
	"Computing normaliser";
	time N := Normaliser( Sym(2*#G), GDouble );
	"Computing Sylow subgroup";
	Sylow2 := SylowSubgroup(N, 2);
	"Computing quotient";
	"Order of quotient:", #Sylow2 / #GDouble;
	time Quoziente, Proiezione := Sylow2/GDouble;
	"Computing section";
	time Sezione := Proiezione^(-1);
	"Computing elements of order 2 in the quotient";
	time Els2 := Subgroups(Quoziente : OrderEqual := 2);
	"Computing subgroups";
	time subs := [ sub< N | GDouble, Sezione(q`subgroup) > : q in Els2 ];
	
    return subs;
end function;


/*
Given a list L of groups, returns a list L' that contains
every isomorphism class of groups in L precisely once
*/
function FilterDuplicates(list)
	CleanList := [];
	for H in list do
		test := true;
		for H2 in CleanList do
			if test then
				test := test and not IsIsomorphic(H, H2);
			end if;
		end for;
		if test then
			CleanList := CleanList cat [H];
		end if;
	end for;
	
    return CleanList;
end function;

function FilterDuplicatesFast(list)
	CleanList := [];
	CleanNames := [];
	for H in list do
		if not GroupName(H) in CleanNames then
			CleanList := CleanList cat [H];
			CleanNames := CleanNames cat [GroupName(H)];
		end if;
	end for;
	
    return CleanList;
end function;


/*
Given a list L of groups [G_i], returns a list L' = [H_j]
where each H_j contains some G_i with index 2, every
group with this property appears in L', and L' contains no
duplicates up to isomorphism.
*/
function DoubleListSlow(L)
	LNew := [];
	for G in L do
		dc := ConstructDoubleCovers(G);
		dc := FilterDuplicates(dc);
		LNew := LNew cat dc;
	end for;
	
    return FilterDuplicates(LNew);
end function;

function DoubleList(L)
	LNew := [];
	for G in L do
		time dc := ConstructDoubleCovers3(G);
		time dc := FilterDuplicatesFast(dc);
		LNew := LNew cat dc;
	end for;
	
    return FilterDuplicatesFast(LNew);
end function;


list2 := [CyclicGroup(2)];
time list4 := DoubleList(list2);
time list8 := DoubleList(list4);
time list16 := DoubleList(list8);
time list32 := DoubleList(list16);
time list64 := DoubleList(list32);

assert #list32 eq NumberOfSmallGroups(32);
assert #list64 eq NumberOfSmallGroups(64);