#############################################################################
##
##  This file is part of GAP, a system for computational discrete algebra.
##  This file's authors include Ákos Seress.
##
##  Copyright of GAP belongs to its developers, whose names are too numerous
##  to list here. Please refer to the COPYRIGHT file for details.
##
##  SPDX-License-Identifier: GPL-2.0-or-later
##


#############################################################################
##
#F  GInverses( <S> )  . . . . . . . . . . . . . . . . . . . . . . . . . local
##
##  <S> must be a stabilizer chain.
##
##  `GInverses' changes `<S>.generators' !
##
BindGlobal( "GInverses", function( S )
    local   inverses,  set,  i;

    set := Set( S.translabels );
    RemoveSet( set, 1 );
    S.generators := S.labels{ set };
    inverses := [  ];
    for i  in [ 1 .. Length( S.generators ) ]  do
        inverses[ i ] := S.generators[ i ] ^ -1;
    od;
    if IsBound( S.stabilizer )  then
        Append( S.generators, S.stabilizer.generators );
    fi;
    return inverses;
end );

#############################################################################
##
#F  DisplayCompositionSeries( <S> ) . . . . . . . . . . . .  display function
##
InstallGlobalFunction( DisplayCompositionSeries, function( S )
    local   f,  i,s;

    # ok, we accept groups too
    if IsGroup( S )  then
        S := CompositionSeries( S );
    fi;

    # if we know the composition series, we know orders of groups, so we may
    # enforce their computation before calling GroupString to display them.
    Perform( S, Size );

    Print( GroupString( S[1], "G" ), "\n" );
    for i  in [2..Length(S)]  do
      f:=Image(NaturalHomomorphismByNormalSubgroup(S[i-1],S[i]));
      s:=IsomorphismTypeInfoFiniteSimpleGroup(f);
      if IsBound(s.shortname) then
        s:=s.shortname;
      else
        s:=s.name;
      fi;
      Print( " | ",s,"\n");
      if i < Length(S)  then
        Print( GroupString( S[i], "S" ), "\n" );
      else
        Print( GroupString( S[i], "1" ), "\n" );
      fi;
    od;
end );

#############################################################################
##
#M  CompositionSeries( <G> )  . . . . composition series of permutation group
##
##  `CompositionSeriesPermGroup' returns the composition series of <G>  as  a
##  list.
##
##  The subgroups in this list have a slightly modified
##  `NaturalHomomorphismByNormalSubgroup' method,
##  which notices if you compute the factor group of one subgroup by the next
##  and return the factor group as a  primitive  permutation  group  in  this
##  case (which is also computed by the function below).  The  factor  groups
##  remember the natural homomorphism since the images of the  generators  of
##  the subgroup are known and the natural  homomorphism can thus be  written
##  as `GroupHomomorphismByImages'.
##
##  The program works for  permutation  groups  of  degree  < 2^20 = 1048576.
##  For higher degrees  `IsSimple'  and  `CasesCSPG'  must  be  extended with
##  longer lists of primitive  groups  from  extensions  in  Kantor's  tables
##  (see JSC. 12(1991), pp. 517-526).  It may also be  necessary  to  modify
##  `FindNormalCSPG'.
##
##  A general reference for the algorithm is:
##  Beals-Seress, 24th Symp. on Theory of Computing 1992.
##
InstallMethod( CompositionSeries,
    "for a permutation group",
    true,
    [ IsPermGroup ], 0,
    function( Gr )
    local   pcgs,
            normals,    # first component of output; normals[i] contains
                        # generators for i^th subgroup in comp. series
            factors,    # second component of output; factors[i] contains
                        # the action of generators in normals[i]
            factorsize, # third component of output; factorsize[i] is the
                        # size of the i^th factor group
            homlist,    # list of homomorphisms applied to input group
            auxiliary,  # if auxiliary[j] is bounded, it contains
                        # a subg. which must be added to kernel of homlist[j]
            index,      # variable recording how many elements of normals are
                        # computed
            workgroup,  # the subnormal factor group we currently work with
            workgrouporbit,
            lastpt,     # degree of workgroup
            tchom,      # transitive constituent homomorphism applied to
                        # intransitive workgroup
            bhom,       # block homomorphism applied to imprimitive workgroup
            fahom,      # factor homomorphism to store
            bl,         # block system in workgroup
            D,          # derived subgroup of workgroup
            top,        # index of D in workgroup
            lenhomlist, # length of homlist
            i, s,  t,   #
            fac,        # factor group as permutation group
            list;       # output of CompositionSeries

    # Solvable groups first.
    pcgs := Pcgs( Gr );
    if pcgs <> fail  then
        list := ShallowCopy( PcSeries( pcgs ) );
        s := list[ 1 ];
        for i  in [ 2 .. Length( list ) ]  do
            t := AsSubgroup( s, list[ i ] );
            fac := CyclicGroup( IsPermGroup,
                           RelativeOrders( pcgs )[ i - 1 ] );
            fahom:=GroupHomomorphismByImagesNC( s, fac,
                 pcgs{ [ i - 1 .. Length( pcgs ) ] },
                 Concatenation( GeneratorsOfGroup( fac ),
                 List( [ i .. Length( pcgs ) ], k -> One( fac ) ) ) );
            Setter( NaturalHomomorphismByNormalSubgroupInParent )( t,fahom);
            AddNaturalHomomorphismsPool(s,t,fahom);
            list[ i ] := t;
            s := t;
        od;
        return list;
    fi;

    # initialize output and work arrays
    normals := [];
    factors := [];
    factorsize := [];
    auxiliary := [];
    homlist := [];
    index := 1;
    workgroup := Gr;

    # workgroup is always a factor group of the input Gr such that a
    # composition series for Gr/workgroup is already computed.
    # Try to get a factor group of workgroup
    while (Size(workgroup) > 1) or (Length(homlist) > 0) do
#Print(List(normals,Length)," ",Size(workgroup),"\n");
        if Size(workgroup) > 1  then
            lastpt := LargestMovedPoint(workgroup);

            # if workgroup is not transitive
            workgrouporbit:= StabChainMutable( workgroup ).orbit;
            if Length(workgrouporbit) < lastpt   then
                tchom :=
                  ActionHomomorphism(workgroup,workgrouporbit,"surjective");
                Add(homlist,tchom);
                workgroup := Image(tchom,workgroup);
            else
                bl := MaximalBlocks(workgroup,[1..lastpt]);

                # if workgroup is not primitive
                if Length(bl) > 1  then
                    bhom:=ActionHomomorphism(workgroup,bl,OnSets,"surjective");
                    workgroup := Image(bhom,workgroup);
                    Add(homlist,bhom);
                else
                    D := DerivedSubgroup(workgroup);
                    top := Size(workgroup)/Size(D);

                    # if workgroup is not perfect
                    if top > 1  then

                        # fill up workgroup/D by cyclic factors
                        index := NonPerfectCSPG(homlist,normals,factors,
                                 auxiliary,factorsize,top,index,D,workgroup);
                        workgroup := D;

                    # otherwise chop off simple factor group from top of
                    # workgroup
                    else
                        workgroup := PerfectCSPG(homlist,normals,factors,
                                       auxiliary,factorsize,index,workgroup);
                        index := index+1;
                    fi;  # nonperfect-perfect

                fi;  # primitive-imprimitive

            fi;  # transitive-intransitive

        # if the workgroup was trivial
        else
            lenhomlist := Length(homlist);

            # pull back natural homs
            PullBackNaturalHomomorphismsPool(homlist[lenhomlist]);

            workgroup := KernelOfMultiplicativeGeneralMapping(
                             homlist[lenhomlist] );

            # if auxiliary[lenhmlist] is bounded, it is faster to augment it
            # by generators of the kernel of `homlist[lenhomlist]'
            if IsBound(auxiliary[lenhomlist])  then
                workgroup := auxiliary[lenhomlist];
                workgroup := ClosureGroup( workgroup, GeneratorsOfGroup(
                                KernelOfMultiplicativeGeneralMapping(
                                    homlist[lenhomlist] ) ) );
            fi;
            Unbind(auxiliary[lenhomlist]);
            Unbind(homlist[lenhomlist]);

        fi; # workgroup is nontrivial-trivial
    od;

    # loop over the subgroups
    #s := SubgroupNC( Gr, normals[1] );
    #SetSize( s, Size( Gr ) );
    s:=Gr;
    list := [ s ];
    for i  in [2..Length(normals)]  do
        t := SubgroupNC( s, normals[i] );
        SetSize( t, Size( s ) / factorsize[i-1] );
        fac := GroupByGenerators( factors[i-1] );
        SetSize( fac, factorsize[i-1] );
        SetIsSimpleGroup( fac, true );
        fahom:=GroupHomomorphismByImagesNC( s, fac,
                        normals[i-1], factors[i-1] );
        #if IsIdenticalObj(Parent(t),s) then
        #  Setter( NaturalHomomorphismByNormalSubgroupInParent )( t,fahom);
        #fi;
        AddNaturalHomomorphismsPool(s, t,fahom);
        Add( list, t );
        s := t;
    od;
    t := TrivialSubgroup( s );
    Assert(1,Size( s )=factorsize[Length(normals)]);
    fac := GroupByGenerators( factors[Length(normals)] );
    SetSize( fac, factorsize[Length(normals)] );
    SetIsSimpleGroup( fac, true );
    fahom:=GroupHomomorphismByImagesNC( s, fac,
                    Last(normals), factors[Length(normals)] );
    if IsIdenticalObj(Parent(t),s) then
      Setter( NaturalHomomorphismByNormalSubgroupInParent )( t,fahom);
    fi;
    AddNaturalHomomorphismsPool(s, t,fahom);
    Add( list, t );

    # return output
    return list;
end );


#############################################################################
##
#F  NonPerfectCSPG()  . . . . . . . .  non perfect case of composition series
##
##  When <workgroup> is not perfect, it fills up the factor group of the
##  commutator subgroup with cyclic factors.
##  Output is the first index in normals which remains undefined
##
InstallGlobalFunction( NonPerfectCSPG,
    function( homlist, normals, factors, auxiliary,
                            factorsize, top, index, D, workgroup )
    local   listlength,   # number of cyclic factors to add to factors
            indexup,      # loop variable for adding the cyclic factors
            oldworkup,    # loop subgroups between
            workup,       # workgroup and derived subgrp
            order,        # index of oldworkup in workup
            orderlist,    # prime factors of order
            g, p,         # generators of workup, oldworkup
            h,            # a power of g
            i;         # loop variables

    # number of primes in factor <workgroup> / <derived subgroup>
    listlength := Length(Factors(Integers,top));
    indexup := index+listlength;
    oldworkup := D;

    # starting with the derived subgroup, add generators g of workgroup
    # each addition produces a cyclic factor group on top of previous;
    # appropriate powers of g will divide the cyclic factor group into
    # factors of prime length
    for g in StabChainMutable( workgroup ).generators do
        if not (g in oldworkup)  then
            # check for error in random computation of derived subgroup
            Assert(1, ForAll ( StabChainMutable( oldworkup ).generators,
                              x->(x^g in oldworkup) ));
            workup := ClosureGroup(oldworkup, g);
            order := Size(workup)/Size(oldworkup);
            orderlist := Factors(Integers,order);
            for i in [1..Length(orderlist)] do

                # h is the power of g which adds prime length factors
                h := g^Product([i+1..Length(orderlist)],x->orderlist[x]);

                # construct entries in factors, normals
                factors[indexup -1] := [];
                normals[indexup -1] := [];
                for p in StabChainMutable( oldworkup ).generators do

                    # p acts trivially in factor group
                    Add(factors[indexup -1],());

                    # preimage of p is a generator in normals
                    Add(normals[indexup -1],PullbackCSPG(p,homlist));

                od;

                # workgroup is a factor group of original input;
                # kernel of homomorphism must be added to gens in normals
                PullbackKernelCSPG(homlist,normals,factors,
                                   auxiliary,indexup-1);

                # add preimage of h to generator list
                Add(normals[indexup-1],PullbackCSPG(h,homlist));

                # add a prime length cycle to factor group action
                Add(factors[indexup-1],
                           PermList(Concatenation([2..orderlist[i]],[1])));
                # size of factor group is a prime

                factorsize[indexup-1] := orderlist[i];
                indexup := indexup -1;

            od;
            oldworkup := workup;
        fi;
    od;

    return index+listlength;
end );


