gap-libs 4r7p5-2 / usr / share / gap / lib / oprt.gi

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


#############################################################################
##
#F  ExternalSet( <arg> )  . . . . . . . . . . . . .  external set constructor
##
InstallMethod( ExternalSet, "G, D, gens, acts, act", true, OrbitsishReq, 0,
    function( G, D, gens, acts, act )
    return ExternalSetByFilterConstructor( IsExternalSet,
                   G, D, gens, acts, act );
end );


#############################################################################
##
#F  ExternalSetByFilterConstructor(<filter>,<G>,<D>,<gens>,<acts>,<act>)
##
InstallGlobalFunction( ExternalSetByFilterConstructor,
    function( filter, G, D, gens, acts, act )
    local   xset;

    xset := rec(  );
    if IsPcgs( gens )  then
        filter := filter and IsExternalSetByPcgs;
    fi;
    if not IsIdenticalObj( gens, acts )  then
        filter := filter and IsExternalSetByActorsRep;
        xset.generators    := gens;
        xset.operators     := acts;
        xset.funcOperation := act;
    else
        filter := filter and IsExternalSetDefaultRep;
    fi;

    # Catch the case that 'D' is an empty list.
    # (Note that an external set shall be a collection and not a list.)
    if IsList( D ) and IsEmpty( D ) then
      D:= EmptyRowVector( CyclotomicsFamily );
    fi;

    Objectify( NewType( FamilyObj( D ), filter ), xset );
    SetActingDomain  ( xset, G );
    SetHomeEnumerator( xset, D );
    if not IsExternalSetByActorsRep( xset )  then
        SetFunctionAction( xset, act );
    fi;
    return xset;
end );


#############################################################################
##
#F  ExternalSetByTypeConstructor(<type>,<G>,<D>,<gens>,<acts>,<act>)
##
# The following function expects the type as first argument,  to avoid a call
# of `NewType'. It is called by `ExternalSubsetOp' and `ExternalOrbitOp' when
# they are called with an external set (which has already stored this type).
#
InstallGlobalFunction( ExternalSetByTypeConstructor,
    function( type, G, D, gens, acts, act )
    local   xset;
    
    xset := Objectify( type, rec(  ) );
    if not IsIdenticalObj( gens, acts )  then
        xset!.generators    := gens;
        xset!.operators     := acts;
        xset!.funcOperation := act;
    fi;
    xset!.ActingDomain   := G;
    xset!.HomeEnumerator := D;
    if not IsExternalSetByActorsRep( xset )  then
        xset!.FunctionAction := act;
    fi;
    return xset;
end );

#############################################################################
##
#M  RestrictedExternalSet
##
InstallMethod(RestrictedExternalSet,"restrict the acting domain",
  true,[IsExternalSet,IsGroup],0,
function(xset,U)
local A,newx;
  A:=ActingDomain(xset);
  if IsSubset(U,A) then
    return xset; # no restriction happens
  fi;
  if IsBound(xset!.gens) then
    # we would have to decompose into generators
    TryNextMethod();
  fi;
  newx:=ExternalSet(U,HomeEnumerator(xset),FunctionAction(xset));
  return newx;
end);

#############################################################################
##
#M  Enumerator( <xset> )  . . . . . . . . . . . . . . . .  the underlying set
##
InstallMethod( Enumerator,"external set -> HomeEnumerator", true,
  [ IsExternalSet ], 0, HomeEnumerator );

#############################################################################
##
#M  FunctionAction( <p>, <g> ) . . . . . . . . . . . .  acting function
##
InstallMethod( FunctionAction,"ExternalSetByActorsRep", true,
  [ IsExternalSetByActorsRep ], 0,
    xset -> function( p, g )
      local pos,actor;
      pos:=Position(xset!.generators,g);
      if pos<>fail then
        actor:=xset!.operators[pos];
      else
	pos:=Position(xset!.generators,g^-1);
	if pos<>fail then
	  actor:=xset!.operators[pos]^-1;
	else
	  Error("need to factor -- not yet implemented");
	fi;
      fi;
      return xset!.funcOperation(p,actor);
#    local   D;
#        D := Enumerator( xset );
#        return D[ PositionCanonical( D, p ) ^
#                  ( g ^ ActionHomomorphismAttr( xset ) ) ];
    end );

#############################################################################
##
#M  PrintObj( <xset> )  . . . . . . . . . . . . . . . . print an external set
##
InstallMethod( PrintObj,"External Set", true, [ IsExternalSet ], 0,
    function( xset )
    Print(HomeEnumerator( xset ));
end );

#############################################################################
##
#M  ViewObj( <xset> )  . . . . . . . . . . . . . . . . print an external set
##
InstallMethod( ViewObj,"External Set", true, [ IsExternalSet ], 0,
function( xset )
local he,i;
  if not HasHomeEnumerator(xset) then
    TryNextMethod();
  fi;
  Print("<xset:");
  he:=HomeEnumerator(xset);
  if Length(he)<15 then
    View(he);
  else
    Print("[");
    for i in [1..15] do
      View(he[i]);
      Print(",");
    od;
    Print(" ...]");
  fi;
  Print(">");
end );

#############################################################################
##
#M  Representative( <xset> )  . . . . . . . . . . first element in enumerator
##
InstallMethod( Representative,"External Set", true, [ IsExternalSet ], 0,
    xset -> Enumerator( xset )[ 1 ] );

#############################################################################
##
#F  ExternalSubset( <arg> ) . . . . . . . . . . . . .  external set on subset
##
InstallMethod( ExternalSubsetOp, "G, D, start, gens, acts, act", true,
        [ IsGroup, IsList, IsList,
          IsList,
          IsList,
          IsFunction ], 0,
    function( G, D, start, gens, acts, act )
    local   xset;
    
    xset := ExternalSetByFilterConstructor( IsExternalSubset,
                    G, D, gens, acts, act );
    xset!.start := Immutable( start );
    return xset;
end );

InstallOtherMethod( ExternalSubsetOp,
        "G, xset, start, gens, acts, act", true,
        [ IsGroup, IsExternalSet, IsList,
          IsList,
          IsList,
          IsFunction ], 0,
    function( G, xset, start, gens, acts, act )
    local   type,  xsset;

    type := TypeObj( xset );

    # The type of an external set can store the type of its external subsets,
    # to avoid repeated calls of `NewType'.
    if not IsBound( type![XSET_XSSETTYPE] )  then
        xsset := ExternalSetByFilterConstructor( IsExternalSubset,
                         G, HomeEnumerator( xset ), gens, acts, act );
        type![XSET_XSSETTYPE] := TypeObj( xsset );
    else
        xsset := ExternalSetByTypeConstructor( type![XSET_XSSETTYPE],
                         G, HomeEnumerator( xset ), gens, acts, act );
    fi;
    
    xsset!.start := Immutable( start );
    return xsset;
end );

InstallOtherMethod( ExternalSubsetOp,
        "G, start, gens, acts, act", true,
        [ IsGroup, IsList,
          IsList,
          IsList,
          IsFunction ], 0,
    function( G, start, gens, acts, act )
    return ExternalSubsetOp( G,
                   Concatenation( Orbits( G, start, gens, acts, act ) ),
                   start, gens, acts, act );
end );


#############################################################################
##
#M  ViewObj( <xset> ) . . . . . . . . . . . . . . . . .  for external subsets
##
InstallMethod( ViewObj, "for external subset", true,
    [ IsExternalSubset ], 0,
    function( xset )
    Print( xset!.start, "^G");
end );


#############################################################################
##
#M  PrintObj( <xset> )  . . . . . . . . . . . . . . . .  for external subsets
##
InstallMethod( PrintObj, "for external subset", true,
    [ IsExternalSubset ], 0,
    function( xset )
    Print( xset!.start, "^G < ", HomeEnumerator( xset ) );
end );
#T It seems to be necessary to distinguish representations
#T for a correct treatment of `PrintObj'.


#############################################################################
##
#M  Enumerator( <xset> )  . . . . . . . . . . . . . . .  for external subsets
##
InstallMethod( Enumerator, "for external subset with home enumerator",
    [ IsExternalSubset and HasHomeEnumerator],
    function( xset )
    local   G,  henum,  gens,  acts,  act,  sublist,  pnt,  pos;
    
    henum := HomeEnumerator( xset );
    if IsPlistRep(henum) and not IsSSortedList(henum) then
      TryNextMethod(); # there is no reason to use the home enumerator
    fi;

    G := ActingDomain( xset );
    if IsExternalSetByActorsRep( xset )  then
        gens := xset!.generators;
        acts := xset!.operators;
        act  := xset!.funcOperation;
    else
        gens := GeneratorsOfGroup( G );
        acts := gens;
        act  := FunctionAction( xset );
    fi;
    sublist := BlistList( [ 1 .. Length( henum ) ], [  ] );
    for pnt  in xset!.start  do
        pos := PositionCanonical( henum, pnt );
        if not sublist[ pos ]  then
            OrbitByPosOp( G, henum, sublist, pos, pnt, gens, acts, act );
        fi;
    od;
    return EnumeratorOfSubset( henum, sublist );
end );

InstallMethod( Enumerator,"for external orbit: compute orbit", true,
  [ IsExternalOrbit ], 0,
function( xset )
  if HasHomeEnumerator(xset) and not IsPlistRep(HomeEnumerator(xset)) then
    TryNextMethod(); # can't do orbit because the home enumerator might
    # imply a different `PositionCanonical' (and thus equivalence of objects)
    # method.
  fi;
  return Orbit(xset,Representative(xset));
end);

InstallMethod( Random,"for external orbit: via acting domain", true,
  [ IsExternalOrbit ], 0,
function( xset )
  if HasHomeEnumerator(xset) and not IsPlistRep(HomeEnumerator(xset)) then
    TryNextMethod(); # can't do orbit because the home enumerator might
    # imply a different `PositionCanonical' (and thus equivalence of objects)
    # method.
  fi;
  return FunctionAction(xset)(Representative(xset),Random(ActingDomain(xset)));
end);

#############################################################################
##
#F  ExternalOrbit( <arg> )  . . . . . . . . . . . . . . external set on orbit
##
InstallMethod( ExternalOrbitOp, "G, D, pnt, gens, acts, act", true,
        OrbitishReq, 0,
    function( G, D, pnt, gens, acts, act )
    local   xorb;
    
    xorb := ExternalSetByFilterConstructor( IsExternalOrbit,
                    G, D, gens, acts, act );
    SetRepresentative( xorb, pnt );
    xorb!.start := Immutable( [ pnt ] );
    return xorb;
end );

InstallOtherMethod( ExternalOrbitOp,
        "G, xset, pnt, gens, acts, act", true,
        [ IsGroup, IsExternalSet, IsObject,
          IsList,
          IsList,
          IsFunction ], 0,
    function( G, xset, pnt, gens, acts, act )
    local   type,  xorb;

    type := TypeObj( xset );
    
    # The type of  an external set  can store the type  of external orbits of
    # its points, to avoid repeated calls of `NewType'.
    if not IsBound( type![XSET_XORBTYPE] )  then
        xorb := ExternalSetByFilterConstructor( IsExternalOrbit,
                        G, HomeEnumerator( xset ), gens, acts, act );
        type![XSET_XORBTYPE] := TypeObj( xorb );
    else
        xorb := ExternalSetByTypeConstructor( type![XSET_XORBTYPE],
                        G, HomeEnumerator( xset ), gens, acts, act );
    fi;
    
    SetRepresentative( xorb, pnt );
    xorb!.start := Immutable( [ pnt ] );
    return xorb;
end );

InstallOtherMethod( ExternalOrbitOp,
        "G, pnt, gens, acts, act", true,
        [ IsGroup, IsObject,
          IsList,
          IsList,
          IsFunction ], 0,
    function( G, pnt, gens, acts, act )
    return ExternalOrbitOp( G, OrbitOp( G, pnt, gens, acts, act ),
                   gens, acts, act );
end );


#############################################################################
##
#M  ViewObj( <xorb> ) . . . . . . . . . . . . . . . . . .  for external orbit
##
InstallMethod( ViewObj, "for external orbit", true,
    [ IsExternalOrbit ], 0,
    function( xorb )
    Print( Representative( xorb ), "^G");
end );


#############################################################################
##
#M  PrintObj( <xorb> )  . . . . . . . . . . . . . . . . .  for external orbit
##
InstallMethod( PrintObj, "for external orbit", true,
    [ IsExternalOrbit ], 0,
    function( xorb )
    Print( Representative( xorb ), "^G < ", HomeEnumerator( xorb ) );
end );
#T It seems to be necessary to distinguish representations
#T for a correct treatment of `PrintObj'.


#############################################################################
##
#M  AsList( <xorb> )  . . . . . . . . . . . . . . . . . .  by orbit algorithm
##
InstallMethod( AsList,"external orbit", true, [ IsExternalOrbit ], 0,
    xorb -> Orbit( xorb, Representative( xorb ) ) );

#############################################################################
##
#M  AsSSortedList( <xorb> )
##
InstallMethod( AsSSortedList,"external orbit", true, [ IsExternalOrbit ], 0,
    xorb -> Set(Orbit( xorb, Representative( xorb ) )) );

#############################################################################
##
#M  <xorb> = <yorb> . . . . . . . . . . . . . . . . . . by ``conjugacy'' test
##
InstallMethod( \=, "xorbs",IsIdenticalObj,
  [ IsExternalOrbit, IsExternalOrbit ], 0,
function( xorb, yorb )
  if not IsIdenticalObj(ActingDomain     (xorb),ActingDomain     (yorb))
  or not IsIdenticalObj(FunctionAction(xorb),FunctionAction(yorb))
      then
      TryNextMethod();
  fi;
  return RepresentativeAction( xorb, Representative( xorb ),
		  Representative( yorb ) ) <> fail;
end );

#############################################################################
##
#M  <xorb> < <yorb>
##
InstallMethod( \<, "xorbs, via AsSSortedList",IsIdenticalObj,
  [ IsExternalOrbit, IsExternalOrbit ], 0,
function( xorb, yorb )
  if not IsIdenticalObj(ActingDomain     (xorb),ActingDomain     (yorb))
  or not IsIdenticalObj(FunctionAction(xorb),FunctionAction(yorb))
      then
      TryNextMethod();
  fi;
  return AsSSortedList(xorb)<AsSSortedList(yorb);
end );

