Quellcodebibliothek Statistik Leitseite products/sources/formale Sprachen/GAP/lib/   (Algebra von RWTH Aachen Version 4.15.1©)  Datei vom 18.9.2025 mit Größe 109 kB image not shown  

Quelle  stbcbckt.gi   Sprache: unbekannt

 
#############################################################################
##
##  This file is part of GAP, a system for computational discrete algebra.
##  This file's authors include Heiko Theißen.
##
##  Copyright of GAP belongs to its developers, whose names are too numerous
##  to list here. Please refer to the COPYRIGHT file for details.
##
##  SPDX-License-Identifier: GPL-2.0-or-later
##
##  This file contains the basic   routines for permutation group   backtrack
##  algorithms that are based  on partitions. These  routines are used in the
##  calculation  of   set   stabilizers, normalizers,    centralizers     and
##  intersections.
##

if not IsBound( LARGE_TASK )  then  LARGE_TASK := false;   fi;

# set some global variables
BindGlobal("STBBCKT_STRING_CENTRALIZER",MakeImmutable("Centralizer"));
BindGlobal("STBBCKT_STRING_REGORB1",MakeImmutable("_RegularOrbit1"));
BindGlobal("STBBCKT_STRING_REGORB2",MakeImmutable("RegularOrbit2"));
BindGlobal("STBBCKT_STRING_REGORB3",MakeImmutable("RegularOrbit3"));
BindGlobal("STBBCKT_STRING_SPLITOFF",MakeImmutable("SplitOffBlock"));
BindGlobal("STBBCKT_STRING_INTERSECTION",MakeImmutable("Intersection"));
BindGlobal("STBBCKT_STRING_PROCESSFIX",MakeImmutable("ProcessFixpoint"));
BindGlobal("STBBCKT_STRING_MAKEBLOX",MakeImmutable("_MakeBlox"));
BindGlobal("STBBCKT_STRING_SUBORBITS0",MakeImmutable("Suborbits0"));
BindGlobal("STBBCKT_STRING_SUBORBITS1",MakeImmutable("Suborbits1"));
BindGlobal("STBBCKT_STRING_SUBORBITS2",MakeImmutable("Suborbits2"));
BindGlobal("STBBCKT_STRING_SUBORBITS3",MakeImmutable("Suborbits3"));
BindGlobal("STBBCKT_STRING_TWOCLOSURE",MakeImmutable("TwoClosure"));


#############################################################################
##
#V  Refinements . . . . . . . . . . . . . . .  record of refinement processes
##
BindGlobal( "Refinements", AtomicRecord() );


#############################################################################
##
#F  IsSlicedPerm( <perm> )  . . . . . . . . . . . . . . . sliced permutations
##
DeclareRepresentation( "IsSlicedPerm", IsPerm and IsComponentObjectRep,
                        [ "length", "word", "lftObj","opr" ] );

#############################################################################
##
#F  UnslicedPerm@( <perm> ) . . . . . . . . . . . . . . . . . . . . . . local
##
InstallGlobalFunction( UnslicedPerm@, function( perm )
    local   prm,  i;

    if IsSlicedPerm( perm )  then
        prm := ();
        for i  in [ 1 .. perm!.length ]  do
            prm := LeftQuotient( perm!.word[ i ], prm );
        od;
        return prm;
    else
        return perm;
    fi;
end );

InstallMethod( \^, "sliced perm",true, [ IsPerm, IsSlicedPerm ], 0,
    function( p, perm )  return p ^ UnslicedPerm@( perm );  end );
InstallMethod( \^, "sliced perm",true, [ IsInt, IsSlicedPerm ], 0,
    function( p, perm )
    local   i;

    for i  in Reversed( [ 1 .. perm!.length ] )  do
        p := p / perm!.word[ i ];
    od;
    return p;
end );

InstallOtherMethod( \/,"sliced perm", true, [ IsObject, IsSlicedPerm ], 0,
    function( p, perm )
    local   i;

    for i  in [ 1 .. perm!.length ]  do
        p := p ^ perm!.word[ i ];
    od;
    return p;
end );

InstallMethod( PrintObj,"sliced perm", true, [ IsSlicedPerm ], 0,
    function( perm )
    Print( "<perm word of length ", perm!.length, ">" );
end );

InstallMethod( ViewObj,"sliced perm", true, [ IsSlicedPerm ], 0,
    function( perm )
    Print( "<perm word of length ", perm!.length, ">" );
end );

DeclareRepresentation( "IsSlicedPermInv", IsPerm and IsComponentObjectRep,
                           [ "length", "word", "lftObj", "opr" ] );

InstallOtherMethod( \^,"sliced perm", true, [ IsObject, IsSlicedPermInv ], 0,
    function( p, perm )
    local   i;

    for i  in [ 1 .. perm!.length ]  do
        p := p ^ perm!.word[ i ];
    od;
    return p;
end );

InstallMethod( PrintObj,"sliced perm", true, [ IsSlicedPermInv ], 0,
    function( perm )
    Print( "<perm word of length ", perm!.length, ">" );
end );

InstallMethod( ViewObj,"sliced perm", true, [ IsSlicedPermInv ], 0,
    function( perm )
    Print( "<perm word of length ", perm!.length, ">" );
end );

#############################################################################
##
#F  PreImageWord( <p>, <word> ) . . . . . . preimage under sliced permutation
##
InstallGlobalFunction( PreImageWord, function( p, word )
    local   i;

    for i  in Reversed( [ 1 .. Length( word ) ] )  do
        p := p / word[ i ];
    od;
    return p;
end );

#############################################################################
##
#F  ExtendedT( <t>, <pnt>, <img>, <G> ) . .  prescribe one more image for <t>
##
InstallGlobalFunction( ExtendedT, function( t, pnt, img, simg, G )
    local   bpt,  len,  edg;

    # Map the image with the part <t> that is already known.
    if simg = 0  then  img := img / t;
                 else  img := simg;     fi;

    # If <G> fixes <pnt>, nothing more can  be changed, so test whether <pnt>
    # = <img>.
    bpt := BasePoint( G );
    if bpt <> pnt  then
        if pnt <> img  then
            return false;
        fi;

    elif not IsBound( G.translabels[ img ] )  then
        return false;
    elif IsSlicedPerm( t )  then
        len := t!.length;
        while img <> bpt  do
            len := len + 1;
            edg := G.transversal[ img ];
            img := img ^ edg;
#            t!.rgtObj := t!.opr( t!.rgtObj, edg );
            t!.word[ len ] := edg;
        od;
        t!.length := len;
    else
        t := LeftQuotient( InverseRepresentative( G, img ), t );
    fi;

    return t;
end );

#############################################################################
##
#F  MeetPartitionStrat( <rbase>,<image>,<S>,<strat> ) .  meet acc. to <strat>
##
InstallGlobalFunction( MeetPartitionStrat, function(rbase,image,S,g,strat )
local  P,  p;

  if Length( strat ) = 0  then
    return false;
  fi;

  P := image.partition;
  for p  in strat  do
    if p[1] =  0  and
      not ProcessFixpoint( image, p[2], FixpointCellNo( P, p[3] ) )
    or p[1] <> 0  and
      SplitCell( P, p[1], S, p[2], g, p[3] ) <> p[3]  then
      return false;
    fi;
  od;
  return true;
end );

