/usr/share/gap/lib/grpperm.gi

############################################################################
##
#W  grpperm.gi                  GAP library                   Heiko Theißen
##
##
#Y  Copyright (C)  1996,  Lehrstuhl D für Mathematik,  RWTH Aachen,  Germany
#Y  (C) 1998 School Math and Comp. Sci., University of St Andrews, Scotland
#Y  Copyright (C) 2002 The GAP Group
##


#############################################################################
##
#M  CanEasilyTestMembership( <permgroup> )
##
InstallTrueMethod(CanEasilyTestMembership,IsPermGroup);


#############################################################################
##
#M  CanEasilyComputePcgs( <permgroup> )
##
##  solvable permgroups can compute a pcgs. (The group may have been told it
##  is solvable without actually doing a solvability test, so we need the
##  implication.)
InstallTrueMethod(CanEasilyComputePcgs,IsPermGroup and IsSolvableGroup);

#M  CanComputeSizeAnySubgroup
InstallTrueMethod(CanComputeSizeAnySubgroup,IsPermGroup);

#############################################################################
##
#M  AsSubgroup( <G>, <U> )  . . . . . . . . . . . .  with stab chain transfer
##
InstallMethod( AsSubgroup,"perm groups",
    IsIdenticalObj, [ IsPermGroup, IsPermGroup ], 0,
    function( G, U )
    local S;
    # test if the parent is already alright
    if HasParent(U) and IsIdenticalObj(Parent(U),G) then
      return U;
    fi;

    if not IsSubset( G, U ) then
      return fail;
    fi;
    S:= SubgroupNC( G, GeneratorsOfGroup( U ) );
    UseIsomorphismRelation( U, S );
    UseSubsetRelation( U, S );
    if HasStabChainMutable( U )  then
        SetStabChainMutable( S, StabChainMutable( U ) );
    fi;
    return S;
end );

#############################################################################
##
#M  CanEasilyComputeWithIndependentGensAbelianGroup( <permgroup> )
##
InstallTrueMethod(CanEasilyComputeWithIndependentGensAbelianGroup,
    IsPermGroup and IsAbelian);

#############################################################################
##
#F  IndependentGeneratorsAbelianPPermGroup( <P>, <p> )  . . . nice generators
##
InstallGlobalFunction( IndependentGeneratorsAbelianPPermGroup,
    function ( P, p )
    local   inds,       # independent generators, result
            pows,       # their powers
            base,       # the base of the vectorspace
            size,       # the size of the group generated by <inds>
            orbs,       # orbits
            trns,       # transversal
            gens,       # remaining generators
            one,        # identity of `P'
            gens2,      # remaining generators for next round
            exp,        # exponent of <P>
            g,          # one generator from <gens>
            h,          # its power
            b,          # basepoint
            c,          # other point in orbit
            i, j, k;    # loop variables

    # initialize the list of independent generators
    inds := [];
    pows := [];
    base := [];
    size := 1;
    orbs := [];
    trns := [];

    # gens are the generators for the remaining group
    gens := GeneratorsOfGroup( P );
    one:= One( P );

    # loop over the exponents
    exp := Maximum( List( gens, g -> LogInt( Order( g ), p ) ) );
    for i  in [exp,exp-1..1]  do

        # loop over the remaining generators
        gens2 := [];
        for j  in [1..Length(gens)]  do
            g := gens[j];
            h := g ^ (p^(i-1));

            # reduce <g> and <h>
            while h <> h^0
              and IsBound(trns[SmallestMovedPoint(h)^h])
            do
                g := g / pows[ trns[SmallestMovedPoint(h)^h] ];
                h := h / base[ trns[SmallestMovedPoint(h)^h] ];
            od;

            # if this is linear indepenent, add it to the generators
            if h <> h^0  then
                Add( inds, g );
                Add( pows, g );
                Add( base, h );
                size := size * p^i;
                b := SmallestMovedPoint(h);
                if not IsBound( orbs[b] )  then
                    orbs[b] := [ b ];
                    trns[b] := [ one ];
                fi;
                for c  in ShallowCopy(orbs[b])  do
                    for k  in [1..p-1]  do
                        Add( orbs[b], c ^ (h^k) );
                        trns[c^(h^k)] := Length(base);
                    od;
                od;

            # otherwise reduce and add to gens2
            else
                Add( gens2, g );

            fi;

        od;

        # prepare for the next round
        gens := gens2;
        pows := List( pows,i->i^ p );

    od;

    # return the indepenent generators
    return inds;
end );

#############################################################################
##
#M  IndependentGeneratorsOfAbelianGroup( <G> )  . . . . . . . nice generators
##
InstallMethod( IndependentGeneratorsOfAbelianGroup, "for perm group",
        [ IsPermGroup and IsAbelian ],
    function ( G )
    local   inds,       # independent generators, result
            one,        # identity of `G'
            p,          # prime factor of group size
            gens,       # generators of <p>-Sylowsubgroup
            g,          # one generator
            o;          # its order

    # loop over all primes
    inds := [];
    one:= One( G );
    for p  in Union( List( GeneratorsOfGroup( G ),
            g -> Factors(Order(g)) ) )  do
      if p <> 1  then

        # compute the generators for the Sylow <p> subgroup
        gens := [];
        for g  in GeneratorsOfGroup( G )  do
            o := Order(g);
            while o mod p = 0  do o := o / p; od;
            if g^o <> g^0  then Add( gens, g^o );  fi;
        od;

        # append the independent generators for the Sylow <p> subgroup
        Append( inds,
                IndependentGeneratorsAbelianPPermGroup(
                    GroupByGenerators( gens, one ), p ) );

      fi;
    od;

    # Sort the independent generators by increasing order
    Sort( inds, function (a,b) return Order(a) < Order(b); end );

    # return the independent generators
    return inds;
end );

#############################################################################
##
#F  OrbitPerms( <gens>, <d> ) . . . . . . . . . . . . . orbit of permutations
##
InstallGlobalFunction( OrbitPerms, function( gens, d )
    local   max,  orb,  new,  pnt,  img,  gen;
  
        # get the largest point <max> moved by the group <G>
        max := LargestMovedPoint( gens );

        # handle fixpoints
        if not d in [1..max]  then
            return [ d ];
        fi;

        # start with the singleton orbit
        orb := [ d ];
        new := BlistList( [1..max], [1..max] );
        new[d] := false;

        # loop over all points found
        for pnt  in orb  do

            # apply all generators <gen>
            for gen  in gens  do
                img := pnt ^ gen;

                # add the image <img> to the orbit if it is new
                if new[img]  then
                    Add( orb, img );
                    new[img] := false;
                fi;

            od;

        od;
        return orb;
end );

#############################################################################
##
#F  OrbitsPerms( <gens>, <D> )  . . . . . . . . . . .  orbits of permutations
##
InstallGlobalFunction( OrbitsPerms, function( gens, D )
    local   max,  dom,  new,  orbs,  fst,  orb,  pnt,  gen,  img;
    
        # get the largest point <max> moved by the group <G>
        max := LargestMovedPoint( gens );
        dom := BlistList( [1..max], D );
        new := BlistList( [1..max], [1..max] );
        orbs := [];

        # repeat until the domain is exhausted
        fst := Position( dom, true );
        while fst <> fail  do

            # start with the singleton orbit
            orb := [ fst ];
            new[fst] := false;
            dom[fst] := false;

            # loop over all points found
            for pnt  in orb  do

                # apply all generators <gen>
                for gen  in gens  do
                    img := pnt ^ gen;

                    # add the image <img> to the orbit if it is new
                    if new[img]  then
                        Add( orb, img );
                        new[img] := false;
                        dom[img] := false;
                    fi;

                od;

            od;

            # add the orbit to the list of orbits and take next point
            Add( orbs, orb );
            fst := Position( dom, true, fst );

        od;

        # add the remaining points of <D>, they are fixed
        for pnt  in D do
	  if pnt>max then
            Add( orbs, [ pnt ] );
	  fi;
        od;
        
        return orbs;
end );


#############################################################################
##
#M  SmallestMovedPoint( <C> ) . . . . . . .  for a collection of permutations
##
SmallestMovedPointPerms := function( C )
    local   min,  m,  gen;
    
    min := infinity;
    for gen  in C  do
        if not IsOne( gen ) then
            m := SmallestMovedPointPerm( gen );
            if m < min  then
                min := m;
            fi;
        fi;
    od;
    return min;
end;

InstallMethod( SmallestMovedPoint,
    "for a collection of permutations",
    true,
    [ IsPermCollection ], 0,
    SmallestMovedPointPerms );

InstallMethod( SmallestMovedPoint,
    "for an empty list",
    true,
    [ IsList and IsEmpty ], 0,
    SmallestMovedPointPerms );


