Math Solutions 1995 October/Math_Solutions_CD-ROM_Walnut_Creek_October_1995.iso/pc/mac/discrete/lib/permag.g

Text File | 1993-05-05 | 30.5 KB | 922 lines

#############################################################################
##
#A  permag.g                    GAP library                    Heiko Thei"sen
##
#A  @(#)$Id: permag.g,v 3.7 1993/02/10 10:11:14 martin Rel $
##
#Y  Copyright 1990-1993,  Lehrstuhl D fuer Mathematik,  RWTH Aachen,  Germany
##
##  This  file  contains  functions  that  calculate  composition  series for
##  solvable permutation groups and convert such groups into ag groups.
##
#H  $Log: permag.g,v $
#H  Revision 3.7  1993/02/10  10:11:14  martin
#H  moved 'PermGroupOps.CompositionSeries' to "permcser"
#H
#H  Revision 3.6  1993/02/09  14:25:55  martin
#H  made undefined globals local
#H
#H  Revision 3.5  1993/01/25  13:02:33  fceller
#H  changed names of abstract gens to "g"
#H
#H  Revision 3.4  1993/01/15  14:46:54  theissen
#H  fixed a minor bug to allow ag conversion of the trivial group
#H
#H  Revision 3.3  1993/01/12  11:49:06  theissen
#H  implemented the use of a bssgs to save memory
#H
#H  Revision 3.2  1992/12/02  09:01:10  fceller
#H  Initial GAP 3.2 revision
#H
##


#############################################################################
##
#F  MaximalBlocksPGroup( <G>, <D> ) . . . . . . . find a maximal block system
##
MaximalBlocksPGroup := function( G, D )
    local   B;

    D := Blocks( G, D );
    if Length(D) = 1  then
        return D;
    fi;
    B := D;
    while not IsPrime(Length(D))  do
        G := Operation( G, D, OnSets );
        D := Blocks( G, [1..Length(G.operationDomain)] );
        B := List( D, d -> Concatenation( Sublist( B, d ) ) );
    od;
    return B;
end;


#############################################################################
##
#F  OrderFactorGroupElement(<G>,<g>)  . . . . . . . . order of G*g in <G,g>/G
##
OrderFactorGroupElement := function( G, g )
    local   ord, div;

    if g = G.identity  then
        return 1;
    fi;
    ord := Order( G, g );
    for div  in Set(FactorsInt(ord))  do
        while ord mod div = 0 and g^(ord/div) in G  do
            ord := ord / div;
        od;
    od;
    return ord;
end;


#############################################################################
##
#F  InsertStabChain( <G>, <base>, <sgs> ) . . . . . insert a stabilizer chain
##
##  This function inserts a stabilizer chain in the group record of <G> along
##  the base <base> using the strong generating set <sgs>.
##
InsertStabChain := function( G, base, sgs )
    local   pt,  orb;

    if 0 < Length(base)  then
        pt  := base[1];
        orb := PermGroupOps.OrbitTransversal( G, pt );
        G.orbit       := orb.orbit;
        G.transversal := orb.transversal;
        G.stabilizer  := rec( identity   := G.identity,
                              generators := Filtered( sgs, s->pt^s=pt ) );
        InsertStabChain( G.stabilizer,
                         Sublist( base, [ 2 .. Length(base) ] ),
                         G.stabilizer.generators );
    fi;
end;