InstallMethod( \=, "xorbs which know their size", IsIdenticalObj,
  [ IsExternalOrbit and HasSize, IsExternalOrbit  and HasSize], 0,
function( xorb, yorb )
  if Size(xorb)<>Size(yorb) then
    return false;
  fi;
  if (Size(xorb)>10  and not HasAsList(yorb))
  or not IsIdenticalObj(ActingDomain     (xorb),ActingDomain     (yorb))
  or not IsIdenticalObj(FunctionAction(xorb),FunctionAction(yorb))
      then
      TryNextMethod();
  fi;
  return Representative( xorb ) in AsList(yorb);
end );

InstallMethod( \=, "xorbs with canonicalRepresentativeDeterminator",
  IsIdenticalObj,
  [IsExternalOrbit and CanEasilyDetermineCanonicalRepresentativeExternalSet,
   IsExternalOrbit and CanEasilyDetermineCanonicalRepresentativeExternalSet],
  0,
function( xorb, yorb )
  if not IsIdenticalObj(ActingDomain     (xorb),ActingDomain     (yorb))
  or not IsIdenticalObj(FunctionAction(xorb),FunctionAction(yorb))
      then
      TryNextMethod();
  fi;
  return CanonicalRepresentativeOfExternalSet( xorb ) =
	  CanonicalRepresentativeOfExternalSet( yorb );
end );

# as this is not necessarily compatible with the global ordering, this
# method is disabled.
# #############################################################################
# ##
# #M  <xorb> < <yorb> . . . . . . . . . . . . . . . . .  by ``canon. rep'' test
# ##
# InstallMethod( \<,"xorbs with canonicalRepresentativeDeterminator",
#   IsIdenticalObj,
#     [ IsExternalOrbit and HasCanonicalRepresentativeDeterminatorOfExternalSet,
#       IsExternalOrbit and HasCanonicalRepresentativeDeterminatorOfExternalSet ],
#         0,
#     function( xorb, yorb )
#     if not IsIdenticalObj(ActingDomain     (xorb),ActingDomain     (yorb))
#     or not IsIdenticalObj(FunctionAction(xorb),FunctionAction(yorb))
#        then
#         TryNextMethod();
#     fi;
#     return CanonicalRepresentativeOfExternalSet( xorb ) <
#            CanonicalRepresentativeOfExternalSet( yorb );
# end );

#############################################################################
##
#M  <pnt> in <xorb> . . . . . . . . . . . . . . . . . . by ``conjugacy'' test
##
InstallMethod( \in,"very small xorbs: test in AsList", IsElmsColls,
  [ IsObject, IsExternalOrbit and HasSize], 0,
function( pnt, xorb )
  if Size(xorb)>10 then
    TryNextMethod();
  fi;
  return pnt in AsList(xorb);
end );

InstallMethod( \in,"xorb: RepresentativeAction", IsElmsColls,
  [ IsObject, IsExternalOrbit ], 0,
function( pnt, xorb )
  return RepresentativeAction( xorb, Representative( xorb ),
                   pnt ) <> fail;
end );

# if we keep a list we will often test the representative
InstallMethod( \in,"xset: Test representative equal", IsElmsColls, [ IsObject,
      IsExternalSet and HasRepresentative ], 
      10, #override even tests in element lists
function( pnt, xset )
  if Representative( xset ) = pnt  then
    return true;
  else
    TryNextMethod();
  fi;
end );

InstallMethod( \in, "xorb: HasEnumerator",IsElmsColls,
  [ IsObject, IsExternalOrbit and HasEnumerator ], 0,
function( pnt, xorb )
local   enum;
    enum := Enumerator( xorb );
    if IsConstantTimeAccessList( enum )  then  return pnt in enum;
                                         else  TryNextMethod();     fi;
end );

InstallMethod(\in,"xorb HasAsList",IsElmsColls,
  [ IsObject,IsExternalOrbit and HasAsList],
  1, # AsList should override Enumerator
function( pnt, xorb )
local  l;
  l := AsList( xorb );
  if IsConstantTimeAccessList( l )  then  return pnt in l;
					else  TryNextMethod();     fi;
end );

InstallMethod(\in,"xorb HasAsSSortedList",IsElmsColls,
  [ IsObject,IsExternalOrbit and HasAsSSortedList],
  2, # AsSSorrtedList should override AsList
function( pnt, xorb )
local  l;
  l := AsSSortedList( xorb );
  if IsConstantTimeAccessList( l )  then  return pnt in l;
					else  TryNextMethod();     fi;
end );

# this method should have a higher priority than the previous to avoid
# searches in vain.
InstallMethod( \in, "by CanonicalRepresentativeDeterminator", 
  IsElmsColls, [ IsObject,
        IsExternalOrbit and
	HasCanonicalRepresentativeDeterminatorOfExternalSet ], 1,
function( pnt, xorb )
local func;
  func:=CanonicalRepresentativeDeterminatorOfExternalSet(xorb);
  return CanonicalRepresentativeOfExternalSet( xorb ) = 
    func(ActingDomain(xorb),pnt)[1];
end );

#############################################################################
##
#M  ActionHomomorphism( <xset> ) . . . . . . . . .  action homomorphism
##
InstallGlobalFunction( ActionHomomorphism, function( arg )
    local   attr,  xset,  p;
    
    if arg[ Length( arg ) ] = "surjective"  or
       arg[ Length( arg ) ] = "onto"  then
        attr := SurjectiveActionHomomorphismAttr;
        Unbind( arg[ Length( arg ) ] );
    else
        attr := ActionHomomorphismAttr;
    fi;
    if Length( arg ) = 1  then
        xset := arg[ 1 ];
    elif     Length( arg ) = 2
         and IsComponentObjectRep( arg[ 2 ] )
         and IsBound( arg[ 2 ]!.actionHomomorphism )
         and IsActionHomomorphism( arg[ 2 ]!.actionHomomorphism )
         and Source( arg[ 2 ]!.actionHomomorphism ) = arg[ 1 ]  then
        return arg[ 2 ]!.actionHomomorphism;  # GAP-3 compatability
    else
        if IsFunction( arg[ Length( arg ) ] )  then  p := 1;
                                               else  p := 0;  fi;
        if Length( arg ) mod 2 = p  then
            xset := CallFuncList( ExternalSet, arg );
        elif IsIdenticalObj( FamilyObj( arg[ 2 ] ),
                          FamilyObj( arg[ 3 ] ) )  then
            xset := CallFuncList( ExternalSubset, arg );
        else
            xset := CallFuncList( ExternalOrbit, arg );
        fi;
    fi;
    return attr( xset );
end );


#############################################################################
##
#M  ActionHomomorphismConstructor( <xset>, <surj> )
##
InstallGlobalFunction( ActionHomomorphismConstructor, function(arg)
local   xset,surj,G,  D,  act,  fam,  filter,  hom,  i;

    xset:=arg[1];surj:=arg[2];
    G := ActingDomain( xset );
    D := HomeEnumerator( xset );
    act := FunctionAction( xset );
    fam := GeneralMappingsFamily( ElementsFamily( FamilyObj( G ) ),
                                  PermutationsFamily );
    if IsExternalSubset( xset )  then
        filter := IsActionHomomorphismSubset;
    else
        filter := IsActionHomomorphism;
    fi;
    if IsPermGroup( G )  then
	filter := filter and IsPermGroupGeneralMapping;
    fi;

    hom := rec(  );
    if Length(arg)>2 then
      filter:=arg[3];
    elif IsExternalSetByActorsRep( xset )  then
        filter := filter and IsActionHomomorphismByActors;
    elif     IsMatrixGroup( G )
         and IsScalarList( D[ 1 ] ) then
      if  act in [ OnPoints, OnRight ]  then
	# we act linearly. This might be used to compute preimages by linear
	# algebra
	# note that we do not test whether the domain actually contains a
	# vector space base. This will be done the first time,
	# `LinearActionBasis' is called (i.e. in the preimages routine).
	filter := filter and IsLinearActionHomomorphism;
      elif act=OnLines then
	filter := filter and IsProjectiveActionHomomorphism;
      fi;

#        if     not IsExternalSubset( xset )
#           and IsDomainEnumerator( D )
#           and IsFreeLeftModule( UnderlyingCollection( D ) )
#           and IsFullRowModule( UnderlyingCollection( D ) )
#           and IsLeftActedOnByDivisionRing( UnderlyingCollection( D ) )  then
#            filter := filter and IsLinearActionHomomorphism;
#        else
#            if IsExternalSubset( xset )  then
#                if HasEnumerator( xset )  then  D := Enumerator( xset );
#                                          else  D := xset!.start;         fi;
#            fi;
#	    Error("hier");
#            if IsSubset( D, IdentityMat
#                       ( Length( D[ 1 ] ), One( D[ 1 ][ 1 ] ) ) )  then
#            fi;
#        fi;
    # test for constituent homomorphism
    elif not IsExternalSubset( xset )
         and IsPermGroup( G )
         and IsList( D ) and IsCyclotomicCollection( D )
         and act = OnPoints  then


        filter := IsConstituentHomomorphism;
        hom.conperm := MappingPermListList( D, [ 1 .. Length( D ) ] );

        # if MappingPermListList took a family/group as an
        # argument then we could patch it instead
        #if IsHomCosetToPerm(One(G)) then
        #    hom.conperm := HomCosetWithImage( Homomorphism(G.1), 
        #                   One(Source(G)), hom.conperm );
        #fi;


    # test for action on disjoint sets of numbers-> blocks homomorphism
    elif not IsExternalSubset( xset )
         and IsPermGroup( G )
         and IsList( D )
         and ForAll( D, IsList and IsSSortedList )
	 and Length(D)>0 and Length(D[1])>0 and IsInt(D[1][1])
         and act = OnSets
	 # disjointness test
	 and Length(Set(Flat(D)))=Sum(List(D,Length))
	 then
        filter := IsBlocksHomomorphism;
        hom.reps := [  ];
        for i  in [ 1 .. Length( D ) ]  do
            hom.reps{ D[ i ] } := i + 0 * D[ i ];
        od;

    # try to find under which circumstances we want to avoid computing
    # images by the action but always use the AsGHBI
    elif 
     # we can decompose into generators 
     (IsPermGroup( G )  or  IsPcGroup( G )) and
     # the action is not harmless
     not (act=OnPoints or act=OnSets or act=OnTuples)

     then
        filter := filter and
	  IsGroupGeneralMappingByAsGroupGeneralMappingByImages;
    # action of fp group
    elif IsSubgroupFpGroup(G) then
      filter:=filter and IsFromFpGroupHomomorphism;
    fi;
    if HasBaseOfGroup( xset )  then
        filter := filter and IsActionHomomorphismByBase;
    fi;
    if surj  then
        filter := filter and IsSurjective;
    fi;
    Objectify( NewType( fam, filter ), hom );
    SetUnderlyingExternalSet( hom, xset );
    return hom;
end );

InstallMethod( ActionHomomorphismAttr,"call OpHomConstructor", true,
  [ IsExternalSet ], 0,
    xset -> ActionHomomorphismConstructor( xset, false ) );

#############################################################################
##
#M  SurjectiveActionHomomorphism( <xset> ) .  surj. action homomorphism
##
InstallMethod( SurjectiveActionHomomorphismAttr,
  "call Ac.Hom.Constructor", true, [ IsExternalSet ], 0,
   xset -> ActionHomomorphismConstructor( xset, true ) );

VPActionHom:=function( hom )
local name;
  name:="homo";
  if HasIsInjective(hom) and IsInjective(hom) then
    name:="mono";
    if HasIsSurjective(hom) and IsSurjective(hom) then
      name:="iso";
    fi;
  elif HasIsSurjective(hom) and IsSurjective(hom) then
    name:="epi";
  fi;
  Print( "<action ",name,"morphism>" );
end;


#############################################################################
##
#F  MultiActionsHomomorphism(G,pnts,ops)
##
InstallGlobalFunction(MultiActionsHomomorphism,function(G,pnts,ops)
  local gens,homs,trans,n,d,gs,i,j,mgi,ran,hom,imgs,c;
  gens:=GeneratorsOfGroup(G);
  homs:=[];
  trans:=[];
  n:=1;

  if Length(pnts)=1 then
    return DoSparseActionHomomorphism(G,[pnts[1]],gens,gens,ops[1],false);
  fi;

  imgs:=List(gens,x->());
  c:=0;
  for i in [1..Length(pnts)] do
    if ForAny(homs,x->FunctionAction(UnderlyingExternalSet(x))=ops[i] and 
                   pnts[i] in HomeEnumerator(UnderlyingExternalSet(x))) then
      Info(InfoGroup,1,"point ",i," already covered");
    else
      hom:=DoSparseActionHomomorphism(G,[pnts[i]],gens,gens,ops[i],false);
      d:=NrMovedPoints(Range(hom));
      if d>0 then 
	c:=c+1;
	homs[c]:=hom;
	trans[c]:=MappingPermListList([1..d],[n..n+d-1]);
	mgi:=MappingGeneratorsImages(hom)[2];
	for j in [1..Length(gens)] do
	  imgs[j]:=imgs[j]*mgi[j]^trans[c];
	od;
	n:=n+d;
      fi;
    fi;
  od;
  ran:=Group(imgs,());
  hom:=GroupHomomorphismByFunction(G,ran,
	function(elm)
	local i,p,q;
	  p:=();
	  for i in [1..Length(homs)] do
	    q:=ImagesRepresentative(homs[i],elm);
	    if q = fail and ValueOption("actioncanfail")=true then
	      return fail;
	    fi;
	    p:=p*(q^trans[i]);
	  od;
	  return p;
	end);

  SetImagesSource(hom,ran);
  SetMappingGeneratorsImages(hom,[gens,imgs]);
  SetAsGroupGeneralMappingByImages( hom, GroupHomomorphismByImagesNC
            ( G, ran, gens, imgs ) );

  return hom;
end);



#############################################################################
##
#M  ViewObj( <hom> )  . . . . . . . . . . . .  view an action homomorphism
##
InstallMethod( ViewObj, "for action homomorphism", true,
    [ IsActionHomomorphism ], 0, VPActionHom);