#############################################################################
##
#M  LargestMovedPoint( <C> )  . . . . . . .  for a collection of permutations 
##
LargestMovedPointPerms := function( C )
    local   max,  m,  gen;
    
    max := 0;
    for gen  in C  do
        if not IsOne( gen )  then
            m := LargestMovedPointPerm( gen );
            if max < m  then
                max := m;
            fi;
        fi;
    od;
    return max;
end;

InstallMethod( LargestMovedPoint,
    "for a collection of permutations",
    true,
    [ IsPermCollection ], 0,
    LargestMovedPointPerms );

InstallMethod( LargestMovedPoint,
    "for an empty list",
    true,
    [ IsList and IsEmpty ], 0,
    LargestMovedPointPerms );


#############################################################################
##
#F  MovedPoints( <C> )  . . . . . . . . . .  for a collection of permutations 
##
InstallGlobalFunction(MovedPointsPerms,function( C )
    local   mov,  gen,  pnt;
    
    mov := [  ];
    for gen  in C  do
        if gen <> One( gen )  then
            for pnt  in [ SmallestMovedPoint( gen ) ..
                           LargestMovedPoint( gen ) ]  do
                if pnt ^ gen <> pnt  then
                    mov[ pnt ] := pnt;
                fi;
            od;
        fi;
    od;
    return Set( mov );
end);

InstallMethod( MovedPoints, "for a collection of permutations", true,
    [ IsPermCollection ], 0,
    MovedPointsPerms );

InstallMethod( MovedPoints, "for an empty list", true,
    [ IsList and IsEmpty ], 0,
    MovedPointsPerms );

InstallMethod( MovedPoints, "for a permutation", true, [ IsPerm ], 0,
function(p)
  return MovedPointsPerms([p]);
end);


#############################################################################
##
#M  NrMovedPoints( <C> )  . . . . . . . . .  for a collection of permutations 
##
NrMovedPointsPerms := function( C )
    local mov, sma, pnt, gen;
    
    mov := 0;
    sma := SmallestMovedPoint( C );
    if sma = infinity  then
        return 0;
    fi;
    for pnt  in [ sma .. LargestMovedPoint( C ) ]  do
      for gen in C do
        if pnt ^ gen <> pnt then
          mov:= mov + 1;
          break;
        fi;
      od;
    od;
    return mov;
end;

InstallMethod( NrMovedPoints,
    "for a collection of permutations",
    [ IsPermCollection ],
    NrMovedPointsPerms );

InstallMethod( NrMovedPoints,
    "for an empty list",
    [ IsList and IsEmpty ],
    NrMovedPointsPerms );


#############################################################################
##
#M  LargestMovedPoint( <G> )  . . . . . . . . . . . . for a permutation group
##
InstallMethod( LargestMovedPoint,
    "for a permutation group",
    true,
    [ IsPermGroup ], 0,
    G -> LargestMovedPoint( GeneratorsOfGroup( G ) ) );


#############################################################################
##
#M  SmallestMovedPoint( <G> ) . . . . . . . . . . . . for a permutation group
##
InstallMethod( SmallestMovedPoint,
    "for a permutation group",
    true,
    [ IsPermGroup ], 0,
    G -> SmallestMovedPoint( GeneratorsOfGroup( G ) ) );


#############################################################################
##
#M  MovedPoints( <G> )  . . . . . . . . . . . . . . . for a permutation group
##
InstallMethod( MovedPoints,
    "for a permutation group",
    true,
    [ IsPermGroup ], 0,
    G -> MovedPoints( GeneratorsOfGroup( G ) ) );


#############################################################################
##
#M  NrMovedPoints( <G> )  . . . . . . . . . . . . . . for a permutation group
##
InstallMethod( NrMovedPoints,
    "for a permutation group",
    true,
    [ IsPermGroup ], 0,
    G -> NrMovedPoints( GeneratorsOfGroup( G ) ) );


#############################################################################
##
#M  BaseOfGroup( <G> )
##
InstallMethod( BaseOfGroup,
    "for a permutation group",
    true,
    [ IsPermGroup ], 0,
    G -> BaseStabChain( StabChainMutable( G ) ) );


#############################################################################
##
#M  Size( <G> ) . . . . . . . . . . . . . . . . . . size of permutation group
##
InstallMethod( Size,
    "for a permutation group",
    true,
    [ IsPermGroup ], 0,
    G -> SizeStabChain( StabChainMutable( G ) ) );


#############################################################################
##
#M  Enumerator( <G> ) . . . . . . . . . . . . enumerator of permutation group
##
BindGlobal( "ElementNumber_PermGroup",
    function( enum, pos )
    local   elm,  S,  len, n;

    S := enum!.stabChain;
    elm := S.identity;
    n:= pos;
    pos := pos - 1;
    while Length( S.genlabels ) <> 0  do
        len := Length( S.orbit );
        elm := LeftQuotient( InverseRepresentative
                       ( S, S.orbit[ pos mod len + 1 ] ), elm );
        pos := QuoInt( pos, len );
        S := S.stabilizer;
    od;

    if pos <> 0 then
      Error( "<enum>[", n, "] must have an assigned value" );
    fi;
    return elm;
end );

BindGlobal( "NumberElement_PermGroup",
    function( enum, elm )
    local   pos,  val,  S,  img;

    if not IsPerm( elm ) then
      return fail;
    fi;

    pos := 1;
    val := 1;
    S := enum!.stabChain;
    while Length( S.genlabels ) <> 0  do
        img := S.orbit[ 1 ] ^ elm;
        #pos := pos + val * ( Position( S.orbit, img ) - 1 );
        pos := pos + val * S.orbitpos[img];
        #val := val * Length( S.orbit );
        val := val * S.ol;
        elm := elm * InverseRepresentative( S, img );
        S := S.stabilizer;
    od;
    if elm <> S.identity  then
        return fail;
    fi;
    return pos;
end );

InstallMethod( Enumerator,
    "for a permutation group",
    [ IsPermGroup ],
function(G)
local e,S,i,p;
    # get a good base
    S:=G;
    p:=[];
    while Size(S)>1 do
      e:=Orbits(S,MovedPoints(S));
      i:=List(e,Length);
      i:=Position(i,Maximum(i));
      Add(p,e[i][1]);
      S:=Stabilizer(S,e[i][1]);
    od;
    Info(InfoGroup,1,"choose base ",p);

    S:=StabChainOp(G,rec(base:=p));
    e:=EnumeratorByFunctions( G, rec(
             ElementNumber := ElementNumber_PermGroup,
             NumberElement := NumberElement_PermGroup,
             Length        := enum -> SizeStabChain( enum!.stabChain ),
             PrintObj      := function( G )
                                  Print( "<enumerator of perm group>" );
                              end,

             stabChain     := S ) );
  while Length(S.genlabels)<>0 do
    S.ol:=Length(S.orbit);
    S.orbitpos:=[];
    for i in [1..Maximum(S.orbit)] do
      p:=Position(S.orbit,i);
      if p<>fail then
	S.orbitpos[i]:=p-1;
      fi;
    od;
    S:=S.stabilizer;
  od;
  return e;
end);


#############################################################################
##
#M  Random( <G> ) . . . . . . . . . . . . . . . . . . . . . .  random element
##
InstallMethod( Random,
    "for a permutation group",
    [ IsPermGroup ], 10,
    function( G )
    local   S,  rnd;

    # go down the stabchain and multiply random representatives
    S := StabChainMutable( G );
    rnd := S.identity;
    while Length( S.genlabels ) <> 0  do
        rnd := LeftQuotient( InverseRepresentative( S,
                       Random( S.orbit ) ), rnd );
        S := S.stabilizer;
    od;

    # return the random element
    return rnd;
end );


#############################################################################
##
#M  <g> in <G>  . . . . . . . . . . . . . . . . . . . . . . . membership test
##
InstallMethod( \in,
    "for a permutation, and a permutation group",
    true,
    [ IsPerm, IsPermGroup and HasGeneratorsOfGroup], 0,
    function( g, G )
    if g = One( G )  or  g in GeneratorsOfGroup( G )  then
        return true;
    else
        G := StabChainMutable( G );
        return SiftedPermutation( G, g ) = G.identity;
    fi;
end );