#############################################################################
##
#F  StratMeetPartition( <rbase>, <P>, <S>, <g> )  . construct a meet strategy
##
##  Entries in <strat> have the following meaning:
##    [p,s,i] (p<>0) means that `0 < |P[p]\cap S[s]/g| = i < |P[p]|',
##            i.e., a new cell with <i> points was appended to <P>
##                  (and these <i> have been taken out of `P[p]'),
##    [0,a,p] means that fixpoint <a> was mapped to fixpoint in `P[p]',
##            i.e., `P[p]' has become a one-point cell.
##
InstallGlobalFunction( StratMeetPartition, function( arg )
    local   P,  S,  # first and second partition
            g,      # permutation such that <P> meet <S> / <g> is constructed
            rbase,  # R-base record, which records processing of fixpoints
            strat,  # meet strategy, the result
            p,  s,  # indices looping over the cells of <P> resp. <S>
            i,      # result of call to `SpliltCell'
            pnt,    # fixpoint to be processed
            cellsP, #\
            blist,  #/ see explanation below
            splits,
            lS,
            cell,  nrcells;

    if not IsPartition( arg[ 1 ] )  then  rbase := arg[ 1 ];  p := 2;
                                    else  rbase := false;     p := 1;  fi;
    P := arg[ p ];
    S := arg[ p + 1 ];
    if Length( arg ) = p + 2  then  g := arg[ p + 2 ];
                              else  g := ();            fi;
    strat := [  ];

    # <cellsP> is a   list whose <a>th entry is   <i> if `a^g  in P[p]'. Then
    # `Set(cellsP{S[s]})'  is  the set of    (numbers of) cells  of <P>  that
    # contain a point from `S[s]/g'. A cell splits iff it contains points for
    # two such values of <s>.
    if IsOne( g )  then
        cellsP := P.cellno;
    else
        cellsP := ListWithIdenticalEntries( Length( P.cellno ), 0 );
        for i  in [ 1 .. NumberCells( P ) ]  do
            cell := Cell( P, i );
            cellsP{ OnTuples( cell, g ) } := i + 0 * cell;
        od;
    fi;

    # If <S> is just a set, it is interpreted as partition ( <S>|<S>^compl ).
    if IsPartition( S )  then
        nrcells := NumberCells( S ) - 1;
        lS:=S;
    else
        nrcells := 1;
        blist := BlistList( [ 1 .. NumberCells( P ) ], cellsP{ S } );
        p := Position( blist, true );
        if p <> fail  then
            IntersectBlist( blist, BlistList( [ 1 .. NumberCells( P ) ],
                cellsP{ Difference( [ 1 .. Length( P.cellno ) ], S ) } ) );
            p := Position( blist, true );
        fi;
        lS:=S;
        S := false;
    fi;

    for s  in [ 1 .. nrcells ]  do
      # now split with cell number s of S.
      if S=false then
        p:=lS;
      else
        p:=Cell(S,s);
      fi;
      p:=cellsP{p}; # the affected P-cells
      p:=Collected(p);
      splits:=[];
      for i in p do
        # a cell will split iff it contains more points than are in the
        # s-cell
        if P.lengths[i[1]]>i[2] then
          Add(splits,i[1]);
        fi;
      od;

      # this code is new, the extensive construction of blists in the old
      # version was awfully slow in larger degrees. ahulpke 11-aug-00
      for p in splits do
            # Last argument `true' means that the cell will split.
            i := SplitCell( P, p, lS, s, g, true );
            if not IsOne( g )  then
                cell := Cell( P, NumberCells( P ) );
                cellsP{ OnTuples( cell, g ) } := NumberCells( P ) + 0 * cell;
            fi;

            if rbase <> false  then
                Add( strat, [ p, s, i ] );

                # If  we have one  or two  new fixpoints, put  them  into the
                # base.
                if i = 1  then
                    pnt := FixpointCellNo( P, NumberCells( P ) );
                    ProcessFixpoint( rbase, pnt );
                    Add( strat, [ 0, pnt, NumberCells( P ) ] );
                    if IsTrivialRBase( rbase )  then
                        return strat;
                    fi;
                fi;
                if P.lengths[ p ] = 1  then
                    pnt := FixpointCellNo( P, p );
                    ProcessFixpoint( rbase, pnt );
                    Add( strat, [ 0, pnt, p ] );
                    if IsTrivialRBase( rbase )  then
                        return strat;
                    fi;
                fi;

            fi;
#            p := Position( blist, true, p );
        od;
    od;
    return strat;
end );

# the following functions are for suborbits given by blists, by element
# lists, or as points (the latter are crucial to save memory)
InstallGlobalFunction(SuboLiBli,function(ran,b)
  if IsInt(b) then
    return [b];
  elif IsBlistRep(b) then
    return ListBlist(ran,b);
  fi;
  return b;
end);

InstallGlobalFunction(SuboSiBli,function(b)
  if IsInt(b) then
    return 1;
  elif IsBlistRep(b) then
    return SizeBlist(b);
  else
    return Length(b);
  fi;
end);

InstallGlobalFunction(SuboTruePos,function(ran,b)
  if IsInt(b) then
    return Position(ran,b);
  elif IsBlistRep(b) then
    return Position(b,true);
  elif HasIsSSortedList(b) and IsSSortedList(b) then
    return Position(ran,MinimumList(b));
  else
    return First([1..Length(ran)],i->ran[i] in b);
  fi;
end);

InstallGlobalFunction(SuboUniteBlist,function(ran,a,b)
  if IsInt(b) then
    a[Position(ran,b)]:=true;
  elif IsBlistRep(b) then
    UniteBlist(a,b);
  else
    #UniteBlist(a,BlistList(ran,b));
    UniteBlistList(ran,a,b);
  fi;
end);

# sb is a list of length 3: [points,subs,blists]. The function returns a
# cell as sorted list of points
InstallGlobalFunction(ConcatSubos,function(ran,sb)
local b,i;
  if Length(sb[3])>0 then
    # blists are used
    b:=ShallowCopy(sb[3][1]);
    for i in [2..Length(sb[3])] do
      UniteBlist(b,sb[3][i]);
    od;
    UniteBlistList(ran,b,sb[1]);
    for i in sb[2] do
      UniteBlistList(ran,b,i);
    od;
    return ListBlist(ran,b);
  elif Length(sb[2])>0 then
    # blists are not used but worth using
    b:=BlistList(ran,sb[1]);
    for i in sb[2] do
      UniteBlistList(ran,b,i);
    od;
    return ListBlist(ran,b);
  else
    b:=ShallowCopy(sb[1]);
    for i in sb[2] do
      UniteSet(b,i);
    od;
    return b;
  fi;
end);