#############################################################################
##
#F  PerfectCSPG() . . . . . . . . . . . .  prefect case of composition series
##
##  Computes maximal normal subgroup of perfect primitive group K and adds
##  its factor group to factors.
##  Output is the maximal normal subgroup NN. In case NN=1 (i.e. K simple),
##  the kernel of homomorphism which produced K is returned
##
InstallGlobalFunction( PerfectCSPG,
    function( homlist, normals, factors, auxiliary,
                         factorsize, index, K )
    local   whichcase,  # var indicating to which case of the O'Nan-Scott
                        # theorem K belongs. When Size(K) and degree do not
                        # determine the case without ambiguity, whichcase
                        # has value as in case of unique nonregular
                        # minimal normal subgroup
            N,          # normal subgroup of K
            prime,      # prime dividing order of degree of K
            stab2,      # stabilizer of first two base points in K
            kerelement, # element of normal subgroup
            ker2,       # conjugate of kerelement
            H,          # normalizer, and then centralizer of stab2
            L,          # set of moved points of stab2
            op,         # operation of H on N
            tchom,      # restriction of H to L
            g,          # generator of subgroups
            lenhomlist, # length of homlist
            kernel,     # output
            ready,      # boolean variable indicating whether normal subgroup
                        # was found
            chainK,
            list;

    while not IsSimpleGroup(K)  do
        whichcase := CasesCSPG(K);

        # becomes true if we find proper normal subgroup by first method
        ready := false;

        # whichcase[1] is true in nonregular minimal normal subgroup case
        if whichcase[1]=1  then
            N := FindNormalCSPG(K, whichcase);

            # check size of result to terminate fast in ambiguous cases
            if 1 < Size(N)  and Size(N) < Size(K)  then
                # K is a factor group with N in the kernel
                K := NinKernelCSPG(K,N,homlist,auxiliary);
                SetDerivedSubgroup( K, K );
#T better set that K is perfect?
                ready := true;
            fi;

       fi;

        # apply regular normal subgroup with nontrivial centralizer method
        if not ready then
            chainK:= StabChainMutable( K );
            stab2 := Stabilizer(K,[ chainK.orbit[1],
                                    chainK.stabilizer.orbit[1]],
                             OnTuples);
            if IsTrivial(stab2) then

               prime := Factors(Integers,whichcase[2])[1];
               N:=Group(One(K));
               repeat
                 kerelement:=Random(K);
                 if NrMovedPoints(kerelement)=LargestMovedPoint(K) and
                     IsOne(kerelement^prime) then
                     ker2:=kerelement^Random(K);
                     if Comm(kerelement,ker2)=One(K) then
                        N := NormalClosure(K, SubgroupNC(K,[kerelement]));
                     fi;
                 fi;
               until Size(N)=whichcase[2];

            else
               list := NormalizerStabCSPG(K);
               H := list[1];
               chainK := list[2];
               if whichcase[1] = 2 then
                  stab2 := Stabilizer( K, [ chainK.orbit[1],
                                           chainK.stabilizer.orbit[1] ],
                                  OnTuples);
                  H := CentralizerNormalCSPG( H, stab2 );
               else
                  L := Orbit( H, StabChainMutable( H ).orbit[1] );
                  tchom := ActionHomomorphism(H,L,"surjective");
                  op := Image( tchom );
                  H := PreImage(tchom,PCore(op,Factors(Integers,whichcase[2])[1]));
                  H := Centre(H);
                  SetIsAbelian( H, true );
              fi;
              N := FindRegularNormalCSPG(K,H,whichcase);
           fi;
           K := NinKernelCSPG(K,N,homlist,auxiliary);
           SetDerivedSubgroup( K, K );
#T better set that K is perfect?
        fi;

    od;

    # add next entry to the CompositionSeries output lists
    factors[index] := [];
    normals[index] := [];
    factorsize[index] := Size(K);
    for g in StabChainMutable( K ).generators do
        Add(factors[index],g); # store generators for image
        Add(normals[index],PullbackCSPG(g,homlist));
    od;

    # add generators for kernel to normals
    PullbackKernelCSPG(homlist,normals,factors,auxiliary,index);
    lenhomlist := Length(homlist);

    # determine output of routine
    if lenhomlist > 0  then
        kernel := KernelOfMultiplicativeGeneralMapping(homlist[lenhomlist]);
        if IsBound(auxiliary[lenhomlist])  then
            kernel := auxiliary[lenhomlist]; # faster to add this way
            kernel := ClosureGroup( kernel,
                KernelOfMultiplicativeGeneralMapping(homlist[lenhomlist]) );
        fi;
        Unbind(homlist[lenhomlist]);
        Unbind(auxiliary[lenhomlist]);

    # case when we found last factor of original group
    else
        kernel := GroupByGenerators( [], () );
    fi;

    return kernel;
end );


#############################################################################
##
#F  CasesCSPG() . . . . . . . . . . . . determine case of O'Nan Scott theorem
##
##  Input: primitive, perfect, nonsimple group G.
##  CasesCSPG determines whether there is a normal subgroup with
##  nontrivial centralizer (output[1] := 2 or 3) or decomposes the
##  degree of G into the form output[2]^output[3], output[1] := 1 (case
##  of nonregular minimal normal subgroup).
##  There are some ambiguous cases, (e.g. degree=2^15) when Size(G)
##  and degree do not determine which case G belongs to. In these cases,
##  the output is as in case of nonregular minimal normal subgroup.
##  This computation duplicates some of what is done in IsSimple.
##
InstallGlobalFunction( CasesCSPG, function(G)
    local   degree,     # degree of G
            g,          # order of G
            primes,     # list of primes in prime decomposition of degree
            output,     # output of routine
            n,m,o,p,  # loop variables
            tab1,       # table of orders of primitive groups
            tab2,       # table of orders of perfect transitive groups
            base;       # prime occurring in order of outer automorphism
                        # group of some group in tab1

    g := Size(G);
    degree := LargestMovedPoint(G);
    if degree>2^20 then
      # see comment before the composition series method
      Error("degree too big");
    fi;

    output := [];

    # case of two regular normal subgroups
    if Size(G)=degree^2  then
        output[1] := 2;
        output[2] := degree;
        return output;
    fi;

    # degree is not prime power
    primes := Factors(Integers,degree);
    if primes[1] < Last(primes) then
        output[1] := 1;
        # only case when index of primitive group in socle is not 2*prime
        if Length(primes)=15  then
            output[2] := 12;
            output[3] := 5;
        else
            output[2] := primes[1]*Last(primes);
            output[3] := Length(primes)/2;
        fi;
        return output;

    # in case of prime power degree, we have to determine the possible
    # orders of G with nonabelian socle. See IsSimple for identification
    # of groups in tab1,tab2
    else
        tab1 := [ ,,,,[60],,[168,2520],[168,20160],[504,181440],,
                  [660,7920,19958400],,[5616,3113510400]];
        tab2 := [ ,,,,[60],[60,360],[168,2520],[168,1344,20160]];
        for n in [5,7,8,9,11,13] do
            for m in [5..8] do
                for o in [1..Length(tab1[n])] do
                    for p in [1..Length(tab2[m])] do
                        if tab1[n][o]=504  then
                            base := 3;
                        else
                            base := 2;
                        fi;
                        if degree=n^m
                          and g mod (tab1[n][o]^m*tab2[m][p]) = 0
                          and (tab1[n][o]^m*tab2[m][p]*base^m) mod g = 0
                        then
                            output[1] := 1;
                            output[2] := n;
                            output[3] := m;
                            return output;
                        fi;
                    od;
                od;
            od;
        od;

        # if the order of G did not satisfy any of the nonabelian socle
        # possibilities, output the abelian socle message
        output[1] := 3;
        output[2] := degree;
        return output;

    fi;

end );


#############################################################################
##
#F  FindNormalCSPG()  . . . . . . . . . . . . . find a proper normal subgroup
##
##  given perfect, primitive G with unique nonregular minimal normal
##  subgroup, the routine returns a proper normal subgroup of G
##
InstallGlobalFunction( FindNormalCSPG, function ( G, whichcase )
    local   n,          # degree of G
            i,          # loop variable
            stabgroup,  # stabilizer subgroup of first point
            orbits,     # list of orbits of stabgroup
            where,      # index of shortest orbit in orbits
            len,        # length of shortest orbit
            tchom,      # trans. constituent homom. of stabgroup
                        # to shortest orbit
            bl,         # blocks in action of stabgroup on shortest orbit
            bhom,       # block homomorphism for the action on bl
            K,          # homomorph image of stabgroup at tchom, bhom
            kernel,     # kernel of bhom
            N;          # output; normal subgroup of G

    # whichcase[1]=1 if G has no normal subgroup with nontrivial
    # centralizer or we cannot determine this fact from Size(G)
    n := LargestMovedPoint(G);
    stabgroup := Stabilizer(G, StabChainMutable( G ).orbit[1],OnPoints);
    orbits := OrbitsDomain(stabgroup,[1..n]);

    # find shortest orbit of stabgroup
    len := n; where := 1;
    for i in [1..Length(orbits)] do
        if (1<Length(orbits[i])) and (Length(orbits[i])< len)  then
            where := i;
            len := Length(orbits[i]);
        fi;
    od;

    # check arith. conditions in order to terminate fast in ambiguous cases
    if len mod whichcase[3] = 0 and len <= whichcase[3]*(whichcase[2]-1) then

        # take action of stabgroup on shortest orbit
        tchom := ActionHomomorphism(stabgroup,orbits[where],"surjective");
        K := Image(tchom,stabgroup);
        bl := MaximalBlocks(K,[1..len]);

        # take action on blocks
        if Length(bl) > 1  then
            bhom := ActionHomomorphism(K,bl,OnSets,"surjective");
            K := Image(bhom,K);
            kernel := KernelOfMultiplicativeGeneralMapping(
                          CompositionMapping(bhom,tchom));
            N := NormalClosure(G,kernel);

            # another check for ambiguous cases
            if Size(N) < Size(G) then
                return N;
            fi;

        fi;
    fi;

    # in ambiguous case, return trivial subgroup
    N := TrivialSubgroup( Parent(G) );
    return N;
end );


#############################################################################
##
#F  FindRegularNormalCSPG()  . . . . . . . . . . find a proper normal subgroup
##
##  given perfect, primitive G with regular minimal normal
##  subgroup(s), the routine returns one
##
InstallGlobalFunction( FindRegularNormalCSPG, function ( G, H, whichcase )

    local core,         # p-core of H
          cosetrep,     # a cosetrep of H.stabilizer
          candidates,   # list of perms; one element is in regular normal sbgrp
          ready,        # boolean to exit loop
          i,            # loop variable
          N,            # regular normal subgroup, output
          chain;

    # case of abelian normal subgroup
    if whichcase[1] <> 2 then
       core := PCore( H, Factors(Integers, whichcase[2])[1] );
       chain:=StabChainOp(core,rec(base:=BaseOfGroup(G),reduced:=false));
       cosetrep := chain.transversal[chain.orbit[2]];
       candidates := AsList(Stabilizer(core,BaseOfGroup(G)[1]))*cosetrep;
       ready := false;
       i:= 0;
       while not ready do
          i := i+1;
          N := NormalClosure(G, SubgroupNC(G, [candidates[i]]) );
          if Size(N) = whichcase[2] then
             ready := true;
          fi;
       od;

     # case of two simple regular normal subgroups
     else
       chain := StabChainOp(H, rec(base := BaseOfGroup(G), reduced := false) );
       cosetrep := chain.transversal[chain.orbit[2]];
       candidates := cosetrep*AsList(Stabilizer(H,BaseOfGroup(G)[1]));
       ready := false;
       i:= 0;
       while not ready do
          i := i+1;
          N := NormalClosure(G, SubgroupNC(G, [candidates[i]]) );
          if Size(N) = whichcase[2] then
             ready := true;
          fi;
       od;
     fi;

     return N;
end );