#############################################################################
##
#M  ClosureGroup( <G>, <gens>, <options> )  . . . . . .  closure with options
##
##  The following function implements the method and may be called even with
##  incomplete (temp) stabilizer chains.
BindGlobal("DoClosurePrmGp",function( G, gens, options )
    local   C,          # closure of < <G>, <obj> >, result
            P, inpar,   # parent of the closure
            g,          # an element of gens
	    o,		# order
	    newgens,    # new generator list
            chain;      # the stabilizer chain created

    options:=ShallowCopy(options); # options will be overwritten
    # if all generators are in <G>, <G> is the closure
    gens := Filtered( gens, gen -> not gen in G );
    if IsEmpty( gens )  then
        return G;
    fi;

    # is the closure normalizing?
    if Length(gens)=1 and 
       ForAll(gens,i->ForAll(GeneratorsOfGroup(G),j->j^i in G)) then
      g:=gens[1];
      if Size(G)>1 then
	o:=First(Difference(DivisorsInt(Order(g)),[1]),j->g^j in G);
	options.limit:=Size(G)*o;
      else
	o:=First(Difference(DivisorsInt(Order(g)),[1]),j->IsOne(g^j));
	options.limit:=o;
      fi;
      options.size:=options.limit;
      #Print(options.limit,"<<\n");
    fi;

    # otherwise decide between random and deterministic methods
    P := Parent( G );
    inpar := IsSubset( P, gens );
    while not inpar  and  not IsIdenticalObj( P, Parent( P ) )  do
        P := Parent( P );
        inpar := IsSubset( P, gens );
    od;
    if inpar  then
        CopyOptionsDefaults( P, options );
    elif not IsBound( options.random )  then
        options.random := DefaultStabChainOptions.random;
    fi;

    # perhaps <G> is normal in <C> with solvable factor group

#AH 5-feb-96: Disabled (see gap-dev discussion).
#    if DefaultStabChainOptions.tryPcgs
#       and ForAll( gens, gen -> ForAll( GeneratorsOfGroup( G ),
#                   g -> g ^ gen in G ) )  then
#        if inpar  then
#            C := SubgroupNC( P,
#                         Concatenation( GeneratorsOfGroup( G ), gens ) );
#        else
#            C := GroupByGenerators
#                 ( Concatenation( GeneratorsOfGroup( G ), gens ) );
#        fi;
#        pcgs:=TryPcgsPermGroup( [ C, G ], false, false, false );
#        if IsPcgs( pcgs )  then
#	  chain:=pcgs!.stabChain;
#	  if inpar  then  C := GroupStabChain( P, chain, true );
#		    else  C := GroupStabChain( chain );           fi;
#	  SetStabChainOptions( C, rec( random := options.random ) );
#
#	  UseSubsetRelation( C, G );
#	  return C;
#        fi;
#    fi;

    # make the base of G compatible with options.base
    chain := StructuralCopy( StabChainMutable( G ) );
    if IsBound( options.base )  then
        ChangeStabChain( chain, options.base,
                IsBound( options.reduced ) and options.reduced );
    fi;

    if LargestMovedPoint( Concatenation( GeneratorsOfGroup( G ),
               gens ) ) <= 100  then
        options := ShallowCopy( options );
        options.base := BaseStabChain( chain );
        for g  in gens  do
            if SiftedPermutation( chain, g ) <> chain.identity  then
                StabChainStrong( chain, [ g ], options );
            fi;
        od;
    else
        chain := ClosureRandomPermGroup( chain, gens, options );
    fi;

    newgens:=Concatenation(GeneratorsOfGroup(G),gens);
    if Length(chain.generators)<=Length(newgens) then
      if inpar  then  C := GroupStabChain( P, chain, true );
		else  C := GroupStabChain( chain );           fi;
    else
      if inpar  then  C := SubgroupNC(P,newgens);
		else  C := Group( newgens,One(G) );           fi;
      SetStabChainMutable(C,chain);
    fi;
    SetStabChainOptions( C, rec( random := options.random ) );

    UseSubsetRelation( C, G );
    return C;
end );

BindGlobal("PG_EMPTY_OPT",rec());

InstallOtherMethod( ClosureGroup,"permgroup, elements, options",
    true,
    [ IsPermGroup, IsList and IsPermCollection, IsRecord ], 0,
    DoClosurePrmGp);

InstallOtherMethod( ClosureGroup, "empty list",true,
        [ IsPermGroup, IsList and IsEmpty ], 0,
function( G, nogens )
    return G;
end );

InstallMethod( ClosureGroup, "permgroup, element",true,
  [ IsPermGroup, IsPerm ], 0,
function( G, g )
  return DoClosurePrmGp( G, [ g ], PG_EMPTY_OPT );
end );

InstallMethod( ClosureGroup, "permgroup, elements",true,
        [ IsPermGroup, IsList and IsPermCollection ], 0,
function( G, gens )
  return DoClosurePrmGp( G, gens, PG_EMPTY_OPT );
end );

InstallOtherMethod( ClosureGroup, "permgroup, element, options",true,
  [ IsPermGroup, IsPerm, IsRecord ], 0,
function( G, g, options )
  return DoClosurePrmGp( G, [ g ], options );
end );

InstallOtherMethod( ClosureGroup, "permgroup, permgroup, options", true,
  [ IsPermGroup, IsPermGroup, IsRecord ], 0,
function( G, H, options )
  return DoClosurePrmGp( G, GeneratorsOfGroup( H ), options );
end );

InstallOtherMethod( ClosureGroup, "empty list and options",true,
        [ IsPermGroup, IsList and IsEmpty, IsRecord ], 0,
    function( G, nogens, options )
    return G;
end );

#############################################################################
##
#M  NormalClosureOp( <G>, <U> ) . . . . . . . . . . . . . . . . in perm group
##
BindGlobal("DoNormalClosurePermGroup",function ( G, U )
    local   N,          # normal closure of <U> in <G>, result
            chain,      # stabilizer chain for the result
            rchain,     # restored version of <chain>, for `VerifySGS'
            options,    # options record for stabilizer chain construction
            gensG,      # generators of the group <G>
            genG,       # one generator of the group <G>
            gensN,      # generators of the group <N>
            genN,       # one generator of the group <N>
            cnj,        # conjugated of a generator of <U>
            random,  k, # values measuring randomness of <chain>
            param,  missing,  correct,  result,  i, one;

    # get a set of monoid generators of <G>
    gensG := GeneratorsOfGroup( G );
    one:= One( G );

    # make a copy of the group to be closed
    N := SubgroupNC( G, GeneratorsOfGroup(U) );
    UseIsomorphismRelation(U,N);
    UseSubsetRelation(U,N);
    SetStabChainMutable( N, StabChainMutable( U ) );
    options := ShallowCopy( StabChainOptions( U ) );
    if IsBound( options.random )  then  random := options.random;
                                  else  random := 1000;            fi;
    options.temp   := true;

    # make list of conjugates to be added to N
    repeat
        gensN := [  ];
        for i  in [ 1 .. 10 ]  do 
            genG := SCRRandomSubproduct( gensG, One( G ) );
            cnj  := SCRRandomSubproduct( Concatenation(
                        GeneratorsOfGroup( N ), gensN ), One( G ) ) ^ genG;
            if not cnj in N  then 
                Add( gensN, cnj );
            fi;
        od;
        if not IsEmpty( gensN )  then
           N := ClosureGroup( N, gensN, options );
        fi;
    until IsEmpty( gensN );
    
    # Guarantee that all conjugates are in the normal  closure: Loop over all
    # generators of N
    gensN := ShallowCopy( GeneratorsOfGroup( N ) );
    for genN  in gensN  do

        # loop over the generators of G
        for genG  in gensG  do

            # make sure that the conjugated element is in the closure
            cnj := genN ^ genG;
            if not cnj in N  then
                N := ClosureGroup( N, [ cnj ], options );
                Add( gensN, cnj );
            fi;

        od;

    od;
    
    # Verify the stabilizer chain.
    chain := StabChainMutable( N );
    if not IsBound( chain.orbits )  then
        if IsBound( chain.aux )  then
            chain := SCRRestoredRecord( chain );
        fi;
        result := chain.identity;
    elif random = 1000  then
        missing := chain.missing;
        correct := chain.correct;
        rchain  := SCRRestoredRecord( chain );
        result  := VerifySGS( rchain, missing, correct );
        if not IsPerm(result) then 
            repeat 
                result := SCRStrongGenTest2(chain,[0,0,1,10/chain.diam,0,0]);
            until result <> one;
        fi;
        chain := rchain;
    else
        k := First([0..14],x->(3/5)^x <= 1-random/1000);
        if IsBound(options.knownBase) then 
            param := [k,4,0,0,0,0];
        else
            param := [QuoInt(k,2),4,QuoInt(k+1,2),4,50,5];
        fi;
        if options.random <= 200 then 
            param[2] := 2;
            param[4] := 2;
        fi;
        result := SCRStrongGenTest( chain, param, chain.orbits,
                          chain.basesize, chain.base,
                          chain.correct, chain.missing );
        if result = chain.identity  and  not chain.correct  then
            result := SCRStrongGenTest2( chain, param );
        fi;
        chain := SCRRestoredRecord( chain );
    fi;
    SetStabChainMutable( N, chain );
    if result <> chain.identity  then
        N := ClosureGroup( N, [ result ] );
    fi;

    # give N the proper randomness
    StabChainOptions(N).random:=random;
    
    # return the normal closure
    UseSubsetRelation( N, U );
    return N;
end );