#############################################################################
##
#F  Suborbits( <G>, <tofix>, <b>, <Omega> ) . . . . . . . . . . . . suborbits
##
##  Returns a record with the following components:
##
##     domain: the set <Omega>
##  stabChainTop: top level of stabilizer chain for  $G_tofix$ (pointwise stabilizer)  with
##             base point <a> (may be different from <b>)
##       conj: an element mapping <b> to <a>
##      which: a list  whose  <p>th entry   is the  number   of the  suborbit
##             containing <p>
##    lengths: a (not strictly) sorted list of suborbit lengths (subdegrees)
##  byLengths: a list whose <i>th entry is the set of numbers of suborbits of
##             the <i>th distinct length appearing in `lengths'
##  partition: the partition into unions of suborbits of equal length
##  The  next three entries  are lists  whose <k>  entry refers  to the <k>th
##  suborbit.
##     blists: the suborbits as boolean lists
##       reps: a transversal  in  <G> s.t.   $a.reps[k]$  lies in  the  <k>th
##             suborbit (reps[k] = `false' if this is impossible)
##  orbitalPartitions:
##             a list to store the `OrbitalPartition' for each suborbit in
##
InstallGlobalFunction( Suborbits, function( arg )
    local   H,  tofix,  b,  Omega,  suborbits,  len,  bylen,
            G,  GG,  a,  conj,  ran,  subs,  all,  k,  pnt,  orb,  gen,
            perm,  omega,  P,  cell,  part,  p,  i, la,bl,
            rep,rep2,te,stabgens;

    # Get the arguments.
    H := arg[ 1 ];
    tofix := arg[ 2 ];
    b     := arg[ 3 ];
    Omega := arg[ 4 ];
    IsRange(Omega);
    if b = 0  then  part := false;  b := Omega[ 1 ];
              else  part := true;   fi;

    G := StabChainMutable( H );
    bl:=Length(BaseStabChain(G));
    conj := One( H );

    # Replace  <H> by  the stabilizer of  all elements  of <tofix> except the
    # last.
    len := Length( tofix );
    for i  in [ 1 .. len ]  do
        conj := conj * InverseRepresentative( G, tofix[ i ] ^ conj );
        G := G.stabilizer;
    od;

    if len <> 0  then
      b := b ^ conj;
      suborbits:=[];
    else
      if not IsBound( H!.suborbits )  then
        H!.suborbits := [  ];
      fi;
      suborbits := H!.suborbits;
    fi;

    # Replace <b> by the minimal element <a> in its <G>-orbit.
    # rep 0 is an element that maps <b> to the orbits base point
    if not IsInBasicOrbit( G, b )  then
      GG := EmptyStabChain( [  ], One( H ), b );
      AddGeneratorsExtendSchreierTree( GG, G.generators );
    else
      GG := G;
    fi;
    a := Minimum( GG.orbit );

    rep:=InverseRepresentative(GG,b);
    rep2:=InverseRepresentative(GG,a)^-1;
    conj := conj * rep*rep2;

    # try whether a and b are in the same path
    #conj := conj * InverseRepresentative( GG, b ) /
    #               InverseRepresentative( GG, a );

    ran := Immutable([ 1 .. Maximum( Omega ) ]);
    IsSSortedList(ran);

    k:=1;
    while k<=Length(suborbits)
      and (suborbits[k][1]<>a or Omega<>suborbits[k][2]) do
      k:=k+1;
    od;
    if k<=Length(suborbits) and suborbits[k][1]=a and Omega=suborbits[k][2] then
      subs := suborbits[ k ][3];
      Info(InfoBckt,2,"Cached suborbits ",a);
    else
      Info(InfoBckt,2,"Enter suborbits ",Size(H),":",a);

        # Construct the suborbits rooted at <a>.
        # GG is a head of a stabilizer chain with base orbit containing
        # b with min elm a
        if not IsIdenticalObj(G,GG) then
          GG:=CopyStabChain( G );
          ChangeStabChain( GG, [ a ], false );
          te:=GG.transversal;
          stabgens:=GG.stabilizer.generators;
          Unbind(GG);
        else
          stabgens:=G.stabilizer.generators;
          # now conjugate with rep, so that we get things based at 'a'
          # rep2 maps the basepoint to a
          te:=ShallowCopy(G.transversal);
          te[G.orbit[1]]:=rep2; # just one mapper further
          te[a]:=G.identity;
          stabgens:=List(stabgens,i->i^rep2);
        fi;

        subs := rec( stabChainTop := rec(orbit:=[a],
                                         transversal:=te,
                                         identity:=G.identity),
                        domain := Omega,
                         which := ListWithIdenticalEntries( Length(ran), 0 ),
                          reps := [ G.identity ],
                          blists:=[],
                       lengths := [ 1 ],
             orbitalPartitions := [  ] );
        subs.blists[1]:=[a];
        subs.which[ a ] := 1;
        if IsRange(Omega) and 1 in Omega then
          all:=BlistList(ran,[]);
        else
          all := BlistList( ran, ran );
          SubtractBlist( all, BlistList( ran, Omega ) );
        fi;
        all[ a ] := true;
        la:=Length(all)-1;

        k := 1;
        pnt := Position( all, false );
        while pnt <> fail  do
          k := k + 1;
          orb := [ pnt ];
          all[ pnt ] := true;
          for p  in orb  do
            for gen  in stabgens  do
              i := p ^ gen;
              if not all[ i ]  then
                Add( orb, i );
                all[ i ] := true;
              fi;
            od;
          od;
          la:=la-Length(orb);
          subs.which{ orb } := k + 0 * orb;
          #if IsInBasicOrbit( G, pnt )  then
          if IsBound(te[pnt]) then
            subs.reps[ k ] := true;
            subs.lengths[ k ] := Length( orb );
          else
            # Suborbits outside the root's orbit get negative length.
            subs.reps[ k ] := false;
            subs.lengths[ k ] := -Length( orb );
          fi;
          #UniteBlist( all, sublique );
          if QuoInt(Length(ran),Length(orb))>100 then
            if Length(orb)=1 then
              subs.blists[ k ] := orb[1];
            else
              subs.blists[ k ] := Immutable(Set(orb));
            fi;
          else
            subs.blists[ k ] := BlistList(ran,orb);
          fi;
          if la=0 then
            pnt:=fail;
          else
            pnt := Position( all, false, pnt );
          fi;
        od;
        subs.sublilen:=Length(subs.blists);

        # store if not too many
        if Length(suborbits)>bl then
          for i in [1..Length(suborbits)-1] do
            suborbits[i]:=suborbits[i+1];
          od;
          suborbits[Length(suborbits)]:=[a,Omega,subs];
        else
          Add(suborbits,[a,Omega,subs]);
        fi;

    fi;

    if part  and  not IsBound( subs.partition )  then
        if not IsBound( subs.lengths )  then
Error("this should not happen 2719");
#            subs.lengths := [  ];
#            for k  in [ 1 .. subs.sublilen ]  do
#                if subs.reps[ k ] = false  then
#                    Add( subs.lengths, -SizeBlist( subs.blists[k] ) );
#                else
#                    Add( subs.lengths, SizeBlist( subs.blists[k] ) );
#                fi;
#            od;
        fi;
        perm := Sortex( subs.lengths ) ^ -1;

        # Determine the partition into unions of suborbits of equal length.
        subs.byLengths := [  ];
        P := [  ];  omega := Set( Omega );  cell := [  ];  bylen := [  ];
        for k  in [ 1 .. Length( subs.lengths ) ]  do
            Append( cell, SuboLiBli( ran, subs.blists[ k ^ perm ] ) );
            AddSet( bylen, k ^ perm );
            if    k = Length( subs.lengths )
               or subs.lengths[ k + 1 ] <> subs.lengths[ k ]  then
                Add( P, cell );  SubtractSet( omega, cell );  cell := [  ];
                Add( subs.byLengths, bylen );  bylen := [  ];
            fi;
        od;
        if Length( omega ) <> 0  then
            Add( P, omega );
        fi;
        subs.partition := Partition( P );
    fi;
    subs := ShallowCopy( subs );
    subs.conj := conj;
    return subs;
end );