#############################################################################
##
#M  PrintObj( <hom> ) . . . . . . . . . . . . print an action homomorphism
##
InstallMethod( PrintObj, "for action homomorphism", true,
    [ IsActionHomomorphism ], 0, VPActionHom);
#T It seems to be difficult to find out what I can use
#T for a correct treatment of `PrintObj'.


#############################################################################
##
#M  Source( <hom> ) . . . . . . . . . . . .  source of action homomorphism
##
InstallMethod( Source, "action homomorphism",true,
  [ IsActionHomomorphism ], 0,
        hom -> ActingDomain( UnderlyingExternalSet( hom ) ) );

#############################################################################
##
#M  Range( <hom> )  . . . . . . . . . . . . . range of action homomorphism
##
InstallMethod( Range,"ophom: S(domain)", true,
  [ IsActionHomomorphism ], 0, hom ->
    SymmetricGroup( Length( HomeEnumerator(UnderlyingExternalSet(hom)) ) ) );

InstallMethod( Range, "surjective action homomorphism",
  [ IsActionHomomorphism and IsSurjective ],
function(hom)
local gens, imgs, ran, i, a, xset,opt;
  gens:=GeneratorsOfGroup( Source( hom ) );
  if false and HasSize(Source(hom)) and Length(gens)>0 then
    imgs:=[ImageElmActionHomomorphism(hom,gens[1])];
    opt:=rec(limit:=Size(Source(hom)));
    if IsBound(hom!.basepos) then
      opt!.knownBase:=hom!.basepos;
    fi;
    ran:=Group(imgs[1]);
    i:=2;
    while i<=Length(gens) and Size(ran)<Size(Source(hom)) do
      a:=ImageElmActionHomomorphism( hom, gens[i]);
      Add(imgs,a);
      ran:=DoClosurePrmGp(ran,[a],opt);
      i:=i+1;
    od;
  else
    imgs:=List(gens,gen->ImageElmActionHomomorphism( hom, gen ) );
    if Length(imgs)=0 then
      ran:= GroupByGenerators( imgs,
		ImageElmActionHomomorphism( hom, One( Source( hom ) ) ) );
    else
      ran:= GroupByGenerators(imgs,One(imgs[1]));
    fi;
  fi;
  # remember a known base
  if HasBaseOfGroup(UnderlyingExternalSet(hom)) then
    xset:=UnderlyingExternalSet(hom);
    if not IsBound( xset!.basePermImage )  then
        xset!.basePermImage:=List(BaseOfGroup( xset ),
				  b->PositionCanonical(Enumerator(xset),b));
    fi;
    SetBaseOfGroup(ran,xset!.basePermImage);
  fi;
  SetMappingGeneratorsImages(hom,[gens{[1..Length(imgs)]},imgs]);
  if HasSize(Source(hom)) then
    StabChainOptions(ran).limit:=Size(Source(hom));
  fi;
  if HasIsInjective(hom) and HasSource(hom) and IsInjective(hom) then
    UseIsomorphismRelation( Source(hom), ran );
  fi;
  return ran;
end);

#############################################################################
##
#M  RestrictedMapping(<ophom>,<U>)
##
InstallMethod(RestrictedMapping,"action homomorphism",
  CollFamSourceEqFamElms,[IsActionHomomorphism,IsGroup],0,
function(hom,U)
local xset,rest;

  xset:=RestrictedExternalSet(UnderlyingExternalSet(hom),U);
  if ValueOption("surjective")=true or (HasIsSurjective(hom) and
    IsSurjective(hom)) then
    rest:=SurjectiveActionHomomorphismAttr( xset );
  else
    rest:=ActionHomomorphismAttr( xset );
  fi;

  if HasIsInjective(hom) and IsInjective(hom) then
    SetIsInjective(rest,true);
  fi;
  if HasIsTotal(hom) and IsTotal(hom) then
    SetIsTotal(rest,true);
  fi;

  return rest;
end);

#############################################################################
##
#F  Action( <arg> )
##
InstallGlobalFunction( Action, function( arg )
    local   hom,  O;
    
    if not IsString(arg[Length(arg)]) then
      Add(arg,"surjective"); # enforce surjective action homomorphism -- we
                             # anyhow compute the image
    fi;
    PushOptions(rec(onlyimage:=true)); # we don't want `ActionHom' to build
                                       # a stabilizer chain.
    hom := CallFuncList( ActionHomomorphism, arg );
    PopOptions();
    O := ImagesSource( hom );
    O!.actionHomomorphism := hom;
    return O;
end );

#############################################################################
##
#F  Orbit( <arg> )  . . . . . . . . . . . . . . . . . . . . . . . . . . orbit
##
InstallMethod( OrbitOp,
        "G, D, pnt, [ 1gen ], [ 1act ], act", true,
        OrbitishReq, 
	20, # we claim this method is very good
    function( G, D, pnt, gens, acts, act )
    if Length( acts ) <> 1  then  TryNextMethod();
                            else  return CycleOp( acts[ 1 ], D, pnt, act );
    fi;
end );

InstallOtherMethod( OrbitOp,
        "G, pnt, [ 1gen ], [ 1act ], act", true,
        [ IsGroup, IsObject,
          IsList,
          IsList,
          IsFunction ], 
	  20, # we claim this method is very good
    function( G, pnt, gens, acts, act )
    if Length( acts ) <> 1  then  TryNextMethod();
                            else  return CycleOp( acts[ 1 ], pnt, act );  fi;
end );

InstallMethod( OrbitOp, "with domain", true, OrbitishReq,0,
function( G, D, pnt, gens, acts, act )
local orb,d,gen,i,p;
  pnt:=Immutable(pnt);
  d:=NewDictionary(pnt,false,D);
  orb := [ pnt ];
  AddDictionary(d,pnt);
  for p in orb do
    for gen in acts do
      i:=act(p,gen);
      MakeImmutable(i);
      if not KnowsDictionary(d,i) then
	Add( orb, i );
	AddDictionary(d,i);
      fi;
    od;
  od;
  return Immutable(orb);
end );


InstallOtherMethod( OrbitOp, "standard orbit algorithm:list", true,
        [ IsGroup, IsObject,
          IsList,
          IsList,
          IsFunction ], 0,
function( G, pnt, gens, acts, act )
local orb,d,gen,i,p,D;
  # try to find a domain
  D:=DomainForAction(pnt,acts,act);
  pnt:=Immutable(pnt);
  d:=NewDictionary(pnt,false,D);
  orb := [ pnt ];
  AddDictionary(d,pnt);
  for p in orb do
    for gen in acts do
      i:=act(p,gen);
      MakeImmutable(i);
      if not KnowsDictionary(d,i) then
	Add( orb, i );
	AddDictionary(d,i);
      fi;
    od;
  od;
  return Immutable(orb);
end );

# all other orbit methods now become obsolete -- the dictionaries do the
# magic.

# InstallMethod( OrbitOp, "with quick position domain", true, [IsGroup,
#   IsList and IsQuickPositionList,IsObject,IsList,IsList,IsFunction],0,
# function( G, D, pnt, gens, acts, act )
#     return OrbitByPosOp( G, D, BlistList( [ 1 .. Length( D ) ], [  ] ),
#                    PositionCanonical( D, pnt ), pnt, gens, acts, act );
# end );

InstallGlobalFunction( OrbitByPosOp,
    function( G, D, blist, pos, pnt, gens, acts, act )
    local   orb,  p,  gen,  img,pofu;
    
    if IsInternalRep(D) then
      pofu:=Position; # avoids one redirection, epsilon faster
    else
      pofu:=PositionCanonical;
    fi;
    blist[ pos ] := true;
    orb := [ pnt ];
    for p  in orb  do
        for gen  in acts  do
            img := act( p, gen );
            pos := pofu( D, img );
            if not blist[ pos ]  then
	      blist[ pos ] := true;
	      #Add( orb, img );
	      Add( orb, D[pos] ); # this way we do not store the new element
	      # but the already existing old one in D. This saves memory.
            fi;
        od;
    od;
    return Immutable( orb );
end );

#############################################################################
##
#M  \^( <p>, <G> ) . . . . . . . orbit of a point under the action of a group
##
##  Returns the orbit of the point <A>p</A> under the action of the group
##  <A>G</A>, with respect to the action <C>OnPoints</C>.
##
InstallOtherMethod( \^, "orbit of a point under the action of a group",
                    ReturnTrue, [ IsObject, IsGroup ], 0,

  function ( p, G )
    return Orbit(G,p,OnPoints);
  end );

#############################################################################
##
#F  OrbitStabilizer( <arg> )  . . . . . . . . . . . . .  orbit and stabilizer
##
InstallMethod( OrbitStabilizerOp, "`OrbitStabilizerAlgorithm' with domain",
        true, OrbitishReq, 0,
function( G, D, pnt, gens, acts, act )
local   orbstab;
  orbstab:=OrbitStabilizerAlgorithm(G,D,false,gens,acts,
                                    rec(pnt:=pnt,act:=act));
  return Immutable( orbstab );
end );
    
InstallOtherMethod( OrbitStabilizerOp, 
        "`OrbitStabilizerAlgorithm' without domain",true,
        [ IsGroup, IsObject, IsList, IsList, IsFunction ], 0,
function( G, pnt, gens, acts, act )
local   orbstab;
  orbstab:=OrbitStabilizerAlgorithm(G,false,false,gens,acts,
                                    rec(pnt:=pnt,act:=act));
  return Immutable( orbstab );
end );

#############################################################################
##
#M OrbitStabilizerAlgorithm
##
InstallMethod( OrbitStabilizerAlgorithm,"use stabilizer size",true,
  [IsGroup and IsFinite and CanComputeSizeAnySubgroup,IsObject,IsObject,
   IsList,IsList,IsRecord],0,
function( G,D,blist,gens,acts, dopr )
local   orb,  stb,  rep,  p,  q,  img,  sch,  i,d,act,
        onlystab, # do we only care about stabilizer?
        getrep, # function to get representative
	actsinv,# inverses of acts
	stopat, # index at which increasal stopped
	notinc, # nr of steps in whiuch we did not increase
	stabsub,# stabilizer seed
	doml,	# maximal orbit length
	dict,	# dictionary
	blico,	# copy of initial blist (to find out the true domain)
	blif,	# flag on whether a blist is given
	crossind,	# index D (via blist) -> orbit position
	ind,	# stabilizer index
	indh,	# 1/2 stabilizer index
	increp,	# do we still want to increase the rep list?
	incstb;	# do we still want to increase the stabilizer?

  d:=Immutable(dopr.pnt);
  if IsBound(dopr.act) then
    act:=dopr.act;
  else
    act:=dopr.opr;
  fi;

  onlystab:=IsBound(dopr.onlystab) and dopr.onlystab=true;

  # try to find a domain
  if IsBool(D) then
    D:=DomainForAction(d,acts,act);
  fi;
  dict:=NewDictionary(d,true,D);

  if IsBound(dopr.stabsub) then
    stabsub:=AsSubgroup(Parent(G),dopr.stabsub);
  else
    stabsub:=TrivialSubgroup(G);
  fi;
  # NC is safe
  stabsub:=ClosureSubgroupNC(stabsub,gens{Filtered([1..Length(acts)],
	    i->act(d,acts[i])=d)});

  if IsBool(D) then
    doml:=Size(G);
  else
    if blist<>false then
      doml:=Size(D)-SizeBlist(blist);
      blico:=ShallowCopy(blist); # the original indices, needed to see what
                                 # a full orbit is
    else
      doml:=Size(D);
    fi;
  fi;

  incstb:=Index(G,stabsub)>1; # do we still include stabilizer elements. If
  # it is `false' the index `ind' must be equal to the orbit size.
  orb := [ d ];

  if incstb=false then
    # do we still need to tick off the orbit in `blist' to
    # please the caller? (see below as well)
    if blist<>false then
      q:=PositionCanonical(D,d);
      blist[q]:=true;
    fi;
    return rec( orbit := orb, stabilizer := G );
  fi;

#  # test for small domains whether the orbit has length 1
#  if doml<10 then
#    if doml=1 or ForAll(acts,i->act( d, i )=d) then
#
#      # do we still need to tick off the orbit in `blist' to
#      # please the caller? (see below as well)
#      if blist<>false then
#	q:=PositionCanonical(D,d);
#	blist[q]:=true;
#      fi;
#
#      return rec( orbit := orb, stabilizer := G );
#    fi;
#    
#  fi;

  AddDictionary(dict,d,1);

  stb := stabsub; # stabilizer seed
  ind:=Size(G);
  indh:=QuoInt(Size(G),2);
  if not IsEmpty( acts )  then

    # using only a factorized transversal can be expensive, in particular if
    # the action is more complicated. We therefore store a certain number of
    # representatives fixed.
    actsinv:=false;

    getrep:=function(pos)
    local a,r;
      a:=rep[pos];
      if not IsInt(a) then
	return a;
      else
	r:=fail;
	while pos>1 and IsInt(a) do
	  if r=fail then
	    r:=gens[a];
	  else
	    r:=gens[a]*r;
	  fi;
	  pos:=LookupDictionary(dict,act(orb[pos],actsinv[a]));
	  a:=rep[pos];
	od;
	if pos>1 then 
	  r:=a*r;
	fi;
	return r;
      fi;
    end;
    notinc:=0;
    increp:=true;

    rep := [ One( gens[ 1 ] ) ];
    p := 1;
    while p <= Length( orb )  do
      for i  in [ 1 .. Length( gens ) ]  do

	img := act( orb[ p ], acts[ i ] );
	MakeImmutable(img);
	q:=LookupDictionary(dict,img);

	if q = fail  then
	  Add( orb, img );
	  AddDictionary(dict,img,Length(orb));

	  if increp then
	    if actsinv=false then
	      Add( rep, rep[ p ] * gens[ i ] );
	    else
	      Add( rep, i );
	    fi;
	    if indh<Length(orb) then
	      # the stabilizer cannot grow any more
	      increp:=false;
	      incstb:=false;
	    fi;
	  fi;

	elif incstb then
	  #sch := rep[ p ] * gens[ i ] / rep[ q ];
	  sch := getrep( p ) * gens[ i ] / getrep( q );
	  if not sch in stb  then
	    notinc:=0;

	    # NC is safe
	    stb:=ClosureSubgroupNC(stb,sch);
	    ind:=Index(G,stb);
	    indh:=QuoInt(ind,2);
	    if indh<Length(orb) then
	      # the stabilizer cannot grow any more
	      increp:=false;
	      incstb:=false;
	    fi;
	  else
	    notinc:=notinc+1;
	    if notinc*50>indh and notinc>1000 then
	      # we have failed often enough -- assume for the moment we have
	      # the right stabilizer
	      #Error("stop stabilizer increase");
	      stopat:=p;
	      incstb:=false; # do not increase the stabilizer, but keep
	                     # representatives
	      actsinv:=List(acts,Inverse);
	    fi;
	  fi;
	fi;

	if increp=false then #we know the stabilizer
	  if onlystab then
	    return rec(stabilizer:=stb);
	  # must the orbit contain the whole domain => extend?
	elif ind=doml and (not IsBool(D)) and Length(orb)<doml then
	    if blist=false then
	      orb:=D;
	    else
	      orb:=D{Filtered([1..Length(blico)],i->blico[i]=false)};
	      # we need to tick off the rest
	      UniteBlist(blist,
		BlistList([1..Length(blist)],[1..Length(blist)]));
	    fi;
	    return rec( orbit := orb, stabilizer := stb );
	  elif  ind=Length(orb) then
	    # we have reached the full orbit. No further tests
	    # neccessary

	    # do we still need to tick off the orbit in `blist' to
	    # please the caller?
	    if blist<>false then
	      # we decided not to use blists for the orbit calculation
	      # but a blist was given in which we have to tick off the
	      # orbit
	      if IsPositionDictionary(dict) then
		UniteBlist(blist,dict!.blist);
	      else
		for img in orb do
		  blist[PositionCanonical(D,img)]:=true;
		od;
	      fi;
	    fi;

	    return rec( orbit := orb, stabilizer := stb );
	  fi;
	fi;
      od;
      p := p + 1;
    od;

    if Size(G)/Size(stb)>Length(orb) then
      #Error("did not really get stabilizer!");
      p:=stopat;
      while p<=Length(orb) do
	img := act( orb[ p ], acts[ i ] );
	MakeImmutable(img);
	q:=LookupDictionary(dict,img);
	if q<=p then
	  sch := getrep( p ) * gens[ i ] / getrep( q );
	  if not sch in stb then
	    stb:=ClosureSubgroupNC(stb,sch);
	    if Size(G)/Size(stb)=Length(orb) then
	      p:=Length(orb); #done
	    fi;
	  fi;
	fi;
	p:=p+1;
      od;
      #Error("after");
    fi;
  fi;

  if blist<>false then
    # we decided not to use blists for the orbit calculation
    # but a blist was given in which we have to tick off the
    # orbit
    if IsPositionDictionary(dict) then
      UniteBlist(blist,dict!.blist);
    else
      for img in orb do
	blist[PositionCanonical(D,img)]:=true;
      od;
    fi;
  fi;

  return rec( orbit := orb, stabilizer := stb );
end );