InstallMethod( NormalClosureOp, "subgroup of perm group",
  IsIdenticalObj, [ IsPermGroup, IsPermGroup ], 0,
function(G,U)
  # test applicability
  if not IsSubset(G,U) then
    TryNextMethod();
  fi;
  return DoNormalClosurePermGroup(G,U);
end);

#############################################################################
##
#M  ConjugateGroup( <G>, <g> )  . . . . . . . . . . . .  of permutation group
##
InstallMethod( ConjugateGroup, "<P>, <g>", true, [ IsPermGroup, IsPerm ], 0,
function( G, g )
local   H,  S;
    
    H := GroupByGenerators( OnTuples( GeneratorsOfGroup( G ), g ), One( G ) );
    if HasStabChainMutable(G) then
      S := EmptyStabChain( [  ], One( H ) );
      ConjugateStabChain( StabChainMutable( G ), S, g, g );
      SetStabChainMutable( H, S );
    elif HasSize(G) then
      SetSize(H,Size(G));
    fi;
    UseIsomorphismRelation( G, H );
    return H;
end );

#############################################################################
##
#M  CommutatorSubgroup( <U>, <V> )  . . . . . . . . . . . . . for perm groups
##
InstallMethod( CommutatorSubgroup, "permgroups", IsIdenticalObj,
        [ IsPermGroup, IsPermGroup ], 0,
    function( U, V )
    local   C,       # the commutator subgroup
            CUV,     # closure of U,V
            doneCUV, # boolean; true if CUV is computed
            u,       # random subproduct of U.generators
            v,       # random subproduct of V.generators 
            comm,    # commutator of u,v
            list,    # list of commutators
            i;       # loop variable

    # [ <U>, <V> ] = normal closure of < [ <u>, <v> ] >.
    C := TrivialSubgroup( U );
    doneCUV := false;
    
    # if there are lot of generators, use random subproducts
    if Length( GeneratorsOfGroup( U ) ) *
       Length( GeneratorsOfGroup( V ) ) > 10  then
        repeat
           list := [];
           for i in [1..10] do
               u := SCRRandomSubproduct( GeneratorsOfGroup( U ), One( U ) );
               v := SCRRandomSubproduct( GeneratorsOfGroup( V ), One( V ) );
               comm := Comm( u, v );
               if not comm in C then
                   Add( list, comm ) ;
               fi;
           od;
           if Length(list) > 0 then
              C := ClosureGroup( C, list, rec( random := 0,
                                                  temp := true ) );
              if not doneCUV then
                  CUV := ClosureGroup( U, V );
                  doneCUV := true;
              fi;
              # force closure test
	      StabChainOptions(C).random:=DefaultStabChainOptions.random;
              C := DoNormalClosurePermGroup( CUV, C );
           fi;
        until IsEmpty( list );
    fi;

    # do the deterministic method; it will also check correctness 
    list := [];
    for u  in GeneratorsOfGroup( U )  do
        for v  in GeneratorsOfGroup( V )  do
           comm := Comm( u, v );
           if not comm in C then 
               Add( list, comm );
           fi;
        od;
    od;
    if not IsEmpty( list )  then
        C := ClosureGroup( C, list, rec( random := 0,
                                           temp := true ) );
        if not doneCUV then
           CUV := ClosureGroup( U, V );
           doneCUV := true;
        fi;
	# force closure test
	StabChainOptions(C).random:=DefaultStabChainOptions.random;
        C := DoNormalClosurePermGroup( CUV, C );
    fi;

    return C;
end );

#############################################################################
##
#M  DerivedSubgroup( <G> )  . . . . . . . . . . . . . .  of permutation group
##
InstallMethod( DerivedSubgroup,"permgrps",true, [ IsPermGroup ], 0,
    function( G )
    local   D,          # derived subgroup of <G>, result
            g, h,       # random subproducts of generators 
            comm,       # their commutator
            list,       # list of commutators
            i,j;  # loop variables

    # find the subgroup generated by the commutators of the generators
    D := TrivialSubgroup( G );

    # if there are >4 generators, use random subproducts
    if Length( GeneratorsOfGroup( G ) ) > 4 then 
        repeat 
            list := [];
            for i in [1..10] do
                g := SCRRandomSubproduct( GeneratorsOfGroup( G ), One( G ) );
                h := SCRRandomSubproduct( GeneratorsOfGroup( G ), One( G ) );
                comm := Comm( g, h );
                if not comm in D then  
                   Add( list, comm );
                fi;
            od;
            if Length(list) > 0 then 
              D := ClosureGroup(D,list,rec( random := 0,
                                               temp := true ) );
              # force closure test
	      if IsBound(StabChainOptions(G).random) then
		StabChainOptions(D).random:=StabChainOptions(G).random;
	      else
		StabChainOptions(D).random:=DefaultStabChainOptions.random;
	      fi;
               D := DoNormalClosurePermGroup( G, D );
            fi;
        until list = [];
    fi;

    # do the deterministic method; it will also check random result
    list := [];
    for i in [ 2 .. Length( GeneratorsOfGroup( G ) ) ]  do
         for j  in [ 1 .. i - 1 ]  do
             comm := Comm( GeneratorsOfGroup( G )[i],
                           GeneratorsOfGroup( G )[j] );
             if not comm in D then 
                 Add( list, comm );
             fi;
         od;
    od;
    if not IsEmpty( list )  then
        D := ClosureGroup(D,list,rec( random := 0,
                                        temp := true ) );
	# give D a proper randomness
	if IsBound(StabChainOptions(G).random) then
	  StabChainOptions(D).random:=StabChainOptions(G).random;
	else
	  StabChainOptions(D).random:=DefaultStabChainOptions.random;
	fi;

        D := DoNormalClosurePermGroup( G, D );
    fi;

    return D;
end );