#############################################################################
##
#F  OrbitalPartition( <subs>, <k> ) . . . . . . . . . . make a nice partition
##
##
## ahulpke, added aug-2-00: If there are only one or two cells, the function
## will return just one cell (the partitions split functions can treat this
## as a special case anyhow).
InstallGlobalFunction( OrbitalPartition, function( subs, k )
local  dom,  # operation domain for the group
        ran,  # range including <dom>, for blist construction
        d,    # number of suborbits, estimate for diameter
        len,  # current path length
        K,    # set of suborbits <k> to process
        Key,  # discriminating information for each suborbit
        key,  # discriminating information for suborbit number <k>
        old,  # farthest distance zone constructed so far
        new,  # new distance zone being constructed
        img,  # new endpoint of path with known predecessor
        o, i, # suborbit of predecessor resp. endpoint
        P,    # points ordered by <key> information, as partition
        typ,  # types of <key> information that occur
        sub,  # suborbit as list of integers
        csiz,
        pos;  # position of cell with given <key> in <P>

  if IsInt( k ) and IsBound( subs.orbitalPartitions[ k ] ) then
    Info(InfoBckt,2,"Orbital partition ",k," cached");
    P:=subs.orbitalPartitions[k];
  else
    ran := Immutable([ 1 .. Length( subs.which ) ]);
    IsSSortedList(ran);
    d   := subs.sublilen;
    if IsRecord( k )  then  K := k.several;
                      else  K := [ k ];      fi;
    Key := 0;
    for k  in K  do
      if IsList( k )  and  Length( k ) = 1  then
        k := k[ 1 ];
      fi;
      key := ListWithIdenticalEntries( d, 0 );

      # Initialize the flooding algorithm for the <k>th suborbit.
      if IsInt( k )  then
        if subs.reps[ k ] = false  then
          sub := 0;
          key[ k ] := -1;
          new := [  ];
        else
          sub := SuboLiBli( ran, subs.blists[ k ] );
          key[ k ] := 1;
          new := [ k ];
        fi;
      else
        #sub := ListBlist( ran, UnionBlist( subs.blists{ k } ) );
        if IsEmpty(k) then
          sub:=[];
        else
          sub:=subs.blists[k[1]];
          if IsInt(sub) then
            sub:=BlistList(ran,[sub]);
          elif not IsBool(sub[1]) then
            sub:=BlistList(ran,sub);
          else
            sub:=ShallowCopy(sub); # don't overwrite
          fi;
          for o in [2..Length(k)] do
            SuboUniteBlist(ran,sub,subs.blists[k[o]]);
          od;
          sub:=ListBlist(ran,sub);
        fi;

        key{ k } := 1 + 0 * k;
        new := Filtered( k, i -> subs.reps[ i ] <> false );
      fi;
      len := 1;

      # If no new points were found in the last round, stop.
      while Length( new ) <> 0  do
        len := len + 1;
        old := new;
        new := [  ];

        # Map the suborbit <sub> with each old representative.
        for o  in old  do
          if subs.reps[o]<>false then
            if subs.reps[ o ] = true  then
              subs.reps[ o ] := InverseRepresentative( subs.stabChainTop,
                  SuboTruePos(ran, subs.blists[ o ] ) ) ^ -1;
            fi;
            for img  in OnTuples( sub, subs.reps[ o ] )  do

              # Find the suborbit <i> of the image.
              i := subs.which[ img ];

              # If this suborbit is encountered for the first time, add
              # it to <new> and store its distance <len>.
              if key[ i ] = 0  then
                Add( new, i );
                key[ i ] := len;
              fi;

              # Store the arrow which starts at suborbit <o>.
              key[ o ] := key[ o ] + d *
                          Length( sub ) ^ ( key[ i ] mod d );
            od;
          else
            Info(InfoWarning,1,"suborbits variant triggered, check!");
          fi;
        od;
      od;

      if sub <> 0  then
        Key := Key * ( d + d * Length( sub ) ^ d ) + key;
      fi;
    od;

    # Partition  <dom> into unions   of  suborbits w.r.t. the  values  of
    # <Key>.
    if Key = 0  then
      P:=[];
      if IsInt( k )  then
        subs.orbitalPartitions[ k ] := P;
      fi;
      return P;
    else

#T1:=Runtime()-T1;
      typ := Set( Key );
      csiz:=ListWithIdenticalEntries(Length(typ),0);
      dom:=List(typ,i->[[],[],[]]);
      for i in [1..Length(Key)] do
        pos := Position( typ, Key[ i ] );
        csiz[pos]:=csiz[pos]+AbsInt(subs.lengths[i]);
        if IsInt(subs.blists[i]) then
          AddSet(dom[pos][1],subs.blists[i]);
        elif IsBlistRep(subs.blists[i]) then
          Add(dom[pos][3],subs.blists[i]);
        else
          Add(dom[pos][2],subs.blists[i]);
        fi;
      od;
      if Sum(csiz)=Length(subs.domain) and Length(typ)=1 then
        P:=[];
        if IsInt( k )  then
          subs.orbitalPartitions[ k ] := P;
        fi;
        return P;
      elif Sum(csiz)=Length(subs.domain) and Length(typ)=2 then
        # only two cells
        # we need to indicate the first cell, the trick to take the sorted
        # one does not work
        P:=ConcatSubos(ran,dom[1]);
        if IsInt( k )  then
          subs.orbitalPartitions[ k ] := P;
        fi;
        return P;
      fi;

      P:=[];
      for pos in [1..Length(typ)] do
        sub := ConcatSubos( ran, dom[pos] );
        Add(P,sub);
      od;
#fi;
#T1:=Runtime()-T1;


      if Sum(List(P,Length)) <> Length(subs.domain)  then
        # there are fixpoints missing
        Add( P, Difference(subs.domain,Union(P)));
      fi;

    fi;

    P := Partition( P );
    if IsInt( k )  then
      subs.orbitalPartitions[ k ] := P;
    fi;

  fi;
  return P;
end );

#############################################################################
##
#F  EmptyRBase( <G>, <Omega>, <P> ) . . . . . . . . . . . . initialize R-base
##
InstallGlobalFunction( EmptyRBase, function( G, Omega, P )
    local   rbase,  pnt;

    rbase := rec( domain := Omega,
                    base := [  ],
                   where := [  ],
                     rfm := [  ],
               partition := StructuralCopy( P ),
                     lev := [  ] );
    if IsList( G )  then
        if IsIdenticalObj( G[ 1 ], G[ 2 ] )  then
            rbase.level2 := true;
        else
            rbase.level2 := CopyStabChain( StabChainImmutable( G[ 2 ] ) );
            rbase.lev2   := [  ];
        fi;
        G := G[ 1 ];
    else
        rbase.level2 := false;
    fi;
#    if IsSymmetricGroupQuick( G )  then
#        Info( InfoBckt, 1, "Searching in symmetric group" );
#        rbase.fix   := [  ];
#        rbase.level := NrMovedPoints( G );
#    else
        rbase.chain := CopyStabChain( StabChainImmutable( G ) );
        rbase.level := rbase.chain;
#    fi;

    # Process all fixpoints in <P>.
    for pnt  in Fixcells( P )  do
        ProcessFixpoint( rbase, pnt );
    od;

    return rbase;
end );

#############################################################################
##
#F  IsTrivialRBase( <rbase> ) . . . . . . . . . . . . . .  is R-base trivial?
##
InstallGlobalFunction( IsTrivialRBase, function( rbase )
    return        IsInt( rbase.level )
              and rbase.level <= 1
           or     IsRecord( rbase.level )
              and Length( rbase.level.genlabels ) = 0;
end );

#############################################################################
##
#F  AddRefinement( <rbase>, <func>, <args> )  . . . . . register R-refinement
##
InstallGlobalFunction( AddRefinement, function( rbase, func, args )
    if    Length( args ) = 0
       or not IsList( Last(args) )
       or Length( Last(args) ) <> 0  then
        Add( Last(rbase.rfm), rec( func := func,
                                                    args := args ) );
        Info( InfoBckt, 1, "Refinement ", func, ": ",
                NumberCells( rbase.partition ), " cells" );
    fi;
end );

#############################################################################
##
#F  ProcessFixpoint( <rbase>|<image>, <pnt> [, <img> ] )  .  process fixpoint
##
##  `ProcessFixpoint( rbase, pnt )' puts in <pnt> as new base point and steps
##  down to the stabilizer, unless <pnt>  is redundant, in which case `false'
##  is returned.
##  `ProcessFixpoint( image, pnt, img )' prescribes <img> as image for <pnt>,
##  extends the permutation and steps down to  the stabilizer. Returns `true'
##  if this was successful and `false' otherwise.
##
InstallGlobalFunction( ProcessFixpoint, function( arg )
    local   rbase,  image,  pnt,  img,  simg,  t;

    if Length( arg ) = 2  then
        rbase := arg[ 1 ];
        pnt   := arg[ 2 ];
        if rbase.level2 <> false  and  rbase.level2 <> true  then
            ChangeStabChain( rbase.level2, [ pnt ] );
            if BasePoint( rbase.level2 ) = pnt  then
                rbase.level2 := rbase.level2.stabilizer;
            fi;
        fi;
        if IsInt( rbase.level )  then
            rbase.level := rbase.level - 1;
        else
            ChangeStabChain( rbase.level, [ pnt ] );
            if BasePoint( rbase.level ) = pnt  then
                rbase.level := rbase.level.stabilizer;
            else
                return false;
            fi;
        fi;
    else
        image := arg[ 1 ];
        pnt   := arg[ 2 ];
        img   := arg[ 3 ];
        if image.perm <> true  then
            if Length( arg ) = 4  then  simg := arg[ 4 ];
                                  else  simg := 0;         fi;
            t := ExtendedT( image.perm, pnt, img, simg, image.level );
            if t = false  then
                return false;
            elif BasePoint( image.level ) = pnt  then
                image.level := image.level.stabilizer;
            fi;
            image.perm := t;
        fi;
        if image.level2 <> false  then
            t := ExtendedT( image.perm2, pnt, img, 0, image.level2 );
            if t = false  then
                return false;
            elif BasePoint( image.level2 ) = pnt  then
                image.level2 := image.level2.stabilizer;
            fi;
            image.perm2 := t;
        fi;
    fi;
    return true;
end );

#############################################################################
##
#F  RegisterRBasePoint( <P>, <rbase>, <pnt> ) . . . . . register R-base point
##
InstallGlobalFunction( RegisterRBasePoint, function( P, rbase, pnt )
    local   O,  strat,  k,  lev;

    if rbase.level2 <> false  and  rbase.level2 <> true  then
        Add( rbase.lev2, rbase.level2 );
    fi;
    Add( rbase.lev, rbase.level );
    Add( rbase.base, pnt );
    k := IsolatePoint( P, pnt );
    Info( InfoBckt, 1, "Level ", Length( rbase.base ), ": ", pnt, ", ",
            P.lengths[ k ] + 1, " possible images" );
    if not ProcessFixpoint( rbase, pnt )  then
        Info(InfoWarning,2,"Warning: R-base point is already fixed" );
    fi;
    Add( rbase.where, k );
    Add( rbase.rfm, [  ] );
    if P.lengths[ k ] = 1  then
        pnt := FixpointCellNo( P, k );
        ProcessFixpoint( rbase, pnt );
        AddRefinement( rbase, STBBCKT_STRING_PROCESSFIX, [ pnt, k ] );
    fi;
    if rbase.level2 <> false  then
        if rbase.level2 = true  then  lev := rbase.level;
                                else  lev := rbase.level2;  fi;
        if not IsInt( lev )  then
            O := OrbitsPartition( lev, rbase.domain );
            strat := StratMeetPartition( rbase, P, O );
            AddRefinement( rbase, STBBCKT_STRING_INTERSECTION, [ O, strat ] );
        fi;
    fi;
end );