#############################################################################
##
#F  ClosureNormalizingElementPermGroup( <H>, <y> )  .  closure of <H> and <y>
##
##  This function uses an idea  of C. Sims  to extend a permutation group <H>
##  with a normalizing element <y>, also extending the stabilizer chain.  The
##  method is faster than the standard algorithm because a generating set for
##  each  stabilizer in the chain can easily be obtained and the extension of
##  the basic  orbits does  not require a full orbit  algorithm.  Besides, if
##  '<H>.bssgs'  is  a base-strong  subnormal generating system upon  call of
##  this procedure the  result will  also have  bound a system '.bssgs'  with
##  this property.
##
##  <y> must be an element of the parent of <H> and must normalize <H>.  This
##  is *not* checked.
##
ClosureNormalizingElementPermGroup := function( H, y )
    local  H,  G,  z,  newgens,  orbit,  pnt,  o,  w,  m,  p;

    # Do not change the original argument.
    H := Copy( H );
    G := H;
    z := y;
    if IsBound( H.bssgs )  then
    newgens := [  ];
    fi;

    while z <> H.identity  do

    # If necessary, extend the base.
    if not IsBound( G.orbit )  then
        G.orbit                     := [ LargestMovedPointPerm( z ) ];
        G.transversal            := [  ];
        G.transversal[ G.orbit[1] ] := H.identity;
        G.stabilizer           := rec( identity   := H.identity,
                        generators := [  ] );
        fi;

    # Extend the orbit with the new generator.
    orbit := ShallowCopy( G.orbit );
    pnt := orbit[ 1 ];
    m := 1;
    w := z;
    while not pnt/w in orbit  do
        for o  in orbit  do
        Add( G.orbit, o/w );
        G.transversal[ o/w ] := z;
        od;
        m := m+1;
        w := w*z;
    od;

          # Put in <z> as generator in the stabilizer chain.
    if not z in G.generators  then
        Add( G.generators, z );
    fi;

    # If a bssgs is  present  and the orbit is  properly extended, put in
    # <z> and its  powers as new polycyclic generators, otherwise let <z>
    # = <w>.
    if m > 1  then
        if IsBound( newgens )  then
        for p  in FactorsInt( m )  do
            Add( newgens, z );
            z := z^p;
        od;
        else
        z := w;
        fi;
    fi;

    # Now <z>  = <w>. Find a cofactor to <z> such that  the product fixes
    # <pnt>.
    while pnt^z <> pnt  do
        z := z * G.transversal[ pnt^z ];
    od;

    # Go down one step in the stabilizer chain.
    G := G.stabilizer;

    od;

    # Extend the bssgs of <H> by <newgens> if it is present.
    if IsBound( newgens )  then
    H.bssgs := Concatenation( newgens, H.bssgs );
    fi;

    return H;

end;

#############################################################################
##

#F  PermGroupOps.OrbitTransversal( <G>, <d> ) . . . . . orbit and transversal
##
PermGroupOps.OrbitTransversal := function( G, d )
    local  max,  orb,  new,  pnt,  gen,  img;

    orb := rec( orbit := [d], transversal := [] );
    orb.transversal[d] := G.identity;
    max := 0;
    for gen  in G.generators  do
        if gen <> G.identity and max < LargestMovedPointPerm(gen)  then
            max := LargestMovedPointPerm(gen);
        fi;
    od;
    if d in [1..max]   then
        new := BlistList( [1..max], [1..max] );
        new[d] := false;
        for pnt  in orb.orbit  do
            for gen  in G.generators  do
                img := pnt/gen;
                if new[img]  then
                    Add( orb.orbit, img );
                    orb.transversal[img] := gen;
                    new[img] := false;
                fi;
            od;
        od;
    fi;
    return orb;

end;