#############################################################################
##
#M  IsSimpleGroup( <G> )  . . . . . . . test if a permutation group is simple
##
##  This  is  a most interesting function.   It  tests whether  a permutation
##  group is  simple  by testing whether the group is  perfect and then  only
##  looking at the size of the group and the degree of a primitive operation.
##  Basically  it uses  the O'Nan--Scott theorem, which gives  a pretty clear
##  description of perfect primitive groups.  This algorithm is described  in
##  William M. Kantor,
##  Finding Composition Factors of Permutation Groups of Degree $n\leq 10^6$,
##  J. Symbolic Computation, 12:517--526, 1991.
##
InstallMethod( IsSimpleGroup,"for permgrp", true, [ IsPermGroup ], 0,
    function ( G )
    local   D,          # operation domain of <G>
            hom,        # transitive constituent or blocks homomorphism
            d,          # degree of <G>
            n, m,       # $d = n^m$
            simple,     # list of orders of simple groups
            transperf,  # list of orders of transitive perfect groups
            s, t;       # loop variables

    # if <G> is the trivial group, it is not simple
    if IsTrivial( G )  then
        return false;
    fi;

    # first find a transitive representation for <G>
    D := Orbit( G, SmallestMovedPoint( G ) );
    if not IsEqualSet( MovedPoints( G ), D )  then
        hom := ActionHomomorphism( G, D,"surjective" );
        if Size( G ) <> Size( Image( hom ) )  then
            return false;
        fi;
        G := Image( hom );
    fi;

    # next find a primitive representation for <G>
    D := Blocks( G, MovedPoints( G ) );
    while Length( D ) <> 1  do
        hom := ActionHomomorphism( G, D, OnSets,"surjective" );
        if Size( G ) <> Size( Image( hom ) )  then
            return false;
        fi;
        G := Image( hom );
        D := Blocks( G, MovedPoints( G ) );
    od;

    # compute the degree $d$ and express it as $d = n^m$
    D := MovedPoints( G );
    d := Length( D );
    n := SmallestRootInt( d );
    m := LogInt( d, n );
    if 10^6 < d  then
        Error("cannot decide whether <G> is simple or not");
    fi;

    # if $G = C_p$, it is simple
    if    IsPrimeInt( Size( G ) )  then
        return true;

    # if $G = A_d$, it is simple (unless $d < 5$)
    elif  IsNaturalAlternatingGroup(G)  then
        return 5 <= d;

    # if $G = S_d$, it is not simple (except $S_2$)
    elif  IsNaturalSymmetricGroup(G) then
        return 2 = d;

    # if $G$ is not perfect, it is not simple (unless $G = C_p$, see above)
    elif  Size( DerivedSubgroup( G ) ) < Size( G )  then
        return false;

    # if $\|G\| = d^2$, it is not simple (Kantor's Lemma 4)
    elif  Size( G ) = d ^ 2  then
        return false;

    # if $d$ is a prime, <G> is simple
    elif  IsPrimeInt( d )  then
        return true;

    # if $G = U(4,2)$, it is simple (operation on 27 points)
    elif  d = 27 and Size( G ) = 25920  then
        return true;

    # if $G = PSL(n,q)$, it is simple (operations on prime power points)
    elif  (  (d =      8 and Size(G) = (7^3-7)/2          )  # PSL(2,7)
          or (d =      9 and Size(G) = (8^3-8)            )  # PSL(2,8)
          or (d =     32 and Size(G) = (31^3-31)/2        )  # PSL(2,31)
          or (d =    121 and Size(G) =        237783237120)  # PSL(5,3)
          or (d =    128 and Size(G) = (127^3-127)/2      )  # PSL(2,127)
          or (d =   8192 and Size(G) = (8191^3-8191)/2    )  # PSL(2,8191)
          or (d = 131072 and Size(G) = (131071^3-131071)/2)  # PSL(2,131071)
          or (d = 524288 and Size(G) = (524287^3-524287)/2)) # PSL(2,524287)
      and IsTransitive( Stabilizer( G, D[1] ), Difference( D, [ D[1] ] ) )
    then
        return true;

    # if $d$ is a prime power, <G> is not simple (except the cases above)
    elif  IsPrimePowerInt( d )  then
        return false;

    # if we don't have at least an $A_5$ acting on the top, <G> is simple
    elif  m < 5  then
        return true;

    # otherwise we must check for some special cases
    else

        # orders of simple subgroups of $S_n$ with primitive normalizer
        simple := [ ,,,,,
          [60,360],,,,                  #  5: A(5), A(6)
          [60,360,1814400],,            # 10: A(5), A(6), A(10)
          [660,7920,95040,239500800],,  # 12: PSL(2,11), M(11), M(12), A(12)
          [1092,43589145600],           # 14: PSL(2,13), A(14)
          [360,2520,20160,653837184000] # 15: A(6), A(7), A(8), A(15)
        ];

        # orders of transitive perfect subgroups of $S_m$
        transperf := [ ,,,,
          [60],                         # 5: A(5)
          [60,360],                     # 6: A(5), A(6)
          [168,2520],                   # 7: PSL(3,2), A(7)
          [168,8168,20160]              # 8: PSL(3,2), AGL(3,2), A(8)
        ];

        # test the special cases (Kantor's Lemma 3)
        for s  in simple[n]  do
            for t  in transperf[m]  do
                if    Size( G ) mod (t * s^m) = 0
                  and (((t * (2*s)^m) mod Size( G ) = 0 and s <> 360)
                    or ((t * (4*s)^m) mod Size( G ) = 0 and s =  360))
                then
                    return false;
                fi;
            od;
        od;

        # otherwise <G> is simple
        return true;

    fi;

end );

#############################################################################
##
#M  IsSolvableGroup( <G> )  . . . . . . . . . . . . . . . .  solvability test
##
InstallMethod( IsSolvableGroup,"for permgrp", true, [ IsPermGroup ], 0,
function(G)
local pcgs;
  pcgs:=TryPcgsPermGroup( G, false, false, true );
  if IsPcgs(pcgs) then
    SetIndicesEANormalSteps(pcgs,pcgs!.permpcgsNormalSteps);
    SetIsPcgsElementaryAbelianSeries(pcgs,true);
    if not HasPcgs(G) then
      SetPcgs(G,pcgs);
    fi;
    if not HasPcgsElementaryAbelianSeries(G) then
      SetPcgsElementaryAbelianSeries(G,pcgs);
    fi;
    return true;
  else
    return false;
  fi;
end);

#############################################################################
##
#M  IsNilpotentGroup( <G> ) . . . . . . . . . . . . . . . . . nilpotency test
##
InstallMethod( IsNilpotentGroup,"for permgrp", true, [ IsPermGroup ], 0,
    G -> IsPcgs(PcgsCentralSeries(G) ) );

#############################################################################
##
#M  PcgsCentralSeries( <G> )
##
InstallOtherMethod( PcgsCentralSeries,"for permgrp", true, [ IsPermGroup ], 0,
function(G)
local pcgs;
  pcgs:=TryPcgsPermGroup( G, true, false, false );
  if IsPcgs(pcgs) then
    SetIndicesCentralNormalSteps(pcgs,pcgs!.permpcgsNormalSteps);
    SetIsPcgsCentralSeries(pcgs,true);
    if not HasPcgs(G) then
      SetPcgs(G,pcgs);
    fi;
  fi;
  return pcgs;
end);

#############################################################################
##
#M  PcgsPCentralSeriesPGroup( <G> )
##
InstallOtherMethod( PcgsPCentralSeriesPGroup,"for permgrp", true, [ IsPermGroup ], 0,
function(G)
local pcgs;
  pcgs:=TryPcgsPermGroup( G, true, true, Factors(Size(G))[1] );
  if IsPcgs(pcgs) then
    SetIndicesPCentralNormalStepsPGroup(pcgs,pcgs!.permpcgsNormalSteps);
    SetIsPcgsPCentralSeriesPGroup(pcgs,true);
  fi;
  if IsPcgs(pcgs) and not HasPcgs(G) then
    SetPcgs(G,pcgs);
  fi;
  return pcgs;
end);

#T AH, 3/26/07: This method triggers a bug with indices. It is not clear
#T form the description why TryPcgs... should also give a drived series. Thus
#T disabled.
# #############################################################################
# ##
# #M  DerivedSeriesOfGroup( <G> ) . . . . . . . . . . . . . . .  derived series
# ##
# InstallMethod( DerivedSeriesOfGroup,"for permgrp", true, [ IsPermGroup ], 0,
#     function( G )
#     local  pcgs,  series;
# 
#     if   (not DefaultStabChainOptions.tryPcgs
#        or ( HasIsSolvableGroup( G )  and  not IsSolvableGroup( G ) ))
#        and not (HasIsSolvableGroup(G) and IsSolvableGroup(G)) then
#         TryNextMethod();
#     fi;
#     
#     pcgs := TryPcgsPermGroup( G, false, true, false );
#     if not IsPcgs( pcgs )  then
#         TryNextMethod();
#     fi;
#     SetIndicesEANormalSteps(pcgs,pcgs!.permpcgsNormalSteps);
#     series := EANormalSeriesByPcgs( pcgs );
#     if not HasDerivedSubgroup( G )  then
#         if Length( series ) > 1  then  SetDerivedSubgroup( G, series[ 2 ] );
#                                  else  SetDerivedSubgroup( G, G );  fi;
#     fi;
#     return series;
# end );

#############################################################################
##
#M  LowerCentralSeriesOfGroup( <G> )  . . . . . . . . .  lower central series
##
InstallMethod( LowerCentralSeriesOfGroup,"for permgrp", true, [ IsPermGroup ], 0,
    function( G )
    local  pcgs,  series;

    if    not DefaultStabChainOptions.tryPcgs
       or HasIsNilpotentGroup( G )  and  not IsNilpotentGroup( G ) 
       and not (HasIsNilpotentGroup(G) and IsNilpotentGroup(G)) then
        TryNextMethod();
    fi;
    
    pcgs := TryPcgsPermGroup( G, true, true, false );
    if not IsPcgs( pcgs )  then
        TryNextMethod();
    fi;
    SetIndicesCentralNormalSteps(pcgs,pcgs!.permpcgsNormalSteps);
    series := CentralNormalSeriesByPcgs( pcgs );
    if not HasDerivedSubgroup( G )  then
        if Length( series ) > 1  then  SetDerivedSubgroup( G, series[ 2 ] );
                                 else  SetDerivedSubgroup( G, G );  fi;
    fi;
    return series;
end );

InstallOtherMethod( ElementaryAbelianSeries, "perm group", true,
 [ IsPermGroup and IsFinite], 
  # we want this to come *before* the method for Pcgs computable groups
  RankFilter(IsPermGroup and CanEasilyComputePcgs and IsFinite)
  -RankFilter(IsPermGroup and IsFinite),
function(G)
local pcgs;
  pcgs:=PcgsElementaryAbelianSeries(G);
  if not IsPcgs(pcgs) then
    TryNextMethod();
  fi;
  return EANormalSeriesByPcgs(pcgs);
end);

