Computation of $\mathcal{S}(G)$

IsContained= libgap.function_factory("\in")

class S_Group:
    
    def __init__(self, G, verbose= False, very_verbose= False):
        
        self.G= G
        self.GG= GG= G.DirectProduct(G)
        self.i1= GG.Embedding(1)
        self.i2= GG.Embedding(2)
        self.p1= GG.Projection(1)
        self.p2= GG.Projection(2)
        
        self.verbose= verbose or very_verbose
        self.very_verbose= very_verbose
        
        
    def compute_generating_pairs(self):
        """Computing pairs of generators up to conjugacy.
        
        This builds the list of all pairs (x, y) (as elements of GG) generating G, up to the diagonal conjugation action of G.
        Saves as self.gens_classes
        """
           
        #diagonal copy of G within G x G
        gensG= self.G.GeneratorsOfGroup()
        diagonal_gens= [ self.i1.Image(g) * self.i2.Image(g)  for g in gensG ]
        diagonal_emb= self.G.GroupHomomorphismByImages(self.GG, gensG, diagonal_gens);
        diagG= diagonal_emb.Image()
        
        #conj classes of pairs
        conj_classes= diagG.OrbitsDomain(self.GG)
        
        if self.verbose:
            print "***", len(conj_classes), "conjugacy classes of pairs"
        
        #keep only the generating pairs
        
        def is_generating_pair(C):
            pair= C[0]
            x= self.p1.Image(pair)
            y= self.p2.Image(pair)
            
            return self.G.Subgroup([x, y]) == self.G
        
        self.gens_classes= [ C for C in conj_classes if is_generating_pair(C) ]

        if self.verbose:
            print "***", len(self.gens_classes), "of them generate G"
    
    
    
    def compute_orbits_on_pairs(self):
        """ Orbits of G x G on pairs
        
        This builds a list of lists. The elements of the inner lists are indices, indicating the elements 
        in gens_classes belonging to the same G x G orbit.
        
        Saves as self.GG_orbits
        """
    
        if self.verbose:
            print "*** computing orbits of GxG"
    
    
        #group them up
        indices= range(len(self.gens_classes))
        partition= []
        
        while len(indices) > 0:
            
            if self.very_verbose:
                print "computing GxG orbits,", len(indices), "indices remaining"
            
            index= indices[0]
            pair= self.gens_classes[index][0] # need [0] to take a representative
            o= self.GG.ConjugacyClass(pair).AsList()
            related= [i for i in indices if self.gens_classes[i][0] in o]
            # clean up the list of indices
            indices= [i for i in indices if i not in related]
            # save
            partition.append( related )
            
        self.GG_orbits= partition # careful! numbered from 0 !
        

    def compute_theta(self):
        """ permutation associated to theta
        
        This computes the permutation induced by the operation (x, y) -> (y, x) on pairs in GG. 
        Saves as self.theta, a list of pairs (=cycles of length 2), numbered from 0
        """
    
    
        if self.verbose:
            print "*** computing theta"
    
        indices= range(len(self.gens_classes))
        pairs= []

        while len(indices) > 0:
            #info
            if self.very_verbose:
                print "computing theta,", len(indices), "remaining"
            
            index= indices[0]
            pair= self.gens_classes[index][0]
            #now let's reverse it
            x= self.p1.Image(pair)
            y= self.p2.Image(pair)
            reverse= self.i1.Image(y) * self.i2.Image(x)    #pair^(p2*i1) * pair^(p1*i2);

            #now where is it in the list?
            found= False
            i= 0
            while not(found):
                if reverse in self.gens_classes[i]:
                    found= True
                else:
                    i= i+1

            # reverse is in position i
            indices.remove(i)
            # have an actual cycle?
            if index != i:
                pairs.append( (index, i) )
                indices.remove(index)

        self.theta= pairs
        

    def compute_delta(self):
        """ permutation associated to delta
        
        This computes the permutation induced by the operation (x, y) -> (y^(-1)x^(-1), y) on pairs in GG. 
        Saves as self.delta, a list of pairs (=cycles of length 2), numbered from 0
        """
    
    
        if self.verbose:
            print "*** computing delta"
    
        indices= range(len(self.gens_classes))
        pairs= []

        while len(indices) > 0:
            #info
            if self.very_verbose:
                print "computing delta,", len(indices), "remaining"
            
            index= indices[0]
            pair= self.gens_classes[index][0]
            #now let's reverse it
            x= self.p1.Image(pair)
            y= self.p2.Image(pair)
            reverse= self.i1.Image(y^(-1)*x^(-1)) * self.i2.Image(y)    #pair^(p2*i1) * pair^(p1*i2);

            #now where is it in the list?
            found= False
            i= 0
            while not(found):
                if reverse in self.gens_classes[i]:
                    found= True
                else:
                    i= i+1

            # reverse is in position i
            indices.remove(i)
            # have an actual cycle?
            if index != i:
                pairs.append( (index, i) )
                indices.remove(index)

        self.delta= pairs

        
        
    def permutation_associated_to(self, f):
        """ permutation associated to the operation f on pairs. 
        
        For example for theta one takes f= lambda x, y : (y, x). In general f(x, y) must be a pair.
        This returns the list of images of the indices in self.gens_classes, under the operation f.
        
        NOTE: as suggested, it is possible to use this to compute theta, for example,
        with f= lambda x, y : (y, x). While this is a good consistency test, it is
        best to call self.compute_theta (or delta), which uses the fact that theta is
        of order two (speeding up things by a factor of 2).
        
        This is normally used for the elements of Out(G).
        """
    
        if self.verbose:
            print "*** computing a permutation"
    
    
        indices= range(len(self.gens_classes))
        images= []

        for k in indices:
            
            if self.very_verbose:
                print len(indices)-k, "indices remaining"
            
            
            pair= self.gens_classes[k][0]
            #now let's apply f
            x= self.p1.Image(pair)
            y= self.p2.Image(pair)
            newpair= f(x, y)
            xprime= newpair[0]
            yprime= newpair[1]
            reverse= self.i1.Image(xprime) * self.i2.Image(yprime) #name 'reverse' thinking about theta

            #now where is it in the list?
            found= False
            i= 0
            while not(found):
                if reverse in self.gens_classes[i]:
                    found= True
                else:
                    i= i+1

            # reverse is in position i
            images.append(i)
        
        return images

    def compute_gens_outG(self):
        """ Generators for the outer automorphism group.

        Saved as self.gens_outG.
        """
        autG= self.G.AutomorphismGroup()
        innG= autG.InnerAutomorphismsAutomorphismGroup()
        gens= [ x for x in autG.GeneratorsOfGroup() if not(IsContained(x, innG)) ]
        self.gens_outG= gens
        
        
    def compute_all(self, reset= False):
        if not(hasattr(self, "gens_classes")) or reset:
            self.compute_generating_pairs()
        if not(hasattr(self, "GG_orbits")) or reset:
            self.compute_orbits_on_pairs()
        if not(hasattr(self, "theta")) or reset:
            self.compute_theta()
        if not(hasattr(self, "delta")) or reset:
            self.compute_delta()
        if not(hasattr(self, "gens_outG")) or reset:
            self.compute_gens_outG()
        if not(hasattr(self, "perms_outG")) or reset:
            self.perms_outG= []
            for a in self.gens_outG:
                f= lambda x, y : (a.Image(x), a.Image(y))
                self.perms_outG.append(self.permutation_associated_to(f))
        
        if self.verbose:
            print "*** preliminaries done"
            
    def setup(self):
        """ Computes all background information, and converts everything into GAP permutations and groups.
        
        self.compute_all() does all the hard work, but the data computed are in raw form:
        permutations are on [0..n] instead of [1..(n+1)], some are given as lists of pairs,
        others as list of images of [0..n]... This does the necessary conversions.
        
        Once the conversions are done, memory is freed by deleting the original data. Self loses
        some attributes, and as a result, if compute_all() or setup() are called again, then
        automatically all computations will start from scratch.        
        """
        
        # make sure everything is ready:
        self.compute_all() # results are cached so this causes no extra work if the computations are already done
        
        #convenient:
        def add_one(L):
            return tuple( x+1 for x in L )
        
        # GxG orbits
        self.blocks= [ add_one(orb) for orb in self.GG_orbits ]
        # theta
        pairs= [ add_one(pair) for pair in self.theta]
        theta= libgap.eval( Permutation(pairs).cycle_string() )
        # delta
        pairs= [ add_one(pair) for pair in self.delta]
        delta= libgap.eval( Permutation(pairs).cycle_string() )
        # outG
        perms= [ libgap.PermList(add_one(p)) for p in self.perms_outG ]
        # group H
        ell= len(self.gens_classes)
        sym= libgap.SymmetricGroup(ell)
        self.H= sym.Subgroup([theta, delta] + perms )
        self.ell= ell
        
        #free some memory
        # only self.H and self.blocks (and self.ell) are needed for factors()
        delattr(self, "gens_classes")
        delattr(self, "GG_orbits")
        delattr(self, "theta")
        delattr(self, "delta")
        delattr(self, "gens_outG")
        delattr(self, "perms_outG")
        
        if self.verbose:
            print "*** setup complete"
        
        
    def factors(self, reset= False):
        """Returns a list of groups, the direct product of which is the S-group of self.G"""
        
        #check if the results were cached earlier:
        if hasattr(self, "factors_list") and not(reset):
            return self.factors_list

        # check whether the 'background information' is available
        if not(hasattr(self, "H")):
            self.setup() 
        
        if self.verbose:
            print "*** computing the factor groups"
        
        #compute the orbits of H
        orbs= self.H.OrbitsDomain([1..self.ell])

        #we gather information on these orbits
        info_orbs= [] # this will be a list of dictionnaries
        conj_cls_H= self.H.ConjugacyClassesSubgroups()
        nb_classes= len(conj_cls_H)

        for orb in orbs:  # for each orbit, record a few things
    
            #compute the stabilizer of an element
            #or rather the conjugacy class of this stabilizer
            stab= self.H.Stabilizer(orb[0]) #now its class
            for i in range(nb_classes):
                if IsContained(stab, conj_cls_H[i]):
                    stab_num= i
                    break
            #ok, stab_num is the number of the conjugacy class of the stabilizer of elements of orb
            
            #now the sizes of the intersections of orb with all blocks:
            # note: intersection of a list of sage integers and a list of GAP integers...
            size_subsets= [len( set(block).intersection(set(orb.sage())) ) for block in self.blocks]
        
            #now store the info before moving on to next orbit
            info_orbs.append({"stab" : stab_num, "sizes" : size_subsets})
            

        # we put together the elements of orbs according to their information in info_orbs
        # orbits with the same info are collected in a single "packet" -- here a list if indices refering to orbs
        # the union of the orbits in one packet is stable under the action of H

        indices= range(len(orbs))
        packets= []
        
        while len(indices) > 0:
            i= indices[0]
            
            packet=  [j for j in indices if info_orbs[j] == info_orbs[i]]
            indices= [j for j in indices if j not in packet]
                
            packets.append( packet )
            
            
        #OK now compute small centralizers
        gensH= self.H.GeneratorsOfGroup()
        progress= 0
        factors= []
    
        for packet in packets:
        
            progress= progress + 1
            if self.very_verbose:
                print "progress:", progress, "out of", len(packets)

    
            #find the union of the orbits in packet
            domain= [x for i in packet for x in orbs[i]]
            #look at the permutation group induced by H,
            #and its centralizer
            newgens= [g.RestrictedPermNC(domain) for g in gensH]
            newH= libgap.Group(newgens)
            S=    libgap.SymmetricGroup(domain)
            Clocal= S.Centralizer(newH)
    
    
            #okay now compute young subgroup 
            gens_young= [ S.Identity() ] 
            domain= [x.sage() for x in domain] # annoying: make sure we use sage integers from now on
            for block in self.blocks:
                subset= list( set(block).intersection(set(domain)) )

                if len(subset) > 1:
                    gens_young.extend( libgap.SymmetricGroup(subset).GeneratorsOfGroup() )

            young_sub= libgap.Group( gens_young )
    
            #compute the intersection
            factor= young_sub.Intersection(Clocal)
            if not(factor.IsTrivial()):
                
                factors.append(factor)
                
                if self.very_verbose:
                    print "found non-trivial factor!"
        
        self.factors_list= factors # cache the results
        return factors

G= libgap.PSL(3, 2)

S= S_Group(G, verbose= True)

S.factors()

L= S.factors() #results were cached

[gp.StructureDescription() for gp in L]