#############################################################################
##
#F  PermGroupOps.TryElementaryAbelianSeries(<G>)  . . . . try to build an eas
##
##  This  function starts with  the series (<K>_1) where <K>_1 is the trivial
##  subgroup     of     <G>     and    extends     this     series    calling
##  'ExtendElementaryAbelianSeriesPermGroup' with  each generator of  <G>  in
##  turn. If <G> is solvable this will result in an elementary abelian series
##  (<G>,...,<K>_1),  otherwise 'ExtendElementaryAbelianSeriesPermGroup' will
##  return  'false' because <G>.upperBoundDerivedLength, which  is calculated
##  according to J.D.  Dixon, is exceeded.  If <G> is solvable, this function
##  binds and  returns  an elementary  abelian  series, otherwise it  returns
##  'false'.  The flag '<G>.isSolvable' is also set by this function.
##
PermGroupOps.TryElementaryAbelianSeries := function( G )
    local   N,  y,  base,  n,  c,  log;

    # if <G> is trivial return
    if G.generators = [  ]  then
        G.elementaryAbelianSeries := [ G ];
        G.isSolvable            := true;
        G.bssgs                   := [  ];
        return G.elementaryAbelianSeries;
    else
        n := PermGroupOps.LargestMovedPoint(G);
    fi;

    # compute an upper bound for the derived length of <G>, assuming that <G>
    # is solvable. According to Dixon (1968) this is (5 log_3(deg(<G>)))/2
    log := 0;
    c   := 1;
    while c < n  do
        log := log + 1;
        c   := c*3;
    od;
    G.upperBoundDerivedLength := 5*log/2;

    # start with the trivial subgroup
    N := TrivialSubgroup(G);
    N.bssgs := [  ];
    InsertStabChain( N, [  ], [  ] );
    G.elementaryAbelianSeries := [ N ];
    for y  in G.generators  do
        N := ExtendElementaryAbelianSeriesPermGroup( G, y, 0 );
        if N = false  then
            G.isSolvable := false;
            Unbind(G.elementaryAbelianSeries);
        Unbind( G.upperBoundDerivedLength );
            return false;
        fi;
    od;

    # copy the information stored in <N> to <G>
    G.bssgs       := N.bssgs;
    G.orbit       := N.orbit;
    G.transversal := N.transversal;
    G.stabilizer  := N.stabilizer;
    G.isSolvable  := true;
    Unbind( G.upperBoundDerivedLength );

    return G.elementaryAbelianSeries;

end;

#############################################################################
##

#F  ExtendElementaryAbelianSeriesPermGroup( <G>, <y>, <der> ) . . . . . local
##
##  This function  assumes that <G>  has bound a  partial  elementary abelian
##  series,  i.e.  a  sequence  (<K>_m,...,<K>_1) of  normal  subgroups  with
##  elementary  abelian  factors.   It   extends  this  series  to  a  series
##  (<K>_n,...,<K>_1) where <K>_n is the  normal closure of <K>_m and <y> and
##  further  normal  subgroups  <K>_{n-1},...,<K>_{m+1} may  be  inserted  to
##  ensure elementary abelian factors.  The method is due to C. Sims and will
##  terminate  if <G> is solvable.  To avoid endless loops the function takes
##  the counter <der>  which must be  0  when this  function is called from a
##  program and which measures the  depth of  the  commutators we are working
##  with.  If <G> is solvable, commutators cannot  have arbitrary depth and a
##  bound on  <der> can  be  calculated  only from  the degree of <G>.   This
##  function    assumes    that    such    a     bound     is    bound     in
##  '<G>.upperBoundDerivedLength' and uses  this  trick to detect whether <G>
##  is not solvable.
##
ExtendElementaryAbelianSeriesPermGroup := function( G, y, der )
    local   N,  M,  U,  Z,  z,  T,  done,  u,  w,  V,  v,  q;

    # if we are too deep in the derived series, then <G> is not solvable
    if der > G.upperBoundDerivedLength  then
        return false;
    fi;

    # try to extend the series
    N := G.elementaryAbelianSeries[1];
    M := N;
    U := [];
    Z := [y];
    for z  in Z  do
        if not z in M  then
            T := U;
            done := false;
            while not done and 0 < Length(T)  do

                # at this point, M=<N,U>
                u := T[1];
                T := Sublist( T, [ 2 .. Length(T) ] );
                w := Comm( u, z );

                # if  <z> and  <u> do  not  commute mod  <N>,  extend <N>  to
                # NormalClosure( <N>,  [<z>,<u>] ),  going  one step  further
                # down the derived series
                if not w in N  then
                    N := ExtendElementaryAbelianSeriesPermGroup(G, w, der+1);
                    if N = false  then
                        return false;
                    fi;

                    # set M := <N,U> for the  new N,  V  <  U  contains those
                    # generators that  are actually  needed in the  extension
                    # process, i.e. M=<N,V>
                    M := N;
                    V := [];
                    for v  in U  do
                        if not v in M  then
                            M := ClosureNormalizingElementPermGroup( M, v );
                            AddSet( V, v );
                        else
                            RemoveSet( T, v );
                        fi;
                    od;

                    # restore U such that M = <N,U>
                    U := V;

                    # if z in M then M is an elementary abelian extension  of
                    # N  containing z, so  no  further conjugates need to  be
                    # considered
                    if z in M  then
                        done := true;
                    fi;
                fi;
            od;

            # extend M with z, if still necessary
            if not done  then
                M := ClosureNormalizingElementPermGroup( M, z );
                AddSet( U, z );
            fi;

            # store the conjugates of z
            Append( Z, List( G.generators, g-> z^g ) );
        fi;
    od;

    if U <> [  ]  then

        # Now M/N is abelian, M=<N,U> and U contains elements of fixed order.
        # Refine M>N to have  elementary abelian factors  by inserting powers
        # of the generators in U, thereby extending N.
        q := OrderFactorGroupElement( N, U[1] );
        while not IsPrime(q)  do
            q := q / FactorsInt(q)[1];
            for u  in U  do
                if not u in N  then
                    N := ClosureNormalizingElementPermGroup( N, u^q );
                fi;
            od;
            G.elementaryAbelianSeries := Concatenation(
                                             [N],
                                             G.elementaryAbelianSeries );
        od;

    # Now M/N is  elementary abelian. Extend the generator list  of M  to
    # contain the generators of N.
        M := N;
        for v  in U  do
        if not v in M  then
            M := ClosureNormalizingElementPermGroup( M, v );
        fi;
        od;
        G.elementaryAbelianSeries := Concatenation(
                         [M],
                         G.elementaryAbelianSeries );

    fi;

    return M;