#############################################################################
##
#F  NinKernelCSPG() . . . . . find homomorphism that contains N in the kernel
##
##  Given a normal subgroup N of G, creates a subgroup H such that the
##  homomorphism to the action on cosets of H contains N in the kernel.
##  Actually, only the image of a subgroup is computed, and we store
##  N in auxiliary to remember that N should be added to kernel of
##  homomorphism.
##  Output is the image at homomorphism
##
InstallGlobalFunction( NinKernelCSPG,
    function ( G, N, homlist, auxiliary )
    local   i,j,        # loop variables
            base,       # base of G
            stab,       # stabilizer of first two base points
            H,HOld,     # subgroups of G
            G1,H1,      # stabilizer chains of G, HOld
            block,      # set of cosets of G1[i]; G1[i] is represented on
                        # images of block
            newrep,     # blocks of imprimitivity in
            bhom,       # block hom. and
            tchom;      # transisitive const. hom. applied to G1[i]

    j := Length(homlist)+1;
    auxiliary[j] := N;

    # find smallest subgroup of G in stabilizer chain which, together with N,
    # generates G
    G1:=StabChainMutable(G);
    base := BaseStabChain(G1);
    G1 := ListStabChain( G1 );
    i := Length(base)+1;
    # first try the stabilizer of first two points
    if Size(N) = LargestMovedPoint(G) then
       stab := AsSubgroup(Parent(G),Stabilizer(G,[base[1],base[2]],OnTuples));
       H := ClosureGroup
             ( stab, GeneratorsOfGroup( N ), rec(size:=Size(N)*Size(stab)) );
    else
       H := ClosureGroup( N, G1[3].generators );
    fi;
    if Size(H) < Size(G) then
       HOld := H;
       i := 2;
    else
    # if did not work, start from bottom of stabilizer chain
       H := N;
       repeat
           HOld := H;
           i := i-1;
           H := ClosureGroup( H, G1[i].generators );
       until Size(H) = Size(G);
    fi;

    # represent G1[i] on cosets of H1[i] := G1[i+1]N \cap G1[i]
    H1 := ListStabChain( StabChainOp( HOld, rec( base := base,
                                              reduced := false ) ) );

    # G1[i] will be represented on images of block
    block := Set( H1[i].orbit );
    G := Stabilizer(G,List([1 .. i-1], x->base[x]),OnTuples);

    # now G is the previous G1[i]
    # find primitive action on images of block
    newrep := MaximalBlocks( G, StabChainMutable( G ).orbit, block );
    if Length(newrep) > 1 then
        bhom := ActionHomomorphism(G,newrep,OnSets,"surjective");
        Add(homlist,bhom);
        G := Image(bhom,G);
    else
        tchom:=ActionHomomorphism(G, StabChainMutable( G ).orbit,"surjective");
        Add(homlist,tchom);
        G := Image(tchom,G);
    fi;

    return G;
end );


#############################################################################
##
#F  RegularNinKernelCSPG()  . . . .  action of G on H and on maximal subgroup
##
##  H is transitive and contains the stabilizer of the first two
##  base points; we want to find the action of G on cosets of H, and
##  then the action of G on cosets of a maximal subgroup K containing H
##  reference: Beals-Seress Lemma 4.3.
##
InstallGlobalFunction( RegularNinKernelCSPG,
    function ( G, H, homlist )
    local   i,j,k,      # loop variables
            base,       # base of G
            chainG,     # stabilizer chain of `G'
            chainH,     # stabilizer chain of `H'
            H1,         # stabilizer chain of G,H
            x,y,        # first two base points of G
            stabgroup,  # stabilizer of x in G
            chainstabgroup,
            Ginverses,  # list of inverses of generators of G
            hgens,      # list of generators of H
            Hinverses,  # list of inverses of generators of H
            stabgens,   # list of generators of stabgroup
            stabinverses, # list of inverses of generators of stabgroup
            block,      # orbit of y in H_x
            orbits,     # images of block in G_x=stabgroup
            a,          # cardinality of orbits
            b,          # cardinality of block
            reprlist,   # for z in stabgroup.orbit, reprlist[z] tells which
                        # element of orbits z belongs to
            reps,       # for representatives z of sets in orbits, reps
                        # contains the cosetrep carrying z to y in stabgroup
                        # (as a word in generators of stabgroup)
            inversereps,# the inverses of words in reps
                        # (as words in stabinverses)
            images,     # list containing the images of generators of G,
                        # acting on cosets of H (there are $a$ cosets,
                        # represented by the elements of orbits)
            v,          # point of permutation domain
            tau,        # the cosetrep of H carrying v to x
                        # (as a word in H.gen's)
            tauinverse, # the inverse of tau (as a word in Hinverses)
            word,       # list of permutations coding a cosetrep of H
            K,          # the factor group of G generated by images
            newrep,     # block system from cosets of K
            c,          # cardinality of newrep
            d,          # size of one block
            newimages,  # list containing the action of generators of G, on
                        # newrep
            hom;        # the homomorphism G->K

    chainG:= StabChainMutable( G );
    base := BaseStabChain(chainG);
    H1 := ListStabChain( StabChainOp( H, rec( base := base,
                                           reduced := false ) ) );
    block := Set( H1[2].orbit );
    x := chainG.orbit[1];
    stabgroup := Stabilizer( G, x, OnPoints );
    orbits := Orbit(stabgroup,block,OnSets);
    chainstabgroup:= StabChainMutable( stabgroup );
    y := chainstabgroup.orbit[1];
    a := Length(orbits);
    b := Length(block);
    reprlist := [];
    for i in [1..a] do
        for k in [1..b] do
            reprlist[orbits[i][k]] := i;
        od;
    od;

    Ginverses := GInverses( chainG );
    chainH:= StabChainMutable( H );
    Hinverses := GInverses( chainH );
    hgens := chainH.generators;

    stabinverses := GInverses( chainstabgroup );
    stabgens := chainstabgroup.generators;

    reps := []; inversereps := [];
    for i in [1..a] do
        reps[i] := CosetRepAsWord( y, orbits[i][1],
                                   chainstabgroup.transversal );
        inversereps[i] := InverseAsWord(reps[i],stabgens,stabinverses);
    od;

    # construct action of G-generators on cosets of H. Each coset of H has a
    # representative in orbits; to find the image of an H coset
    # at multiplication by G.generators[i], take element of H coset such that
    # the product with G.generators[i] fixes x. Then the image of the coset
    # can be read from the position in orbits (cf. Lemma 4.3)
    images := [];
    for i in [1..Length( chainG.generators )] do
        images[i] := [];
        for j in [1..a] do
            v := ImageInWord(x^Ginverses[i],reps[j]);
            tau := CosetRepAsWord( x, v, chainH.transversal );
            tauinverse := InverseAsWord(tau,hgens,Hinverses);
            word := Concatenation(tauinverse,inversereps[j],
                                  [ chainG.generators[i] ]);
            images[i][j] := reprlist[ImageInWord(y,word)];
        od;
        images[i] := PermList(images[i]);
    od;
    K := GroupByGenerators(images,());

    # check whether new representation is primitive. If not, construct action
    # on maximal block system
    newrep := MaximalBlocks(K,[1..a]);
    if Length(newrep) > 1  then
        c := Length(newrep);
        d := Length(newrep[1]);
        reprlist := [];
        for i in [1..c] do
            for k in [1..d] do
                reprlist[newrep[i][k]] := i;
            od;
        od;
        newimages := [];
        for i in [1..Length( chainG.generators )] do
            newimages[i] := [];
            for k in [1..c] do
                newimages[i][k] := reprlist[newrep[k][1]^images[i]];
            od;
            newimages[i] := PermList(newimages[i]);
        od;
        K := GroupByGenerators(newimages,());
        hom := GroupHomomorphismByImagesNC( G, K,
                   chainG.generators, newimages );
    else
        hom := GroupHomomorphismByImagesNC( G, K,
                   chainG.generators, images );
    fi;
    j := Length(homlist)+1;
    homlist[j] := hom;
    K := Image(homlist[j],G);
    SetDerivedSubgroup( K, K );
#T better set that K is perfect?

    return K;
end );


#############################################################################
##
#F  NormalizerStabCSPG( <G> ) . . . . . . .  normalizer of 2 point stabilizer
##
##  Given a primitive, perfect group <G> which has a regular normal subgroup
##  with nontrivial centralizer,
##  the output is a list of length two, the first entry being N_G(G_{xy})
##  and the second entry being a stabilizer chain of <G>.
##
InstallGlobalFunction( NormalizerStabCSPG, function(G)
    local   n,          # degree of G
            chainG,     # stabilizer chain of `G'
            chainstab,  # stabilizer chain of a point stabilizer in `G'
            orbits,     # orbits of stabgroup
            len,        # minimal length of stabgroup orbits
            where,      # index of minimal length orbit
            i,          # loop variable
            chainstab2, # chain of stabilizer of first two base points in G
            x,y,        # first two base points
            normalizer, # output group N_G(G_{xy})
            L,          # fixed points of stabgroup2
            yL,         # intersection of L and y-orbit in stabgroup
            orbity,     # orbit of y in normalizer_x;
                        # eventually, orbity must be yL
            orbitx,     # orbit of x in normalizer;
                        # eventually, orbitx must be L
            u,v,        # points in permutation domain
            tau,sigma,p,# cosetreps of G, stabgroup
            Ltau;       # image of L under tau

    n := LargestMovedPoint(G);
    chainG:= StabChainMutable( G );
    chainstab := chainG.stabilizer;

    # If necessary, make base change to achieve that second base point is
    # in smallest orbit of stabilizer.
    orbits := OrbitsPerms( chainstab.generators, [1..n] );
    len := n; where := 1;
    for i in [1..Length(orbits)] do
        if (1<Length(orbits[i])) and (Length(orbits[i])< len)  then
            where := i;
            len := Length(orbits[i]);
        fi;
    od;
    if Length( chainstab.orbit ) > len  then
      chainG:= StabChainOp( G, [ chainG.orbit[1], orbits[where][1] ] );
      chainstab:= chainG.stabilizer;
    fi;
    x := chainG.orbit[1];
    y := chainstab.orbit[1];
    chainstab2 := chainstab.stabilizer;

    # compute normalizer. Method: Beals-Seress, Lemma 7.1
    L := Difference( [1..n], MovedPoints( chainstab2.generators ) );
    yL := Intersection( L, chainstab.orbit );

    # initialize normalizer to G_{xy}
    normalizer := rec( generators := ShallowCopy( chainstab2.generators) );
    orbity := OrbitPerms(normalizer.generators,y);
    while Length(orbity) < Length(yL) do
        v := Difference(yL,orbity)[1];
        p := Product( CosetRepAsWord( y, v, chainstab.transversal ) );
        Add(normalizer.generators,p);
        orbity := OrbitPerms(normalizer.generators,y);
    od;
    normalizer.stabChain2 := EmptyStabChain( [  ], (), y );
    AddGeneratorsExtendSchreierTree(normalizer.stabChain2,normalizer.generators);
    normalizer.stabChain2.stabilizer:= chainstab2;

    orbitx := OrbitPerms(normalizer.generators,x);
    while Length(orbitx) < Length(L) do
        v := Difference(L,orbitx)[1];
        tau := Product( CosetRepAsWord( x, v, chainG.transversal ) );
        Ltau := OnSets(L,tau);
        u := Intersection( Ltau, chainstab.orbit )[1];
        sigma := Product( CosetRepAsWord( y, u, chainstab.transversal ) );
        Add(normalizer.generators,tau*sigma);
        orbitx := OrbitPerms(normalizer.generators,x);
    od;
    normalizer.stabChain := EmptyStabChain( [  ], (), x );
    AddGeneratorsExtendSchreierTree(normalizer.stabChain,normalizer.generators);
    normalizer.stabChain.stabilizer:=normalizer.stabChain2;

    normalizer := GroupStabChain( Parent( G ), normalizer.stabChain, true );
    return [normalizer, chainG];
end );