InstallMethod( OrbitStabilizerAlgorithm,"collect stabilizer generators",true,
  [IsGroup,IsObject,IsObject, IsList,IsList,IsRecord],0,
function( G,D,blist,gens, acts, dopr )
local   orb,  stb,  rep,  p,  q,  img,  sch,  i,d,act,
	stabsub,	# stabilizer seed
	dict,  		# dictionary
	crossind;	# index D (via blist) -> orbit position

  d:=Immutable(dopr.pnt);
  if IsBound(dopr.act) then
    act:=dopr.act;
  else
    act:=dopr.opr;
  fi;

  # try to find a domain
  if IsBool(D) then
    D:=DomainForAction(d,acts,act);
  fi;

  if IsBound(dopr.stabsub) then
    stabsub:=AsSubgroup(Parent(G),dopr.stabsub);
  else
    stabsub:=TrivialSubgroup(G);
  fi;

  dict:=NewDictionary(d,true,D);

  # `false' the index `ind' must be equal to the orbit size.
  orb := [ d ];
  AddDictionary(dict,d,1);

  stb := stabsub; # stabilizer seed
  if not IsEmpty( acts )  then
    rep := [ One( gens[ 1 ] ) ];
    p := 1;
    while p <= Length( orb )  do
      for i  in [ 1 .. Length( gens ) ]  do

	img := act( orb[ p ], acts[ i ] );
	MakeImmutable(img);

	q:=LookupDictionary(dict,img);

	if q = fail  then
	  Add( orb, img );
	  AddDictionary(dict,img,Length(orb));
	  Add( rep, rep[ p ] * gens[ i ] );
	else
	  sch := rep[ p ] * gens[ i ] / rep[ q ];
	  # NC is safe
	  stb:=ClosureSubgroupNC(stb,sch);
	fi;

      od;
      p := p + 1;
    od;

  fi;

  # can we compute the index from the orbit length?
  if HasSize(G) then
    if IsFinite(G) then
      SetSize(stb,Size(G)/Length(orb));
    else
      SetSize(stb,infinity);
    fi;
  fi;

  # do we care about a blist?
  if blist<>false then
    if IsPositionDictionary(dict) then
      # just copy over
      UniteBlist(blist,dict!.blist);
    else
      # tick off by hand
      for i in orb do
        blist[PositionCanonical(D,i)]:=true;
      od;
    fi;
  fi;

  return rec( orbit := orb, stabilizer := stb );
end );

#############################################################################
##
#F  Orbits( <arg> ) . . . . . . . . . . . . . . . . . . . . . . . . .  orbits
##

BindGlobal("OrbitsByPosOp",function( G, D, gens, acts, act )
    local   blist,  orbs,  next,  orb;
    
    blist := BlistList( [ 1 .. Length( D ) ], [  ] );
    orbs := [  ];
    for next in [1..Length(D)] do
      if blist[next]=false then
	# by calling `OrbitByPosOp' we avoid testing for positions twice.
	orb:=OrbitByPosOp(G,D,blist,next,D[next],gens,acts,act);
        Add( orbs, orb );
      fi;
    od;
    return Immutable( orbs );
end );

InstallMethod( OrbitsDomain, "for quick position domains", true,
  [ IsGroup, IsList and IsQuickPositionList, IsList, IsList, IsFunction ], 0,
  OrbitsByPosOp);

InstallMethod( OrbitsDomain, "for arbitrary domains", true,
    OrbitsishReq, 0,
function( G, D, gens, acts, act )
local   orbs, orb,sort,plist,pos,use,o;
  
  if Length(D)>0 and not IsMutable(D) and HasIsSSortedList(D) and IsSSortedList(D) 
    and CanEasilySortElements(D[1]) then
    return OrbitsByPosOp( G, D, gens, acts, act );
  fi;
  sort:=Length(D)>0 and CanEasilySortElements(D[1]);
  plist:=IsPlistRep(D);
  if not plist then
    use:=BlistList([1..Length(D)],[]);
  fi;
  orbs := [  ];
  pos:=1;
  while Length(D)>0  and pos<=Length(D) do
    orb := OrbitOp( G,D, D[pos], gens, acts, act );
    Add( orbs, orb );
    if plist then
      if sort then
	D:=Difference(D,orb);
	MakeImmutable(D); # to remember sortedness
      else
	D:=Filtered(D,i-> not i in orb);
      fi;
    else
      for o in orb do
        use[PositionCanonical(D,o)]:=true;
      od;
      # not plist -- do not take difference as there may be special
      # `PositionCanonical' method.
      while pos<=Length(D) and use[pos] do 
	pos:=pos+1;
      od;
    fi;
  od;
  return Immutable( orbs );
end );

InstallMethod( OrbitsDomain, "empty domain", true,
    [ IsGroup, IsList and IsEmpty, IsList, IsList, IsFunction ], 0,
function( G, D, gens, acts, act )
    return Immutable( [  ] );
end );

InstallOtherMethod(OrbitsDomain,"group without domain",true,[ IsGroup ], 0,
function( G )
  Error("You must give a domain on which the group acts");
end );

InstallMethod( Orbits, "for arbitrary domains", true, OrbitsishReq, 0,
function( G, D, gens, acts, act )
local   orbs, orb,sort,plist,pos,use,o,nc,ld,ld1;
    
  sort:=Length(D)>0 and CanEasilySortElements(D[1]);
  plist:=IsPlistRep(D);
  if not plist then
    use:=BlistList([1..Length(D)],[]);
  fi;
  nc:=true;
  ld1:=Length(D);
  orbs := [  ];
  pos:=1;
  while Length(D)>0  and pos<=Length(D) do
    orb := OrbitOp( G,D[pos], gens, acts, act );
    Add( orbs, orb );
    if plist then
      ld:=Length(D);
      if sort then
	D:=Difference(D,orb);
	MakeImmutable(D); # to remember sortedness
      else
	D:=Filtered(D,i-> not i in orb);
      fi;
      if Length(D)+Length(orb)>ld then
	nc:=false; # there are elements in `orb' not in D
      fi;
    else
      for o in orb do
        use[PositionCanonical(D,o)]:=true;
      od;
      # not plist -- do not take difference as there may be special
      # `PositionCanonical' method.
      while pos<=Length(D) and use[pos] do 
	pos:=pos+1;
      od;
    fi;
  od;
  if nc and ld1>10000 then
    Info(InfoPerformance,1,
    "You are calculating `Orbits' with a large set of seeds.\n",
      "#I  If you gave a domain and not seeds consider `OrbitsDomain' instead.");
  fi;
  return Immutable( orbs );
end );

InstallMethod( OrbitsDomain, "empty domain", true,
    [ IsGroup, IsList and IsEmpty, IsList, IsList, IsFunction ], 0,
function( G, D, gens, acts, act )
    return Immutable( [  ] );
end );

InstallOtherMethod( Orbits, "group without domain", true, [ IsGroup ], 0,
function( G )
  Error("You must give a domain on which the group acts");
end );

#############################################################################
##
#F  SparseActionHomomorphism( <arg> )   action homomorphism on `[1..n]'
##
InstallMethod( SparseActionHomomorphismOp,
        "domain given", true,
        [ IsGroup, IsList, IsList, IsList, IsList, IsFunction ], 0,
function( G, D, start, gens, acts, act )
local   list,  ps,  p,  i,  gen,  img,  pos,  imgs,  hom,orb,ran,xset;

  orb := List( start, p -> PositionCanonical( D, p ) );
  list := List( gens, gen -> [  ] );
  ps := 1;
  while ps <= Length( orb )  do
      p := D[ orb[ ps ] ];
      for i  in [ 1 .. Length( gens ) ]  do
	  gen := acts[ i ];
	  img := PositionCanonical( D, act( p, gen ) );
	  pos := Position( orb, img );
	  if pos = fail  then
	      Add( orb, img );
	      pos := Length( orb );
	  fi;
	  list[ i ][ ps ] := pos;
      od;
      ps := ps + 1;
  od;
  imgs := List( list, PermList );
  xset := ExternalSet( G, D{orb}, gens, acts, act);
  SetBaseOfGroup( xset, start );
  p:=RUN_IN_GGMBI; # no niceomorphism translation here
  RUN_IN_GGMBI:=true;
  hom := ActionHomomorphism(xset,"surjective" );
    ran:= Group( imgs, () );  # `imgs' has been created with `PermList'
  SetRange(hom,ran);
  SetImagesSource(hom,ran);
  SetAsGroupGeneralMappingByImages( hom, GroupHomomorphismByImagesNC
	  ( G, ran, gens, imgs ) );

  # We know that the points corresponding to `start' give a base. We can use
  # this to get images quickly, using a stabilizer chain in the permutation
  # group
  SetFilterObj( hom, IsActionHomomorphismByBase );
  RUN_IN_GGMBI:=p;
  return hom;
end );

#############################################################################
##
#F  DoSparseActionHomomorphism( <arg> )
##
InstallGlobalFunction(DoSparseActionHomomorphism,
function(G,start,gens,acts,act,sort)
local dict,p,i,img,imgs,hom,permimg,orb,imgn,ran,D,xset;

  # get a dictionary

  if IsMatrix(start) and Length(start)>0 and Length(start)=Length(start[1]) then
    # if we have matrices, we need to give a domain as well, to ensure the
    # right field
    D:=DomainForAction(start[1],acts,act);
  else # just base on the start values
    D:=fail;
  fi;
  dict:=NewDictionary(start[1],true,D);

  orb:=List(start,x->x); # do force list rep.
  for i in [1..Length(orb)] do
    AddDictionary(dict,orb[i],i);
  od;

  permimg:=List(acts,i->[]);

  # orbit algorithm with image keeper
  p:=1;
  while p<=Length(orb) do
    for i in [1..Length(gens)] do
      img := act(orb[p],acts[i]);
      imgn:=LookupDictionary(dict,img);
      if imgn=fail then
	Add(orb,img);
	AddDictionary(dict,img,Length(orb));
	permimg[i][p]:=Length(orb);
      else
	permimg[i][p]:=imgn;
      fi;
    od;
    p:=p+1;
  od;

  # any asymptotic argument is pointless here: In practice sorting is much
  # quicker than image computation.
  if sort then
    imgs:=Sortex(orb); # permutation we must apply to the points to be sorted.
    # was: permimg:=List(permimg,i->OnTuples(Permuted(i,imgs),imgs));
    # run in loop to save memory
    for i in [1..Length(permimg)] do
      permimg[i]:=Permuted(permimg[i],imgs);
      permimg[i]:=OnTuples(permimg[i],imgs);
    od;
  fi;
  
  for i in [1..Length(permimg)] do
    permimg[i]:=PermList(permimg[i]);
  od;

  # We know that the points corresponding to `start' give a base. We can use
  # this to get images quickly, using a stabilizer chain in the permutation
  # group
  if fail in permimg then
    Error("not permutations");
  fi;

  imgs:=permimg;
  ran:= Group( imgs, () );  # `imgs' has been created with `PermList'

  xset := ExternalSet( G, orb, gens, acts, act);
  if IsMatrix(start) and (act=OnPoints or act=OnRight or act=OnLines) then
    # act on vectors -- if we have a basis we have a base
    p:=RankMat(start);
    if p=Length(start[1]) then
      SetBaseOfGroup( xset, start );
    elif RankMat(orb{[1..Minimum(Length(orb),200)]})=Length(start[1]) then
      start:=ShallowCopy(start);
      i:=0;
      # we know we will be successful
      while p<Length(start[1]) do
	i:=i+1;
	if RankMat(Concatenation(start,[orb[i]]))>p then
	  Add(start,orb[i]);
	  p:=p+1;
	fi;
      od;
      SetBaseOfGroup( xset, start );
    fi;
  fi;

  p:=RUN_IN_GGMBI; # no niceomorphism translation here
  RUN_IN_GGMBI:=true;
  hom := ActionHomomorphism( xset,"surjective" );
  SetRange(hom,ran);
  SetImagesSource(hom,ran);
  SetMappingGeneratorsImages(hom,[gens,imgs]);
  SetAsGroupGeneralMappingByImages( hom, GroupHomomorphismByImagesNC
            ( G, ran, gens, imgs ) );

  if HasBaseOfGroup(xset) then
    SetFilterObj( hom, IsActionHomomorphismByBase );
  fi;
  RUN_IN_GGMBI:=p;

  return hom;
end);