#InstallOtherMethod( ElementaryAbelianSeries, fam -> IsIdenticalObj
#        ( fam, CollectionsFamily( CollectionsFamily( PermutationsFamily ) ) ),
#        [ IsList ], 0,
#    function( list )
#    local  pcgs;
#    
#    if    IsEmpty( list )
#       or not DefaultStabChainOptions.tryPcgs
#       or (      HasIsSolvableGroup( list[ 1 ] )
#            and not IsSolvableGroup( list[ 1 ] ) )  then
#        TryNextMethod();
##    fi;
#    
##    pcgs := TryPcgsPermGroup( list, false, false, true );
#    if not IsPcgs( pcgs )  then  return fail;
#                           else  return ElementaryAbelianSeries( pcgs );  fi;
#end );

#############################################################################
##
#M  PCentralSeries( <G>, <p> )  . . . . . . . . . . . . . .  p-central series
##
InstallMethod( PCentralSeriesOp,"for permgrp", true, [ IsPermGroup, IsPosInt ], 0,
function( G, p )
    local  pcgs;

  if Length(Factors(Size(G)))>1 then
    TryNextMethod();
  fi;
  pcgs:=PcgsPCentralSeriesPGroup(G);
  if not IsPcgs(pcgs) then
    TryNextMethod();
  fi;
  return PCentralNormalSeriesByPcgsPGroup(pcgs);
end);

#############################################################################
##
#M  SylowSubgroup( <G>, <p> ) . . . . . . . . . . . . . .  Sylow $p$-subgroup
##
InstallMethod( SylowSubgroupOp,"permutation groups", true,
  [ IsPermGroup, IsPosInt ], 0,
function( G, p )
local   S;
  if not HasIsSolvableGroup(G) and IsSolvableGroup(G) then
    # enforce solvable methods if we detected anew that the group is
    # solvable.
    return SylowSubgroupOp(G,p);
  fi;
  S := SylowSubgroupPermGroup( G, p );
  S := GroupStabChain( G, StabChainMutable( S ) );
  SetIsNilpotentGroup( S, true );
  if Size(S) > 1 then
    SetIsPGroup( S, true );
    SetPrimePGroup( S, p );
  fi;
  return S;
end );

InstallGlobalFunction( SylowSubgroupPermGroup, function( G, p )
    local   S,          # <p>-Sylow subgroup of <G>, result
            q,          # largest power of <p> dividing the size of <G>
            D,          # domain of operation of <G>
            O,          # one orbit of <G> in this domain
            B,          # blocks of the operation of <G> on <D>
            f,          # action homomorphism of <G> on <O> or <B>
            T,          # <p>-Sylow subgroup in the image of <f>
            g, g2,      # one <p> element of <G>
            C, C2;      # centralizer of <g> in <G>
    
    # get the size of the <p>-Sylow subgroup
    q := 1;  while Size( G ) mod (q * p) = 0  do q := q * p;  od;

    # handle trivial subgroup
    if   q = 1  then
        return TrivialSubgroup( G );
    fi;
    # go down in stabilizers as long as possible
    S := StabChainMutable( G );
    while Length( S.orbit ) mod p <> 0  do
        S := S.stabilizer;
    od;
    G := GroupStabChain( G, S, true );

    # handle full group
    if q = Size( G )  then
        return G;
    fi;

    # handle <p>-Sylow subgroups of size <p>
    if q = p  then
        repeat g := Random( G );  until Order( g ) mod p = 0;
        g := g ^ (Order( g ) / p);
        return SubgroupNC( G, [ g ] );
    fi;

    # if the group is not transitive work with the transitive constituents
    D := MovedPoints( G );
    if not IsTransitive( G, D )  then
	Info(InfoGroup,1,"PermSylow: transitive");
        S := G;
        D := ShallowCopy( D );
        while q < Size( S )  do
            O := Orbit( S, D[1] );
            f := ActionHomomorphism( S, O,"surjective" );
            T := SylowSubgroupPermGroup( Range( f ), p );
            S := PreImagesSet( f, T );
            SubtractSet( D, O );
        od;
        return S;
    fi;

    # if the group is not primitive work in the image first
    B := Blocks( G, D );
    if Length( B ) <> 1  then
	Info(InfoGroup,1,"PermSylow: blocks");
        f := ActionHomomorphism( G, B, OnSets,"surjective" );
        T := SylowSubgroupPermGroup( Range( f ), p );
        if Size( T ) < Size( Range( f ) )  then
            T := PreImagesSet( f, T );
            S := SylowSubgroupPermGroup( T , p) ;
            return S;
        fi;
    fi;

    # find a <p> element whose centralizer contains a full <p>-Sylow subgroup
    Info(InfoGroup,1,"PermSylow: element");
    repeat g := Random( G );  until Order( g ) mod p = 0;
    g := g ^ (Order( g ) / p);
    C := Centralizer( G, g );
    while GcdInt( q, Size( C ) ) < q  do
        repeat g2 := Random( C );  until Order( g2 ) mod p = 0;
        g2 := g2 ^ (Order( g2 ) / p);
        C2 := Centralizer( G, g2 );
        if GcdInt( q, Size( C ) ) < GcdInt( q, Size( C2 ) )  then
            C := C2;  g := g2;
        fi;
    od;

    # the centralizer acts on the cycles of the <p> element
    B := List( Cycles( g, D ), Set );
    Info(InfoGroup,1,"PermSylow: cycleaction");
    f := ActionHomomorphism( C, B, OnSets,"surjective" );
    T := SylowSubgroupPermGroup( Range( f ), p );
    S := PreImagesSet( f, T );
    return S;

end );

#############################################################################
##
#M  Socle( <G> )  . . . . . . . . . . .  socle of primitive permutation group
##
InstallMethod( Socle,"for permgrp", true, [ IsPermGroup ], 0,
    function( G )
    local   Omega,  deg,  shortcut,  coll,  d,  m,  c,  ds,  L,  z,  ord,
            p,  i;
    
    Omega := MovedPoints( G );
    
    if not IsPrimitive( G, Omega )  then
        TryNextMethod();
    fi;
    
    # Affine groups first.
    L := Earns( G, Omega );
    if L <> fail  then
        return L;
    fi;
    
    deg := Length( Omega );
    if deg >= 6103515625  then
        TryNextMethod();
    elif deg < 12960000  then
        shortcut := true;
        if deg >= 3125  then
            coll := Collected( FactorsInt( deg ) );
            d := Gcd( List( coll, c -> c[ 2 ] ) );
            if d mod 5 = 0  then
                m := 1;
                for c  in coll  do
                    m := m * c[ 1 ] ^ ( c[ 2 ] / d );
                od;
                if m >= 5  then
                    shortcut := false;
                fi;
            fi;
        fi;
        if shortcut  then
            ds := DerivedSeriesOfGroup( G );
            return ds[ Length( ds ) ];
        fi;
    fi;
    
    coll := Collected( FactorsInt( Size( G ) ) );
    if deg < 78125  then
        p := coll[ Length( coll ) ][ 1 ];
    else
        i := Length( coll );
        while coll[ i ][ 2 ] = 1  do
            i := i - 1;
        od;
        p := coll[ i ][ 1 ];
    fi;
    repeat
        repeat
            z := Random( G );
            ord := Order( z );
        until ord mod p = 0;
        z := z ^ ( ord / p );
    until Index( G, Centralizer( G, z ) ) mod p <> 0;
    L := NormalClosure( G, SubgroupNC( G, [ z ] ) );
    if deg >= 78125  then
        ds := DerivedSeriesOfGroup( L );
        L := ds[ Length( ds ) ];
    fi;
    if IsSemiRegular( L, Omega )  then
        L := ClosureSubgroup( L, Centralizer( G, L ) );
    fi;
    return L;
end );
    
#############################################################################
##
#M  FrattiniSubgroup( <G> ) . . . . . . . . . . . .  for permutation p-groups
##
InstallMethod( FrattiniSubgroup,"for permgrp", true, [ IsPermGroup ], 0,
    function( G )
    local   fac,  p,  l,  k,  i,  j;

    fac := Set( FactorsInt( Size( G ) ) );
    if Length( fac ) > 1  then
        TryNextMethod();
    elif fac[1]=1 then
      return G;
    fi;
    p := fac[ 1 ];
    l := GeneratorsOfGroup( G );
    k := [ l[1]^p ];
    for i  in [ 2 .. Length(l) ]  do
        for j  in [ 1 .. i-1 ]  do
            Add(k, Comm(l[i], l[j]));
        od;
        Add(k, l[i]^p);
    od;
    return SolvableNormalClosurePermGroup( G, k );
end );