#############################################################################
##
#F  TransStabCSPG() . . . embed a 2 point stabilizer in a transitive subgroup
##
##  given a subgroup H of G which contains G_{xy}, the stabilizer of the
##  first two points in G, and a theoretical guarantee that there is a
##  proper transitive subgroup K containing H, the routine finds such K
##
InstallGlobalFunction( TransStabCSPG, function(G,H)
    local   n,          # degree of G
            chainG,     # stabilizer chain of `G'
            chainH,     # stabilizer chain of `H'
            x,y,        # first two points of the base of G
            stabgroup,  # stabilizer of x in G
            chainstabgroup,
            hstabgroup, # stabilizer of x in H
            chainhstabgroup,
            u,v,        # indices of points in G.orbit, stabgroup.orbit
            g,          # list of permutations whose product is
                        # (semi)random element of G
            notinH,     # boolean; true if g is not in H
            word,       # list of permutations whose product is
                        # (semi)random element of <H,g>
            len,        # length of word
            hword,      # list of permutations giving random element of H
            tau,sigma,  # lists of permutations whose
            tau1,sigma1,# products are coset representatives
            i,j,k,      # loop variables
            K;          # K=<H,g>

    #Print(Size(G),",",Size(H));
    n := LargestMovedPoint(G);
    chainG:= StabChainMutable( G );
    x := chainG.orbit[1];
    stabgroup := Stabilizer(G,x,OnPoints);
    chainstabgroup := StabChainMutable( stabgroup );
    y := chainstabgroup.orbit[1];
    hstabgroup := Stabilizer(H,x,OnPoints);
    chainhstabgroup:= StabChainMutable( hstabgroup );
    chainH:= StabChainMutable( H );
    ExtendStabChain( chainH, BaseStabChain(chainG) );
    ExtendStabChain( chainhstabgroup, BaseStabChain( chainstabgroup ) );

    # try to embed H into bigger subgroups; stop when result is transitive
    repeat
        # Print("brum");

        # first, take random element of G\H
        repeat
            v := Random(1, Length( chainG.orbit ));
            g := CosetRepAsWord( x, chainG.orbit[v], chainG.transversal );
            u := Random(1, Length( chainstabgroup.orbit ));
            Append(g,CosetRepAsWord( y, chainstabgroup.orbit[u],
                                        chainstabgroup.transversal ));
            notinH := false;
            v := ImageInWord(x,g);
            if not IsBound( chainH.transversal[v] ) then
                notinH := true;
            else
                u := ImageInWord(y,g);
                u := ImageInWord( u, CosetRepAsWord( x, v,
                                         chainH.transversal ) );
                if not IsBound( chainhstabgroup.transversal[u] ) then
                    notinH := true;
                fi;
            fi;
        until  notinH;

        for i in [1..n] do

            # construct semirandom element of <H,g>
            word := [];
            for j in [1..5] do
                len := Length(word);
                for k in [1..Length(g)] do
                    word[len+k] := g[k];
                od;
                len := Length(word);
                hword := RandomElmAsWord(H);
                for k in [1..Length(hword)] do
                    word[len+k] := hword[k];
                od;
            od;

            # check whether word is in H;
            # if not, then let g=cosetrep of word in G_{xy}
            v := ImageInWord(x,word);
            tau := CosetRepAsWord( x, v, chainH.transversal );
            if tau = []  then
                tau1 := CosetRepAsWord( x, v, chainG.transversal );
                u := ImageInWord(y,word);
                u := ImageInWord(u,tau1);
                sigma1 := CosetRepAsWord( y, u, chainstabgroup.transversal );
                g := Concatenation(tau1,sigma1);
            else
                u := ImageInWord(y,word);
                u := ImageInWord(u,tau);
                sigma := CosetRepAsWord( y, u, chainhstabgroup.transversal );
                if sigma = []  then
                    tau1 := CosetRepAsWord( x, v, chainG.transversal );
                    u := ImageInWord(y,word);
                    u := ImageInWord(u,tau1);
                    sigma1 := CosetRepAsWord( y, u,
                                  chainstabgroup.transversal );
                    g := Concatenation(tau1,sigma1);
                fi;
            fi;
        od;

        # check whether H,g generate a proper subgroup of G
        K := ClosureGroup(H,Product(g));
        if 1 < Size(G)/Size(K)  then
            H := K;
            chainH:= StabChainMutable( H );
            #Print(Size(H));
            hstabgroup := Stabilizer(H,x,OnPoints);
            chainhstabgroup:= StabChainMutable( hstabgroup );
            ExtendStabChain( chainhstabgroup, BaseStabChain(chainstabgroup) );
            ExtendStabChain( chainH, BaseStabChain(chainG) );
        fi;

    until Length( chainH.orbit ) = n;

    return H;
end );


#############################################################################
##
#F  PullbackKernelCSPG()  . . . . . . . . . . . . . . . pull back the kernels
##
InstallGlobalFunction( PullbackKernelCSPG,
    function( homlist, normals, factors, auxiliary, index )
    local   lenhomlist, # length of homlist
            i, j,       # loop variables
            gens,       # list of generators in kernels
                        # of homomorphisms in homlist
            k,          # kernel
            kg,         # kernel generators
            g;          # a member of gens

    # for each kernel, compute preimages of the kernel generators in the
    # input group add these to generators of the current subnormal subgroup
    # in the composition series
    lenhomlist := Length(homlist);
    for i in [1..lenhomlist] do
       k:=KernelOfMultiplicativeGeneralMapping(homlist[i]);
       kg:=GeneratorsOfGroup(k);
       if IsBound(auxiliary[i])  then
           gens := Union( GeneratorsOfGroup( k ),
                         StabChainMutable( auxiliary[i] ).generators);
           if Length(gens)>6 then
             g:=Group(gens,());
             if IsSubset(auxiliary[i],k) then
               SetSize(g,Size(auxiliary[i]));
             else
               StabChainOptions(g).limit:=Size(k)*Size(auxiliary[i]);
             fi;
             gens:=SmallGeneratingSet(g);
           fi;
       else
         if Length(kg)>5 then
           gens:=SmallGeneratingSet(k);
         else
           gens := kg;
         fi;
       fi;
       for g in gens do
           for j in [1..i-1] do
               g := PreImagesRepresentative(homlist[i-j],g);
           od;
           Add(normals[index],g);
           Add(factors[index],());
       od;
    od;
end );


#############################################################################
##
#F  PullbackCSPG()  . . . . . . . . . . . . . . . . . . . . . . . . pull back
##
InstallGlobalFunction( PullbackCSPG, function(p,homlist)
    local   i,          # loop variable
            lenhomlist; # length of homlist

    # compute a preimage of the permutation p in the input group
    lenhomlist := Length(homlist);
    for i in [1..lenhomlist] do
        p := PreImagesRepresentative(homlist[lenhomlist+1-i],p);
    od;
    return p;
end );


#############################################################################
##
#F  CosetRepAsWord()  . . . . . . . . .  write a coset representative as word
##
##  returns the cosetrep carrying y to the base point x as a word in the
##  generators. If y is not in the orbit of x, returns []
##
InstallGlobalFunction( CosetRepAsWord, function(x,y,transversal)
    local   word,       # list of permutations
            point;      # element of permutation domain

    word := [];
    if IsBound(transversal[y])  then
        point := y;
        repeat
            word[Length(word)+1] := transversal[point];
            point := point^transversal[point];
        until point = x;
    fi;
    return word;
end );


#############################################################################
##
#F  ImageInWord() . . .  image of a point under a permutation written as word
##
##  computes the image of x when the list of permutations word is applied
##
InstallGlobalFunction( ImageInWord, function(x,word)
    local   i,          # loop variable
            value;      # element of permutation domain

    value := x;
    for i in [1..Length(word)] do
        value := value^word[i];
    od;
    return value;
end );


#############################################################################
##
#F  SiftAsWord( <chain>, <perm> ) . . . .  sift a permutation written as word
##
##  given a list <perm> of permutations and a stabilizer chain <chain> for
##  the group $G$, the routine computes the residue at the sifting of perm
##  through the SGS of $G$.
##  The output is a list of length 2: the first component is the siftee,
##  as a word, the second component is 0 if perm in $G$, and i if the siftee
##  on the i^th level could not be computed.
##
#T <perm> is changed!
##
InstallGlobalFunction( SiftAsWord, function( chain, perm )
    local   i,          # loop variable
            y,          # element of permutation domain
            word,       # the list collecting the siftee of perm
            len,        # length of word
            coset,      # word representing a coset in a stabilizer
            index,      # the level where the siftee cannot be computed
            stb;        # the stabilizer group we currently work with

    # perm must be a list of permutations itself!
    stb :=  chain;
    word := perm;
    index := 0;
    while IsBound(stb.stabilizer) do
       index:=index+1;
       y:=ImageInWord(stb.orbit[1],word);
       if IsBound(stb.transversal[y]) then
          coset :=  CosetRepAsWord(stb.orbit[1],y,stb.transversal);
          len := Length(word);
          for i in [1..Length(coset)] do
              word[len+i] := coset[i];
          od;
          stb:=stb.stabilizer;
       else
          return([word,index]);
       fi;
    od;

    index := 0;
    return [word,index];
end );


#############################################################################
##
#F  InverseAsWord() . . . . . . . . . .  invert a permutation written as word
##
##  given a list of permutations "list", the inverses of these permutations
##  in inverselist, and a list of permutations "word" with elements from
##  list, returns the inverse of word as a list of inverses from inverselist
##
InstallGlobalFunction( InverseAsWord, function(word,list,inverselist)
    local   i,          # loop variable
            p,          # position
            inverse;    # the inverse of word

    if word = [ () ]  then
        return word;
    fi;
    inverse := [];
    for i in [1..Length(word)] do
      # identity tests are cheaper if the degree gets bigger.
      p:=PositionProperty(list,j->IsIdenticalObj(j,word[Length(word)+1-i]));
      if p=fail then
        # this is very unlikely to happen.
        p:=Position(list,word[Length(word)+1-i]);
      fi;
      inverse[i] := inverselist[p];
    od;
    return inverse;
end );


#############################################################################
##
#F  RandomElmAsWord( <chain> )  . . . . . . .  random element written as word
##
##  given an stabilizer chain <chain> for the group $G$, returns a uniformly
##  distributed random element of $G$,
##  as a word in the strong generators
##
InstallGlobalFunction( RandomElmAsWord, function( chain )
    local  i,       # loop variable
           word,    # the random element
           len,     # length of word
           stb,     # the stabilizer group we currently work with
           v,       # index of random element of stb.orbit
           coset;   # word representing a coset
    word:=[];
    stb:= chain;
    while IsBound(stb.stabilizer) do
       v := Random(1,Length(stb.orbit));
       coset := CosetRepAsWord(stb.orbit[1],stb.orbit[v],stb.transversal);
       len := Length(word);
       for i in [1..Length(coset)] do
           word[len+i] := coset[i];
       od;
       stb:=stb.stabilizer;
    od;
    return  word;

end );