#############################################################################
##
#M  SparseActionHomomorphism( <arg> )
##
InstallOtherMethod( SparseActionHomomorphismOp,
  "no domain given", true,
  [ IsGroup, IsList, IsList, IsList, IsFunction ], 0,
function( G, start, gens, acts, act )
  return DoSparseActionHomomorphism(G,start,gens,acts,act,false);
end); 

#############################################################################
##
#M  SortedSparseActionHomomorphism( <arg> )
##
InstallOtherMethod( SortedSparseActionHomomorphismOp,
  "no domain given", true,
  [ IsGroup, IsList, IsList, IsList, IsFunction ], 0,
function( G, start, gens, acts, act )
  return DoSparseActionHomomorphism(G,start,gens,acts,act,true);
end );

#############################################################################
##
#F  ExternalOrbits( <arg> ) . . . . . . . . . . . .  list of transitive xsets
##
InstallMethod( ExternalOrbits,
    "G, D, gens, acts, act",
    true,
    OrbitsishReq, 0,
    function( G, D, gens, acts, act )
    local   blist,  orbs,  next,  pnt,  orb;

    blist := BlistList( [ 1 .. Length( D ) ], [  ] );
    orbs := [  ];
    for next in [1..Length(D)] do
      if blist[next]=false then
        pnt := D[ next ];
        orb := ExternalOrbitOp( G, D, pnt, gens, acts, act );
        #SetCanonicalRepresentativeOfExternalSet( orb, pnt );
        SetEnumerator( orb, OrbitByPosOp( G, D, blist, next, pnt,
                gens, acts, act ) );
        Add( orbs, orb );
      fi;
    od;
    return Immutable( orbs );
end );

InstallOtherMethod( ExternalOrbits,
    "G, xset, gens, acts, act",
    true,
    [ IsGroup, IsExternalSet,
      IsList,
      IsList,
      IsFunction ], 0,
    function( G, xset, gens, acts, act )
    local   D,  blist,  orbs,  next,  pnt,  orb;

    D := Enumerator( xset );
    blist := BlistList( [ 1 .. Length( D ) ], [  ] );
    orbs := [  ];
    for next in [1..Length(D)] do
      if blist[next]=false then
        pnt := D[ next ];
        orb := ExternalOrbitOp( G, xset, pnt, gens, acts, act );
        #SetCanonicalRepresentativeOfExternalSet( orb, pnt );
        SetEnumerator( orb, OrbitByPosOp( G, D, blist, next, pnt,
                gens, acts, act ) );
        Add( orbs, orb );
      fi;
    od;
    return Immutable( orbs );
end );

#############################################################################
##
#F  ExternalOrbitsStabilizers( <arg> )  . . . . . .  list of transitive xsets
##
BindGlobal("ExtOrbStabDom",function( G, xsetD,D, gens, acts, act )
local   blist,  orbs,  next,  pnt,  orb,  orbstab,actrec;

    orbs := [  ];
    if IsEmpty( D ) then
      return Immutable( orbs );
    else
      blist:= BlistList( [ 1 .. Length( D ) ], [  ] );
    fi;
    for next in [1..Length(D)] do
      if blist[next]=false then
        pnt := D[ next ];
        orb := ExternalOrbitOp( G, xsetD, pnt, gens, acts, act );
        # was orbstab := OrbitStabilizer( G, D, pnt, gens, acts, act );
	actrec:=rec(pnt:=pnt, act:=act );
	# Does the external set give a kernel? Use it!
	if IsExternalSet(xsetD) and HasActionKernelExternalSet(xsetD) then
	  actrec.stabsub:=ActionKernelExternalSet(xsetD);
	fi;
	orbstab := OrbitStabilizerAlgorithm( G, D, blist, gens, acts, actrec);
        #SetCanonicalRepresentativeOfExternalSet( orb, pnt );
	if IsSSortedList(orbstab.orbit) then
	  SetAsSSortedList( orb, orbstab.orbit );
	else
	  SetAsList( orb, orbstab.orbit );
	fi;
        SetEnumerator( orb, orbstab.orbit );
        SetStabilizerOfExternalSet( orb, orbstab.stabilizer );
        Add( orbs, orb );
      fi;
    od;
    return Immutable( orbs );
end );

InstallMethod( ExternalOrbitsStabilizers,
    "arbitrary domain",
    true,
    OrbitsishReq, 0,
function( G, D, gens, acts, act )
  return ExtOrbStabDom(G,D,D,gens,acts,act);
end );

InstallOtherMethod( ExternalOrbitsStabilizers,
    "external set",
    true,
    [ IsGroup, IsExternalSet, IsList, IsList, IsFunction ], 0,
function( G, xset, gens, acts, act )
  return ExtOrbStabDom(G,xset,Enumerator(xset),gens,acts,act);
end );

#############################################################################
##
#F  Permutation( <arg> )  . . . . . . . . . . . . . . . . . . . . permutation
##
InstallGlobalFunction( Permutation, function( arg )
    local   g,  D,  gens,  acts,  act,  xset,  hom;

    # Get the arguments.
    g := arg[ 1 ];
    if Length( arg ) = 2  and  IsExternalSet( arg[ 2 ] )  then
        xset := arg[ 2 ];
        D := Enumerator( xset );
        if IsExternalSetByActorsRep( xset )  then
            gens := xset!.generators;
            acts := xset!.operators;
            act  := xset!.funcOperation;
        else
            act := FunctionAction( xset );
        fi;
    else
        D := arg[ 2 ];
        if IsDomain( D )  then
            D := Enumerator( D );
        fi;
        if IsFunction( arg[ Length( arg ) ] )  then
            act := arg[ Length( arg ) ];
        else
            act := OnPoints;
        fi;
        if Length( arg ) > 3  then
            gens := arg[ 3 ];
            acts := arg[ 4 ];
        fi;
    fi;
    
    if IsBound( gens )  and  not IsIdenticalObj( gens, acts )  then
        hom := ActionHomomorphismAttr( ExternalSetByFilterConstructor
                       ( IsExternalSet,
                         GroupByGenerators( gens ), D, gens, acts, act ) );
        return ImagesRepresentative( hom, g );
    else
        return PermutationOp( g, D, act );
    fi;
end );
                                
InstallMethod( PermutationOp, "object on list", true,
  [ IsObject, IsList, IsFunction ], 0,
    function( g, D, act )
    local   list,  blist,  fst,  old,  new,  pnt,perm;
    
    perm:=();
    if IsPlistRep(D) and Length(D)>2 
       and CanEasilySortElements(D[1]) then
      if not IsSSortedList(D) then
	D:=ShallowCopy(D);
        perm:=Sortex(D);
	D:=Immutable(D);
	SetIsSSortedList(D,true); # ought to be unnecessary, just be safe
      fi;
    fi;
    list := [  ];
    blist := BlistList( [ 1 .. Length( D ) ], [  ] );
    fst := Position( blist, false );
    while fst <> fail  do
        pnt := D[ fst ];
        new := fst;
        repeat
            old := new;
            pnt := act( pnt, g );
            new := PositionCanonical( D, pnt );
	    if new=fail then
	      Info(InfoWarning,2,"PermutationOp: mapping does not leave the domain invariant");
	      return fail;
	    elif blist[new] then
	      Info(InfoWarning,2,"PermutationOp: mapping is not injective");
	      return fail;
	    fi;
            blist[ new ] := true;
            list[ old ] := new;
        until new = fst;
        fst := Position( blist, false, fst );
    od;
    new:=PermList( list );
    if not IsOne(perm) then
      perm:=perm^-1;
      new:=new^perm;
    fi;
    return new;
end );

#############################################################################
##
#F  PermutationCycle( <arg> ) . . . . . . . . . . . . . . . cycle permutation
##
InstallGlobalFunction( PermutationCycle, function( arg )
    local   g,  D,  pnt,  gens,  acts,  act,  xset,  hom;

    # Get the arguments.
    g := arg[ 1 ];
    if Length( arg ) = 3  and  IsExternalSet( arg[ 2 ] )  then
        xset := arg[ 2 ];
        pnt  := arg[ 3 ];
        D := Enumerator( xset );
        if IsExternalSetByActorsRep( xset )  then
            gens := xset!.generators;
            acts := xset!.operators;
            act  := xset!.funcOperation;
        else
            act := FunctionAction( xset );
        fi;
    else
        D := arg[ 2 ];
        if IsDomain( D )  then
            D := Enumerator( D );
        fi;
        pnt := arg[ 3 ];
        if IsFunction( arg[ Length( arg ) ] )  then
            act := arg[ Length( arg ) ];
        else
            act := OnPoints;
        fi;
        if Length( arg ) > 4  then
            gens := arg[ 4 ];
            acts := arg[ 5 ];
        fi;
    fi;
    
    if IsBound( gens )  and  not IsIdenticalObj( gens, acts )  then
        hom := ActionHomomorphismAttr( ExternalSetByFilterConstructor
                       ( IsExternalSet,
                         GroupByGenerators( gens ), D, gens, acts, act ) );
        g := ImagesRepresentative( hom, g );
        return PermutationOp( g, CycleOp( g, PositionCanonical( D, pnt ),
                       OnPoints ), OnPoints );
    else
        return PermutationCycleOp( g, D, pnt, act );
    fi;
end );
                                
InstallMethod( PermutationCycleOp,"of object in list", true,
        [ IsObject, IsList, IsObject, IsFunction ], 0,
    function( g, D, pnt, act )
    local   list,  old,  new,  fst;
    
    list := [1..Size(D)];
    fst := PositionCanonical( D, pnt );
    if fst = fail  then
        return ();
    fi;
    new := fst;
    repeat
        old := new;
        pnt := act( pnt, g );
        new := PositionCanonical( D, pnt );
        if new = fail  then
            return fail;
        fi;
        list[ old ] := new;
    until new = fst;
    return PermList( list );
end );

#############################################################################
##
#F  Cycle( <arg> )  . . . . . . . . . . . . . . . . . . . . . . . . . . cycle
##
InstallGlobalFunction( Cycle, function( arg )
    local   g,  D,  pnt,  gens,  acts,  act,  xset,  hom,  p;
    
    # Get the arguments.
    g := arg[ 1 ];
    if Length( arg ) = 3  and  IsExternalSet( arg[ 2 ] )  then
        xset := arg[ 2 ];
        pnt  := arg[ 3 ];
        D := Enumerator( xset );
        if IsExternalSetByActorsRep( xset )  then
            gens := xset!.generators;
            acts := xset!.operators;
            act  := xset!.funcOperation;
        else
            act := FunctionAction( xset );
        fi;
    else
        if Length( arg ) > 2  and
           IsIdenticalObj( FamilyObj( arg[ 2 ] ),
                        CollectionsFamily( FamilyObj( arg[ 3 ] ) ) )  then
            D := arg[ 2 ];
            if IsDomain( D )  then
                D := Enumerator( D );
            fi;
            p := 3;
        else
            p := 2;
        fi;
        pnt := arg[ p ];
        if IsFunction( arg[ Length( arg ) ] )  then
            act := arg[ Length( arg ) ];
        else
            act := OnPoints;
        fi;
        if Length( arg ) > p + 1  then
            gens := arg[ p + 1 ];
            acts := arg[ p + 2 ];
        fi;
    fi;
    
    if IsBound( gens )  and  not IsIdenticalObj( gens, acts )  then
        hom := ActionHomomorphismAttr( ExternalOrbitOp
               ( GroupByGenerators( gens ), D, pnt, gens, acts, act ) );
        return D{ CycleOp( ImagesRepresentative( hom, g ),
                       PositionCanonical( D, pnt ), OnPoints ) };
    elif IsBound( D )  then
        return CycleOp( g, D, pnt, act );
    else
        return CycleOp( g, pnt, act );
    fi;
end );

InstallMethod( CycleOp,"of object in list", true,
        [ IsObject, IsList, IsObject, IsFunction ], 0,
    function( g, D, pnt, act )
    return CycleOp( g, pnt, act );
end );

CycleByPosOp := function( g, D, blist, fst, pnt, act )
    local   cyc,  new;
    
    cyc := [  ];
    new := fst;
    repeat
        Add( cyc, pnt );
        pnt := act( pnt, g );
        new := PositionCanonical( D, pnt );
        blist[ new ] := true;
    until new = fst;
    return Immutable( cyc );
end;

InstallOtherMethod( CycleOp, true, [ IsObject, IsObject, IsFunction ], 0,
    function( g, pnt, act )
    local   orb,  img;
    
    orb := [ pnt ];
    img := act( pnt, g );
    while img <> pnt  do
        Add( orb, img );
        img := act( img, g );
    od;
    return Immutable( orb );
end );