#############################################################################
##
#F  NextRBasePoint( <P>, <rbase> [, <order> ] ) . . .  find next R-base point
##
InstallGlobalFunction( NextRBasePoint, function( arg )
    local  rbase,    # R-base to be extended
           P,        # partition of <Omega> to be refined
           order,    # order in which to try the cells of <Omega>
           lens,     # sequence of cell lengths of <P>
           p,        # the next point chosen
           k,  l;    # loop variables

    # Get the arguments.
    P     := arg[ 1 ];
    rbase := arg[ 2 ];
    if Length( arg ) > 2  then  order := arg[ 3 ];
                          else  order := false;     fi;

    # When  this is called,   there is  a point  that   is neither  fixed  by
    # <rbase.level> nor in <P>.
    lens := P.lengths;
    p := fail;
    if order <> false  then
        if IsInt( rbase.level )  then
            p := PositionProperty( order, p ->
                         lens[ CellNoPoint(P,p ) ] <> 1 );
        else
            p := PositionProperty( order, p ->
                         lens[ CellNoPoint(P,p) ] <> 1
                     and not IsFixedStabilizer( rbase.level, p ) );
        fi;
    fi;

    if p <> fail  then
        p := order[ p ];
    else
        lens := ShallowCopy( lens );
        order := [ 1 .. NumberCells( P ) ];
        SortParallel( lens, order );
        k := PositionProperty( lens, x -> x <> 1 );
        l := fail;
        while l = fail  do
            if IsInt( rbase.level )  then
                l := 1;
            else
                l := PositionProperty
                     ( P.firsts[ order[ k ] ] - 1 + [ 1 .. lens[ k ] ],
                       i -> not IsFixedStabilizer( rbase.level,
                               P.points[ i ] ) );
            fi;
            k := k + 1;
        od;
        p := P.points[ P.firsts[ order[ k - 1 ] ] - 1 + l ];
    fi;

    RegisterRBasePoint( P, rbase, p );
end );

#############################################################################
##
#F  RRefine( <rbase>, <image>, <uscore> ) . . . . . . . . . apply refinements
##
InstallGlobalFunction( RRefine, function( rbase, image, uscore )
local  Rf,  t;

  if not uscore then
    for Rf  in rbase.rfm[ image.depth ]  do
      t := CallFuncList( Refinements.( Rf.func ), Concatenation
                    ( [ rbase, image ], Rf.args ) );
      if   t = false  then  return fail;
      elif t <> true  then  return t;     fi;
    od;
    return true;
  else
    for Rf  in rbase.rfm[ image.depth ]  do
      if Rf.func[ 1 ] = '_'  then
        t := CallFuncList( Refinements.( Rf.func ), Concatenation
                      ( [ rbase, image ], Rf.args ) );
        if   t = false  then  return fail;
        elif t <> true  then  return t;     fi;
      fi;
    od;
    return true;
  fi;

  #old code
  for Rf  in rbase.rfm[ image.depth ]  do
      if not uscore  or  Rf.func[ 1 ] = '_'  then
          t := CallFuncList( Refinements.( Rf.func ), Concatenation
                        ( [ rbase, image ], Rf.args ) );
          if   t = false  then  return fail;
          elif t <> true  then  return t;     fi;
      fi;
  od;
  return true;

end );

#############################################################################
##
#F  PBIsMinimal( <range>, <a>, <b>, <S> ) . . . . . . . . . . minimality test
##
InstallGlobalFunction( PBIsMinimal, function( range, a, b, S )
    local   orb,  old,  pnt,  l,  img;

    if IsInBasicOrbit( S, b )  then
        return ForAll( S.orbit, p -> a <= p );
    elif b < a                      then  return false;
    elif IsFixedStabilizer( S, b )  then  return true;   fi;

    orb := [ b ];
    old := BlistList( range, orb );
    for pnt  in orb  do
        for l  in S.genlabels  do
            img := pnt ^ S.labels[ l ];
            if not old[ img ]  then
                if img < a  then
                    return false;
                fi;
                old[ img ] := true;
                Add( orb, img );
            fi;
        od;
    od;
    return true;
end );

#############################################################################
##
#F  SubtractBlistOrbitStabChain( <blist>, <R>, <pnt> )  remove orbit as blist
##
InstallGlobalFunction( SubtractBlistOrbitStabChain, function( blist, R, pnt )
    local   orb,  gen,  img;

    orb := [ pnt ];
    blist[ pnt ] := false;
    for pnt  in orb  do
        for gen  in R.generators  do
            img := pnt ^ gen;
            if blist[ img ]  then
                blist[ img ] := false;
                Add( orb, img );
            fi;
        od;
    od;
end );