#############################################################################
##
#M  PCore() . . . . . . . . . . . . . . . . . . p core of a permutation group
##
##  O_p(G), the p-core of G, is the maximal normal p-subgroup
##  Output of routine: the subgroup O_p(workgroup)
##  reference: Luks-Seress
##
InstallMethod( PCoreOp,
    "for a permutation group, and a positive integer",
    true,
    [ IsPermGroup, IsPosInt ], 0,
    function(workgroup,p)
    local   n,          # degree of workgroup
            G,          # a factor group of workgroup
            list,       # the record workgroup.compositionSeries
            normals,    # gens for the subgroups in the composition series
            factorsize, # the sizes of factor groups in composition series
            index,      # loop variable running through the indices of
                        # subgroups in the composition series
            primes,     # list of primes in the factorization of numbers
            ppart,      # p-part of Size(G)
            homlist,    # list of homomorphisms applied to workgroup
            lenhomlist, # length of homlist
            K, N,       # subnormal subgroups of G from composition series
            g,          # generator of K
            C,          # centralizer of N in K
            D,          # the p-part of C
            order,      # order of a generator of C
            H,          # first solvable, then also
                        # abelian normal p-subgroup of G
            series,     # the derived series of H; H becomes abelian when it
                        # is redefined as last nontrivial term of series
            actionlist, # record of G action on transitive
                        # constituent pieces of H
            Ggens,      # generators of stab. chain of `G'
            i, j,       # loop variables
            image,      # list of images of generators of G
                        # acting on pieces of H
            GG,         # the image of G at this action
            hom,        # the homomorphism from G to GG
            pgenlist;   # list of generators for the p-core

    # handle trivial cases
    if not IsPrime(p)  then
        return TrivialSubgroup(workgroup);
    fi;
    if IsTrivial(workgroup)  then
        return TrivialSubgroup(workgroup);
    fi;
    if Size(workgroup) mod p <> 0 then
       # p does not divide Size(workgroup)
       return TrivialSubgroup(workgroup);
    fi;

    #handle nilpotent case directly
    if IsNilpotentGroup( workgroup ) then
           # compute the p-part of generators of workgroup
           primes := Collected( Factors( Size(workgroup) ) );
           ppart := p^primes[PositionProperty( primes, x->x[1]=p )][2];
           pgenlist := [];
           for g in StabChainMutable( workgroup ).generators do
               Add( pgenlist, g^( Size(workgroup)/( ppart ) ) );
           od;
           D := SubgroupNC( workgroup, pgenlist );
           if ppart > 1 then
               SetIsPGroup( D, true );
               SetPrimePGroup( D, p );
               SetSylowSubgroup( workgroup, p, D );
               SetHallSubgroup( workgroup, [p], D );
           fi;
           return D;
    fi;

    n := LargestMovedPoint(workgroup);
    G := workgroup;
    list := CompositionSeries(G);
    # normals := Copy(list[1]);
    # factorsize := list[3];
    normals := List( [1..Length(list)-1],
                     i->ShallowCopy(StabChainMutable(list[i]).generators));
    factorsize := List([1..Length(list)-1],i->Size(list[i])/Size(list[i+1]));
    Add(normals, [()]);
    homlist := [];
    index := Length(factorsize);

    # try to find smallest subgroup in composition series with nontrivial
    # p-core. The normal closure of this p-core is a solvable normal
    # p-subgroup of G; taking commutator subgroups, find abelian normal
    # p-subgroup of G.
    # represent G acting on transitive constituent pieces of abelian normal
    # p-subgroup; kernel is abelian p-group. Take image at this action, and
    # repeat
    while index > 0 do
        if factorsize[index] <> p  then
            index := index-1;
        else
            N := SubgroupNC(Parent(G),normals[index+1]);

            # define K := SubGroup(Parent(G),normals[index]);
            # N has trivial p-core; check whether K has nontrivial one
            # K=N is possible when we work in homomorphic images of original
            if ForAll(normals[index], x -> x in N)  then
                index := index-1;
            else
                K := ClosureGroup( N,normals[index],
                                          rec( size:=p*Size(N) ) );
                C := CentralizerNormalCSPG(K,N);
                # O_p(K) is cyclic or trivial; it must show up in C
                # C is always abelian; check whether it has p-part
                D := [];
                C:= GeneratorsOfGroup( C );
                for i in [1..Length( C )] do
                    order := Order(C[i]);
                    if order mod p = 0  then
                        D[i] := C[i]^(order/p);
                    else
                        D[i] := ();
                    fi;
                od;

                # redefine C as the p-core of C
                C := SubgroupNC(Parent(K),D);
                if IsTrivial(C)  then
                    index := index-1;
                else
                    H := NormalClosure(G,C);
                    series := DerivedSeriesOfGroup(H);
                    H := series[Length(series)-1];

                    # at that moment, H is abelian normal in G
                    # define new action of G with H in the kernel
                    actionlist := ActionAbelianCSPG(H,n);

                    Ggens:= StabChainMutable( G ).generators;
                    image:= List( Ggens,
                                g -> ImageOnAbelianCSPG( g, actionlist ) );

                    # take homomorphic image of G
                    GG := GroupByGenerators(image,());
                    hom:=GroupHomomorphismByImagesNC(G,GG,Ggens,image);
                    Add(homlist,hom);
                    #force makemapping
                    SetSize(GG,Size(G)/Size(
                      KernelOfMultiplicativeGeneralMapping( hom )));
                    # find new action of subgroups in composition series
                    for i in [1..index] do
                        for j in [1..Length(normals[i])] do
                            normals[i][j] :=
#                                ImageOnAbelianCSPG(normals[i][j],actionlist);
                            Image(hom,normals[i][j]);
                        od;
                    od;

                    G := GG;
                    index := index-1;

                fi;         # IsTrivial(C)

            fi;             # K = N

        fi;                 # factorsize[index] <> p

    od;

    # create output;
    # the p-core is the kernel of homomorphisms applied to workgroup
    lenhomlist := Length(homlist);
    if lenhomlist = 0  then
        pgenlist := [()];
    else
        pgenlist := [];
        for i in [1..lenhomlist] do
            for g in GeneratorsOfGroup( KernelOfMultiplicativeGeneralMapping(
                                            homlist[i] ) ) do
                for j in [1..i-1] do
                    g := PreImagesRepresentative(homlist[i-j],g);
                od;
                Add(pgenlist,g);
            od;
        od;
    fi;
    D := SubgroupNC(workgroup,pgenlist);
    if not ForAll(pgenlist,IsOne) then
        SetIsPGroup( D, true );
        SetPrimePGroup( D, p );
    fi;
    return D;
end );


#############################################################################
##
#M  SolvableRadical( <G> )  . . . . . solvable radical of a permutation group
##
##  the radical is the maximal solvable normal subgroup
##  output of routine: the subgroup radical of workgroup
##  reference: Luks-Seress
##
InstallMethod( SolvableRadical,
    "for a permutation group",
    [ IsPermGroup ],
    function(workgroup)
    local   n,          # degree of workgroup
            G,          # a factor group of workgroup
            list,       # the record workgroup.compositionSeries
            normals,    # gens for the subgroups in the composition series
            factorsize, # the sizes of factor groups in composition series
            index,      # loop variable running through the indices of
                        # subgroups in the composition series
            primes,     # list of primes in the factorization of numbers
            homlist,    # list of homomorphisms applied to workgroup
            lenhomlist, # length of homlist
            K, N,       # subnormal subgroups of G from composition series
            g,          # generator of K
            C,          # centralizer of N in K
            H,          # first solvable,
                        # then also abelian normal subgroup of G
            series,     # the derived series of H; H becomes abelian when it
                        # is redefined as last nontrivial term of series
            actionlist, # record of G action on transitive
                        # constituent pieces of H
            Ggens,      # generators of stab. chain of `G'
            i, j,       # loop variables
            image,      # list of images of generators of G
                        # acting on pieces of H
            GG,         # the image of G at this action
            hom,        # the homomorphism from G to GG
            map,        # natural homomorphism for radical.
            solvable,   # list of generators for the radical
            o,          # orbits of G
            b,          # blocks
            TryReduction;# function to test whether a hom. can reduce

    if IsTrivial(workgroup)  then
        return TrivialSubgroup(workgroup);
    fi;

    if IsSolvableGroup(workgroup) then
        return workgroup;
    fi;

    n := LargestMovedPoint(workgroup);
    G := workgroup;

    # if the degree is big, try to reduce it in a first step
    if n>1000 then

      TryReduction:=function(hom)
      local s,f,k,map;
        s:=Size(G)/Size(Image(hom)); # kernel size
        # is the kernel solvable? If yes we can go to the image
        f:=Collected(Factors(s));
        # at most 2 primes or all primes to power 1 -> Solvable
        if Length(f)<3 or ForAll(f,i->i[2]=1) then
          Info(InfoGroup,1,"solvable kernel size ",f);
          # OK, transfer result back
          k:=SolvableRadical(Image(hom));
          solvable:=PreImage(hom,k);
          map:=hom*NaturalHomomorphismByNormalSubgroup(Image(hom),k);
          SetKernelOfMultiplicativeGeneralMapping(map,solvable);
          AddNaturalHomomorphismsPool(G,solvable,map);
          return solvable;
        fi;
        return fail;
      end;

      # try orbits
      o:=ShallowCopy(Orbits(G,MovedPoints(G)));
      if Length(o)>1 then
        SortBy(o, Length);
        for i in o do
          Info(InfoGroup,1,"trying orbit length ",Length(o));
          hom:=ActionHomomorphism(G,i,"surjective");
          K:=TryReduction(hom);
          if K<>fail then
            return K;
          fi;
        od;
      fi;
      # try blocks on orbits
      for i in o do
        b:=Blocks(G,i);
        if Length(b)>1 then
          Info(InfoGroup,1,"trying blocks length ",Length(b));
          hom:=ActionHomomorphism(G,b,OnSets,"surjective");
          K:=TryReduction(hom);
          if K<>fail then
            return K;
          fi;
        fi;
      od;
    fi;

    list := CompositionSeries(G);
    # normals := Copy(list[1]);
    # factorsize := list[3];

    #was:
    #normals := List( [1..Length(list)-1],
    #                 i->ShallowCopy(StabChainMutable(list[i]).generators));
    # but not all subgroups in the comp.ser have their own stabchain.
    normals:=[];
    for i in [1..Length(list)-1] do
      if HasStabChainMutable(list[i]) then
        normals[i]:=ShallowCopy(StabChainMutable(list[i]).generators);
      else
        normals[i]:=ShallowCopy(GeneratorsOfGroup(list[i]));
      fi;
    od;

    factorsize := List([1..Length(list)-1],i->Size(list[i])/Size(list[i+1]));
    Add(normals, [()]);
    homlist := [];
    index := Length(factorsize);

    # try to find smallest subgroup in composition series with nontrivial
    # radical. The normal closure of this radical is a solvable normal
    # subgroup of G; taking commutator subgroups, find abelian normal
    # subgroup of G.
    # represent G acting on transitive constituent pieces of abelian normal
    # subgroup; kernel is abelian normal.
    # Take image at this action, and repeat
    while index > 0 do
        primes := Factors(Integers,factorsize[index]);

        # if the factor group is not cyclic, no chance for nontrivial radical
        if Length(primes) > 1  then
            index := index-1;
        else
            N := SubgroupNC(Parent(G),normals[index+1]);

            # define K := SubGroup(Parent(G),normals[index]);
            # N has trivial radical; check whether K has nontrivial one
            # K=N is possible when we work in homomorphic images of original
            if ForAll(normals[index], x -> x in N)  then
                index := index-1;
            else
                K := ClosureGroup( N,normals[index],
                                     rec( size:=factorsize[index]*Size(N) ) );
                Size(K);
                C := CentralizerNormalCSPG(K,N);

                # radical of K is cyclic or trivial; it has to show up in C
                if IsTrivial(C)  then
                    index := index-1;
                else
                    H := NormalClosure(G,C);
                    series := DerivedSeriesOfGroup(H);
                    H := series[Length(series)-1];

                    # at that moment, H is abelian normal in G
                    # define new action of G with H in the kernel
                    actionlist := ActionAbelianCSPG(H,n);

                    Ggens:= StabChainMutable( G ).generators;
                    if Length(Ggens)>5*Length(GeneratorsOfGroup(G)) then
                      Ggens:=GeneratorsOfGroup(G);
                    fi;
                    image:= List( Ggens,
                                g -> ImageOnAbelianCSPG( g, actionlist ) );

                    # take homomorphic image of G
                    GG := GroupByGenerators(image,());
                    hom := GroupHomomorphismByImagesNC(G,GG,
                                                     Ggens,image);
                    Add(homlist,hom);
                    #force makemapping and deduce image order
                    SetSize(GG,Size(G)/Size(
                      KernelOfMultiplicativeGeneralMapping( hom )));
                    # find new action of subgroups in composition series
                    for i in [1..index] do
                        for j in [1..Length(normals[i])] do
                            normals[i][j] :=
#                                ImageOnAbelianCSPG(normals[i][j],actionlist);
                            Image(hom,normals[i][j]);
                        od;
                    od;
                    Unbind(actionlist); # big object that is not needed later
                    G := GG;
                    index := index-1;

                fi;         # IsTrivial(C)

            fi;             # K = N

        fi;                 # Length(primes)>1

    od;

    # create output;
    # the radical is the kernel of homomorphisms applied to workgroup
    lenhomlist := Length(homlist);
    if lenhomlist = 0  then
        return TrivialSubgroup(workgroup);
    else
        solvable := [];
        for i in [1..lenhomlist] do
            for g in GeneratorsOfGroup( KernelOfMultiplicativeGeneralMapping(
                                            homlist[i] ) ) do
                for j in [1..i-1] do
                    g := PreImagesRepresentative(homlist[i-j],g);
                od;
                Add(solvable,g);
            od;
        od;
    fi;

    # construct the natural hom.
    map:=[];
    for i in GeneratorsOfGroup(workgroup) do
      g:=i;
      for j in [1..lenhomlist] do
        g:=ImageElm(homlist[j],g);
      od;
      Add(map,g);
    od;

    solvable:=SubgroupNC(workgroup,solvable);
    g:=Group(map,());
    SetSize(g,Index(workgroup,solvable));
    SetSolvableRadical(g,TrivialSubgroup(g));
    map:=GroupHomomorphismByImagesNC(workgroup,g,
                                     GeneratorsOfGroup(workgroup),map);
    SetKernelOfMultiplicativeGeneralMapping(map,solvable);
    AddNaturalHomomorphismsPool(workgroup,solvable,map);
    return solvable;
end );


