// ******************************* Code **************************************** 
 
// The function "BraidMove" applies the i-th braid move on "S"

BraidMove:=function(S,i)
	c:=S[i]*S[i+1]*S[i]^-1;
return Insert(Remove(S,i+1),i,c); end function;

// The function "BraidOrbit" determines the orbit of a given generating 
// triple "S" under the braid group. 

BraidOrbit:=function(S)
    orb:={ }; Trash:={S}; 
    repeat
        ExtractRep(~Trash,~gens); Include(~orb, gens);
		for i in [1,2] do newgens:=BraidMove(gens,i);
	    if newgens notin orb then 
			    Include(~Trash, newgens);
		    end if; 
	    end for;
    until IsEmpty(Trash);
    return orb;
end function;

// The function "STriplesC" determines the spherical generating triples "S:=[g1,g2,g3]" for "G" of type 
// "[d,d,l]", where "g3" is contained in the normal subgroup "A" and "g1, g2" are outside of "A".

STriplesC:=function(G,A,d,l)
    AC:={}; SphTripC:={};
    n:=#G;
    for g in G do 
        if g notin A and Order(g) eq d then 
            AC:=Include(AC,g);
        end if;
    end for;

    for g1 in AC do 
    for g2 in AC do 
        g3:=(g1*g2)^-1;
        if g3 in A and Order(g3) eq l and #sub<G|g1,g2> eq n then 
            SphTripC:=Include(SphTripC,[g1,g2,g3]);
        end if;
    end for;
    end for;
    
    return SphTripC;
end function;

// The function "ElsOfOrd" determines all elements of order "n" in a given group "G". 

ElsOfOrd:=func<G, d | {g: g in G| Order(g) eq d}>;

// The function "STriplesE" determines all (ordered) spherical generating triples "S:=[g1,g2,g3]" 
// of "G" of type "T=[n1,n2,n3]". 

STriplesE:=function(G,T)
    n1:=T[1]; n2:=T[2]; n3:=T[3];
    SphTripE:={}; n:=#G;

    for g1 in ElsOfOrd(G,n1) do 
    for g2 in ElsOfOrd(G,n2) do 
       g3:=(g1*g2)^-1;
            if Order(g3) eq n3 and #sub<G|g1,g2> eq n then 
                SphTripE:=Include(SphTripE,[g1,g2,g3]);
            end if;
    end for; end for;

    return SphTripE;
end function;

// The function "BraidOrbitReps" returns for a given set of spherical generating triples 
// precisely one representative for each orbit under the action of the braid group. 

BraidOrbitReps:=function(SetOfSpericalGens)
    Repres:=[];
	while not IsEmpty(SetOfSpericalGens) do 
		V:=Rep(SetOfSpericalGens); 
		Append(~Repres,V); 
		orb:=BraidOrbit(V);
		for W in orb do 
			Exclude(~SetOfSpericalGens, W);
		end for; 
	end while; 

    return Repres;
end function;

// The function "AutTup" performs the simultaneous conjugation of a generating triple "S"  by 
// an automorphism of the group. 

AutTup:=function(phi,S);
    return [phi(S[1]),phi(S[2]),phi(S[3])];
end function;

// The function "BraidAndAutOrbitReps" returns for a given set of spherical generating triples 
// precisely one representative for each orbit under the action of "Aut(G) x B_3". 

BraidAndAutOrbitReps:=function(SetOfSpericalGens,G)
    Aut:=AutomorphismGroup(G); 
    f,P:= PermutationRepresentation(Aut);   
    Repres:=[];
	while not IsEmpty(SetOfSpericalGens) do 
		S:=Rep(SetOfSpericalGens); 
		Append(~Repres,S); 
		orb:=BraidOrbit(S);
			for T in orb do 
                        for p in P do 
                        phi:=p @@ f; 
				Exclude(~SetOfSpericalGens, AutTup(phi,T));
			end for; 
                        end for;
	end while; 

    return Repres;
end function;

// The function "BraidAutPairReps" checks if two given pairs of spherical generating triples belong 
// to the same orbit  of "Aut(G) x B_3 x B_3", where "Aut(G)" simultaneously and "B_3" separately on
// each generating triple. 

BraidAutPairReps:=function(Pair1,Pair2,G)
    Aut:=AutomorphismGroup(G); 
    f,P:= PermutationRepresentation(Aut);
    for S in BraidOrbit(Pair1[1]) do
    for T in BraidOrbit(Pair1[2]) do 
    for p in P do 
        phi:=p @@ f; 
        Snew:=AutTup(phi,S); Tnew:=AutTup(phi,T);
        if [Snew,Tnew] eq Pair2 then 
            return true;
        end if;
    end for; end for; end for;
    
    return false;
end function;

// The function "ConjTup" performs the "complex conjugation" of a generating triple "S". 

ConjTup:=function(S);
    return [S[1]^-1,S[1]*S[3],S[3]^-1];
end function;

// With the procedure "ClassificationProcedure" we verify the claims in the proof of Theorem 6.6.

ClassificationProcedure:=procedure(G,A,d,l,T)
    SetE:=STriplesE(G,T);
    SetC:=STriplesC(G,A,d,l);
    
    numberC:=#BraidAndAutOrbitReps(SetC,G);
    print "Cardinality of S_C:", numberC;
    numberE:=#BraidOrbitReps(SetE);
    print "Cardinality of S_E:", numberE;
    
    S:=BraidAndAutOrbitReps(SetC,G)[1];
    S1:=BraidOrbitReps(SetE)[1];    
    S2:=BraidOrbitReps(SetE)[2];
    
    test:=S1[3]*S2[3] in A;
    print "The two generating triples S1, S2 in the set S_E induce different canonical representations:", test;
    
    sameorbit:=BraidAutPairReps([S1,S],[S2,S],G);
    print "The pairs [S1,S] and [S2,S] belong to the same orbit:", sameorbit;
    
    conjsameorbit:=BraidAutPairReps([S1,S],[ConjTup(S2),ConjTup(S)],G);
    print "The [S1,S] and the complex conjugate of [S2,S] belong to the same orbit:", conjsameorbit;
end procedure;


// Finally, we define the 3 groups "Gd=A \rtimes Z_d" together with the branching signature "T" correspronding to 
// the elliptic curve "E" and the numbers "d" and "l" for the branching signature "[d,d,l]" corresponding to C_d

G3:=SmallGroup(21,1);
t3:=G3.2; A3:=sub<G3|t3>;
d3:=3; l3:=7;
T3:=[3,3,3];

G4:=SmallGroup(20,3);
t4:=G4.3; A4:=sub<G4|t4>;
d4:=4; l4:=5;
T4:=[2,4,4];

G6:=SmallGroup(18,3);
t6:=G6.3; A6:=sub<G6|t6>;
d6:=6; l6:=3;
T6:=[2,3,6];

// Here we run the classification procedure: 

ClassificationProcedure(G3,A3,d3,l3,T3);