#############################################################################
##
#F  Cycles( <arg> ) . . . . . . . . . . . . . . . . . . . . . . . . .  cycles
##
InstallGlobalFunction( Cycles, function( arg )
    local   g,  D,  gens,  acts,  act,  xset,  hom;
    
    # Get the arguments.
    g := arg[ 1 ];
    if Length( arg ) = 2  and  IsExternalSet( arg[ 2 ] )  then
        xset := arg[ 2 ];
        D := Enumerator( xset );
        if IsExternalSetByActorsRep( xset )  then
            gens := xset!.generators;
            acts := xset!.operators;
            act  := xset!.funcOperation;
        else
            act := FunctionAction( xset );
        fi;
        D := Enumerator( xset );
    else
        D := arg[ 2 ];
        if IsDomain( D )  then
            D := Enumerator( D );
        fi;
        if IsFunction( arg[ Length( arg ) ] )  then
            act := arg[ Length( arg ) ];
        else
            act := OnPoints;
        fi;
        if Length( arg ) > 3  then
            gens := arg[ 3 ];
            acts := arg[ 4 ];
        fi;
    fi;
    
    if IsBound( gens )  and  not IsIdenticalObj( gens, acts )  then
        hom := ActionHomomorphismAttr( ExternalSetByFilterConstructor
                       ( IsExternalSet,
                         GroupByGenerators( gens ), D, gens, acts, act ) );
        return List( CyclesOp( ImagesRepresentative( hom, g ),
                       [ 1 .. Length( D ) ], OnPoints ), cyc -> D{ cyc } );
    else
        return CyclesOp( g, D, act );
    fi;
end );

InstallMethod( CyclesOp, true, [ IsObject, IsList, IsFunction ], 1,
    function( g, D, act )
    local   blist,  orbs,  next,  pnt,  pos,  orb;
    
    IsSSortedList(D);
    blist := BlistList( [ 1 .. Length( D ) ], [  ] );
    orbs := [  ];
    next := 1;
    while next <> fail do
        pnt := D[ next ];
        orb := CycleOp( g, D[ next ], act );
        Add( orbs, orb );
        for pnt  in orb  do
            pos := PositionCanonical( D, pnt );
            if pos <> fail  then
                blist[ pos ] := true;
            fi;
        od;
        next := Position( blist, false, next );
    od;
    return Immutable( orbs );
end );

#############################################################################
##
#F  Blocks( <arg> ) . . . . . . . . . . . . . . . . . . . . . . . . .  blocks
##
InstallOtherMethod( BlocksOp,
        "G, D, gens, acts, act", true,
        [ IsGroup, IsList,
          IsList,
          IsList,
          IsFunction ], 0,
    function( G, D, gens, acts, act )
    return BlocksOp( G, D, [  ], gens, acts, act );
end );

InstallMethod( BlocksOp,
        "via action homomorphism", true,
        [ IsGroup, IsList, IsList,
          IsList,
          IsList,
          IsFunction ], 0,
    function( G, D, seed, gens, acts, act )
    local   hom,  B;
    
    if Length(D)=1 then return Immutable([D]);fi;
    hom := ActionHomomorphism( G, D, gens, acts, act );
    B := Blocks( ImagesSource( hom ), [ 1 .. Length( D ) ] );
    B:=List( B, b -> D{ b } );
    # force sortedness
    if Length(B[1])>0 and CanEasilySortElements(B[1][1]) then
      B:=AsSSortedList(List(B,i->Immutable(Set(i))));
      IsSSortedList(B);
    fi;
    return B;
end );

InstallMethod( BlocksOp,
        "G, [  ], seed, gens, acts, act", true,
        [ IsGroup, IsList and IsEmpty, IsList,
          IsList,
          IsList,
          IsFunction ], 
	  20, # we claim this method is very good
    function( G, D, seed, gens, acts, act )
    return Immutable( [  ] );
end );

#############################################################################
##
#F  MaximalBlocks( <arg> )  . . . . . . . . . . . . . . . . .  maximal blocks
##
InstallOtherMethod( MaximalBlocksOp,
        "G, D, gens, acts, act", true,
        [ IsGroup, IsList,
          IsList,
          IsList,
          IsFunction ], 0,
    function( G, D, gens, acts, act )
    return MaximalBlocksOp( G, D, [  ], gens, acts, act );
end );

InstallMethod( MaximalBlocksOp,
        "G, D, seed, gens, acts, act", true,
        [ IsGroup, IsList, IsList,
          IsList,
          IsList,
          IsFunction ], 0,
    function ( G, D, seed, gens, acts, act )
    local   blks,       # blocks, result
            H,          # image of <G>
            blksH,      # blocks of <H>
	    onsetact;   # induces set action

    blks := BlocksOp( G, D, seed, gens, acts, act );

    # iterate until the action becomes primitive
    H := G;
    blksH := blks;
    onsetact:=function(l,g)
      return Set(List(l,i->act(i,g)));
    end;

    while Length( blksH ) <> 1  do
        H     := Action( H, blksH, onsetact );
        blksH := Blocks( H, [1..Length(blksH)] );
        if Length( blksH ) <> 1  then
            blks := List( blksH, bl -> Union( blks{ bl } ) );
        fi;
    od;

    # return the blocks <blks>
    return Immutable( blks );
end );

#############################################################################
##
#F  OrbitLength( <arg> )  . . . . . . . . . . . . . . . . . . .  orbit length
##
InstallMethod( OrbitLengthOp,"compute orbit", true, OrbitishReq, 0,
    function( G, D, pnt, gens, acts, act )
    return Length( OrbitOp( G, D, pnt, gens, acts, act ) );
end );

InstallOtherMethod( OrbitLengthOp,"compute orbit", true,
        [ IsGroup, IsObject,
          IsList,
          IsList,
          IsFunction ], 0,
    function( G, pnt, gens, acts, act )
    return Length( OrbitOp( G, pnt, gens, acts, act ) );
end );


#############################################################################
##
#F  OrbitLengths( <arg> )
##
InstallMethod( OrbitLengths,"compute orbits", true,
    OrbitsishReq, 0,
    function( G, D, gens, acts, act )
    return Immutable( List( Orbits( G, D, gens, acts, act ), Length ) );
end );

#############################################################################
##
#F  OrbitLengthsDomain( <arg> )
##
InstallMethod( OrbitLengthsDomain,"compute orbits", true,
    OrbitsishReq, 0,
    function( G, D, gens, acts, act )
    return Immutable( List( OrbitsDomain( G, D, gens, acts, act ), Length ) );
end );


#############################################################################
##
#F  CycleLength( <arg> )  . . . . . . . . . . . . . . . . . . .  cycle length
##
InstallGlobalFunction( CycleLength, function( arg )
    local   g,  D,  pnt,  gens,  acts,  act,  xset,  hom,  p;
    
    # test arguments
    if Length(arg)<2 or not IsMultiplicativeElementWithInverse(arg[1]) then
      Error("usage: CycleLength(<g>,<D>,<pnt>[,<act>])");
    fi;

    # Get the arguments.
    g := arg[ 1 ];
    if IsExternalSet( arg[ 2 ] )  then
        xset := arg[ 2 ];
        pnt := arg[ 3 ];
        if HasHomeEnumerator( xset )  then
            D := HomeEnumerator( xset );
        fi;
        act := FunctionAction( xset );
    else
      if Length( arg ) > 2  and
	  IsIdenticalObj( FamilyObj( arg[ 2 ] ),
		      CollectionsFamily( FamilyObj( arg[ 3 ] ) ) )  then
	  D := arg[ 2 ];
	  if IsDomain( D )  then
	      D := Enumerator( D );
	  fi;
	  p := 3;
      else
	  p := 2;
      fi;
      pnt := arg[ p ];
      if IsFunction( arg[ Length( arg ) ] )  then
	  act := arg[ Length( arg ) ];
      else
	  act := OnPoints;
      fi;
      if Length( arg ) > p + 1  then
	gens := arg[ p + 1 ];
	acts := arg[ p + 2 ];
	if not IsIdenticalObj( gens, acts )  then
	  xset:=ExternalOrbitOp(GroupByGenerators(gens),D,pnt,gens,acts,act);
	fi;
      fi;
    fi;
    
    if IsBound(xset) and IsExternalSetByActorsRep(xset) then
      # in all other cases the homomorphism was ignored anyhow
      hom := ActionHomomorphismAttr(xset);
      return CycleLengthOp( ImagesRepresentative( hom, g ),
		      PositionCanonical( D, pnt ), OnPoints );
    elif IsBound( D )  then
      return CycleLengthOp( g, D, pnt, act );
    else
      return CycleLengthOp( g, pnt, act );
    fi;
end );

InstallMethod( CycleLengthOp, true,
        [ IsObject, IsList, IsObject, IsFunction ], 0,
    function( g, D, pnt, act )
    return Length( CycleOp( g, D, pnt, act ) );
end );

InstallOtherMethod( CycleLengthOp, true,
        [ IsObject, IsObject, IsFunction ], 0,
    function( g, pnt, act )
    return Length( CycleOp( g, pnt, act ) );
end );

#############################################################################
##
#F  CycleLengths( <arg> ) . . . . . . . . . . . . . . . . . . . cycle lengths
##
InstallGlobalFunction( CycleLengths, function( arg )
    local   g,  D,  gens,  acts,  act,  xset,  hom;
    
    # test arguments
    if Length(arg)<2 or not IsMultiplicativeElementWithInverse(arg[1]) then
      Error("usage: CycleLengths(<g>,<D>[,<act>])");
    fi;
    
    # Get the arguments.
    g := arg[ 1 ];
    if IsExternalSet( arg[ 2 ] )  then
        xset := arg[ 2 ];
        D := Enumerator( xset );
        act := FunctionAction( xset );
        hom := ActionHomomorphismAttr( xset );
    else
        D := arg[ 2 ];
        if IsDomain( D )  then
            D := Enumerator( D );
        fi;
        if IsFunction( arg[ Length( arg ) ] )  then
            act := arg[ Length( arg ) ];
        else
            act := OnPoints;
        fi;
        if Length( arg ) > 3  then
            gens := arg[ 3 ];
            acts := arg[ 4 ];
            if not IsIdenticalObj( gens, acts )  then
                hom := ActionHomomorphismAttr
                       ( ExternalSetByFilterConstructor( IsExternalSet,
                         GroupByGenerators( gens ), D, gens, acts, act ) );
            fi;
        fi;
    fi;
    
    if IsBound( hom )  and  IsActionHomomorphismByActors( hom )  then
        return CycleLengthsOp( ImagesRepresentative( hom, g ),
                       [ 1 .. Length( D ) ], OnPoints );
    else
        return CycleLengthsOp( g, D, act );
    fi;
end );

InstallMethod( CycleLengthsOp, true, [ IsObject, IsList, IsFunction ], 0,
    function( g, D, act )
    return Immutable( List( CyclesOp( g, D, act ), Length ) );
end );

#############################################################################
##
#F  CycleIndex( <arg> ) . . . . . . . . . . . . . . . . . . . cycle lengths
##
InstallGlobalFunction( CycleIndex, function( arg )
local cs, g, dom, op;
  
  # get/test arguments
  cs:=Length(arg)>0 
	and (IsMultiplicativeElementWithInverse(arg[1]) or IsGroup(arg[1]));
  if cs then
    g:=arg[1];
    if Length(arg)<2 then
      cs:= IsPerm(g) or IsPermGroup(g);
      if cs then
	dom:=MovedPoints(g);
      fi;
    else
      dom:=arg[2];
    fi;

    if Length(arg)<3 then
      op:=OnPoints;
    else
      op:=arg[3];
      cs:=cs and IsFunction(op);
    fi;
  fi;
  if not cs then
    Error("usage: CycleIndex(<g>,<Omega>[,<act>])");
  fi;
  return CycleIndexOp( g, dom, op );
end );

InstallOtherMethod(CycleIndexOp,"element",true,
  [IsMultiplicativeElementWithInverse,IsListOrCollection,IsFunction ],0,
function( g, dom, act )
local c, i;
  c:=Indeterminate(Rationals,1)^0;
  for i in CycleLengthsOp(g,dom,act) do
    c:=c*Indeterminate(Rationals,i);
  od;
  return c;
end);

InstallMethod(CycleIndexOp,"finite group",true,
  [IsGroup and IsFinite,IsListOrCollection,IsFunction ],0,
function( g, dom, act )
  return 1/Size(g)*
  Sum(ConjugacyClasses(g),i->Size(i)*CycleIndexOp(Representative(i),dom,act));
end);

#############################################################################
##
#F  IsTransitive( <G>, <D>, <gens>, <acts>, <act> ) . . . . transitivity test
##
InstallMethod( IsTransitive,
    "compare with orbit of element",
    true,
    OrbitsishReq, 0,
function( G, D, gens, acts, act )
    return Length(D)=0 or IsSubset( OrbitOp( G, D[ 1 ], gens, acts, act ), D );
end );


#############################################################################
##
#F  Transitivity( <arg> ) . . . . . . . . . . . . . . . . transitivity degree
##
InstallMethod( Transitivity,"of the image of an ophom",
    true,
    OrbitsishReq, 0,
    function( G, D, gens, acts, act )
    local   hom;
    
    hom := ActionHomomorphism( G, D, gens, acts, act );
    return Transitivity( ImagesSource( hom ), [ 1 .. Length( D ) ] );
end );

InstallMethod( Transitivity,
    "G, [  ], gens, perms, act",
    true,
    [ IsGroup, IsList and IsEmpty,
      IsList,
      IsList,
      IsFunction ], 
      20, # we claim this method is very good
    function( G, D, gens, acts, act )
    return 0;
end );


#############################################################################
##
#F  IsPrimitive( <G>, <D>, <gens>, <acts>, <act> )  . . . .  primitivity test
##
InstallMethod( IsPrimitive,"transitive and no blocks",
    true,
    OrbitsishReq, 0,
    function( G, D, gens, acts, act )
    return     IsTransitive( G, D, gens, acts, act )
           and Length( Blocks( G, D, gens, acts, act ) ) <= 1;
end );


#############################################################################
##
#M  SetEarns( <G>, fail ) . . . . . . . . . . . . . . . . .  never set `fail'
##
InstallOtherMethod( SetEarns,"never set fail",
    true, [ IsGroup, IsBool ], 0,
function( G, failval )
    Setter( IsPrimitiveAffine )( G, false );
end );


#############################################################################
##
#F  IsPrimitiveAffine( <arg> )
##
InstallMethod( IsPrimitiveAffine,"primitive and earns",
    true,
    OrbitsishReq, 0,
    function( G, D, gens, acts, act )
    return     IsPrimitive( G, D, gens, acts, act )
           and Earns( G, D, gens, acts, act ) <> fail;
end );