#############################################################################
##
#M  Centre( <G> ) . . . . . . . . . . . . . . . center of a permutation group
##
##  constructs the center of G.
##  Reference: Beals-Seress, 24th Symp. on Theory of Computing 1992, sect. 9
##
InstallMethod( Centre,
    "for a permutation group",
    [ IsPermGroup ],
    function(G)
    local   n,          # degree of G
            orbits,     # list of orbits of G
            base,       # lexicographically smallest (in list) base of G
            i,j,        # loop variables
            reps,       # array recording which orbit of G the points in
                        # perm. domain belong to
            domain,     # union of G orbits which contain base points
            significant,# indices of orbits in "orbits" that belong to domain
            max,        # loop variable, used at definition of significant
            len,        # length of domain
            tchom,      # trans. const. homom, restricting G to domain
            GG,         # the image of tchom
            chainGG,    # stabilizer chain of `GG'
            chainGGG,   # stabilizer chain of `GGG'
            orbit,      # an orbit of GG
            tchom2,     # trans. const. homom, restricting GG to orbit
            GGG,        # the image of GG at tchom2
            hgens,      # list of generators for the direct product of
                        # centralizers of GG in Sym(orbit), for orbits of GG
            order,      # order of `GroupByGenerators( hgens, () )'
            centr,      # the centralizer of GG in Sym(orbit)
            inverse2,   # inverse of the conjugating permutation of tchom2
            g,          # generator of centr
            cent;       # center of GG

    if IsTrivial(G)  then
       return TrivialSubgroup(G);
    fi;

    base := BaseStabChain(StabChainMutable(G));
    n := Maximum( Maximum( base ), LargestMovedPoint(G) );
    orbits := OrbitsDomain(G,[1..n]);
    # orbits := List( orbits, Set );

    # handle case of transitive G directly
    if Length(orbits) = 1  then
        centr := CentralizerTransSymmCSPG( G, StabChainMutable( G ) );
        if IsEmpty( GeneratorsOfGroup( centr ) ) then
           return TrivialSubgroup( G );
        else
           order := Size(centr);
           cent := IntersectionNormalClosurePermGroup(G,centr,order*Size(G));
           Assert( 1, IsAbelian( cent ) );
           SetIsAbelian( cent, true );
           return cent;
        fi;
    fi;

    # for intransitive G, find which orbit contains which
    # points of permutation domain
    reps := [];
    for i in [1..Length(orbits)] do
        for j in [1..Length(orbits[i])] do
            reps[orbits[i][j]] := i;
        od;
    od;

    # take union of significant orbits which contain base points
    max := reps[base[1]];
    significant := [max];
    domain := ShallowCopy(orbits[max]);
    for i in [2..Length(base)] do
        if not (reps[base[i]] in significant)  then
            max := reps[base[i]];
            Append(domain,orbits[max]);
            Add(significant,max);
        fi;
    od;
    len := Length(domain);

    # restrict G to significant orbits
    if n = len then
       GG := G;
    else
       tchom := ActionHomomorphism(G,domain,"surjective");
       GG := Image(tchom,G);
    fi;

    # handle case of transitive GG directly
    if Length(significant) = 1  then
        centr := CentralizerTransSymmCSPG( GG, StabChainMutable( GG ) );
        if IsEmpty( GeneratorsOfGroup( centr ) ) then
           return TrivialSubgroup( G );
        else
           order := Size( centr );
           cent := IntersectionNormalClosurePermGroup(GG,centr,order*Size(GG));
           cent:= PreImages(tchom,cent);
           Assert( 1, IsAbelian( cent ) );
           SetIsAbelian( cent, true );
           return cent;
        fi;
    fi;

    # case of intransitive GG
    # for each orbit of GG, construct generators of centralizer of GG in
    # Sym(orbit).  hgens is a list of generators for the direct product of
    # these centralizers.
    # the group generated by hgens contains the center of GG
    hgens := [];
    order := 1;
    for i in significant do
        if n = len then
           orbit := orbits[i];
        else
           orbit := OnTuples(orbits[i],tchom!.conperm);
        fi;
        tchom2 := ActionHomomorphism(GG,orbit,"surjective");
        GGG := Image(tchom2,GG);
        chainGG:= StabChainOp( GG, [ orbit[1] ] );
        chainGGG:= StabChainMutable( GGG );
        chainGGG.stabFxdPnts:=[ orbit[1]^tchom2!.conperm,
            OnTuples( Difference(orbit,
                      MovedPoints( chainGG.stabilizer.generators ) ),
                      tchom2!.conperm ) ];
        centr := CentralizerTransSymmCSPG( GGG, chainGGG );
        if not IsEmpty( GeneratorsOfGroup( centr ) ) then
           order := order * Size( centr );
           inverse2 := tchom2!.conperm^(-1);
           for g in StabChainMutable( centr ).generators do
               Add(hgens,g^inverse2);
           od;
        fi;
    od;

    if order = 1 then
        return TrivialSubgroup( G );
    else
        cent := IntersectionNormalClosurePermGroup
                 ( GG, GroupByGenerators(hgens,()), order*Size(GG) );
        if n <> len then
          cent:= PreImages( tchom, cent );
        fi;
        Assert( 1, IsAbelian( cent ) );
        SetIsAbelian( cent, true );
        return cent;
    fi;
end );


#############################################################################
##
#F  CentralizerNormalCSPG() . . . . . . . .  centralizer of a normal subgroup
##
##  computes the centralizer of a NORMAL subgroup N in G.
##  Reference: Luks-Seress
##
InstallGlobalFunction( CentralizerNormalCSPG, function(G,N)
    local   n,          # degree of G
            orbits,     # list of orbits of G
            base,       # lexicographically smallest (in list) base of G
            i,j,        # loop variables
            reps,       # array recording which orbit of G the points in
                        # perm. domain belong to
            domain,     # union of G orbits which contain base points
            significant,# indices of orbits in "orbits" that belong to domain
            max,        # loop variable, used at definition of significant
            len,        # length of domain
            tchom,      # trans. const. homom, restricting G to domain
            GG,         # the image of G at tchom
            NN,         # the image of N at tchom
            orbit,      # an orbit of GG
            tchom2,     # trans. const. homom, restricting GG to orbit
            GGG,        # the image of GG at tchom2
            NNN,        # the image of NN at tchom2
            hgens,      # list of generators for the direct product of
                        # centralizers of NN in GG restricted to Sym(orbit),
                        # for orbits of GG
            order,      # order of Group(hgens,())
            centrnorm,  # centralizer of NN in GG restricted to Sym(orbit)
            inverse2,   # inverse of the conjugating permutation of tchom2
            g,          # loop variable for generators
            central;    # centralizer of NN in GG

    if IsTrivial(N)  then
        return G;
    fi;

    n := LargestMovedPoint(G);
    orbits := OrbitsDomain(G,[1..n]);
    #orbits := List( orbits, Set );

    # handle case of transitive G directly
    if Length(orbits) = 1  then
        centrnorm := CentralizerNormalTransCSPG(G,N);
        return centrnorm;
    fi;

    # for intransitive G, find which orbit contains which
    #points of permutation domain
    reps := [];
    for i in [1..Length(orbits)] do
        for j in [1..Length(orbits[i])] do
            reps[orbits[i][j]] := i;
        od;
    od;
    #list := Concatenation(orbits);
    #MakeStabChain(G,list);

    # take union of significant orbits which contain base points
    base := BaseStabChain(StabChainMutable(G));
    max := reps[base[1]];
    significant := [max];
    domain := ShallowCopy(orbits[max]);
    for i in [2..Length(base)] do
        if not (reps[base[i]] in significant)  then
            max := reps[base[i]];
            Append(domain,orbits[max]);
            Add(significant,max);
        fi;
    od;
    len := Length(domain);

    # restrict G,N to significant orbits
    if n = len then
       GG := G;
       NN := N;
    else
       tchom := ActionHomomorphism(G,domain,"surjective");
       GG := Image(tchom,G);
       NN := Image(tchom,N);
    fi;

    # handle case of transitive GG directly
    if Length(significant) = 1  then
        centrnorm := CentralizerNormalTransCSPG(GG,NN);
        return PreImages(tchom,centrnorm);
    fi;

    # case of intransitive GG
    # for each GG orbit, compute the centralizer of NN in GG, restricted to
    # the orbit. hgens contains generators for the direct product of these
    # centralizers; the group generated by hgens contains the centralizer of
    # NN in GG
    hgens := [];
    order := 1;
    for i in significant do
        if n = len then
            orbit := orbits[i];
        else
            orbit := OnTuples(orbits[i],tchom!.conperm);
        fi;
        # restrict GG, NN to orbit
        tchom2 := ActionHomomorphism(GG,orbit,"surjective");
        GGG := Image(tchom2,GG);
        NNN := Image(tchom2,NN);

        # compute centralizer of NNN in GGG
        centrnorm := CentralizerNormalTransCSPG(GGG,NNN);
        inverse2 := tchom2!.conperm^(-1);
        order := order * Size(centrnorm);

        # determine how the centralizer acts on domain
        for g in StabChainMutable( centrnorm ).generators do
            Add(hgens,g^inverse2);
        od;
    od;

    if order = 1 then
       return TrivialSubgroup( Parent(G) );
    else
       central := IntersectionNormalClosurePermGroup
                    ( GG, GroupByGenerators(hgens,()), order*Size(GG) );
    fi;

    if n = len then
       return central;
    else
       return PreImages(tchom,central);
    fi;
end );