end;


#############################################################################
##
#F  BaseStrongSubnormalGeneratingSetPPermGroup(<G>)  base and sgs for p-group
##
##  For intransitive <G>, this function  uses a  constituent  homomorphism to
##  get  a  base and  sgs  from  base  and  sgs for  kernel and  image.   For
##  imprimitive <G>, it likewise uses a blocks homomorphism to get a base and
##  sgs  from those of kernel  and  image. Finally,  if <G> is transitive and
##  primitive,  it must be isomorphic to <Z>_<p> so base and sgs are obvious.
##  The result  is bound to '<G>.polycyclicGenerators'  and  the  base can be
##  retrieved by 'Base(<G>)'.
##
BaseStrongSubnormalGeneratingSetPPermGroup := function( G )
    local   degree, pi, K, I, Delta, blocks, i, reps, rep, s, t, done;

    # if <G> is trivial return
    if IsTrivial( G ) then
    G.polycyclicGenerators := [  ];
        InsertStabChain( G, [  ], [  ] );
        return;
    fi;

    # find a non-trivial orbit
    degree := PermGroupOps.LargestMovedPoint( G );
    i := 0;
    repeat
        i := i+1;
        Delta := PermGroupOps.OrbitTransversal( G, i );
    until 1 < Length(Delta.orbit);

    # <G> is intransitive
    if Length(Delta.orbit)<degree then
        pi := OperationHomomorphism( G, Operation(G,Delta.orbit) );
        K  := Kernel(pi);
        BaseStrongSubnormalGeneratingSetPPermGroup(K);
        I := Image(pi);
        BaseStrongSubnormalGeneratingSetPPermGroup(I);
        G.polycyclicGenerators := Concatenation(
            List( I.polycyclicGenerators,
                  s -> PreImagesRepresentative(pi,s) ),
            K.polycyclicGenerators );
        InsertStabChain(
            G,
            Concatenation(List(Base(I),i->i^(pi.conperm^-1)),Base(K)),
            G.polycyclicGenerators );

    # <G> is transitive
    else
        blocks := MaximalBlocks( G, [1..degree] );

        # <G> is imprimitive
        if 1 < Length(blocks)  then
            K := Stabilizer( G, blocks[1], OnSets );
            BaseStrongSubnormalGeneratingSetPPermGroup(K);
            i := First( G.generators, gen->not gen in K );
            G.polycyclicGenerators := Concatenation(
        [i],
        K.polycyclicGenerators );
            InsertStabChain( G, Base(K), G.polycyclicGenerators );

        # <G> is primitive, thus isomorphic to <Z>_<p>
        else
            G.polycyclicGenerators := [ G.generators[1] ];
            InsertStabChain( G, [degree], G.polycyclicGenerators );
        fi;
    fi;