#############################################################################
##
#F  PartitionBacktrack( <G>, <Pr>, <repr>, <rbase>, <data>, <L>, <R> )  . . .
##
InstallGlobalFunction( PartitionBacktrack,
    function( G, Pr, repr, rbase, data, L, R )
    local  PBEnumerate,
           blen,         # length of R-base
           rep,          # representative or `false', the result
           branch,       # level where $Lstab\ne Rstab$ starts
           image,        # image information running through the tree
           oldcel,       # old value of <image.partition.cellno>
           orb,  org,    # intersected (mapped) basic orbits of <G>
           orB,          # backup of <orb>
           range,        # range for construction of <orb>
           fix,  fixP,   # fixpoints of partitions at root of search tree
           obj,  prm,    # temporary variables for constructed permutation
           nrback,       # backtrack counter
           bail,         # do we want to bail out quickly?
           i,  dd,  p;   # loop variables

#############################################################################
##
#F      PBEnumerate( ... )  . . . . . . . recursive enumeration of a subgroup
##
    PBEnumerate := function( d, wasTriv )
        local  undoto,   # number of cells of <P> wanted after undoing
               oldprm,   #\
               oldprm2,  #/  old values of <image>
               a,        # current R-base point
               m,        # initial number of candidates in <orb>
               max,      # maximal number of candidates still needed
               b,        # image of base point currently being considered
               t;        # group element constructed, to be handed upwards

        if image.perm = false  then
            return fail;
        fi;
        image.depth := d;

        # Store the original values of <image.*>.
        undoto := NumberCells( image.partition );
        if image.perm = true  then
            oldcel := image.partition;
        else
            oldcel := image.partition.cellno;
            if IsSlicedPerm( image.perm ) then  oldprm := image.perm!.length;
                                          else  oldprm := image.perm;
                                          fi;
        fi;
        if image.level2 <> false  then  oldprm2 := image.perm2;
                                  else  oldprm2 := false;        fi;

        # Recursion comes to an end  if all base  points have been prescribed
        # images.
        if d > Length( rbase.base )  then
            if IsTrivialRBase( rbase )  then
                blen := Length( rbase.base );

                # Do     not  add the   identity    element  in the  subgroup
                # construction.
                if wasTriv  then

                    # In the subgroup case, assign to  <L> and <R> stabilizer
                    # chains when the R-base is complete.
                    L := ListStabChain( CopyStabChain( StabChainOp( L,
                                 rec( base := rbase.base,
                                   reduced := false ) ) ) );
                    R := ShallowCopy( L );

                    if image.perm <> true  then
                        Info( InfoBckt, 1, "Stabilizer chain with depths ",
                                DepthSchreierTrees( rbase.chain ) );
                    fi;
                    Info( InfoBckt, 1, "Indices: ",
                          IndicesStabChain( L[ 1 ] ) );
                    return fail;

                else
                    if image.perm = true  then
                        prm := MappingPermListList
                               ( rbase.fix[ Length( rbase.base ) ],
                                 Fixcells( image.partition ) );
                    else
                        prm := image.perm;
                    fi;
                    if image.level2 <> false  then
                        prm := UnslicedPerm@( prm );
                        if SiftedPermutation( image.level2,
                                   prm / UnslicedPerm@( image.perm2 ) )
                           = image.level2.identity  then
                            return prm;
                        fi;
                    elif Pr( prm )  then
                        return UnslicedPerm@( prm );
                    fi;
                    return fail;
                fi;

            # Construct the   next refinement  level. This  also  initializes
            # <image.partition> for the case ``image = base point''.
            else
                if not repr  then
                    oldcel := StructuralCopy( oldcel );
                fi;
                rbase.nextLevel( rbase.partition, rbase );
                if image.perm = true  then
                    Add( rbase.fix, Fixcells( rbase.partition ) );
                fi;
                Add( org, ListWithIdenticalEntries( Length( range ), 0 ) );
                if repr  then

                    # In  the representative  case,  change  the   stabilizer
                    # chains of <L> and <R>.
                    ChangeStabChain( L[ d ], [ rbase.base[ d ] ], false );
                    L[ d + 1 ] := L[ d ].stabilizer;
                    ChangeStabChain( R[ d ], [ rbase.base[ d ] ], false );
                    R[ d + 1 ] := R[ d ].stabilizer;

                fi;
            fi;

        fi;
        a := rbase.base[ d ];
        Info(InfoBckt,3,Ordinal(d)," basepoint: ",a);

        # Intersect  the current cell of <P>  with  the mapped basic orbit of
        # <G> (and also with the one of <H> in the intersection case).
        if image.perm = true  then
            orb[ d ] := BlistList( range, Cell( oldcel, rbase.where[ d ] ) );
            if image.level2 <> false  then
                b := Position( orb[ d ], true );
                while b <> fail  do
                    if not IsInBasicOrbit( rbase.lev2[ d ], b / image.perm2 )
                       then
                        orb[ d ][ b ] := false;
                    fi;
                    b := Position( orb[ d ], true, b );
                od;
            fi;
        else
            orb[ d ] := BlistList( range, [  ] );
            for p  in rbase.lev[ d ].orbit  do
                b := p ^ image.perm;
                if oldcel[ b ] = rbase.where[ d ]
               and ( image.level2 = false
                  or IsInBasicOrbit( rbase.lev2[d], b/image.perm2 ) )  then
                    orb[ d ][ b ] := true;
                    org[ d ][ b ] := p;
                fi;
            od;
        fi;
        if d=1 and ForAll(GeneratorsOfGroup(G),x->a^x=a) then
          orb[d][a]:=true; # ensure a is a possible image (can happen if
                           # acting on permutations with more points)
        fi;

        orB[ d ] := StructuralCopy( orb[ d ] );

        # Loop  over the candidate images  for the  current base point. First
        # the special case ``image = base'' up to current level.
        if wasTriv  then
            image.bimg[ d ] := a;

            # Refinements that start with '_' must be executed even when base
            # = image since they modify `image.data' etc.
            RRefine( rbase, image, true );

            # Recursion.
            PBEnumerate( d + 1, true );
            image.depth := d;

            # Now we  can  remove  the  entire   <R>-orbit of <a>  from   the
            # candidate list.
            SubtractBlist( orb[ d ], BlistList( range, L[ d ].orbit ) );

        fi;

        # Only the early points of the orbit have to be considered.
        m := SizeBlist( orB[ d ] );
        if m < Length( L[ d ].orbit )  then
            return fail;
        fi;
        max := PositionNthTrueBlist( orB[ d ],
                       m - Length( L[ d ].orbit ) + 1 );
        if wasTriv  and  a > max  then
            m := m - 1;
            if m < Length( L[ d ].orbit )  then
                return fail;
            fi;
            max := PositionNthTrueBlist( orB[ d ],
                           m - Length( L[ d ].orbit ) + 1 );
        fi;

        # Now the other possible images.
        b := Position( orb[ d ], true );
        if b <> fail  and  b > max  then
            b := fail;
        fi;
        while b <> fail  do

            # Try to prune the node with prop 8(ii) of Leon's paper.
            if not repr  and  not wasTriv  and  IsBound( R[ d ].orbit )  then
                dd := branch;
                while dd < d  do
                    if IsInBasicOrbit( L[ dd ], a )  and  not PBIsMinimal
                       ( range, R[ dd ].orbit[ 1 ], b, R[ d ] )  then
                        Info( InfoBckt, 3, d, ": point ", b,
                                " pruned by minimality condition" );
                        dd := d + 1;
                    else
                        dd := dd + 1;
                    fi;
                od;
            else
                dd := d;
            fi;

            if dd = d  then

                # Undo the  changes made to  <image.partition>, <image.level>
                # and <image.perm>.
                for i  in [ undoto+1 .. NumberCells( image.partition ) ]  do
                    UndoRefinement( image.partition );
                od;
                if image.perm <> true  then
                    image.level := rbase.lev[ d ];
                    if IsSlicedPerm( image.perm )  then
                        image.perm!.length := oldprm;
                         # Here and below the code that refers to `rgtObj` was used to avoid multiplication
                         # of permutations. It has been commented out for a long time, but accidentally remained
                         # documented in `doc/ref/stbchain.xml` until its withdrawal in 2018.
#                        image.perm!.rgtObj := oldrgt;
                    else
                        image.perm := oldprm;
                    fi;
                fi;
                if image.level2 <> false  then
                    image.level2 := rbase.lev2[ d ];
                    image.perm2  := oldprm2;
                fi;

                # If <b> could not be prescribed as image for  <a>, or if the
                # refinement was impossible, give up for this image.
                image.bimg[ d ] := b;
                IsolatePoint( image.partition, b );

                if ProcessFixpoint( image, a, b, org[ d ][ b ] )  then
#Error(a," ",b," ",Cells(rbase.partition),Cells(image.partition));
                    t := RRefine( rbase, image, false );
                else
                    t := fail;
                fi;

                if t <> fail  then

                    # Subgroup case, base <> image   at current level:   <R>,
                    #   which until now is identical to  <L>, must be changed
                    #   without affecting <L>, so take a copy.
                    if wasTriv  and  IsIdenticalObj( L[ d ], R[ d ] )  then
                        R{ [ d .. Length( rbase.base ) ] } := List(
                        L{ [ d .. Length( rbase.base ) ] }, CopyStabChain );
                        branch := d;
                    fi;

                    if 2 * d <= blen  then
                        ChangeStabChain( R[ d ], [ b ], false );
                        R[ d + 1 ] := R[ d ].stabilizer;
                    else
                        if IsBound( R[ d ].stabilizer )  then
                            R[ d + 1 ] := StrongGeneratorsStabChain( R[ d ] );
                        else
                            R[ d + 1 ] := R[ d ].generators;
                        fi;
                        R[ d + 1 ] := rec( generators := Filtered
                            ( R[ d + 1 ], gen -> b ^ gen = b ) );
                    fi;

                else
                    Info( InfoBckt, 5, d, ": point ", b,
                            " pruned by partition condition" );
                fi;

                # Recursion.
                if t = true  then
                    t := PBEnumerate( d + 1, false );
                    nrback:=nrback+1;
                    if bail and nrback>500 then
                      return infinity; # bail out, this will bail out
                                       # recursively
                    fi;
                    image.depth := d;
                fi;

                # If   <t>   =   `fail', either   the   recursive   call  was
                #   unsuccessful,  or all new  elements   have been added  to
                #   levels  below  the current one   (this happens if  base =
                #   image up to current level).
                if t <> fail  then

                    # Representative case, element found: Return it.
                    # Subgroup case, base <> image  before current level:  We
                    #   need  only find  a representative  because we already
                    #   know the stabilizer of <L> at an earlier level.
                    if repr  or  not wasTriv  then
                        return t;

                    # Subgroup case, base  <> image at current level: Enlarge
                    #   <L>    with  <t>. Decrease <max>     according to the
                    #   enlarged <L>. Reset <R> to the enlarged <L>.
                    else
                        for dd  in [ 1 .. d ]  do
                            AddGeneratorsExtendSchreierTree( L[ dd ], [ t ] );
                        od;
                        Info( InfoBckt, 1, "Level ", d,
                                ": ", IndicesStabChain( L[ 1 ] ) );
                        if m < Length( L[ d ].orbit )  then
                            return fail;
                        fi;
                        max := PositionNthTrueBlist( orB[ d ],
                                       m - Length( L[ d ].orbit ) + 1 );
                        R{ [ d .. Length( rbase.base ) ] } := List(
                        L{ [ d .. Length( rbase.base ) ] }, CopyStabChain );
                    fi;

                fi;

                # Now  we can remove the   entire <R>-orbit  of <b> from  the
                # candidate list.
                if      IsBound( R[ d ].translabels )
                    and IsBound( R[ d ].translabels[ b ] )  then
                    SubtractBlist( orb[ d ],
                            BlistList( range, R[ d ].orbit ) );
                else
                    SubtractBlistOrbitStabChain( orb[ d ], R[ d ], b );
                fi;

            fi;

            b := Position( orb[ d ], true, b );
            if b <> fail  and  b > max  then
                b := fail;
            fi;
        od;

        return fail;
    end;

##
#F      main function . . . . . . . . . . . . . . . . . . . . . . . . . . . .
##

    nrback:=0; # count the number of times we jumped up
    bail:=repr and ValueOption("bailout")=true;

    # If necessary, convert <Pr> from a list to a function.
    if     IsList( Pr )
       and (    IsTrivial( G )
             #or IsSymmetricGroupQuick( G )
             ) then
        obj := rec( lftObj := Pr[ 1 ],
#                    rgtObj := Pr[ 2 ],
                       opr := Pr[ 3 ],
                      prop := Pr[ 4 ] );
        Pr := gen -> obj.prop
              ( rec( lftObj := obj.lftObj
#             ,
#                     rgtObj := obj.opr( obj.rgtObj, gen ^ -1 )
            ) );
    fi;

    # Trivial cases first.
    if IsTrivial( G )  then
        if   not repr        then  return G;
        elif Pr( One( G ) )  then  return One( G );
                             else  return fail;      fi;
    fi;

    # Construct the <image>.
    image := rec( data := data,
                  bimg := [  ],
                 depth := 1 );
    if repr  then  image.partition := data[ 1 ];
             else  image.partition := rbase.partition;  fi;
    if IsBool( rbase.level2 )  then
        image.level2 := false;
    else
        image.level2 := rbase.level2;
        image.perm2  := rbase.level2.identity;
    fi;

    # If  <Pr> is  function,   multiply  permutations. Otherwise, keep   them
    # factorized.
#    if IsSymmetricGroupQuick( G )  then
#        image.perm := true;
#    else
        if IsList( Pr )  then
            image.perm := Objectify
                ( NewType( PermutationsFamily, IsSlicedPerm ),
                  rec( length := 0, word := [  ] ) );
            image.perm!.lftObj := Pr[ 1 ];
#            image.perm!.rgtObj := Pr[ 2 ];
            image.perm!.opr    := Pr[ 3 ];
            Pr                 := Pr[ 4 ];
        else
            image.perm := One( G );
        fi;
        image.level := rbase.chain;
#    fi;

    if repr  then

        # In the representative case, map the  fixpoints of the partitions at
        # the root of the search tree.
        if rbase.partition.lengths <> image.partition.lengths  then
            image.perm := false;
        else
            fix  := Fixcells( rbase.partition );
            fixP := Fixcells( image.partition );
            for i  in [ 1 .. Length( fix ) ]  do
                ProcessFixpoint( image, fix[ i ], fixP[ i ] );
            od;
        fi;

        # In   the representative case,   assign  to <L>  and <R>  stabilizer
        # chains.
        L := ListStabChain( CopyStabChain( StabChainImmutable( L ) ) );
        R := ListStabChain( CopyStabChain( StabChainImmutable( R ) ) );

    fi;

    org := [  ];  orb := [  ];  orB := [  ];
    range := [ 1 .. Last(rbase.domain) ];
    blen := infinity;
    rep := PBEnumerate( 1, not repr );
    if not repr  then
        ReduceStabChain( L[ 1 ] );
        return GroupStabChain( G, L[ 1 ], true );
    else
        return rep;
    fi;
end );