#############################################################################
##
#F  IsSemiRegular( <arg> )  . . . . . . . . . . . . . . . semiregularity test
##
InstallMethod( IsSemiRegular,"via ophom",
    true,
    OrbitsishReq, 0,
    function( G, D, gens, acts, act )
    local   hom;
    
    hom := ActionHomomorphism( G, D, gens, acts, act );
    return IsSemiRegular( ImagesSource( hom ), [ 1 .. Length( D ) ] );
end );

InstallMethod( IsSemiRegular,
    "G, [  ], gens, perms, act",
    true,
    [ IsGroup, IsList and IsEmpty,
      IsList,
      IsList,
      IsFunction ], 
      20, # we claim this method is very good
    function( G, D, gens, acts, act )
    return true;
end );

InstallMethod( IsSemiRegular,
    "G, D, gens, [  ], act",
    true,
    [ IsGroup, IsList,
      IsList,
      IsList and IsEmpty,
      IsFunction ], 
      20, # we claim this method is very good
    function( G, D, gens, acts, act )
    return IsTrivial( G );
end );


#############################################################################
##
#F  IsRegular( <arg> )  . . . . . . . . . . . . . . . . . . . regularity test
##
InstallMethod( IsRegular,"transitive and semiregular",
    true,
    OrbitsishReq, 0,
    function( G, D, gens, acts, act )
    return     IsTransitive( G, D, gens, acts, act )
           and IsSemiRegular( G, D, gens, acts, act );
end );


#############################################################################
##
#F  RepresentativeAction( <arg> )  . . . . . . . .  representative element
##
InstallGlobalFunction( RepresentativeAction, function( arg )
local   G,  D,  d,  e,  gens,  acts,  act,  xset,  hom,  p,  rep;
    
    if IsExternalSet( arg[ 1 ] )  then
        xset := arg[ 1 ];
        d := arg[ 2 ];
        e := arg[ 3 ];
        G := ActingDomain( xset );

	# catch a trivial case (that is called from some operations often)
	if d=e then
	  return One(G);
	fi;

        if HasHomeEnumerator( xset )  then
            D := HomeEnumerator( xset );
        fi;
        if IsExternalSetByActorsRep( xset )  then
            gens := xset!.generators;
            acts := xset!.operators;
            act  := xset!.funcOperation;
        else
            act := FunctionAction( xset );
        fi;
    else
        G := arg[ 1 ];
        if Length( arg ) > 2  and
	  IsIdenticalObj( FamilyObj( arg[ 2 ] ),
		      CollectionsFamily( FamilyObj( arg[ 3 ] ) ) )  then
	  D := arg[ 2 ];
	  if IsDomain( D )  then
	      D := Enumerator( D );
	  fi;
	  p := 3;
	else
	  p := 2;
        fi;
        d := arg[ p     ];
        e := arg[ p + 1 ];

	# catch a trivial case (that is called from some operations often)
	if d=e then
	  return One(G);
	fi;

        if IsFunction( arg[ Length( arg ) ] )  then
            act := arg[ Length( arg ) ];
        else
            act := OnPoints;
        fi;
        if Length( arg ) > p + 2  then
	  gens := arg[ p + 2 ];
	  acts := arg[ p + 3 ];
	  if not IsPcgs( gens )  and  not IsIdenticalObj( gens, acts )  then
	    if not IsBound( D )  then
	      D := OrbitOp( G, d, gens, acts, act );
	      # don't make it a subset!
	      xset:=ExternalSet(G,D,gens,acts,act);
	    else
	      D := OrbitOp( G, D,d, gens, acts, act );
	      # don't make it a subset!
	      xset:=ExternalSet(G,D,gens,acts,act);
	      #xset:=ExternalOrbitOp( G, D, d, gens, acts, act );
	    fi;
	  fi;
        fi;
    fi;
    
    if IsBound(xset) and IsExternalSetByActorsRep(xset) then
      # in all other cases the homomorphism was ignored anyhow
      hom := ActionHomomorphismAttr(xset);

      d := PositionCanonical( D, d );  e := PositionCanonical( D, e );
      rep := RepresentativeActionOp( ImagesSource( hom ), d, e,
		      OnPoints );
      if rep <> fail  then
	rep := PreImagesRepresentative( hom, rep );
      fi;
      return rep;
    elif IsBound( D )  then
      if IsBound( gens )  and  IsPcgs( gens )  then
	return RepresentativeAction( G, D, d, e, gens, acts, act );
      else
	return RepresentativeActionOp( G, D, d, e, act );
      fi;
    else
        return RepresentativeActionOp( G, d, e, act );
    fi;
end );

InstallMethod( RepresentativeActionOp,"ignore domain",
    true,
    [ IsGroup, IsList, IsObject, IsObject, IsFunction ], 0,
    function( G, D, d, e, act )        
    return RepresentativeActionOp( G, d, e, act );
end );

InstallOtherMethod( RepresentativeActionOp,
  "orbit algorithm: trace transversal", true,
        [ IsGroup, IsObject, IsObject, IsFunction ], 0,
    function( G, d, e, act )
    local   rep,        # representative, result
            orb,        # orbit
            gen,        # generator of the group <G>
            pnt,        # point in the orbit <orb>
            img,        # image of the point <pnt> under the generator <gen>
            by,         # <by>[<pnt>] is a gen taking <frm>[<pnt>] to <pnt>
	    dict,	# dictionary
	    pos,	# position
            frm;        # where <frm>[<pnt>] lies earlier in <orb> than <pnt>

    d:=Immutable(d);
    e:=Immutable(e);
    dict:=NewDictionary(d,true);
    orb := [ d ];
    AddDictionary(dict,d,1);

    if act=OnPairs or act=OnTuples and CanComputeSizeAnySubgroup(G) then

      if Length( d ) <> Length( e ) then
        return fail;
      fi;

      # a well-behaving group acts on tuples. We compute the representative
      # iteratively, by mapping element for element
      rep:=One(G);
      d:=ShallowCopy(d);
      for pnt in [1..Length(d)] do
        img:=RepresentativeAction(G,d[pnt],e[pnt],OnPoints);
	if img=fail then
	  return fail;
	fi;
	rep:=rep*img;
	for pos in [pnt+1..Length(d)] do
	  d[pos]:=OnPoints(d[pos],img);
	od;
	G:=Stabilizer(G,e[pnt],OnPoints);
      od;
      return rep;
    else
      # standard action. If act is OnPoints, it should be as fast as pnt^gen.
      # So there should be no reason to split cases.
      if d = e  then return One( G );  fi;
      by  := [ One( G ) ];
      frm := [ 1 ];
      for pnt  in orb  do
	  for gen  in GeneratorsOfGroup( G )  do
	      img := act(pnt,gen);
	      MakeImmutable(img);
	      if img = e  then
		  rep := gen;
		  while pnt <> d  do
		    pos:=LookupDictionary(dict,pnt);
		    rep := by[ pos ] * rep;
		    pnt := frm[ pos ];
		  od;
		  Assert(2,act(d,rep)=e);
		  return rep;
	      elif not KnowsDictionary(dict,img) then
		  Add( orb, img );
		  AddDictionary( dict, img, Length(orb) );
		  Add( frm, pnt );
		  Add( by,  gen );
	      fi;
	  od;
      od;
      return fail;
    fi;

#    # other action
#    else
#        if d = e  then return One( G );  fi;
#        by  := [ One( G ) ];
#        frm := [ 1 ];
#        for pnt  in orb  do
#            for gen  in GeneratorsOfGroup( G )  do
#                img := act( pnt, gen );
#                if img = e  then
#                    rep := gen;
#                    while pnt <> d  do
#                        rep := by[ Position(orb,pnt) ] * rep;
#                        pnt := frm[ Position(orb,pnt) ];
#                    od;
#                    return rep;
#                elif not img in set  then
#                    Add( orb, img );
#		    if cansort then
#		      AddSet( set, img );
#		    fi;
#                    Add( frm, pnt );
#                    Add( by,  gen );
#                fi;
#            od;
#        od;
#        return fail;
#
#    fi;
#
#    # special case for action on pairs
#    elif act = OnPairs  then
#        if d = e  then return One( G );  fi;
#        by  := [ One( G ) ];
#        frm := [ 1 ];
#        for pnt  in orb  do
#            for gen  in GeneratorsOfGroup( G )  do
#                img := [ pnt[1]^gen, pnt[2]^gen ];
#                if img = e  then
#                    rep := gen;
#                    while pnt <> d  do
#                        rep := by[ Position(orb,pnt) ] * rep;
#                        pnt := frm[ Position(orb,pnt) ];
#                    od;
#                    return rep;
#                elif not img in set  then
#                    Add( orb, img );
#		    if cansort then
#		      AddSet( set, img );
#		    fi;
#                    Add( frm, pnt );
#                    Add( by,  gen );
#                fi;
#            od;
#        od;
#        return fail;
end );

#############################################################################
##
#F  Stabilizer( <arg> ) . . . . . . . . . . . . . . . . . . . . .  stabilizer
##
InstallGlobalFunction( Stabilizer, function( arg )
    if Length( arg ) = 1  then
        return StabilizerOfExternalSet( arg[ 1 ] );
    else
        return CallFuncList( StabilizerFunc, arg );
    fi;
end );

InstallMethod( StabilizerOp,
        "`OrbitStabilizerAlgorithm' with domain",true,
        OrbitishReq, 0,
function( G, D, d, gens, acts, act )
local   orbstab;
  
  orbstab:=OrbitStabilizerAlgorithm(G,D,false,gens,acts,
	      rec(pnt:=d,act:=act,onlystab:=true));
  return orbstab.stabilizer;

end );

InstallOtherMethod( StabilizerOp,
        "`OrbitStabilizerAlgorithm' without domain",true,
        [ IsGroup, IsObject, IsList, IsList, IsFunction ], 0,
function( G, d, gens, acts, act )
local   stb,  p,  orbstab;
  
  if     IsIdenticalObj( gens, acts )
    and act = OnTuples  or  act = OnPairs  then
    # for tuples compute the stabilizer iteratively
    stb := G;
    for p  in d  do
	stb := StabilizerOp( stb, p, GeneratorsOfGroup( stb ),
			GeneratorsOfGroup( stb ), OnPoints );
    od;
  else
    orbstab:=OrbitStabilizerAlgorithm(G,false,false,gens,acts,
                                      rec(pnt:=d,act:=act,onlystab:=true));
    stb := orbstab.stabilizer;
  fi;
  return stb;
end );

#############################################################################
##
#F  RankAction( <arg> ) . . . . . . . . . . . . . . . number of suborbits
##
InstallMethod( RankAction,"via ophom",
    true,
    OrbitsishReq, 0,
    function( G, D, gens, acts, act )
    local   hom;

    hom := ActionHomomorphism( G, D, gens, acts, act );
    return RankAction( Image( hom ), [ 1 .. Length( D ) ] );
end );

InstallMethod( RankAction,
    "G, ints, gens, perms, act",
    true,
    [ IsGroup, IsList and IsCyclotomicCollection,
      IsList,
      IsList,
      IsFunction ], 0,
    function( G, D, gens, acts, act )
    if    act <> OnPoints
       or not IsIdenticalObj( gens, acts )  then
        TryNextMethod();
    fi;
    return Length( OrbitsDomain( Stabilizer( G, D, D[ 1 ], act ),
                   D, act ) );
end );

InstallMethod( RankAction,
    "G, [  ], gens, perms, act",
    true,
    [ IsGroup, IsList and IsEmpty,
      IsList,
      IsList,
      IsFunction ], 
      20, # we claim this method is very good
    function( G, D, gens, acts, act )
    return 0;
end );


#############################################################################
##
#M  CanonicalRepresentativeOfExternalSet( <xset> )  . . . . . . . . . . . . .
##
InstallMethod( CanonicalRepresentativeOfExternalSet,"smallest element", true,
        [ IsExternalSet ], 0,
    function( xset )
    local   aslist;
    
    aslist := AsList( xset );
    return First( HomeEnumerator( xset ), p -> p in aslist );
end );

# for external sets that know how to get the canonical representative
InstallMethod( CanonicalRepresentativeOfExternalSet, 
      "by CanonicalRepresentativeDeterminator",
      true,
      [ IsExternalSet
        and HasCanonicalRepresentativeDeterminatorOfExternalSet ],
      0,
function( xset )
local func,can;

  func:=CanonicalRepresentativeDeterminatorOfExternalSet(xset);
  can:=func(ActingDomain(xset),Representative(xset));
  # note the stabilizer we got for free
  if not HasStabilizerOfExternalSet(xset) and IsBound(can[2]) then
    SetStabilizerOfExternalSet(xset,can[2]^(can[3]^-1));
  fi;
  return can[1];
end ) ;

#############################################################################
##
#M  ActorOfExternalSet( <xset> ) . . . . . . . . . . . . . . . . . . . . .
##
InstallMethod( ActorOfExternalSet, true, [ IsExternalSet ], 0,
    xset -> RepresentativeAction( xset, Representative( xset ),
            CanonicalRepresentativeOfExternalSet( xset ) ) );

#############################################################################
##
#M  StabilizerOfExternalSet( <xset> ) . . . . . . . . . . . . . . . . . . . .
##
InstallMethod( StabilizerOfExternalSet,"stabilizer of the represenattive",
  true, [ IsExternalSet ], 0,
        xset -> Stabilizer( xset, Representative( xset ) ) );

#############################################################################
##
#M  ImageElmActionHomomorphism( <hom>, <elm> )
##
InstallGlobalFunction(ImageElmActionHomomorphism,function( hom, elm )
local   xset,he,gp,p,canfail,bas,fun;
  canfail:=ValueOption("actioncanfail")=true;
  xset := UnderlyingExternalSet( hom );
  he:=HomeEnumerator(xset);
  fun:=FunctionAction(xset);
  # can we compute the image cheaper using stabilizer chain methods?
  if not canfail and HasImagesSource(hom) then
    gp:=ImagesSource(hom);
    if HasStabChainMutable(gp) or HasStabChainImmutable(gp) then
      bas:=BaseStabChain(StabChain(gp));
      if Length(bas)*50<Length(he) then
	p:=RepresentativeActionOp(gp,bas,
	    List(bas,x->PositionCanonical(he,fun(he[x],elm))),
	    OnTuples);
        return p;
      fi;
    fi;
  fi;

  p:=Permutation(elm,he,fun);
  if p=fail then
    if canfail then
      return fail;
    fi;
    Error("Action not well-defined. See the manual section\n",
	  "``Action on canonical representatives''.");
  fi;
  return p;
end );