#############################################################################
##
#F  CentralizerNormalTransCSPG()  . . . centralizer of normal in transitive G
##
##  computes C_G(N) with G transitive, N normal in G
##  reference: Luks-Seress
##
InstallGlobalFunction( CentralizerNormalTransCSPG, function(G,N)
    local   chainG,     # stabilizer chain of `G'
            chainN,     # stabilizer chain of `N'
            n,          # degree of G
            x,          # the first base point of G
            stabgroup,  # stabilizer of x in N
            U,          # an orbit of centralizer of N in S_n
            orbits,     # list of orbits of centralizer of N is S_n
            bhom,       # block homomorphism from G to action on orbits
            GG,         # the kernel of bhom
            GGgens,     # generators of a stabilizer chain of `GG'
            Ginverses,  # list of inverses of generators of G
            Ninverses,  # list of inverses of generators of N
            norbits,    # list of orbits of N
            orbitlength,# the length of the N orbits
                        # (all are of the same size)
            reprlist,   # list recording which orbit of N contains a point of
                        # permutation domain
            positionlist,
                        # list recording the position of a point within its
                        # N orbit
            positiongenlist,
                        # list of length orbitlength; i^th entry records
                        # the position of the generator in N.generators
                        # which occurs in N.transversal at N.orbit[i]
            len,        # number of N orbits intersecting U
            diff,       # loop variable denoting a subset of U
                        # used at creation of N orbits which intersect U
            new,        # an orbit of N intersecting U
            i,j,k,m,    # loop variables
            y,u,s,      # points of permutation domain
            set,        # loop variable denoting subset of [1..n], used at
                        # creation of covering of [1..n] by orbits of N
            newlen,     # loop variable counting the total length of N orbits
                        # at the covering of [1..n]
            word,       # a coset representative of G or N, as a word
            tchom,      # transitive constituent homomorphism restricting
                        # N to N.orbit
            inverse,    # the inverse of tchom.conperm
            img,        # image of `N' under `tchom'
            centr,      # the centralizer of N in Sym(N.orbit)
            chaincentr, # stabilizer chain of `centr'
            hom,        # homomorphism of GG whose kernel is C_G(N)
            images,     # list of images of generators of GG at hom
            top,bottom,g,
                        # permutations used at the creation of images
            K;          # image of GG at hom

    if IsTrivial(N)  then
        return G;
    fi;

    chainG:= StabChainMutable( G );
    x := chainG.orbit[1];
    chainN:= StabChainOp( N, [x] );

    # handle transitive N directly
    if Length( chainN.orbit ) = Length( chainG.orbit ) then
        centr := CentralizerTransSymmCSPG( N, chainN );
        if Size(centr) > 1 then
        return IntersectionNormalClosurePermGroup( G, centr,
                   Size( centr ) * Size( G ) );
        else
           return TrivialSubgroup( Parent(G) );
        fi;
    fi;

    n := LargestMovedPoint(G);
    stabgroup := Stabilizer(N,x,OnPoints);
    U := Difference([1..n],MovedPoints(stabgroup));
    if Length(U) = 1 then
        return TrivialSubgroup( Parent(G) );
    fi;
    orbits:=Blocks(G,[1..n],U);

    # orbits contains the orbits of the centralizer of N in S_n;
    # so C_G(N) must fix setwise the elements of orbits
    bhom := ActionHomomorphism(G,orbits,OnSets,"surjective");
    GG := KernelOfMultiplicativeGeneralMapping( bhom );
    if IsTrivial(GG)  then
        return TrivialSubgroup( Parent(G) );
    fi;
    Size(GG);

    Ginverses := GInverses( chainG );
    Ninverses := GInverses( chainN );

    # we partition [1..n] into the orbits of N, and compute the
    # identification between equivalent orbits (equivalent in the sense
    # that the centralizer of N in S_n exchanges them). After that, we
    # conjugate the union of equivalent orbits to cover [1..n]
    norbits := [ chainN.orbit ];
    orbitlength := Length( chainN.orbit );
    positionlist := [];
    reprlist := [];
    positiongenlist := [];
    for i in [1..orbitlength] do
        positionlist[ chainN.orbit[i] ] := i;
        reprlist[ chainN.orbit[i] ] := 1;
        positiongenlist[i]:= Position( chainN.generators,
                                 chainN.transversal[ chainN.orbit[i] ] );
    od;
    diff := Difference(U,norbits[1]);
    len := 1;

    # create the orbits of N equivalent to the first one
    while diff <> [] do
        len := len+1;
        y := diff[1];
        new := [y];
        positionlist[y] := 1;
        reprlist[y] := len;
        for i in [2..orbitlength] do
            u := chainN.orbit[i] ^ chainN.generators[ positiongenlist[i] ];
            new[i] := new[positionlist[u]]^Ninverses[positiongenlist[i]];
            positionlist[new[i]] := i;
            reprlist[new[i]] := len;
        od;
        Add(norbits,new);
        diff := Difference(diff,new);
    od;

    # if the domain is not covered, create further orbits of N
    if len*orbitlength < n  then
        set := Difference([1..n],Union(norbits));
        for k in [2..n/(len*orbitlength)] do
            newlen := (k-1)*len;
            y := set[1];
            word := CosetRepAsWord( x, y, chainG.transversal );
            word := InverseAsWord( word, chainG.generators, Ginverses );
            for i in [1..len] do
                norbits[newlen+i] := [];
                for j in [1..orbitlength] do
                    norbits[newlen+i][j] := ImageInWord(norbits[i][j],word);
                    positionlist[norbits[newlen+i][j]] := j;
                    reprlist[norbits[newlen+i][j]] := newlen+i;
                od;
                set := Difference(set,norbits[newlen+i]);
            od;
        od;
    fi;

    # compute centralizer of N in first orbit; centralizer in other orbits
    # is obtained from identification between orbits
    tchom := ActionHomomorphism( N, chainN.orbit,"surjective" );
    inverse := tchom!.conperm^(-1);
    img:= Image( tchom, N );
    centr := CentralizerTransSymmCSPG( img, StabChainMutable( img ) );

    # compute (and store) transversal of centr
    chaincentr:= EmptyStabChain( [  ], (), x^tchom!.conperm );
    AddGeneratorsExtendSchreierTree( chaincentr, GeneratorsOfGroup(centr));

    # compute images at homomorphism of GG, g -> g c_g^{-1} (cf. Luks-Seress)
    # the kernel of this homomorphism is C_G(N)
    images := [];
    GGgens := StabChainMutable( GG ).generators;
    for i in [1..Length( GGgens)] do
        images[i] := [];

        # top is the permutation in the wreath product which pulls back g to
        # orbits of N
        top := [];
        for j in [1..Length(norbits)] do
            k := reprlist[norbits[j][1]^GGgens[i]];
            for m in [1..orbitlength] do
                top[norbits[k][m]] := norbits[j][m];
            od;
        od;
        top := PermList(top);
        g := GGgens[i]*top;

        # pull back each leading point in norbits by centralizer of N
        bottom := [];
        for j in [1..Length(norbits)] do
            k := positionlist[norbits[j][1]^g];
            word := CosetRepAsWord( x^tchom!.conperm,
                                    chainN.orbit[k]^tchom!.conperm,
                                 chaincentr.transversal);
            for m in [1..orbitlength] do
                s := (ImageInWord( chainN.orbit[m]^tchom!.conperm,
                                   word ))^inverse;
                bottom[norbits[j][m]] := norbits[j][positionlist[s]];
            od;
         od;
         bottom := PermList(bottom);
         images[i] := g*bottom;
    od;

    K := GroupByGenerators(images,());
    hom := GroupHomomorphismByImagesNC(GG,K,GGgens,images);
    return KernelOfMultiplicativeGeneralMapping( hom );
end );


#############################################################################
##
#F  CentralizerTransSymmCSPG()  . . . . .  centralizer of transitive G in S_n
##
##  computes the centralizer of a transitive group G in S_n
##
InstallGlobalFunction( CentralizerTransSymmCSPG, function( G, chainG )
    local   n,          # the degree of G
            x,          # the first base point
            L,          # the set of fixed points of stabgroup
            orbitx,     # the orbit of x in the centralizer;
                        # eventually, orbitx=L
            y,          # a point in L
            z,          # loop variable running through permutation domain
            h,          # a coset representative of G, written as word in the
                        # generators
            gens,       # list of generators for the centralizer
            gen,        # an element of gens
            Ggens,      # generators of G
            Ginverses,  # list of inverses for the generators of G
            H;          # output group
    if IsTrivial(G)  then
        return TrivialSubgroup( Parent(G) );
    fi;

    if IsBound( chainG.stabFxdPnts ) then
       x := chainG.stabFxdPnts[1];
       L := chainG.stabFxdPnts[2];
       n := LargestMovedPoint(G);
       if not IsBound( chainG.orbit ) or chainG.orbit[1] <> x then
          chainG := EmptyStabChain( [  ], (), x );
          AddGeneratorsExtendSchreierTree( chainG, GeneratorsOfGroup(G) );
       fi;
    else
       n := LargestMovedPoint(G);
       x := chainG.orbit[1];
       L := Difference( [ 1 .. n ],
                        MovedPoints( chainG.stabilizer.generators ) );
    fi;

    Ginverses := GInverses( chainG );
    Ggens := chainG.generators;

    # the centralizer of G is semiregular, acting transitively on L
    orbitx := [x];
    gens := [];
    while Length(orbitx) < Length(L) do

        # construct element of centralizer which carries x to new point in L
        gen := [];
        y := Difference(L,orbitx)[1];
        for z in [1..n] do
            h := CosetRepAsWord( x, z, chainG.transversal );
            h := InverseAsWord(h,Ggens,Ginverses);
            gen[z] := ImageInWord(y,h);
        od;
        Add(gens,PermList(gen));
        orbitx := OrbitPerms(gens,x);
    od;

    H := SubgroupNC( G, gens );
    SetSize( H, Length( L ) );
    return H;
end );


#############################################################################
##
#F  IntersectionNormalClosurePermGroup(<G>,<H>[,order]) . . . intersection of
#F                                   normal closure of <H> under <G> with <G>
##
##  computes $H^G \cap G$ as subgroup of Parent(G)
##
InstallGlobalFunction( IntersectionNormalClosurePermGroup,
    function(arg)
    local   G,H,        # the groups to be handled
            n,          # maximum of degrees of G,H
            i,          # loop variables
            conperm,    # perm exchanging first and second n points
            newgens,    # set of extended generators
            options,    # options record for stabilizer computation
            group;      # the group generated by newgens
                        # stabilizing the second n points, we get H^G \cap G

    G := arg[1];
    H := arg[2];

    if IsTrivial(G) or IsTrivial(H)  then
        return TrivialSubgroup( Parent(G) );
    fi;

    n := Maximum(LargestMovedPoint(G),
                 LargestMovedPoint(H));
    conperm := PermList( Concatenation( [n+1 .. 2*n] , [1 .. n] ) );
    # extend the generators of G acting on [n+1..2n] exactly as on [1..n]
    newgens := List( StabChainMutable( G ).generators,
                     g -> g * ( g^conperm ) );

    # from the generators of H, create permutations which act on [n+1..2n]
    # as the original generator on [1..n] and which act trivially on [1..n]
    for i in StabChainMutable( H ).generators do
      Add( newgens, i^conperm );
    od;

    group := GroupByGenerators(newgens,());

    # create options record for stabilizer chain computation
    options := rec( base := [n+1..2*n] );
    #if size of group is part of input, use it
    if Length(arg) = 3 then
       options.size := arg[3];
       # if H is normalized by G and G,H already have stabilizer chains
       # then compute base for group
       #if ( IsBound(G.size) or IsBound(G.stabChain) ) and
       #   ( IsBound(H.size) or IsBound(H.stabChain) )  then
       #   if Size(G) * Size(H) = arg[3] then
       #      options.knownBase :=
       #      Concatenation( List( Base(H), x -> n + x ), Base(G) ) ;
       #   fi;
       #fi;
    fi;
    StabChain(group,options);
#T is this meaningful ??
    group := Stabilizer(group,[n+1 .. 2*n],OnTuples);
    return AsSubgroup( Parent(G),group);
end );


#############################################################################
##
#F  ActionAbelianCSPG() . . . . . . . . . action of abelian permutation group
##
##  given an abelian subgroup H of S_n, the routine codes the action of
##  H on its orbits. The output is an array of length 7, describing this
##  action; the components of this array are described at the local variable
##  section
##
InstallGlobalFunction( ActionAbelianCSPG, function(H,n)
    local   i,j,k,      # loop variables
            orbits,     # list of orbits of H; 6th element of output
            action,     # list; the i^th element contains a list of
                        # generators for the action of H on i^th orbit
            inverse,    # inverse[i][k] is the inverse of action[i][k]
                        # 1st element of output
            Hgens,      # generators of `H'
            C,          # C[i] is the stabilizer chain of the group
                        # generated by action[i]
                        # 2nd element of output
            chainC,     # one stabilizer chain in `C'
            positionlist,
                        # for i in [1..n], positionlist[i] gives the position
                        # of i in its H orbit. 3rd element of output
            reprlist,   # for i in [1..n], reprlist[i] gives the position of
                        # the H orbit of i in orbits. 4th element of output
            cpositiongenlist,
                        # cpositionlength[i][k] gives the position in
                        # action[i] of the C[i] generator which occurs in
                        # C[i].transversal[k]. 5th element of output
            cumulativelength;
                        # cumulativelength[i] is the sum of lengths of first
                        # i-1 elements of orbits. 7th element of output

    orbits := OrbitsDomain(H,[1..n]);
    cumulativelength := [0];
    for i in [1..Length(orbits)-1] do
        cumulativelength[i+1] := cumulativelength[i]+Length(orbits[i]);
    od;

    positionlist := [];
    reprlist := [];
    for i in [1..Length(orbits)] do
        for j in [1..Length(orbits[i])] do
            positionlist[orbits[i][j]] := j;
            reprlist[orbits[i][j]] := i;
        od;
    od;

    # action[i][k] is the action of H.generators[k] on the i^th orbit of H,
    # viewed as a permutation on [1..Length(orbits[i])]
    action := [];
    inverse := [];
    Hgens:= StabChainMutable( H ).generators;
    for i in [1..Length(orbits)] do
        action[i] := [];
        inverse[i] := [];
        for k in [1..Length(Hgens)] do
            action[i][k] := [];
            for j in [1..Length(orbits[i])] do
                action[i][k][j]:=positionlist[orbits[i][j]^Hgens[k]];
            od;
            action[i][k] := PermList(action[i][k]);
            inverse[i][k] := action[i][k]^(-1);
        od;
    od;

    C := [];
    cpositiongenlist := [];
    for i in [1..Length(orbits)] do
        cpositiongenlist[i] := [];

        # create stabilizer chain C[i]
        chainC := EmptyStabChain( [  ], (), 1 );
        AddGeneratorsExtendSchreierTree( chainC, action[i] );
        C[i]:= chainC;
        Add(action[i],());
        Add(inverse[i],());

        # determine position of generators occurring in transversal
        for j in [1..Length( chainC.orbit )] do
            cpositiongenlist[i][j]:=Position(action[i],chainC.transversal[j]);
        od;

    od;

    return [inverse,C,positionlist,reprlist,
            cpositiongenlist,orbits,cumulativelength];
end );