#############################################################################
##
#F  Refinements.ProcessFixpoint( <pnt>, <cellnum> )  . . .  process a fixpoint
##
InstallGlobalFunction(Refinements_ProcessFixpoint,
function( rbase, image, pnt, cellnum )
    local   img;

    img := FixpointCellNo( image.partition, cellnum );
    return ProcessFixpoint( image, pnt, img );
end);
Refinements.(STBBCKT_STRING_PROCESSFIX) := Refinements_ProcessFixpoint;

#############################################################################
##
#F  Refinements.Intersection( <O>, <strat> )  . . . . . . . . . . second type
##
InstallGlobalFunction(Refinements_Intersection,
function( rbase, image, Q, strat )
    local   t;

    if image.level2 = false  then  t := image.perm;
                             else  t := image.perm2;  fi;
    if IsSlicedPerm( t )  then
        t := ShallowCopy( t );
        SET_TYPE_COMOBJ( t, NewType( PermutationsFamily, IsSlicedPermInv ) );
    else
        t := t ^ -1;
    fi;
    return MeetPartitionStrat( rbase, image, Q, t, strat );
end);
Refinements.(STBBCKT_STRING_INTERSECTION) := Refinements_Intersection;

#############################################################################
##
#F  Refinements.Centralizer(<no>,<g>,<pnt>,<strat>) . P meet Pz for one point
##
InstallGlobalFunction(Refinements_Centralizer,
function( rbase, image, cellnum, g, pnt, strat )
    local   P,  img;

    P := image.partition;
    img := FixpointCellNo( P, cellnum ) ^ image.data[ g + 1 ];
    return     IsolatePoint( P, img ) = strat
           and ProcessFixpoint( image, pnt, img );
end);
Refinements.(STBBCKT_STRING_CENTRALIZER) := Refinements_Centralizer;

#############################################################################
##
#F  Refinements._MakeBlox( <rbase>, <image>, <len> )  . . . . . . . make blox
##
InstallGlobalFunction(Refinements__MakeBlox,
function( rbase, image, len )
    local   F;

    F := image.data[ 2 ];
    image.data[ 4 ] := Partition( Blocks( F, rbase.domain,
                               image.bimg{ [ 1, len ] } ) );
    return Collected( rbase.blox.lengths ) =
           Collected( image.data[ 4 ].lengths );
end);
Refinements.(STBBCKT_STRING_MAKEBLOX) := Refinements__MakeBlox;

#############################################################################
##
#F  Refinements.SplitOffBlock( <k>, <strat> ) . . . . . . . . split off block
##
InstallGlobalFunction(Refinements_SplitOffBlock,
function( rbase, image, k, strat )
    local   B,  a,  orb;

    B   := image.data[ 4 ];
    a   := FixpointCellNo( image.partition, k );
    orb := Cell( B, CellNoPoint(B,a) );
    if Length( orb ) = Length( rbase.domain )  then
        return false;
    else
        return MeetPartitionStrat( rbase, image, orb, (),strat );
    fi;
end);
Refinements.(STBBCKT_STRING_SPLITOFF) := Refinements_SplitOffBlock;

#############################################################################
##
#F  Refinements._RegularOrbit1( <d>, <len> )  . . . . . . extend mapped orbit
##
##  Computes orbit and transversal `bF' for group <F>  = `data[6]' regular on
##  that orbit.
##
InstallGlobalFunction(Refinements__RegularOrbit1,
function( rbase, image, d, len )
    local   F,  trees;

    trees := image.data[ 5 ];
    if d = 1  then
        F := image.data[ 6 ];
        image.regorb := EmptyStabChain( [  ], One( F ), image.bimg[ d ] );
        AddGeneratorsExtendSchreierTree( image.regorb,
                GeneratorsOfGroup( F ) );
        if Length( image.regorb.orbit ) <> Length( rbase.regorb.orbit )  then
            return false;
        fi;
        trees[ d ] := EmptyStabChain( [  ], One( F ),
                              image.regorb.orbit[ 1 ] );
    else
        trees[ d ] := StructuralCopy( trees[ d - 1 ] );
        AddGeneratorsExtendSchreierTree( trees[ d ],
          [ QuickInverseRepresentative
            ( image.regorb, image.bimg[ d ] ) ^ -1 ] );
        if Length( trees[ d ].orbit ) <> len  then
            return false;
        fi;
    fi;
    return true;
end);
Refinements.(STBBCKT_STRING_REGORB1) := Refinements__RegularOrbit1;