end;


#############################################################################
##
#F  ExponentsPermSolvablePermGroup( <G>, <g> [, <start> ] ) . . . . . . . . .
#F  . . . . . .  exponents of <g> as normal word in the pag generators of <G>
##
##  For  a  solvable  permutation  group <G> with bssgs  '<G>.bssgs' and  and
##  element  <g> in <G>, this function determines the  exponent vector  <exp>
##  such     that      '<g>     =     <G>.bssgs[<start>]^exp[<start>]     ...
##  <G>.bssgs[<n>]^exp[<n>]'.  If  a value  for  <start> is  known, it can be
##  given as third argument, otherwise it is defaulted to 1.
##
ExponentsPermSolvablePermGroup := function( arg )
    local  G,  g,  start,  exp,  eps,  h,  H,  pnt,  i;

    # Get the arguments.
    G := arg[ 1 ];
    g := arg[ 2 ];
    if Length( arg ) > 2  then
    start := arg[ 3 ];
    else
    start := 1;
    fi;

    exp := [  ];

    # Mind the offset <start>.
    for i  in [ start .. Length( G.bssgs ) ]  do

        # Find the base level of the <i>-th generator, remove the part of <g>
    # not fixing the earlier basepoints.
    h := g;
         H := G;
        while IsBound( H.orbit )  and  H.orbit[ 1 ] ^ G.bssgs[ i ] =
                       H.orbit[ 1 ]  do
        pnt := H.orbit[ 1 ];
        while pnt ^ h <> pnt  do
        h := h * H.transversal[ pnt^h ];
        od;
            H := H.stabilizer;
      od;

        # Determine the <i>-th exponent.
    eps := 0;
        pnt := H.orbit[ 1 ] ^ h;
    while H.transversal[ pnt ] = G.bssgs[ i ]  do
        eps := eps + 1;
        pnt := pnt / G.bssgs[ i ];
    od;
        exp[ i ] := eps;

    # Remove next factor of <g>.
    g := G.bssgs[ i ] ^ eps mod g;

    od;

    # Return the result.
    return exp;

end;