#############################################################################
##
#M  ApproximateSuborbitsStabilizerPermGroup(<G>,<pnt>) . . . approximation of
##  the orbits of Stab_G(pnt) on all points of the orbit pnt^G. (As not
##  all schreier generators are used, the results may be the orbits of a
##  subgroup.)
##
InstallGlobalFunction( ApproximateSuborbitsStabilizerPermGroup,
    function(G,punkt)
    local one, orbit, trans, stab, gen, pnt, img, norb, no, i, j, changed,
          processStabGen, StabGenAvailable, stgp, orblen, gens;

    one:= One( G );
  if HasStabChainMutable( G ) and punkt in StabChainMutable(G).orbit then
# if we already have a stab chain and bother computing the proper
# stabilizer, a trivial orbit algorithm seems best.
    return Concatenation([[punkt]],
              OrbitsDomain(Stabilizer(G,punkt),
	                   Difference(MovedPoints(G),[punkt])));
#    orbit:=Difference(StabChainMutable(G).orbit,[punkt]);
#    stab:=Stabilizer(G,punkt));
#    # catch trivial case
#    if Length(stab)=0 then
#      stab:=[ one ];
#    fi;
#    stgp:=1;
#
#    processStabGen:=function()
#      stgp:=stgp+1;
#      if stgp>Length(stab) then
#	StabGenAvailable:=false;
#      fi;
#      return stab[stgp-1];
#    end;

  else
    # compute orbit and transversal
    orbit := [ punkt ];
    trans := [];
    trans[ punkt ] := one;
    stab:=[];
    for pnt  in orbit  do
      for gen  in GeneratorsOfGroup(G)  do
        img:=pnt/gen;
	if not IsBound( trans[ img ] )  then
	  Add( orbit, img );
	  trans[ img ] := gen;
	fi;
      od;
    od;
    orblen:=Length(orbit);
    orbit:=orbit{[2..Length(orbit)]};

    processStabGen:=function()
      i:=orblen; 
      gen:= one;
      while 1<=i do
	gen:=gen*Random(GeneratorsOfGroup(G));
	i:=QuoInt(i,2);
      od;
      while punkt^gen<>punkt do
        gen:=gen*trans[punkt^gen];
      od;
      return gen;
    end;

  fi;

  StabGenAvailable:=true;
  # as the first approximation take a subgroup generated by three nontrivial
  # Schreier generators
  changed:=0;
  gens:=[];
  while changed<=3 and StabGenAvailable do
    gen:=processStabGen();
    if gen = one then
      if Random([1..10])>1 then
        changed:=changed+1;
      fi;
    else
      changed:=changed+1;
      Add(gens,gen);
    fi;
  od;
  if Length(gens)=0 then
    orbit:=List(orbit,i->[i]);
  else
    orbit:=OrbitsPerms(gens,orbit);
  fi;
  orbit:=List(orbit,i->Immutable(Set(i)));

  #trivial case
  if Length(orbit)=0 then
    changed:=infinity;
  else
    changed:=0;
  fi;

  # compute suborbits, fusing iteratively under generators
  while changed<=2 and StabGenAvailable do

    gen:=processStabGen();

    if gen = one then
      if Random([1..10])>1 then
        changed:=changed+1;
      fi;
    else
      changed:=changed+1;

      # fuse under gen-action
      norb:=[];
      for i in [1..Length(orbit)] do
	if IsBound(orbit[i]) then
	  no:=ShallowCopy(orbit[i]);
	  Unbind(orbit[i]);
	  i:=i+1;
	  for pnt in no do
	    img:=pnt^gen;
	    if not img in no then
	      img:=First([i..Length(orbit)],j->IsBound(orbit[j])
			 and img in orbit[j]);
	      if img<>fail then
		changed:=0;
		UniteSet(no,orbit[img]);
		Unbind(orbit[img]);
	      fi;
	    fi;
	  od;
	  Add(norb,Immutable(no));
	fi;
      od;
      orbit:=norb;
    fi;
  od;

  return Concatenation([[punkt]],orbit);
end );

#############################################################################
##
#M  AllBlocks . . . Representatives of all block systems
##
InstallMethod( AllBlocks, "generic", [ IsPermGroup ],
function(g)
local gens, one, dom,DoBlocks,pool,subo;

  DoBlocks:=function(b)
  local bl,bld,i,t,n;
    # find representatives of all potential suborbits
    bld:=List(Filtered(subo,i->not ForAny(b,j->j in i)),i->i[1]);
    bl:=[];
    if not IsPrime(Length(dom)/Length(b)) then
      for i in bld do
        t:=Union(b,[i]);
        n:=Blocks(g,dom,t);
        if Length(n)>1 and #ok durch pool:ForAll(Difference(n[1],t),j->j>i) and
           not n[1] in pool then
          t:=n[1];
          Add(pool,t);
          bl:=Concatenation(bl,[t],DoBlocks(t));
        fi;
      od;
    fi;
    return bl;
  end;

    one:= One( g );
  gens:= Filtered( GeneratorsOfGroup(g), i -> i <> one );
  if Length(gens)=0 then return []; fi;
  subo:=ApproximateSuborbitsStabilizerPermGroup(g,
          Minimum(List(gens,SmallestMovedPoint)));
  dom:=Union(subo);
  pool:=[];
  return DoBlocks(subo[1]);
end);

#############################################################################
##
#F  SignPermGroup
##
InstallGlobalFunction( SignPermGroup, function(g)
  if ForAll(GeneratorsOfGroup(g),i->SignPerm(i)=1) then
    return 1;
  else
    return -1;
  fi;
end );

CreateAllCycleStructures := function(n)
local i,j,l,m;
  l:=[];
  for i in Partitions(n) do
    m:=[];
    for j in i do
      if j>1 then
        if IsBound(m[j-1]) then
          m[j-1]:=m[j-1]+1;
        else 
          m[j-1]:=1;
        fi;
      fi;
    od;
    Add(l,m);
  od;
  return l;
end;

#############################################################################
##
#F  CycleStructuresGroup(G)  list of all cyclestructures occ. in a perm group
##
InstallGlobalFunction( CycleStructuresGroup, function(g)
local l,m,i, one;
  l:=CreateAllCycleStructures(Length(MovedPoints(g)));
  m:=List([1..Length(l)-1],i->0);
    one:= One( g );
  for i in ConjugacyClasses(g) do
    if Representative(i) <> one then
      m[Position(l,CycleStructurePerm(Representative(i)))-1]:=1;
    fi;
  od;
  return m;
end );

#############################################################################
##
#M  SmallGeneratingSet(<G>) 
##
InstallMethod(SmallGeneratingSet,"random and generators subset, randsims",true,
  [IsPermGroup],0,
function (G)
local  i, j, U, gens,o,v,a,sel;

  # remove obvious redundancies
  gens := ShallowCopy(Set(GeneratorsOfGroup(G)));
  o:=List(gens,Order);
  SortParallel(o,gens,function(a,b) return a>b;end);
  sel:=Filtered([1..Length(gens)],x->o[x]>1);
  for i in [1..Length(gens)] do
    if i in sel then
      U:=Filtered(sel,x->x>i and IsInt(o[i]/o[x])); # potential powers
      v:=[];
      for j in U do
	a:=SmallestMovedPoint(gens[j]);
	if IsSubset(OrbitPerms([gens[i]],a),OrbitPerms([gens[j]],a)) then
	  Add(v,j);
	fi;
      od;
      # v are the possible powers
      for j in v do
	a:=gens[i]^(o[i]/o[j]);
	if ForAny(Filtered([1..o[j]],z->Gcd(z,o[j])=1),x->a^x=gens[j]) then
	  RemoveSet(sel,j);
	fi;
      od;
    fi;
  od;
  gens:=gens{sel};

  # try pc methods first
  if Length(MovedPoints(G))<1000 and HasIsSolvableGroup(G) 
     and IsSolvableGroup(G) and Length(gens)>3 then
    return MinimalGeneratingSet(G);
  fi;

  if Length(gens)>2 then
    i:=2;
    while i<=3 and i<Length(gens) do
      # try to find a small generating system by random search
      j:=1;
      while j<=5 and i<Length(gens) do
        U:=Subgroup(G,List([1..i],j->Random(G)));
        StabChain(U,rec(random:=1));
        if Size(U)=Size(G) then
          gens:=Set(GeneratorsOfGroup(U));
        fi;
        j:=j+1;
      od;
      i:=i+1;
    od;
  fi;
  i := 1;
  if not IsAbelian(G) then
    i:=i+1;
  fi;
  while i < Length(gens)  do
    # random did not improve much, try subsets
    U:=Subgroup(G,gens{Difference([1..Length(gens)],[i])});
    if Size(U)<Size(G) then
      i:=i+1;
    else
      gens:=Set(GeneratorsOfGroup(U));
    fi;
  od;
  return gens;
end);

#############################################################################
##
#M  GeneratorsSmallest(<G>) . . . . . . . . . . . . . for permutation groups
##
DeclareGlobalFunction( "GeneratorsSmallestStab" );