#############################################################################
##
#F  Refinements.RegularOrbit2( <d>, <orb>, <strat> )  . . . meet mapped orbit
##
##  Compute images `bhg' of `bh' under  $g$ in `trees[<d>].orbit = bE$ ($h\in
##  E$).
##  Entries in <strat> have the following meaning:
##    [i,j] means  that the image `bhg\in  P[j]' of  `bh  = orb[<i>]'  can be
##          calculated from `bg'.
##   [-p,j] means that fixpoint <p> was mapped to fixpoint in `P[j]',
##          i.e., `P[j]' has become a one-point cell.
##
InstallGlobalFunction(Refinements_RegularOrbit2,
function( rbase, image, d, orbit, strat )
    local   P,  trees,  orb,  i;

    P     := image.partition;
    trees := image.data[ 5 ];
    orb   := trees[ d ].orbit;
    for i  in strat  do
        if (   i[ 1 ] < 0
           and not ProcessFixpoint( image, -i[1], FixpointCellNo(P,i[2]) ) )
        or (   i[ 1 ] > 0
           and (    IsolatePoint( P, orb[ i[ 1 ] ] ) <> i[ 2 ]
                 or not ProcessFixpoint( image, orbit[i[1]], orb[i[1]] ) ) )
           then  return false;
        fi;
    od;
    return true;
end);
Refinements.(STBBCKT_STRING_REGORB2) := Refinements_RegularOrbit2;

#############################################################################
##
#F  Refinements.RegularOrbit3( <f>, <strat> ) . . . . .  find images of orbit
##
##  Register images `yhg' of `yh' under $g$ in an arbitrary orbit `yE' ($h\in
##  E$). `yg\in P[f]' is a one-point cell.
##  Entries in <strat> have the following meaning:
##    [yh,i,j] means that  the image `yhg\in P[j]' of  `yh' can be calculated
##             from `yg' and `bhg\in P[i]' (a one-point cell).
##      [-p,j] means that fixpoint <p> was mapped to fixpoint in `P[j]',
##             i.e., `P[j]' has become a one-point cell.
##
InstallGlobalFunction(Refinements_RegularOrbit3,
function( rbase, image, f, strat )
    local   P,  yg,  bhg,  hg,  yhg,  i;

    P   := image.partition;
    yg  := FixpointCellNo( P, f );
    for i  in strat  do
        if i[ 1 ] < 0  then
            if not ProcessFixpoint( image, -i[1], FixpointCellNo(P,i[2]) )
               then
                return false;
            fi;
        else
            bhg := FixpointCellNo( P, i[ 2 ] );
            hg  := InverseRepresentativeWord( image.regorb, bhg );
            yhg := PreImageWord( yg, hg );
            if    IsolatePoint( P, yhg ) <> i[ 3 ]
               or not ProcessFixpoint( image, i[ 1 ], yhg )  then
                return false;
            fi;
        fi;
    od;
    return true;
end);
Refinements.(STBBCKT_STRING_REGORB3) := Refinements_RegularOrbit3;

#############################################################################
##
#F  Refinements.Suborbits0( <tra>, <f>, <lens>, <byLen>, <strat> ) subdegrees
##
##  Computes   suborbits of the stabilizer in   <F> =  `image.data[2]' of the
##  fixpoint in cell no. <f>.  (If <F> is multiply  transitive, replace it by
##  the stabilizer of the first <tra>-1 images of R-base points.)
##
##  Returns `true' if (1)~the  list  of suborbit lengths (subdegrees)  equals
##  <lens>, (2)~the list of subdegree  frequencies equals <byLen> and (3)~the
##  meet  with  the  partition  into unions   of   suborbits of equal  length
##  succeeds.
##
InstallGlobalFunction(Refinements_Suborbits0,
function( rbase, image, tra, f, lens, byLen, strat )
    local   F,  pnt,  subs;

    F    := image.data[ 2 ];
    pnt  := FixpointCellNo( image.partition, f );
    subs := Suborbits( F, image.bimg{ [ 1 .. tra - 1 ] }, pnt,
                    rbase.domain );
    if    subs.lengths <> lens
       or List( subs.byLengths, Length ) <> byLen  then
        return false;
    else
        return MeetPartitionStrat( rbase, image, subs.partition, subs.conj,
                       strat );
    fi;
end);
Refinements.(STBBCKT_STRING_SUBORBITS0):=Refinements_Suborbits0;

#############################################################################
##
#F  Refinements.Suborbits1( <rbase>, <image>, <tra>, <f>, <k>, <strat> )  . .
##
##  Meets  the  image partition with the  orbital  partition of the  union of
##  orbital graphs of suborbits of length `subs.byLengths[ <k> ]'. (<tra> and
##  <f> as in `Suborbits0'.)
##
InstallGlobalFunction(Refinements_Suborbits1,
function( rbase, image, tra, f, k, strat )
    local   F,  pnt,  subs,  Q;

    F    := image.data[ 2 ];
    pnt  := FixpointCellNo( image.partition, f );
    subs := Suborbits( F, image.bimg{ [ 1 .. tra - 1 ] }, pnt,
                    rbase.domain );
    Q := OrbitalPartition( subs, subs.byLengths[ k ] );
    return MeetPartitionStrat( rbase, image, Q, subs.conj, strat );
end);
Refinements.(STBBCKT_STRING_SUBORBITS1) := Refinements_Suborbits1;

#############################################################################
##
#F  Refinements.Suborbits2( <rbase>, <image>, <tra>, <f>, <start>, <coll> ) .
##
##  Computes  for each suborbit the  intersection sizes with cells <start> or
##  more in the image partition. Stores the  result in `data[3]' (needed only
##  on this level,  hence no  '_'). Returns  `true'  if the collected  result
##  equals <coll>.
##
InstallGlobalFunction(Refinements_Suborbits2,
function( rbase, image, tra, f, start, coll )
    local   F,  types,  pnt,  subs,  i, k;

    F    := image.data[ 2 ];
    pnt  := FixpointCellNo( image.partition, f );
    subs := Suborbits( F, image.bimg{ [ 1 .. tra - 1 ] }, pnt,
                    rbase.domain );
    if start = 1  then
        image.data[ 3 ] := List( subs.blists, o -> [ -SuboSiBli( o ) ] );
    fi;
    types := image.data[ 3 ];
    for i  in [ start .. NumberCells( image.partition ) ]  do
        for k  in Set( subs.which
          { OnTuples( Cell( image.partition, i ), subs.conj ) } )  do
            AddSet( types[ k ], i );
        od;
    od;
    return Collected( types ) = coll;
end);
Refinements.(STBBCKT_STRING_SUBORBITS2) := Refinements_Suborbits2;

#############################################################################
##
#F  Refinements.Suborbits3( <rbase>, <image>, <tra>, <f>, <typ>, <strat> )  .
##
##  Meets  the image  partition with  the orbital partition   of the union of
##  orbital  graphs of suborbits of type  <typ>. Returns `false' if there are
##  not <many> of them. (<tra> and <f> as in `Suborbits0'.)
##
InstallGlobalFunction(Refinements_Suborbits3,
function( rbase, image, tra, f, typ, many, strat )
    local   F,  types,  pnt,  subs,  k,  Q;

    F     := image.data[ 2 ];
    types := image.data[ 3 ];
    pnt   := FixpointCellNo( image.partition, f );
    subs  := Suborbits( F, image.bimg{ [ 1 .. tra - 1 ] }, pnt,
                     rbase.domain );
    k := Filtered( [ 1 .. subs.sublilen ], k -> types[ k ] = typ );
    if Length( k ) <> many  then
        return false;
    else
        Q := OrbitalPartition( subs, k );
        return MeetPartitionStrat( rbase, image, Q, subs.conj, strat );
    fi;
end);
Refinements.(STBBCKT_STRING_SUBORBITS3) := Refinements_Suborbits3;

#############################################################################
##
#F  Refinements.TwoClosure( <G>, <Q>, <d>, <strat> )  . . . . . . two-closure
##
InstallGlobalFunction(Refinements_TwoClosure,
--> --------------------

--> maximum size reached

--> --------------------

[ Dauer der Verarbeitung: 0.54 Sekunden  (vorverarbeitet)  ]