#############################################################################
##
#F  PcPresentationPermGroup( <G>, <series>, <pgens>, <index>, <isNilp> )  . .
##
##  This function calculates a  pc presentation  for <G> along  a composition
##  series that  refines <series>.   <pgens> must contain the generator  list
##  according  to  the  composition  series  of  <G>.   <series>  can  be  an
##  elementary  abelian series, in which  case  <index> is a  list such  that
##  '<pgens>[<index[i]>]'  ...  '<pgens>[<index[i+1]>-1]' are the  generators
##  of the <i>-th elementary abelian factor in  <series>.  Otherwise <series>
##  must    equal    '<G>.compositionSeries'    and    <index>   must   equal
##  '[1..Length(<G>.compositionSeries)+1]'.
##
##  If <isNilp>  is true,  <series> must equal '<G>.compositionSeries'  which
##  must be a central series of <G>.  In this case, the composition series of
##  the resulting ag group is also a central series.
##
##  This     function    is    called    by    'PermGroupOps.AgGroup'     and
##  'PermGroupOps.PgGroup'.
##
PcPresentationPermGroup := function( G, series, pgens, index, isNilp )
    local   PC,  m,  p,  i,  i2,  n,  n2,  start,  k,  exp,  word,  rel;

    # <PC> will hold the presentation
    m  := Length( pgens );
    PC := rec( relations  := [],
               generators := WordList( m, "g" ) );

    # Find the relations of the p-th powers. Use  the  vector space structure
    # of the elementary abelian factors.
    for i  in [ 1 .. Length(series)-1 ]  do
        p := OrderFactorGroupElement(
                 series[i+1],
                 pgens[index[i]] );
        for n  in [ index[i] .. index[i+1]-1 ]  do
            word := PC.generators[n] ^ p;
        rel  := word^0;
        exp  := ExponentsPermSolvablePermGroup
            ( G, pgens[n]^p, index[i+1] );
            for k  in [ index[i+1] .. m ]  do
        rel := rel * PC.generators[k]^exp[k];
            od;
            Add( PC.relations, word/rel );
        od;
    od;

    # Find the relations of the commutators.
    for i  in [ 1 .. Length(series)-1 ]  do
        for n  in [ index[i] .. index[i+1]-1 ]  do
            for i2  in [ 1 .. i-1 ]  do
                if isNilp then
                    start := n+1;
                else
                    start := index[i2+1];
                fi;
                for n2  in [ index[i2] .. index[i2+1]-1 ]  do
                    word := Comm( PC.generators[n], PC.generators[n2] );
            rel  := word^0;
              exp  := ExponentsPermSolvablePermGroup
                ( G, Comm( pgens[n], pgens[n2] ), start );
                    for k  in [ start .. m ]  do
            rel := rel * PC.generators[k]^exp[k];
                    od;
                    Add( PC.relations, word/rel );
                od;
            od;
            for n2  in [ index[i] .. n-1 ]  do
                word := Comm( PC.generators[n], PC.generators[n2] );
        rel  := word^0;
        exp  := ExponentsPermSolvablePermGroup
                ( G, Comm( pgens[n], pgens[n2] ), index[i+1] );
                for k  in [ index[i+1] .. m ]  do
            rel := rel * PC.generators[k]^exp[k];
                od;
                Add( PC.relations, word/rel );
            od;
        od;
    od;
    return PC;

end;

#############################################################################
##
#F  CompositionSeriesSolvablePermGroup( <G> ) . . . . . .  composition series
##
CompositionSeriesSolvablePermGroup := function( G )
    local  compositionSeries,  N,  elabstep,  i,  pgens,  g,  ord,  cof,
           fac,  newgens, gen;

    # If  a subnormal series (with abelian factors) and polycyclic generators
    # are  given, refine  the  series  to have <Z>_<p>-factors  by adding the
    # polycyclic generators and their powers one by one.
    if IsBound(G.subnormalSeries)  then
        pgens             := G.polycyclicGenerators;
    newgens           := [  ];
        compositionSeries := [  ];
        for i  in [ 1 .. Length(pgens) ]  do
            N := G.subnormalSeries[i+1];
            g := pgens[i];
            ord := OrderFactorGroupElement( N, g );
            cof := 1;
            for fac  in FactorsInt(ord)  do
                if cof = 1  then
                    Add( newgens, g );
                    Add( compositionSeries, G.subnormalSeries[i] );
                    g := g ^ fac;
                elif cof < ord  then
                    Add( newgens, g );
                    Add( compositionSeries,
                         ClosureNormalizingElementPermGroup( N, g ) );
                    g := g ^ fac;
                fi;
                cof := cof * fac;
            od;
        od;
        Add( compositionSeries, TrivialSubgroup(G) );
    G.polycyclicGenerators := newgens;

    else

    # Otherwise: First determine an elementary abelian series for  <G> in
    # order to get a bssgs.
    ElementaryAbelianSeries( G );

    # Start with the trivial subgroup.
    N := TrivialSubgroup( G );
    compositionSeries := [ N ];

    # Loop over the elementary abelian series.
    for elabstep  in Reversed
        ( [ 1 .. Length( G.elementaryAbelianSeries ) - 1 ] )  do

        # For  each elementary abelian factor,  add  in the  intermediate
        # composition factors given by the pag system of <G>.
        for gen  in Reversed( Sublist
        ( G.elementaryAbelianSeries[ elabstep ].bssgs, [ 2 ..
          Length( G.elementaryAbelianSeries[ elabstep   ].bssgs ) -
              Length( G.elementaryAbelianSeries[ elabstep+1 ].bssgs ) ]
                ) )  do
        N := ClosureNormalizingElementPermGroup( N, gen );
        compositionSeries := Concatenation( [ N ],
                                 compositionSeries );
        od;

        # Finally add in the subgroup from the elementary abelian series.
        N := G.elementaryAbelianSeries[ elabstep ];
        compositionSeries := Concatenation( [ N ], compositionSeries );

    od;

    fi;

    # Return the result.
    return compositionSeries;