#############################################################################
##
#M  ImagesRepresentative( <hom>, <elm> )  . . . . . . . . . for action hom
##
InstallMethod( ImagesRepresentative,"for action hom", FamSourceEqFamElm,
        [ IsActionHomomorphism, IsMultiplicativeElementWithInverse ], 0,
  ImageElmActionHomomorphism);

InstallMethod( ImagesRepresentative, "for action hom that is `ByAsGroup'",
  FamSourceEqFamElm,
  [ IsGroupGeneralMappingByAsGroupGeneralMappingByImages 
    and IsActionHomomorphism, IsMultiplicativeElementWithInverse ], 0,
function( hom, elm )
  return ImagesRepresentative( AsGroupGeneralMappingByImages( hom ), elm );
end );


#############################################################################
##
#M  MappingGeneratorsImages( <map> )  . . . . .  for group homomorphism
##
InstallMethod( MappingGeneratorsImages, "for action hom that is `ByAsGroup'",
    true, [ IsGroupGeneralMappingByAsGroupGeneralMappingByImages and
	    IsActionHomomorphism ], 0,
function( map )
local gens;
  gens:= GeneratorsOfGroup( PreImagesRange( map ) );
  return [gens, List( gens, g -> ImageElmActionHomomorphism( map, g ) ) ];
end );

#############################################################################
##
#M  KernelOfMultiplicativeGeneralMapping( <ophom>, <elm> )
##
InstallMethod( KernelOfMultiplicativeGeneralMapping,
  "for action homomorphism", true, [ IsActionHomomorphism ], 0,
function( hom )
local map,mapi;
  if HasIsSurjective(hom) and IsSurjective(hom) then
    return KernelOfMultiplicativeGeneralMapping(
	    AsGroupGeneralMappingByImages( hom ) );
  else
    Range( hom );
    mapi := MappingGeneratorsImages( hom );
    map := GroupHomomorphismByImagesNC( Source( hom ), ImagesSource( hom ),
              mapi[1], mapi[2] );
    CopyMappingAttributes( hom,map );
    return KernelOfMultiplicativeGeneralMapping(map);
  fi;
end );

#############################################################################
##
#M  ImagesRepresentative( <hom>, <elm> )  . . . . . . . .  if a base is known
##
InstallMethod( ImagesRepresentative, "using `RepresentativeAction'",
  FamSourceEqFamElm, [ IsActionHomomorphismByBase and HasImagesSource,
          IsMultiplicativeElementWithInverse ], 0,
    function( hom, elm )
    local   xset,  D,  act,  imgs;

    xset := UnderlyingExternalSet( hom );
    D := HomeEnumerator( xset );
    act := FunctionAction( xset );
    if not IsBound( xset!.basePermImage )  then
        xset!.basePermImage := List( BaseOfGroup( xset ),
                                    b -> PositionCanonical( D, b ) );
    fi;
    imgs:=List(BaseOfGroup(xset),b->PositionCanonical(D,act(b,elm)));
    return RepresentativeActionOp( ImagesSource( hom ),
                   xset!.basePermImage, imgs, OnTuples );
end );

#############################################################################
##
#M  ImagesSource( <hom> ) . . . . . . . . . . . . . . . . . set base in image
##
InstallMethod( ImagesSource,"actionHomomorphismByBase", true,
        [ IsActionHomomorphismByBase ], 0,
    function( hom )
    local   xset,  img,  D;
    
    xset := UnderlyingExternalSet( hom );
    img := ImagesSet( hom, Source( hom ) );
    if not HasStabChainMutable( img )  and  not HasBaseOfGroup( img )  then
        if not IsBound( xset!.basePermImage )  then
            D := HomeEnumerator( xset );
            xset!.basePermImage := List( BaseOfGroup( xset ),
                                        b -> PositionCanonical( D, b ) );
        fi;
        SetBaseOfGroup( img, xset!.basePermImage );
    fi;
    return img;
end );
    
#############################################################################
##
#M  ImagesRepresentative( <hom>, <elm> )  . . . . .  restricted `Permutation'
##
InstallMethod( ImagesRepresentative,"restricted perm", FamSourceEqFamElm,
        [ IsActionHomomorphismSubset,
          IsMultiplicativeElementWithInverse ], 0,
    function( hom, elm )
    local   xset;
    
    xset := UnderlyingExternalSet( hom );
    return RestrictedPermNC( Permutation( elm, HomeEnumerator( xset ),
        FunctionAction( xset ) ),
        MovedPoints( ImagesSource( AsGroupGeneralMappingByImages( hom ) ) ) );
end );

#############################################################################
##
#M  PreImagesRepresentative( <hom>, <elm> ) . . . . . . . . . .  build matrix
##
InstallMethod( PreImagesRepresentative,"IsLinearActionHomomorphism",
  FamRangeEqFamElm, [ IsLinearActionHomomorphism, IsPerm ], 0,
function( hom, elm )
  local   V,  base,  mat,  b,xset,lab,f;
  
  # is this method applicable? Test whether the domain contains a vector
  # space basis (respectively just get this basis).
  lab:=LinearActionBasis(hom);
  if lab=fail then
    TryNextMethod();
  fi;

  # PreImagesRepresentative does not test membership
  #if not elm in Image( hom )  then return fail; fi;
  xset:=UnderlyingExternalSet(hom);
  V := HomeEnumerator(xset);
  f:=DefaultFieldOfMatrixGroup(Source(hom));

  if not IsBound(hom!.linActBasisPositions) then
    hom!.linActBasisPositions:=List(lab,i->PositionCanonical(V,i));
  fi;
  if not IsBound(hom!.linActInverse) then
    lab:=ImmutableMatrix(f,lab);
    hom!.linActInverse:=Inverse(lab);
  fi;

  elm:=OnTuples(hom!.linActBasisPositions,elm); # image points
  elm:=V{elm}; # the corresponding vectors
  f:=DefaultFieldOfMatrixGroup(Source(hom));
  elm:=ImmutableMatrix(f,elm);

  return hom!.linActInverse*elm;
end );

#############################################################################
##
#M  PreImagesRepresentative( <hom>, <elm> ) . . . . . . . . . .  build matrix
##
InstallMethod( PreImagesRepresentative,"IsProjectiveActionHomomorphism",
  FamRangeEqFamElm, [ IsProjectiveActionHomomorphism, IsPerm ], 0,
function( hom, elm )
  local   V,  base,  mat,  b,xset,lab,f,dim,start,time,sol,i;
  
  # is this method applicable? Test whether the domain contains a vector
  # space basis (respectively just get this basis).
  lab:=LinearActionBasis(hom);
  if lab=fail then
    TryNextMethod();
  fi;

  # PreImagesRepresentative does not test membership
  #if not elm in Image( hom )  then return fail; fi;
  xset:=UnderlyingExternalSet(hom);
  V := HomeEnumerator(xset);
  f:=DefaultFieldOfMatrixGroup(Source(hom));
  dim:=DimensionOfMatrixGroup(Source(hom));

  elm:=OnTuples(hom!.projActBasisPositions,elm); # image points
  elm:=V{elm}; # the corresponding vectors

  mat:=elm{[1..dim]};
  sol:=SolutionMat(mat,elm[dim+1]);
  for i in [1..dim] do
    mat[i]:=sol[i]*mat[i];
  od;
  mat:=hom!.projActInverse*ImmutableMatrix(f,mat);

  # correct scalar using determinant if needed
  if hom!.correctionFactors[1]<>fail then
    V:=DeterminantMat(mat);
    if not IsOne(V) then
      mat:=mat*hom!.correctionFactors[2][
	      PositionSorted(hom!.correctionFactors[1],V)];
    fi;
  fi;

  return mat;
end);

#############################################################################
##
#A  LinearActionBasis(<hom>)
##
InstallMethod(LinearActionBasis,"find basis in domain",true,
  [IsLinearActionHomomorphism],0,
function(hom)
local xset,dom,D,b,t,i,r,pos;
  xset:=UnderlyingExternalSet(hom);
  if Size(xset)=0 then
    return fail;
  fi;
  pos:=[];
  # if there is a base, check whether it's full rank, if yes, take it
  if HasBaseOfGroup(xset) 
     and RankMat(BaseOfGroup(xset))=Length(BaseOfGroup(xset)[1]) then
    # this implies injectivity
    SetIsInjective(hom,true);
    return BaseOfGroup(xset);
  fi;
  # otherwise we've to find a basis from the domain.
  D:=HomeEnumerator(xset);
  b:=[];
  t:=[];
  r:=Length(D[1]);
  i:=1;
  while Length(b)<r and i<=Length(D) do
    if RankMat(Concatenation(t,[D[i]]))>Length(t) then
      # new indep. vector
      Add(b,D[i]);
      Add(pos,i);
      Add(t,ShallowCopy(D[i]));
      TriangulizeMat(t); # for faster rank tests
    fi;
    i:=i+1;
  od;
  if Length(b)=r then
    # this implies injectivity
    hom!.linActBasisPositions:=pos;
    SetIsInjective(hom,true);
    return b;
  else
    return fail; # does not span
  fi;
end);

#############################################################################
##
#A  LinearActionBasis(<hom>)
##
InstallOtherMethod(LinearActionBasis,"projective with extra vector",true,
  [IsProjectiveActionHomomorphism],0,
function(hom)
local xset,dom,D,b,t,i,r,binv,pos,kero,dets,roots,dim,f;
  xset:=UnderlyingExternalSet(hom);
  if Size(xset)=0 then
    return fail;
  fi;

  # will the determinants suffice to get suitable scalars?
  dim:=DimensionOfMatrixGroup(Source(hom));
  f:=DefaultFieldOfMatrixGroup(Source(hom));

  roots:=Set(RootsOfUPol(f,X(f)^dim-1));

  D:=List(GeneratorsOfGroup(Source(hom)),DeterminantMat);
  D:=Elements(Group(D));

  if Length(roots)<=1 then
    # 1 will always be root
    kero:=[One(f)];
  elif IsNaturalGL(Source(hom)) then
    # the full GL clearly will contain the kernel
    kero:=roots; # to skip test
  elif not IsSubset(D,roots) then
    # even the kernel determinants are not reached, so clearly kernel not in
    return fail;
  else
    kero:=List(Elements(KernelOfMultiplicativeGeneralMapping(hom)),x->x[1][1]^dim);
  fi;
  
  if not IsSubset(kero,roots) then
    # we cannot fix the scalar with the determinant
    return fail;
  fi;

  dets:=[];
  roots:=[];
  for i in Filtered(Elements(f),x->not IsZero(x)) do
    b:=i^dim;
    if not b in dets then
      Add(dets,b);
      Add(roots,i^-1); # the factor by which we must correct
    fi;
  od;
  SortParallel(dets,roots);

  if IsSubset(D,dets) then
    dets:=fail; # not that we do not need to correct with determinant as all
                # values are fine
  fi;

  # find a basis from the domain.
  D:=HomeEnumerator(xset);
  b:=[];
  t:=[];
  r:=Length(D[1]);
  i:=1;
  pos:=[];
  while Length(b)<r and i<=Length(D) do
    if RankMat(Concatenation(t,[D[i]]))>Length(t) then
      # new indep. vector
      Add(b,D[i]);
      Add(pos,i);
      Add(t,ShallowCopy(D[i]));
      TriangulizeMat(t); # for faster rank tests
    fi;
    i:=i+1;
  od;
  if Length(b)<r then
    return fail;
  fi;

  # try to find a vector that has nonzero coefficients for all b
  binv:=Inverse(ImmutableMatrix(f,b));
  while i<=Length(D) do
    if ForAll(D[i]*binv,x->not IsZero(x)) then
      Add(b,D[i]);
      Add(pos,i);
      hom!.projActBasisPositions:=pos;
      hom!.projActInverse:=ImmutableMatrix(f,binv*Inverse(DiagonalMat(D[i]*binv)));
      hom!.correctionFactors:=[dets,roots];
      return ImmutableMatrix(f,b);
    fi;
    i:=i+1;
  od;

  return fail; # no extra vector found
end);

#############################################################################
##
#M  DomainForAction( <pnt>, <acts>,<act> )
##
InstallMethod(DomainForAction,"default: fail",true,
  [IsObject,IsListOrCollection,IsFunction],0,
function(pnt,acts,act)
  return fail;
end);

#############################################################################
##
#M  AbelianSubfactorAction(<G>,<M>,<N>)
##
InstallMethod(AbelianSubfactorAction,"generic:use modulo pcgs",true,
  [IsGroup,IsGroup,IsGroup],0,
function(G,M,N)
local p,n,f,o,v,ran,exp,H,phi,alpha;
  p:=ModuloPcgs(M,N);
  n:=Length(p);
  f:=GF(RelativeOrders(p)[1]);
  o:=One(f);
  v:=f^n;
  ran:=[1..n];
  exp:=ListWithIdenticalEntries(n,0);
  f:=Size(f);
  phi:=LinearActionLayer(G,p);
  H:=Group(phi);
  UseFactorRelation(G,fail,H);
  phi:=GroupHomomorphismByImagesNC(G,H,GeneratorsOfGroup(G),phi);

  alpha:=GroupToAdditiveGroupHomomorphismByFunction(M,v,function(e)
    e:=ExponentsOfPcElement(p,e)*o;
    MakeImmutable(e);
    ConvertToVectorRep(e,f);
    return e;
  end,
  function(r)
  local i,l;
    l:=exp;
    for i in ran do
      l[i]:=Int(r[i]);
    od;
    return PcElementByExponentsNC(p,l);
  end);
  return [phi,alpha,p];
end);


#############################################################################
##
#M  IsInjective( <acthom> )
##
##  This is triggered by a fallback method of PreImageElm if it is not
##  yet known whether the action homomorphism is injective or not.
##  If there exists a LinearActionBasis, then the hom is injective
##  and the better method for PreImageElm is taken.
##
InstallMethod( IsInjective, "for a linear action homomorphism",
  [IsLinearActionHomomorphism],
  function( a )
    local b;   
    b := LinearActionBasis(a);
    if b = fail then          
        TryNextMethod();
    fi;                 
    return true;
  end );


#############################################################################
##
#E