InstallGlobalFunction( GeneratorsSmallestStab, function ( S )
    local   gens,       # smallest generating system of <S>, result
            gen,        # one generator in <gens>
            orb,        # basic orbit of <S>
            pnt,        # one point in <orb>
            T,          # stabilizer in <S>
	    o2,img,oo,og; # private orbit algorithm

    # handle the anchor case
    if Length(S.generators) = 0  then
        return [];
    fi;

    # now get the smallest generating system of the stabilizer
    gens := GeneratorsSmallestStab( S.stabilizer );

    # get the sorted orbit (the basepoint will be the first point)
    orb := Set( S.orbit );
    SubtractSet( orb, [S.orbit[1]] );

    # handle the other points in the orbit
    while Length(orb) <> 0  do

        # take the smallest point (coset) and one representative
        pnt := orb[1];
        gen := S.identity;
        while S.orbit[1] ^ gen <> pnt  do
           gen := LeftQuotient( S.transversal[ pnt / gen ], gen );
        od;

        # the next generator is the smallest element in this coset
        T := S.stabilizer;
        while Length(T.generators) <> 0  do
            pnt := Minimum( OnTuples( T.orbit, gen ) );
            while T.orbit[1] ^ gen <> pnt  do
                gen := LeftQuotient( T.transversal[ pnt / gen ], gen );
            od;
            T := T.stabilizer;
        od;

        # add this generator to the generators list and reduce orbit
        Add( gens, gen );

        #SubtractSet( orb,
        #    Orbit( GroupByGenerators( gens, S.identity ), S.orbit[1] ) );

	# here we want to be really fast, so use a private orbit algorithm
	o2:=[S.orbit[1]];
	RemoveSet(orb,S.orbit[1]);
	for oo in o2 do
	  for og in gens do
	    img:=oo^og;
	    if img in orb then
	      Add(o2,img);
	      RemoveSet(orb,img);
	    fi;
	  od;
	od;

    od;

    # return the smallest generating system
    return gens;
end );

InstallMethod(GeneratorsSmallest,"perm group via minimal stab chain",
  [IsPermGroup],
function(G)
  # call the recursive function to do the work
  return GeneratorsSmallestStab(MinimalStabChain(G));
end);

InstallMethod(LargestElementGroup,"perm group via minimal stab chain",
  [IsPermGroup],
  # call the recursive function to do the work
  G -> LargestElementStabChain( MinimalStabChain( G ), One( G ) ) );


#############################################################################
##
#M  KnowsHowToDecompose( <G>, <gens> )
##
InstallMethod(KnowsHowToDecompose,
    "perm group and generators: always true",
    IsIdenticalObj,
    [ IsPermGroup, IsList ],
    ReturnTrue);


#############################################################################
##
#M  ViewObj( <G> )
##
InstallMethod( ViewObj,
    "for a permutation group",
    true,
    [ IsPermGroup and HasGeneratorsOfGroup ], 0,
function(G)
  local gens;
  gens:= GeneratorsOfGroup( G );
  if 30 < Length( gens ) * LargestMovedPoint( G ) / GAPInfo.ViewLength then
    Print("<permutation group");
    if HasSize(G) then
      Print(" of size ",Size(G));
    fi;
    Print( " with ", Length( gens ), " generators>" );
  else
    Print("Group(");
    if IsEmpty( gens ) or ForAll(gens,i->Order(i)=1) then
      ViewObj( One( G ) );
    else
      ViewObj( gens );
    fi;
    Print(")");
  fi;
end);


#############################################################################
##
#M  AsList( <G> ) elements of perm group
##
InstallMethod( AsList,"permgp: AsSSortedList", true,
  [ IsPermGroup ], 0, AsSSortedList );

#############################################################################
##
#M  AsSSortedList( <G> ) elements of perm group
##
InstallMethod( AsSSortedList,"via stabchain",true, [ IsPermGroup ], 0,
function( G )
  return ElementsStabChain( StabChainMutable(G));
end );

InstallMethod( AsSSortedListNonstored,"via stabchain", true, [ IsPermGroup ], 0,
function( G )
  return ElementsStabChain( StabChainMutable(G));
end );

#############################################################################
##
#M  ONanScottType( <G> )
##
InstallMethod(ONanScottType,"primitive permgroups",true,[IsPermGroup],0,
function(G)
local dom,s,cs,t,ts,o;
  dom:=MovedPoints(G);
  if not IsPrimitive(G,dom) then
    Error("<G> must be primitive");
  fi;
  s:=Socle(G);
  if IsAbelian(s) then
    return "1";
  elif IsSimpleGroup(s) then
    return "2";
  elif Length(dom)=Size(s) then
    return "5";
  else
    # now get one simple factor of the socle
    cs:=CompositionSeries(s);
    # if the group is diagonal, the diagonal together with a maximal socle
    # normal subgroup generates the whole socle, so this normal subgroup
    # acts transitively. For product type this is not the case.
    t:=cs[Length(cs)-1];
    if IsTransitive(t,dom) then
      # type 3
      if Length(cs)=3 and IsNormal(G,t) then
	# intransitive on 2 components:
	return "3a";
      else
	return "3b";
      fi;
    else
      # type 4
      ts:=Orbit(G,t);
      if Length(ts)=2 and NormalClosure(G,t)<>s then
        return "4a";
      fi;
      # find the block on the T's first: Those t which keep the orbit are
      # glued together diagonally
      o:=Orbit(t,dom[1]);
      ts:=Filtered(ts,i->i=t or IsSubset(o,Orbit(i,dom[1])));
      if Length(ts)=1 then
        return "4c"; 
      else
	return "4b";
      fi;
    fi;
  fi;
end);

#############################################################################
##
#M  SocleTypePrimitiveGroup( <G> )
##
InstallMethod(SocleTypePrimitiveGroup,"primitive permgroups",true,
  [IsPermGroup],0,function(G)
local s,cs,t,id,r;
  s:=Socle(G);
  cs:=CompositionSeries(s);
  t:=cs[Length(cs)-1];
  id:=IsomorphismTypeInfoFiniteSimpleGroup(t);
  r:=rec(series:=id.series,
	     width:=Length(cs)-1,
	     name:=id.name);
  if IsBound(id.parameter) then
    r.parameter:=id.parameter;
  else
    r.parameter:=id.name;
  fi;
  return r;
end);

# this function is needed for the transitive and primitive groups libraries
BindGlobal("STGSelFunc",function(a,b)
  if IsFunction(b) then
    return b(a);
  elif IsList(b) and not IsString(b) then
   if IsList(a) and a=b then 
     return true;
   fi;
   return a in b;
 else
   return a=b;
 fi;
end);

InstallGlobalFunction(DiagonalSocleAction,function(g,nr)
local ran,d,emb,u,act;
  ran:=[1..nr];
  d:=DirectProduct(List(ran,i->g));
  emb:=List(ran,i->Embedding(d,i));
  u:=List(GeneratorsOfGroup(g),i->Product(ran,j->Image(emb[j],i)));
  u:=SubgroupNC(d,u);
  act:=ActionHomomorphism(d,RightTransversal(d,u),OnRight,"surjective");
  return Image(act,d);
end);

InstallGlobalFunction(ReducedPermdegree,function(g)
  local dom, deg, orb, gdeg, p, ndom, s,hom;
  dom:=MovedPoints(g);
  deg:=Length(dom);
  orb:=ShallowCopy(Orbits(g,dom));
  if Length(orb)=1 then return fail;fi;
  gdeg:=20*LogInt(Size(g),2); # good degree
  Sort(orb,function(a,b) return Length(a)<Length(b);end);
  p:=PositionProperty(orb,i->Length(i)>=gdeg);
  if p=fail then p:=Length(orb);fi;
  ndom:=orb[p];
  if Length(ndom)*2>Length(dom) then return fail;fi;
  s:=Stabilizer(g,orb[p],OnTuples);
  p:=p+1;
  if p>Length(orb) then p:=1;fi;
  while Size(s)>1 do
    while not ForAny(orb[p],j->ForAny(GeneratorsOfGroup(s),k->j^k<>j)) do
      p:=p+1;
      if p>Length(orb) then p:=1;fi;
    od;
    ndom:=Union(ndom,orb[p]);
    s:=Stabilizer(s,orb[p],OnTuples);
  od;
  if Length(ndom)*2>Length(dom) then return fail;fi;
  hom:=ActionHomomorphism(g,ndom,"surjective");
  SetIsInjective(hom,true);
  SetSize(Image(hom),Size(g));
  return hom;
end);

#############################################################################
##
#E
gap-libs 4r6p5-3 / usr / share / gap / lib / grpperm.gi