end;


#############################################################################
##
#F  SubnormalSeriesPPermGroup(<G>)  . . . . .  subnormal series for <p>-group
##
##  This function returns a subnormal series for <G>. A polycyclic generating
##  system for this series  that also contains strong generating sets for all
##  members of the subnormal series is bound to '<G>.polycyclicGenerators'.
##
SubnormalSeriesPPermGroup := function( G )
    local   subnormalSeries,  H,  gen;

    BaseStrongSubnormalGeneratingSetPPermGroup( G );
    H := TrivialSubgroup(G);
    InsertStabChain( H, Base(G), [  ] );
    subnormalSeries := [ H ];
    for gen in Reversed( Sublist( G.polycyclicGenerators,
        [ 2 .. Length(G.polycyclicGenerators) ] ) )  do
        H := ClosureNormalizingElementPermGroup( H, gen );
        subnormalSeries := Concatenation( [H], subnormalSeries );
    od;
    subnormalSeries := Concatenation( [G], subnormalSeries );
    return subnormalSeries;

end;


#############################################################################
##
#F  CentralCompositionSeriesPPermGroup( <G> ) . .  central composition series
##
##  This function calls 'SubnormalSeriesPPermGroup' to  calculate a subnormal
##  series with  strong generating  set and then improves it to be a  central
##  series maintaining the good properties of the generating set.  The method
##  used is known as Holt's algorithm.
##
CentralCompositionSeriesPPermGroup := function( G )
    local   words,  H,  gens,  base,  a,  f,  n,  i,  j,  h,  k,  l,  eps,
        old_gen_k, list,  pgens,  tmp,  s,  p,  compositionSeries;

    # <G> must be a <p>-group
    tmp := Set( FactorsInt( Size(G) ) );
    if 1 < Length(tmp)  then
        Error( "<G> must be a p-group" );
    fi;
    p := tmp[1];

    # calculate  a  base-adapted  subnormal  series   and   refine  it  to  a
    # composition series
    G.subnormalSeries := SubnormalSeriesPPermGroup(G);
    compositionSeries := CompositionSeriesSolvablePermGroup(G);
    base  := Base(G);
    pgens := G.polycyclicGenerators;
    n     := Length(pgens);

    # improve the composition series to a  p-central series  maintaining  the
    # base-adaption (Holt's algorithm)
    for i  in Reversed([1..n-2])  do
        for j  in Reversed([i+2..n])  do
            h := Comm( pgens[j], pgens[i] );
            while not h in compositionSeries[j+1]  do
                f := First([i+1..n], x -> not h in compositionSeries[x+1]);
        eps := 1;
        h := pgens[f] mod h;
        while not h in compositionSeries[ f+1 ]  do
            eps := eps + 1;
            h := pgens[f] mod h;
        od;
                if eps > 1  then
                    h := h^((1/eps) mod p);
                fi;
                k := f;
                for l  in [f+1..j]  do
            eps := 0;
                    while not h in compositionSeries[ l+1 ]  do
            eps := eps + 1;
            h := pgens[l] mod h;
            od;
                    if eps > 0  then
                        old_gen_k := pgens[k];
                        a := PositionProperty( base, b->b ^ pgens[l] <> b );
                        tmp := pgens[k] * pgens[l]^eps;
                        for s  in [ k+1 .. l ]  do
                            pgens[s-1] := pgens[s];
                        od;
                        pgens[l] := tmp;
                        if a <> false
                           and ForAll(Sublist(base,[1..a]),b->b^old_gen_k=b)
                        then
                            pgens[l-1] := old_gen_k;
                        fi;
                        k := l;
                    fi;
                od;
                if k <> j then
                    tmp := pgens[k];
                    for s  in [ k+1 .. j ]  do
                        pgens[s-1] := pgens[s];
                    od;
                    pgens[j] := tmp;
                fi;
                for l  in [ f+1 .. j ]  do
                    gens := Sublist( pgens, [l..n] );
                    compositionSeries[l] := Subgroup( Parent(G), gens );
                    InsertStabChain( compositionSeries[l], base, gens );
                od;
                h := Comm( pgens[j], pgens[i] );
            od;
        od;
    od;

    return compositionSeries;