#############################################################################
##
#F  ImageOnAbelianCSPG( <g>, <actionlist> ) . .  image of normalizing element
#F  . . . . . . . . . . . . . . . . . . . . . . . . . .  on orbits of abelian
##
##  Given the action of an abelian group $H$ encoded in <actionlist> by the
##  subroutine `ActionAbelianCSPG', and a permutation <g> normalizing H,
##  this subroutine computes the conjugation action of <g> on the transitive
##  constituent pieces of $H$.
##
InstallGlobalFunction( ImageOnAbelianCSPG, function(g,actionlist)
    local   i,s,        # loop variables
            orbits,     # list of orbits of H
      # let action denote the list with the i^th element containing a list of
      # generators for the action of H on i^th orbit
            inverse,    # inverse[i][k] is the inverse of action[i][k]
            C,          # C[i] is a stabilizer chain of the group generated
                        # by action[i]
            positionlist,
                        # for i in [1..n], positionlist[i] gives the position
                        # of i in its H orbit
            reprlist,   # for i in [1..n], reprlist[i] gives the position of
                        # the H orbit of i in orbits
            cpositiongenlist,
                        # cpositionlength[i][k] gives the position in
                        # action[i] of the C[i] generator which occurs in
                        # C[i].transversal[k]
            cumulativelength,
                        # cumulativelength[i] is the sum of lengths of first
                        # i-1 elements of orbits
            j,          # index of H-orbit in orbits which is the image of
                        # the i^th H-orbit
            x,          # position of element of i^th orbit which is mapped
                        # by g to first element of j^th orbit
            inv,        # the inverse of g
            gimage,     # output of the routine; conjugation action of g
            image,t;    # see explanation in body of routine

    inv := g^(-1);
    gimage := [];
    inverse := actionlist[1];
    C := actionlist[2];
    positionlist := actionlist[3];
    reprlist := actionlist[4];
    cpositiongenlist := actionlist[5];
    orbits := actionlist[6];
    cumulativelength := actionlist[7];

    # the transitive constituent pieces of H are regarded as a list of
    # length n; the (unique) piece carrying the first point of i^th orbit
    # to k^th point of i^th orbit is in the position cumulativelength[i]+k.
    # gimage will contain the conjugation action of g on the elements of
    # this list
    for i in [1..Length(orbits)] do

        # determine which orbit contains the images of pieces from i^th orbit
        j := reprlist[orbits[i][1]^g];

        # for each piece h from i^th orbit, we have to determine the image of
        # orbits[j][1] at the permutation g^(-1)*h*g
        # from regularity of action on orbits, this image determines the
        # conjugate first, compute the images of orbits[j][1] in g^(-1)*h,
        # and store the result in the array image. Then determine the
        # g-image of the result and store it in gimage.
        # This way, elements in "image" can be used more times,
        # and the running time is linear (no hidden log factors).
        x := positionlist[orbits[j][1]^inv];
        image := [x];
        gimage[cumulativelength[i]+1] := cumulativelength[j]+1;
        for s in [2..Length(C[i].orbit)] do

            # t is the predecessor in Schreier tree of C[i].orbit[s]
            t := C[i].orbit[s]^C[i].transversal[C[i].orbit[s]];
            image[C[i].orbit[s]] :=
                     image[t]^inverse[i][cpositiongenlist[i][C[i].orbit[s]]];
            gimage[cumulativelength[i]+C[i].orbit[s]] :=
                     cumulativelength[j]
                    +positionlist[orbits[i][image[C[i].orbit[s]]]^g];
        od;
    od;

    gimage := PermList(gimage);
    return gimage;
end );


# ser is descending subnormal series, nt a descending series of normal subs
BindGlobal("ChangeSeriesThrough",function(ser,nt)
local new,start,n,i,tail, up,u,v;

  new:=Reversed(ser); # we step up (closure is easier than intersection)

  # make nt also increasing
  nt:=ShallowCopy(nt);
  SortBy(nt, Size);
#Print(List(nt,Size),"\n");

  start:=1;
  while Length(nt)>0 do
    n:=nt[1];
    nt:=nt{[2..Length(nt)]};
    ser:=new;
    new:=ser{[1..start-1]};
    i:=start;
    while i<=Length(ser) and IsSubset(n,ser[i]) do
      Add(new,ser[i]);
      i:=i+1;
    od;
    # now n does not contain ser[i]

    # was n actually in the series?
    if Last(new)=n then
      # yes, go on and add the rest of the series
      start:=i; # next time start from next step
    else
      # no generate/intersect

      # in each step either we ascend in intersection with n or in closure with
      # n
      tail:=[];
      up:=[n];
      u:=ClosureGroup(n,ser[i]);
      Add(up,u);
      i:=i+1;
      while not IsSubset(ser[i],n) do
        v:=ClosureGroup(n,ser[i]);
        if Size(v)=Size(u) then
          # no increase, need for tail
          Add(tail,i);
        else
          Add(up,v);
          u:=v;
        fi;
        i:=i+1;
      od;
#Print("A",List([1..Length(new)-1],x->Size(new[x+1])/Size(new[x])),"\n");

      # now ser[i] contains n.
      for i in tail do
        Add(new,NormalIntersection(n,ser[i]));
      od;
#Print("B",List([1..Length(new)-1],x->Size(new[x+1])/Size(new[x])),"\n");
      start:=Length(new)+1;
      Append(new,up);
#Print("C",List([1..Length(new)-1],x->Size(new[x+1])/Size(new[x])),"\n");
      i:=i+1;
      while i<=Length(ser) and Size(Last(new))>=Size(ser[i]) do
        i:=i+1;
      od;

    fi;

    # add the rest
    while i<=Length(ser) do
      Add(new,ser[i]);
      i:=i+1;
    od;
#Print("D",List([1..Length(new)-1],x->Size(new[x+1])/Size(new[x])),"\n");

  od;
  return Reversed(new);
end);


#############################################################################
##
#F  ChiefSeriesOfGroup( [<H>, ]<G>[, <through>] )
##
InstallGlobalFunction( ChiefSeriesOfGroup, function(arg)
local G,H,nser,U,i,j,k,cs,n,mat,mats,row,p,one,m,v,ser,gens,r,dim,im,
      through,ocs;
  G:=arg[1];
  H:=G;
  through:=[];
  if Length(arg)=2 then
    if IsGroup(arg[2]) then
      H:=arg[1];
      G:=arg[2];
    else
      through:=arg[2];
    fi;
  elif Length(arg)>2 then
    H:=arg[1];
    G:=arg[2];
    through:=arg[3];
  fi;

  if Length(through)>0 then
    nser:=ChiefSeriesOfGroup(G,H);
    nser:=ChangeSeriesThrough(nser,through);
    return nser;
  fi;

  nser:=[G];
  U:=G;
  while Size(U)>1 do
    # get maximal normal subgroup
    if Size(U)<Size(G) and Size(ocs[1])/Size(U)<1000 then
      n:=List(ocs,i->Intersection(U,i));
      cs:=[U];
      for i in [2..Length(n)] do
        if Size(Last(cs))>Size(n[i]) then
          Add(cs,n[i]);
        fi;
      od;
    else
      cs:=CompositionSeries(U);
    fi;
    ocs:=cs;
    # add composition factors which are normal
    n:=2;
    while n<=Length(cs) and Length(through)=0 and
      # IsNormal(H,cs[n]) do
      ForAll(GeneratorsOfGroup(H),x->ForAll(GeneratorsOfGroup(cs[n]),
                                       y->y^x in cs[n])) do
      U:=cs[n];
      Add(nser,U);
      n:=n+1;
    od;

    if n<=Length(cs) then
      cs:=cs[n];

      if Length(through)>0 then
        if Size(U)=Size(through[1]) then
          through:=through{[2..Length(through)]};
        fi;
        if Length(through)>0 and not IsSubgroup(cs,through[1]) then
          # enforce way through
          Info(InfoGroup,1,"force");
          n:=NaturalHomomorphismByNormalSubgroup(U,through[1]);
          cs:=CompositionSeries(Image(n));
          cs:=cs[2];
          cs:=PreImage(n,cs);
        fi;
      fi;

      #n:=Core(H,cs);
      n:=cs;
      i:=1;
      gens:=GeneratorsOfGroup(H);
      while i<=Length(gens) do
        if not ForAll(GeneratorsOfGroup(n), x->x^gens[i] in n) then
          if IsIdenticalObj(FamilyObj(One(n)),FamilyObj(gens[i])) then
            n:=Intersection(n,n^gens[i]);
          else
            n:=Intersection(n,Image(gens[i],n));
          fi;
          i:=1;
        else
          i:=i+1;
        fi;
      od;

      #o:=GroupOnSubgroupsOrbit(H,cs);
      #Info(InfoGroup,1,"orblen=",Length(o));
      #n:=Intersection(o);
      #n:=o[1];
      #for i in o{[2..Length(o)]} do
        #n:=IntersectionNormalClosurePermGroup(n,i);
      #od;
      if HasAbelianFactorGroup(U,cs) then
        # abelian case, utilize MeatAxe to chop

        p:=Index(U,cs);
        one:=One(GF(p));

        # first get series
        v:=n;
        ser:=[n];
        gens:=[];
        while Size(v)<Size(U) do
          repeat
            r:=Random(U);
          until not r in v;
          Add(gens,r);
          v:=ClosureGroup(v,r);
          Add(ser,v);
        od;

        dim:=Length(gens);
        ser:=Reversed(ser);
        gens:=Reversed(gens);

        # now construct matrices for operation
        mats:=[];
        for i in GeneratorsOfGroup(H) do
          mat:=[];
          for j in gens do
            im:=j^i;
            row:=[];
            for k in [1..dim] do
              if not im in ser[k+1] then
                # test power which does
                # Ug^l=U im
                r:=First([1..p],l->im/gens[k]^l in ser[k+1]);
                Add(row,r);
                im:=im/gens[k]^r;
              else
                Add(row,0);
              fi;
            od;
            row:=row*one;
            Add(mat,row);
          od;
          Add(mats,mat);
        od;

        m:=GModuleByMats(mats,GF(p));
        r:=MTX.BasesCompositionSeries(m);
        v:=[];
        for i in r do
          im:=n;
          for j in i do
            im:=ClosureGroup(im,Product([1..dim],k->gens[k]^IntFFE(j[k])));
          od;
          # only intermediates
          if Size(im)<Size(U) then
            Add(v,im);
          fi;
        od;
        v:=Reversed(v); # MTX sorts already
        #Sort(v,function(a,b) return Size(a)>Size(b);end);
        Info(InfoGroup,2,"i:",List(v,Size));

        #note the intermediates
        nser:=Concatenation(nser,v);

      else
        # nonabelian, as transitive operation on the components no proper
        # intermediate normal subgroup possible
        Add(nser,n);
      fi;

    else
      n:=cs[n-1];
    fi;
    Info(InfoGroup,1,"Step ",Index(U,n));
    U:=n;
  od;
  return nser;
end );


#############################################################################
##
#M  ChiefSeries( <G> )
##
InstallMethod( ChiefSeries,
    "generic method for a group",
    true,
    [ IsGroup ], 0,
    ChiefSeriesOfGroup );


#############################################################################
##
#M  ChiefSeriesUnderAction( <G>, <H> )
##
InstallMethod( ChiefSeriesUnderAction,
    "generic method for two groups",
    true,
    [ IsGroup, IsGroup ], 0,
    ChiefSeriesOfGroup );


#############################################################################
##
#M  ChiefSeriesThrough( <G>, <list> )
#M  ChiefSeriesThrough( <G>, <H>, <list> )
##
InstallMethod( ChiefSeriesThrough,
    "generic method for a group and a list",
    true,
    [ IsGroup, IsList ], 0,
    ChiefSeriesOfGroup );

InstallOtherMethod( ChiefSeriesThrough,
    "generic method for two groups and a list",
    true,
    [ IsGroup, IsGroup, IsList ], 0,
    ChiefSeriesOfGroup );