end;


#############################################################################
##

#F  PermGroupOps.AgGroup(<G>) . . . . .  make an ag group out of a perm group
##
PermGroupOps.AgGroup := function( G )
    local   PC,  A,  series,  index;

    # Get an elementary abelian series and a bssgs of <G>. Find the positions
    # in the bssgs where the generators of each factor in the series start.
    series := ElementaryAbelianSeries( G );
    index  := List( series, x -> Length( G.bssgs ) - Length( x.bssgs )+1 );

    # Calculate a pc presentation for <G>.
    PC := PcPresentationPermGroup( G, series, G.bssgs, index, false );

    # Construct the ag group <A> and the bijection between <A> and <G>.
    A := AgGroupFpGroup(PC);
    A.bijection := GroupHomomorphismByImages(A,G,A.generators,G.bssgs);
    A.bijection.isMapping           := true;
    A.bijection.isGroupHomomorphism := true;
    A.bijection.isInjective         := true;
    A.bijection.isMonomorphism      := true;
    A.bijection.isSurjective        := true;
    A.bijection.isEpimorphism       := true;
    A.bijection.isBijection         := true;
    A.bijection.isIsomorphism       := true;
    return A;

end;

#############################################################################
##
#F  PermGroupOps.PgGroup(<G>) . . . . . . . .  pg group out of <p>-perm group
##
PermGroupOps.PgGroup := function( G )
    local   PC,  A,  series,  index;

    # Find a central series of <G>.
    series := CentralCompositionSeriesPPermGroup( G );
    index  := [ 1 .. Length( series )+1 ];

    # Calculate a pc presentation for <G>.
    PC := PcPresentationPermGroup( G, series, G.polycyclicGenerators,
                       index, true );

    # and construct the ag group <A> and the bijection between <A> and <G>
    A := AgGroupFpGroup(PC);
    A.bijection := GroupHomomorphismByImages(A,G,A.generators,
                         G.polycyclicGenerators);
    A.bijection.isMapping           := true;
    A.bijection.isGroupHomomorphism := true;
    A.bijection.isInjective         := true;
    A.bijection.isMonomorphism      := true;
    A.bijection.isSurjective        := true;
    A.bijection.isEpimorphism       := true;
    A.bijection.isBijection         := true;
    A.bijection.isIsomorphism       := true;
    return A;

end;


#############################################################################
##
#F  PermGroupOps.ElementaryAbelianSeries( <G> ) . . elementary abelian series
##
PermGroupOps.ElementaryAbelianSeries := function( G )
    if IsList(G)  then
        return GroupOps.ElementaryAbelianSeries(G);
    else
        PermGroupOps.TryElementaryAbelianSeries(G);
        return G.elementaryAbelianSeries;
    fi;
end;


#############################################################################
##
#F  PermGroupOps.IsSolvable( <G> )  . . . . . . . . . .  test for solvability
##
##  'PermGroupOps.IsSolvable' calls 'PermGroupOps.TryElementaryAbelianSeries'
##  and returns the flag '<G>.isSolvable' which is set by that function.
##
PermGroupOps.IsSolvable := function( G )
    PermGroupOps.TryElementaryAbelianSeries(G);
    return G.isSolvable;
end;


#############################################################################
##

#E  Emacs . . . . . . . . . . . . . . . . . . . . . . . local emacs variables
##
##  Local Variables:
##  mode:           outline
##  outline-regexp: "#F\\|#V\\|#E"
##  fill-column:    77
##  fill-prefix:    "##  "
##  eval:           (hide-body)
##  End:
##