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

Quelle  morpheus.gi   Sprache: unbekannt

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

#############################################################################
##
#V  MORPHEUSELMS . . . .  limit up to which size to store element lists
##
MORPHEUSELMS := 50000;

# this method calculates a chief series invariant under `hom` and calculates
# orders of group elements in factors of this series under action of `hom`.
# Every time an orbit length is found, `hom` is replaced by the appropriate
# power. Initially small chief factors are preferred. In the end all
# generators are used while stepping through the series descendingly, thus
# ensuring the proper order is found.
InstallMethod(Order,"for automorphisms",true,[IsGroupHomomorphism],0,
function(hom)
local map,phi,o,lo,i,j,start,img,d,nat,ser,jord,first;
  d:=Source(hom);
  if not (HasIsFinite(d) and IsFinite(d)) then
    TryNextMethod();
  fi;
  if Size(d)<=10000 then
    ser:=[d,TrivialSubgroup(d)]; # no need to be clever if small
  else
    if HasAutomorphismGroup(d) then
      if IsBound(d!.characteristicSeries) then
        ser:=d!.characteristicSeries;
      else
        ser:=ChiefSeries(d); # could try to be more clever, introduce attribute
        # `CharacteristicSeries`.
        ser:=Filtered(ser,x->ForAll(GeneratorsOfGroup(AutomorphismGroup(d)),
          a->ForAll(GeneratorsOfGroup(x),y->ImageElm(a,y) in x)));
        d!.characteristicSeries:=ser;
      fi;
    else
      ser:=ChiefSeries(d); # could try to be more clever, introduce attribute
      # `CharacteristicSeries`.
      ser:=Filtered(ser,
            x->ForAll(GeneratorsOfGroup(x),y->ImageElm(hom,y) in x));
    fi;
  fi;

  # try to do factors in ascending order in the hope to get short orbits
  # first
  jord:=[2..Length(ser)]; # order in which we go through factors
  if Length(ser)>2 then
    i:=List(jord,x->Size(ser[x-1])/Size(ser[x]));
    SortParallel(i,jord);
  fi;

  o:=1;
  phi:=hom;
  map:=MappingGeneratorsImages(phi);

  first:=true;
  while map[1]<>map[2] do
    for j in jord do
      i:=1;
      while i<=Length(map[1]) do
        # the first time, do only the generators from prior layer
        if (not first)
           or (map[1][i] in ser[j-1] and not map[1][i] in ser[j]) then

          lo:=1;
          if j<Length(ser) then
            nat:=NaturalHomomorphismByNormalSubgroup(d,ser[j]);
            start:=ImagesRepresentative(nat,map[1][i]);
            img:=map[2][i];
            while ImagesRepresentative(nat,img)<>start do
              img:=ImagesRepresentative(phi,img);
              lo:=lo+1;

              # do the bijectivity test only if high local order, then it
              # does not matter. IsBijective is cached, so second test is
              # cheap.
              if lo=1000 and not IsBijective(hom) then
                Error("<hom> must be bijective");
              fi;

            od;

          else
            start:=map[1][i];
            img:=map[2][i];
            while img<>start do
              img:=ImagesRepresentative(phi,img);
              lo:=lo+1;

              # do the bijectivity test only if high local order, then it
              # does not matter. IsBijective is cached, so second test is
              # cheap.
              if lo=1000 and not IsBijective(hom) then
                Error("<hom> must be bijective");
              fi;

            od;
          fi;

          if lo>1 then
            o:=o*lo;
            #if i<Length(map[1]) then
              phi:=phi^lo;
              map:=MappingGeneratorsImages(phi);
              i:=0; # restart search, as generator set may have changed.
            #fi;
          fi;
        fi;
        i:=i+1;
      od;
    od;

    # if iterating make `jord` standard to we don't skip generators
    jord:=[2..Length(ser)];
    first:=false;
  od;

  return o;
end);

#############################################################################
##
#M  AutomorphismDomain(<G>)
##
##  If <G> consists of automorphisms of <H>, this attribute returns <H>.
InstallMethod( AutomorphismDomain, "use source of one",true,
  [IsGroupOfAutomorphisms],0,
function(G)
  return Source(One(G));
end);

DeclareRepresentation("IsActionHomomorphismAutomGroup",
  IsActionHomomorphismByBase,["basepos"]);

#############################################################################
##
#M  IsGroupOfAutomorphisms(<G>)
##
InstallMethod( IsGroupOfAutomorphisms, "test generators and one",true,
  [IsGroup],0,
function(G)
local s;
  if IsGeneralMapping(One(G)) then
    s:=Source(One(G));
    if Range(One(G))=s and ForAll(GeneratorsOfGroup(G),
      g->IsGroupGeneralMapping(g) and IsSPGeneralMapping(g) and IsMapping(g)
         and IsInjective(g) and IsSurjective(g) and Source(g)=s
         and Range(g)=s) then
      SetAutomorphismDomain(G,s);
      # imply finiteness
      if IsFinite(s) then
        SetIsFinite(G,true);
      fi;
      return true;
    fi;
  fi;
  return false;
end);

#############################################################################
##
#M  IsGroupOfAutomorphismsFiniteGroup(<G>)
##
InstallMethod( IsGroupOfAutomorphismsFiniteGroup,"default",true,
  [IsGroup],0,
  G->IsGroupOfAutomorphisms(G) and IsFinite(AutomorphismDomain(G)));

# Try to embed automorphisms into wreath product.
BindGlobal("AutomorphismWreathEmbedding",function(au,g)
local gens, inn,out, nonperm, syno, orb, orbi, perms, free, rep, i, maxl, gen,
      img, j, conj, sm, cen, n, w, emb, ge, no,reps,synom,ginn,oemb,act,
      ogens,genimgs,oo,op;

  gens:=GeneratorsOfGroup(g);
  if Size(Centre(g))>1 then
    return fail;
  fi;
  #sym:=SymmetricGroup(MovedPoints(g));
  #syno:=Normalizer(sym,g);

  inn:=Filtered(GeneratorsOfGroup(au),IsInnerAutomorphism);
  out:=Filtered(GeneratorsOfGroup(au),i->not IsInnerAutomorphism(i));
  nonperm:=Filtered(out,i->not IsConjugatorAutomorphism(i));
  syno:=g;
  #syno:=Group(List(Filtered(GeneratorsOfGroup(au),IsInnerAutomorphism),
#        ConjugatorOfConjugatorIsomorphism),One(g));
  for i in Filtered(out,IsConjugatorAutomorphism) do
    syno:=ClosureGroup(syno,ConjugatorOfConjugatorIsomorphism(i));
  od;
  #nonperm:=Filtered(out,i->not IsInnerAutomorphism(i));
  # enumerate cosets of subgroup of conjugator isomorphisms
  orb:=[IdentityMapping(g)];
  orbi:=[IdentityMapping(g)];
  perms:=List(nonperm,i->[]);
  free:=FreeGroup(Length(nonperm));
  rep:=[One(free)];
  i:=1;
  maxl:=NrMovedPoints(g);
  while i<=Length(orb) and Length(orb)<maxl do
    for w in [1..Length(nonperm)] do
      gen:=nonperm[w];
      img:=orb[i]*gen;
      j:=1;
      conj:=fail;
      while conj=fail and j<=Length(orb) do
        sm:=img*orbi[j];
        if IsConjugatorAutomorphism(sm) then
          conj:=ConjugatorOfConjugatorIsomorphism(sm);
        else
          j:=j+1;
        fi;
      od;
      #j:=First([1..Length(orb)],k->IsConjugatorAutomorphism(img*orbi[k]));
      if conj=fail then
        Add(orb,img);
        Add(orbi,InverseGeneralMapping(img));
        Add(rep,rep[i]*GeneratorsOfGroup(free)[w]);
        perms[w][i]:=Length(orb);
      else
        perms[w][i]:=j;
        if not conj in syno then
          syno:=ClosureGroup(syno,conj);
        fi;
      fi;
    od;
    i:=i+1;
  od;

  if not IsTransitive(syno,MovedPoints(syno)) then
    return fail;
    # # characteristic product?
    # w:=Orbits(syno,MovedPoints(syno));
    # n:=List(w,x->Stabilizer(g,Difference(MovedPoints(g),x),OnTuples));
    # if ForAll(n,x->ForAll(GeneratorsOfGroup(au),y->Image(y,x)=x)) then
    #   Error("WR5");
    # fi;
  fi;
  if Length(GeneratorsOfGroup(syno))>5 then
    syno:=Group(SmallGeneratingSet(syno));
  fi;

  cen:=Centralizer(syno,g);
  Info(InfoMorph,2,"|syno|=",Size(syno)," |cen|=",Size(cen));
  if Size(cen)>1 then
    w:=syno;
    syno:=ComplementClassesRepresentatives(syno,cen);
    if Length(syno)=0 then
      return fail; # not unique permauts
    fi;
    syno:=syno[1];
    synom:=GroupHomomorphismByImagesNC(w,syno,
            Concatenation(GeneratorsOfGroup(syno),GeneratorsOfGroup(cen)),
            Concatenation(GeneratorsOfGroup(syno),List(GeneratorsOfGroup(cen),x->One(syno))));
  else
    synom:=IdentityMapping(syno);
  fi;

  # try wreath embedding
  if Length(orb)<maxl then
    Info(InfoMorph,1,Length(orb)," copies");
    perms:=List(perms,PermList);
    Info(InfoMorph,2,List(rep,i->MappedWord(i,GeneratorsOfGroup(free),perms)));
    n:=Length(orb);
    w:=WreathProduct(syno,SymmetricGroup(n));
    no:=KernelOfMultiplicativeGeneralMapping(Projection(w));
    gens:=SmallGeneratingSet(g);
    emb:=List(gens,
          i->Product(List([1..n],j->Image(Embedding(w,j),Image(synom,Image(orbi[j],i))))));
    ge:=SubgroupNC(w,emb);
    emb:=GroupHomomorphismByImagesNC(g,ge,gens,emb);
    reps:=[];
    oo:=List(Orbits(no,MovedPoints(no)),Set);
    ForAll(oo,IsRange);
    act:=[];
    for i in out do

      # how does it permute the components?
      conj:=List(orb,x->x*i);
      conj:=List(conj,x->PositionProperty(orb,
        y->IsConjugatorAutomorphism(x/y)));
      conj:=PermList(conj);
      if conj=fail then return fail;fi;
      conj:=ImagesRepresentative(Embedding(w,n+1),conj);

      #gen:=RepresentativeAction(no,GeneratorsOfGroup(ge),
      #  List(gens,j->Image(emb,ImagesRepresentative(i,j))^(conj^-1)),OnTuples);
      ogens:=GeneratorsOfGroup(ge);
      genimgs:=List(gens,j->Image(emb,ImagesRepresentative(i,j))^(conj^-1));
      gen:=One(no);
      for op in [1..Length(oo)] do
        if not IsBound(act[op]) then
          Info(InfoMorph,2,"oo=",op,oo[op]);
          act[op]:=Group(List(GeneratorsOfGroup(no),x->RestrictedPerm(x,oo[op])));
          SetSize(act,Size(syno));
        fi;

        sm:=RepresentativeAction(act[op],
          List(ogens,x->RestrictedPerm(x,oo[op])),
          List(genimgs,x->RestrictedPerm(x,oo[op])),OnTuples);
        if not IsPerm(sm) then return fail;fi;
        #ogens:=List(ogens,x->x^sm);
        gen:=gen*sm;
      od;

      if not OnTuples(ogens,gen)=genimgs then
        Error("conjugation error!");
      fi;

      #if not IsPerm(gen) then return fail;fi;
      gen:=gen*conj;
      Add(reps,gen);
    od;

    #no:=Normalizer(w,ge);
    #no:=ClosureGroup(ge,reps);
    ginn:=List(inn,ConjugatorOfConjugatorIsomorphism);
    no:=Group(List(ginn,i->Image(emb,i)), One(w));
    oemb:=emb;
    if Size(no)<Size(ge) then
      emb:=RestrictedMapping(emb,Group(ginn,()));
    fi;
    no:=ClosureGroup(no,reps);
    cen:=Centralizer(no,ge);
    if Size(no)/Size(cen)<Length(orb) then
      return fail;
    fi;

    if Size(cen)>1 then
      no:=ComplementClassesRepresentatives(no,cen);
      if Length(no)>0 then
        no:=no[1];
      else
        return fail;
      fi;
    fi;
    #
    #if Size(no)/Size(syno)<>Length(orb) then
    #  Error("wreath embedding failed");
    #fi;
    sm:=SmallerDegreePermutationRepresentation(ClosureGroup(ge,no):cheap);
    no:=Image(sm,no);
    if IsIdenticalObj(emb,oemb) then
      emb:=emb*sm;
      return [no,emb,emb,Image(emb,ginn)];
    else
      emb:=emb*sm;
      oemb:=oemb*sm;
      return [no,emb,oemb,Group(Image(oemb,ginn),One(w))];
    fi;
  fi;
  return fail;
end);

#############################################################################
##
#F  AssignNiceMonomorphismAutomorphismGroup(<autgrp>,<g>)
##
# try to find a small faithful action for an automorphism group
InstallGlobalFunction(AssignNiceMonomorphismAutomorphismGroup,function(au,g)
local hom, gens, c, ran, r, cen, img, u, orbs,
      ser,pos, o, i, j, best,actbase,action,finish,
      bestdeg,deg,baddegree,auo,of,cl,store, offset,claselms,fix,new,
      preproc,postproc;

  finish:=function(hom)
    SetIsGroupHomomorphism(hom,true);
    SetIsBijective(hom,true);
    SetFilterObj(hom,IsNiceMonomorphism);
    SetNiceMonomorphism(au,hom);
    SetIsHandledByNiceMonomorphism(au,true);
  end;

  # avoid storing list of class elements anew
  claselms:=function(clas)
    if HasAsSSortedList(clas) then
      return AsSSortedList(clas);
    else
      return MakeImmutable(Set(Orbit(g,Representative(clas))));
    fi;
  end;

  hom:=fail;

  if not IsFinite(g) then
    Error("can't do!");
  else
    SetIsFinite(au,true);
  fi;

  if IsFpGroup(g) then
    # no sane person should work with automorphism groups of fp groups, as
    # this is doubly inefficient, but if someone really does, this is a
    # shortcut that avoids canonization issues.
    c:=Filtered(Elements(g),x->not IsOne(x));
    hom:=ActionHomomorphism(au,c,
           function(e,a) return Image(a,e);end,"surjective");
    finish(hom); return;
  elif IsAbelian(g) or (Size(g)<50000 and Size(DerivedSubgroup(g))^2<Size(g)) then
    # for close to abelian groups, just act on orbits of generators
    if IsAbelian(g) then
      gens:=IndependentGeneratorsOfAbelianGroup(g);
    else
      gens:=SmallGeneratingSet(g);
    fi;
    c:=[];
    for i in gens do
      c:=Union(c,Orbit(au,i));
    od;
    hom:=NiceMonomorphismAutomGroup(au,c,gens);
    finish(hom); return;
  fi;

  # if no centre and all automorphism conjugator, try to extend exiting permrep
  if Size(Centre(g))=1 and IsPermGroup(g) and
     ForAll(GeneratorsOfGroup(au),IsConjugatorAutomorphism) then
    ran:= Group( List( GeneratorsOfGroup( au ),
                      ConjugatorOfConjugatorIsomorphism ),
                One( g ) );
    Info(InfoMorph,1,"All automorphisms are conjugator");
    Size(ran); #enforce size calculation

    # if `ran' has a centralizing bit, we're still out of luck.
    # TODO: try whether there is a centralizer complement into which we
    # could go.

    if Size(Centralizer(ran,g))=1 then
      r:=ran; # the group of conjugating elements so far
      cen:=TrivialSubgroup(r);

      hom:=GroupHomomorphismByFunction(au,ran,
        function(auto)
          if not IsConjugatorAutomorphism(auto) then
            return fail;
          fi;
          img:=ConjugatorOfConjugatorIsomorphism( auto );
          if not img in ran then
            # There is still something centralizing left.
            if not img in r then
              # get the cenralizing bit
              r:=ClosureGroup(r,img);
              cen:=Centralizer(r,g);
            fi;
            # get the right coset element
            img:=First(List(Enumerator(cen),i->i*img),i->i in ran);
          fi;
          return img;
        end,
        function(elm)
          return ConjugatorAutomorphismNC( g, elm );
        end);
      SetRange( hom,ran );
      finish(hom); return;
    fi;
  fi;

  actbase:=ValueOption("autactbase");

  bestdeg:=infinity;
  # what degree would we consider immediately acceptable?
  if actbase<>fail then
    if ForAny(actbase,x->not IsSubset(g,x)) then
      Error("illegal actbase given!");
    fi;
    baddegree:=RootInt(Sum(actbase,Size)^2,3);
  else
    baddegree:=RootInt(Size(g)^3,4);
  fi;
  if IsPermGroup(g) then baddegree:=Minimum(baddegree,Maximum(2000,NrMovedPoints(g)*10));fi;
    Info(InfoMorph,4,"degree limit ",baddegree);

  ser:=StructuralSeriesOfGroup(g);
  # eliminate very small subgroups -- there are just very few elements in
  r:=RootInt(Size(ser),4);
  ser:=Filtered(ser,x->Size(x)=1 or Size(x)>r);
  ser:=List(ser,x->SubgroupNC(g,SmallGeneratingSet(x)));

  pos:=2;

  # try action on elements, through orbits on classes
  auo:=Group(Filtered(GeneratorsOfGroup(au),x->not IsInnerAutomorphism(x)),One(au));

  if IsPermGroup(g) and Size(SolvableRadical(g))^2>Size(g) then
    FittingFreeLiftSetup(g);
    preproc:=x->TFCanonicalClassRepresentative(g,[Representative(x)])[1][2];
    postproc:=rep->ConjugacyClass(g,rep);
    action:=function(crep,aut)
      return TFCanonicalClassRepresentative(g,[ImagesRepresentative(aut,crep)])[1][2];
    end;
  else
    action:=function(class,aut)
      local c,s;
      c:=ConjugacyClass(g,ImagesRepresentative(aut,Representative(class)));
      if HasSize(class) then
        SetSize(c,Size(class));
      fi;
      if HasStabilizerOfExternalSet(class) then
        Size(StabilizerOfExternalSet(class)); # force size known to transfer
        s:=Image(aut,StabilizerOfExternalSet(class));
        SetStabilizerOfExternalSet(c,s);
      fi;
      return c;
    end;
    preproc:=fail;
    postproc:=fail;
  fi;

  # first try classes that generate, starting with small ones
  if actbase<>fail then
    c:=Concatenation(List(actbase,GeneratorsOfGroup));
  else
    c:=GeneratorsOfGroup(g);
  fi;
  # split c into prime power parts
  c:=Concatenation(List(c,PrimePowerComponents));

  c:=Unique(List(c,x->ConjugacyClass(g,x)));
  SortBy(c,Size);
  r:=Filtered([1..Length(c)],x->Representative(c[x]) in ser[pos]);
  c:=c{Concatenation(Difference([1..Length(c)],r),r)};
  o:=[];
  fix:=Filtered(List(c,Representative),
    x->ForAll(GeneratorsOfGroup(au),z->ImagesRepresentative(z,x)=x));
  u:=SubgroupNC(g,fix);

  while Size(u)<Size(g) do
    if Size(u)>1 then SortBy(c,Size);fi; # once an outside class is known, sort
    of:=First(c,x->not Representative(x) in u);

    if of=fail then
      # rarely reps are in u
      of:=First(c,
        x->ForAny(GeneratorsOfGroup(g),y->not Representative(x)^y in u));
    fi;
    if MemoryUsage(Representative(of))*Size(of)
      # if a single class only requires
      <10000
      # for element storage, just store its elements
      then
      of:=Orbit(auo,claselms(of),OnSets);
    elif preproc<>fail then
      of:=List(Orbit(auo,preproc(of),action),postproc);
    else
      of:=Orbit(auo,of,action);
    fi;
    Info(InfoMorph,4,"Genclass orbit ",Length(of)," size ",Sum(of,Size));
    Append(o,of);
    u:=ClosureGroup(u,List(of,Representative));
    u:=NormalClosure(g,u);
  od;

  bestdeg:=Sum(o,Size);
  Info(InfoMorph,4,"Generators degree ",bestdeg);
  store:=o;
  best:=[function() return Union(List(store,claselms));end,OnPoints];
  if bestdeg<baddegree then
    hom:=ActionHomomorphism(au,Union(List(o,claselms)),"surjective");
    finish(hom); return;
  fi;

  # fuse classes under g
  if actbase<>fail then
    if HasConjugacyClasses(g) then
      cl:=Filtered(ConjugacyClasses(g),x->ForAny(actbase,y->Representative(x) in y));
    else
      cl:=Concatenation(List(actbase,x->Filtered(ConjugacyClasses(x),y->Order(Representative(y))>1)));

      o:=cl;
      cl:=[];
      for i in o do
# test to avoid duplicates is more expensive than it is worth
#        if not ForAny(cl,x->Order(Representative(i))=Order(Representative(x))
#          and (Size(x) mod Size(i)=0) and Representative(i) in x) then
          r:=ConjugacyClass(g,Representative(i));
          Add(cl,r);
#        fi;
      od;

      Info(InfoMorph,2,"actbase ",Length(cl), " classes");
    fi;
  else
    cl:=ShallowCopy(ConjugacyClasses(g));
    Info(InfoMorph,2,Length(cl)," classes");
  fi;
  SortBy(cl,Size);

  # split classes in patterns
  o:=[];
  offset:=0; # degree that comes for minimum generating just factor
  i:=First([1..Length(cl)],x->not Representative(cl[x]) in ser[pos]);
  while i<>fail do
    r:=[Order(Representative(cl[i])),Size(cl[i])];
    u:=Filtered([1..Length(cl)],
      x->[Order(Representative(cl[x])),Size(cl[x])]=r and not Representative(cl[x]) in ser[pos]);
    c:=cl{u};
    cl:=cl{Difference([1..Length(cl)],u)};
    if Length(c)>0
      # no need to process classes that will not improve
      and Size(c[1])<bestdeg-offset then
      if Length(c)>1 then
        if Size(c[1])<250 and
          MemoryUsage(Representative(c[1]))*Size(c[1])*Length(c)
          # at most
          <10^7
          # bytes storage
          then
          c:=List(c,claselms);
          Sort(c);
          Info(InfoMorph,4,"Element sets");
          if actbase<>fail then
            of:=Orbits(auo,c,OnSets); # not necc. domain
          else
            of:=OrbitsDomain(auo,c,OnSets);
          fi;
        else
          if actbase<>fail then
            if preproc<>fail then
              of:=List(Orbits(auo,List(c,preproc),action),
                x->List(x,postproc));
            else
              of:=Orbits(auo,c,action); # not necc. domain
            fi;
          else
            if preproc<>fail then
              of:=List(OrbitsDomain(auo,List(c,preproc),action),
                x->List(x,postproc));
            else
              of:=OrbitsDomain(auo,c,action);
            fi;
          fi;
        fi;
      else
        of:=[c];
      fi;
      u:=List(of,x->rec(siz:=Size(x[1])*Length(x),clas:=x,reps:=[Representative(x[1])],
        closure:=NormalClosure(g,SubgroupNC(g,
          Concatenation(List(x,Representative),fix))) ));
      of:=[];
      for j in u do
        if j.siz>1 and j.siz<bestdeg
            and ForAll(o,x->x.siz>j.siz or x.closure<>j.closure)
            and ForAll(of,x->x.siz>j.siz or x.closure<>j.closure) then
          Add(of,j);
        fi;
      od;


      of:=Filtered(of,x->x.siz<bestdeg and x.siz>1);

      Info(InfoMorph,4,"Local ",Length(of)," orbits from ",Length(c),
        " sizes=",Collected(List(of,x->x.siz)));

      # combine with previous
      u:=[];
      for j in o do
        for r in of do
          if j.siz+r.siz<bestdeg
            and not ForAny(j.reps,x-> x in r.closure)
            and not ForAny(r.reps,x-> x in j.closure)
              then
                new:=rec(siz:=j.siz+r.siz,clas:=Concatenation(j.clas,r.clas),
                  reps:=Concatenation(j.reps,r.reps),
                  closure:=ClosureGroup(j.closure,r.closure));
                # don't add if known closure but no better price
                if ForAll(o,x->x.siz>new.siz or x.closure<>new.closure)
                    and ForAll(of,x->x.siz>new.siz or x.closure<>new.closure)
                    and ForAll(u,x->x.siz>new.siz or x.closure<>new.closure)
                     then
                  Add(u,new);
                fi;
          fi;
        od;
      od;
      Append(of,u);

      SortBy(of,x->x.siz);

      u:=First(of,x->Size(x.closure)=Size(g));

      if u<>fail and u.siz<bestdeg then
        if u.siz<baddegree then
          hom:=ActionHomomorphism(au,Union(List(u.clas,claselms)),"surjective");
          finish(hom); return;
        else
          store:=u;
          best:=[function() return Union(List(store.clas,claselms));end,
            OnPoints];
          bestdeg:=u.siz;
          Info(InfoMorph,4,"Improved bestdeg to ",bestdeg);
        fi;
      fi;

      Append(o,of);
      SortBy(o,x->x.siz);
    fi;

    i:=First([1..Length(cl)],x->not Representative(cl[x]) in ser[pos]);
    while i=fail and bestdeg>10*baddegree and pos<Length(ser) do
      u:=First(o,x->Size(ClosureGroup(ser[pos],x.closure))=Size(g));
      if u<>fail then
        offset:=u.siz;
      fi;
      pos:=pos+1;
      i:=First([1..Length(cl)],x->not Representative(cl[x]) in ser[pos]);
      if (bestdeg-offset)/bestdeg<
      # don't bother if all but
        1/4
      # of the points are needed anyhow for the factor. In that case abort
      then  i:=fail;
      fi;

    od;
  od;

  Info(InfoMorph,3,Length(o)," class orbits");

  if Size(SolvableRadical(g))=1 and IsNonabelianSimpleGroup(Socle(g)) then
    Info(InfoMorph,2,"Try ARG");
    img:=AutomorphismRepresentingGroup(g,GeneratorsOfGroup(au));
    # make a hom from auts to perm group
    ran:=Image(img[2],g);
    deg:=NrMovedPoints(ran);
    if deg<bestdeg then
      bestdeg:=deg;
      r:=List(GeneratorsOfGroup(g),i->Image(img[2],i));
      store:=[ran,r];
      hom:=GroupHomomorphismByFunction(au,img[1],
        function(auto)
          if IsInnerAutomorphism(auto) then
            return Image(img[2],ConjugatorOfConjugatorIsomorphism(auto));
          fi;
          return RepresentativeAction(img[1],store[2],
                    List(GeneratorsOfGroup(g),i->Image(img[2],Image(auto,i))),
                    OnTuples);
        end,
        function(perm)
          if perm in store[1] then
            return ConjugatorAutomorphismNC(g,
                      PreImagesRepresentative(img[2],perm));
          fi;
          return GroupHomomorphismByImagesNC(g,g,GeneratorsOfGroup(g),
                    List(store[2],i->PreImagesRepresentative(img[2],i^perm)));
        end);
      if bestdeg<baddegree then
        finish(hom); return;
      fi;
      best:=hom;
    fi;
  fi;

  # try to embed into wreath according to non-conjugator
  if IsPermGroup(g) then
    img:=AutomorphismWreathEmbedding(au,g);
  else
    img:=fail;
  fi;
  if img<>fail and NrMovedPoints(img[4])<bestdeg then
    Info(InfoMorph,2,"AWE succeeds");
    # make a hom from auts to perm group
    ran:=img[4];
    deg:=NrMovedPoints(ran);
    bestdeg:=deg;
    r:=List(GeneratorsOfGroup(g),i->Image(img[3],i));
    store:=[img[1],img[2],img[3],r,ran];
    hom:=GroupHomomorphismByFunction(au,img[1],
      function(auto)
        if IsConjugatorAutomorphism(auto) and
          ConjugatorOfConjugatorIsomorphism(auto) in Source(store[2]) then
          return Image(store[2],ConjugatorOfConjugatorIsomorphism(auto));
        fi;
        return RepresentativeAction(store[1],store[4],
                  List(GeneratorsOfGroup(g),i->Image(img[3],Image(auto,i))),OnTuples);
      end,
      function(perm)
        if perm in store[5] then
          return ConjugatorAutomorphismNC(g,
                    PreImagesRepresentative(store[2],perm));
        fi;
        return GroupHomomorphismByImagesNC(g,g,GeneratorsOfGroup(g),
                  List(store[4],i->PreImagesRepresentative(store[3],i^perm)));
      end);
    if bestdeg<baddegree then
      finish(hom); return;
    fi;
    best:=hom;
  fi;

  # if no centre and permgroup, try to blow up perm rep
  if Size(Centre(g))=1 and IsPermGroup(g) then

    orbs:=List(Orbits(g,MovedPoints(g)),Set);
    SortBy(orbs,Length);

    # we act on lists of [orbits,group]. This way comparison becomes
    # cheap, as typically the orbits suffice to identify

    if IsPermGroup(g) then
      action:=function(obj,hom)
              local img;
                img:=Image(hom,obj[2]);
                SetSize(img,Size(obj[2]));
                return Immutable([OrbitsMovedPoints(img),img]);
              end;
    else
      action:=function(sub,autom)
        local img;
        img:=Image(autom,sub);
        SetSize(img,Size(sub));
        return img;
      end;
    fi;

    for o in orbs do
      r:=Stabilizer(g,o,OnTuples);
      if hom=fail and #have not yet found any
        not ForAll(GeneratorsOfGroup(au),x->Image(x,r)=r or
        Difference(MovedPoints(g),MovedPoints(Image(x,r))) in orbs) then
        u:=Size(r);
        r:=Subgroup(Parent(r),SmallGeneratingSet(r));
        SetSize(r,u);
        if IsPermGroup(g) then
          r:=Orbit(au,Immutable([OrbitsMovedPoints(r),r]),action);
          r:=List(r,x->x[2]);
        else
          r:=Orbit(au,r,action);
        fi;

        u:=Stabilizer(g,o[1]);
        if Size(Intersection(r))=1 and
          # this orbit and automorphism images represent all of g, so can
          # be used to represent automorphisms
          u=Normalizer(g,u) then
          # point stabilizer is self-normalizing, so action on
          # cosets=action by conjugation
          u:=Group(SmallGeneratingSet(u));
          if IsPermGroup(g) then
            u:=Orbit(au,Immutable([OrbitsMovedPoints(u),u]),action);
          else
            u:=Orbit(au,u,action);
          fi;
          deg:=Length(u);
          if deg<bestdeg then
            bestdeg:=deg;
            if bestdeg<baddegree then
              hom:=ActionHomomorphism(au,u,action,"surjective");
              finish(hom); return;
            fi;
            store:=u;
            best:=[function() return store;end,action];
          fi;

        fi;
      fi;
    od;
  fi;

#if bestdeg>10*baddegree then Error("bad degree");fi;

  Info(InfoMorph,1,"Go back to best rep so far ",bestdeg);
  # go back to what we had before
  if IsList(best) then
    hom:=ActionHomomorphism(au,best[1](),best[2],"surjective");
  else
    hom:=best;
  fi;
  finish(hom);return;

#  # fall back on element sets
#  Info(InfoMorph,1,"General Case");
#
#
#  # general case: compute small domain
#  gens:=[];
#  dom:=[];
#  u:=TrivialSubgroup(g);
#  subs:=[];
#  orbs:=[];
#  while Size(u)<Size(g) do
#    # find a reasonable element
#    cnt:=0;
#    br:=false;
#    bv:=0;
#
#    # use classes from before --
#    #if HasConjugacyClasses(g) then
#    for r in cl do
#      if IsPrimePowerInt(Order(Representative(r))) and
#          not Representative(r) in  u then
#        v:=ClosureGroup(u,Representative(r));
#        if allinner then
#          val:=Size(Centralizer(r))*Size(NormalClosure(g,v));
#        else
#          val:=Size(Centralizer(r))*Size(v);
#        fi;
#        if val>bv then
#          br:=Representative(r);
#          bv:=val;
#        fi;
#      fi;
#    od;
#
##    else
##      if actbase=fail then
###        actbase:=[g];
##      fi;
##      repeat
##        cnt:=cnt+1;
##        repeat
##          r:=Random(Random(actbase));
##        until not r in u;
##        # force small prime power order
##        if not IsPrimePowerInt(Order(r)) then
##          v:=List(Collected(Factors(Order(r))),x->r^(x[1]^x[2]));
##          v:=Filtered(v,x->not x in u);;
##          v:=List(v,x->x^First(Reversed(DivisorsInt(Order(x))),
##            y->not x^y in u)); # small order power not in u
##          SortBy(v,Order);
##          r:=First(v,x->not x in u); # if all are in u, r would be as well
##        fi;
##
##        v:=ClosureGroup(u,r);
##        #if allinner then
##        #  val:=Size(Centralizer(g,r))*Size(NormalClosure(g,v));
##        #else
##          val:=Size(Centralizer(g,r));
##        #fi;
##        if val>bv then
##          br:=r;
##          bv:=val;
##        fi;
##      until bv>Size(g)/2^Int(cnt/50);
##    fi;
#
#    r:=br;
#    Info(InfoMorph,4,"Element class length ",Index(g,Centralizer(g,r))," after ",cnt);
#
#    #if Index(g,Centralizer(g,r))>2*10^4 then Error("big degree debug");fi;
#
#    if allinner then
#      u:=NormalClosure(g,ClosureGroup(u,r));
#    else
#      u:=ClosureGroup(u,r);
#    fi;
#
#    #calculate orbit and closure
#    o:=Orbit(au,r);
#    v:=TrivialSubgroup(g);
#    i:=1;
#    while i<=Length(o) do
#      if not o[i] in v then
#        if allinner then
#          v:=NormalClosure(g,ClosureGroup(v,o[i]));
#        else
#          v:=ClosureGroup(v,o[i]);
#        fi;
#        if Size(v)=Size(g) then
#          i:=Length(o);
#        fi;
#      fi;
#      i:=i+1;
#    od;
#    u:=ClosureGroup(u,v);
#
#    i:=1;
#    while Length(o)>0 and i<=Length(subs) do
#      if IsSubset(subs[i],v) then
#        o:=[];
#      elif IsSubset(v,subs[i]) then
#        subs[i]:=v;
#        orbs[i]:=o;
#        gens[i]:=r;
#        o:=[];
#      fi;
#      i:=i+1;
#    od;
#    if Length(o)>0 then
#      Add(subs,v);
#      Add(orbs,o);
#      Add(gens,r);
#    fi;
#  od;
#
#  # now find the smallest subset of domains
#  comb:=Filtered(Combinations([1..Length(subs)]),i->Length(i)>0);
#  bv:=infinity;
#  for i in comb do
#    val:=Sum(List(orbs{i},Length));
#    if val<bv then
#      v:=subs[i[1]];
#      for r in [2..Length(i)] do
#        v:=ClosureGroup(v,subs[i[r]]);
#      od;
#      if Size(v)=Size(g) then
#        best:=i;
#        bv:=val;
#      fi;
#    fi;
#  od;
#  gens:=gens{best};
#  dom:=Union(orbs{best});
#  Unbind(orbs);
#
#  if Length(dom)>bestdeg then
#    # go back to what we had before
#    if IsList(best) then
#      hom:=ActionHomomorphism(au,best[1](),best[2],"surjective");
#    else
#      hom:=best;
#    fi;
#    finish(hom);return;
#  fi;
#
#  u:=SubgroupNC(g,gens);
#  while Size(u)<Size(g) do
#    repeat
#      r:=Random(dom);
#    until not r in u;
#    Add(gens,r);
#    u:=ClosureSubgroupNC(u,r);
#  od;
#  Info(InfoMorph,1,"Found generating set of ",Length(gens)," elements",
#        List(gens,Order));
#  hom:=NiceMonomorphismAutomGroup(au,dom,gens);
#
#  finish(hom);
#
end);

#############################################################################
##
#F  NiceMonomorphismAutomGroup
##
InstallGlobalFunction(NiceMonomorphismAutomGroup,
function(aut,elms,elmsgens)
local xset,fam,hom;
  One(aut); # to avoid infinite recursion once the niceo is set

  elmsgens:=Filtered(elmsgens,i->i in elms); # safety feature
  #if Size(Group(elmsgens))<>Size(Source(One(aut))) then Error("holler1"); fi;
  xset:=ExternalSet(aut,elms);
  SetBaseOfGroup(xset,elmsgens);
  fam := GeneralMappingsFamily( ElementsFamily( FamilyObj( aut ) ),
                                PermutationsFamily );
  hom := rec(  );
  hom:=Objectify(NewType(fam,
                IsActionHomomorphismAutomGroup and IsSurjective ),hom);
  SetIsInjective(hom,true);
  SetUnderlyingExternalSet( hom, xset );
  hom!.basepos:=List(elmsgens,i->Position(elms,i));
  SetRange( hom, Image( hom ) );
  Setter(SurjectiveActionHomomorphismAttr)(xset,hom);
  Setter(IsomorphismPermGroup)(aut,ActionHomomorphism(xset,"surjective"));
  hom:=ActionHomomorphism(xset,"surjective");
  SetFilterObj(hom,IsNiceMonomorphism);
  return hom;

end);

#############################################################################
##
#M  PreImagesRepresentative   for OpHomAutomGrp
##
InstallMethod(PreImagesRepresentative,"AutomGroup Niceomorphism",
  FamRangeEqFamElm,[IsActionHomomorphismAutomGroup,IsPerm],0,
function(hom,elm)
local xset,g,imgs;
  xset:= UnderlyingExternalSet( hom );
  g:=Source(One(ActingDomain(xset)));
  imgs:=OnTuples(hom!.basepos,elm);
  imgs:=Enumerator(xset){imgs};
  #if g<>Group(BaseOfGroup(xset)) then Error("holler"); fi;
  elm:=GroupHomomorphismByImagesNC(g,g,BaseOfGroup(xset),imgs);
  SetIsBijective(elm,true);
  return elm;
end);


#############################################################################
##
#F  MorFroWords(<gens>) . . . . . . create some pseudo-random words in <gens>
##                                                featuring the MeatAxe's FRO
InstallGlobalFunction(MorFroWords,function(gens)
local list,a,b,ab,i;
  list:=[];
  ab:=gens[1];
  for i in [2..Length(gens)] do
    a:=ab;
    b:=gens[i];
    ab:=a*b;
    list:=Concatenation(list,
         [ab,ab^2*b,ab^3*b,ab^4*b,ab^2*b*ab^3*b,ab^5*b,ab^2*b*ab^3*b*ab*b,
         ab*(ab*b)^2*ab^3*b,a*b^4*a,ab*a^3*b]);
  od;
  return list;
end);


#############################################################################
##
#F  MorRatClasses(<G>) . . . . . . . . . . . local rationalization of classes
##
InstallGlobalFunction(MorRatClasses,function(GR)
local r,c,u,j,i;
  Info(InfoMorph,2,"RationalizeClasses");
  r:=[];
  for c in RationalClasses(GR) do
    u:=Subgroup(GR,[Representative(c)]);
    j:=DecomposedRationalClass(c);
    Add(r,rec(representative:=u,
                class:=j[1],
                classes:=j,
                size:=Size(c)));
  od;

  for i in r do
    i.size:=Sum(i.classes,Size);
  od;
  return r;
end);

#############################################################################
##
#F  MorMaxFusClasses(<l>) . .  maximal possible morphism fusion of classlists
##
InstallGlobalFunction(MorMaxFusClasses,function(r)
local i,j,flag,cl;
  # cl is the maximal fusion among the rational classes.
  cl:=[];
  for i in r do
    j:=0;
    flag:=true;
    while flag and j<Length(cl) do
      j:=j+1;
      flag:=not(Size(i.class)=Size(cl[j][1].class) and
                  i.size=cl[j][1].size and
                  Size(i.representative)=Size(cl[j][1].representative));
    od;
    if flag then
      Add(cl,[i]);
    else
      Add(cl[j],i);
    fi;
  od;

  # sort classes by size
  SortBy(cl,a->Sum(a,i->i.size));
  return cl;
end);

#############################################################################
##
#F  SomeVerbalSubgroups
##
## correspond simultaneously some verbal subgroups in g and h
BindGlobal("SomeVerbalSubgroups",function(g,h)
local l,m,i,j,cg,ch,pg;
  l:=[g];
  m:=[h];
  i:=1;
  while i<=Length(l) do
    for j in [1..i] do
      cg:=CommutatorSubgroup(l[i],l[j]);
      ch:=CommutatorSubgroup(m[i],m[j]);
      pg:=Position(l,cg);
      if pg=fail then
        Add(l,cg);
        Add(m,ch);
      else
        while m[pg]<>ch do
          pg:=Position(l,cg,pg+1);
          if pg=fail then
            Add(l,cg);
            Add(m,ch);
            pg:=Length(m);
          fi;
        od;
      fi;
    od;
    i:=i+1;
  od;
  return [l,m];
end);

#############################################################################
##
#F  MorClassLoop(<range>,<classes>,<params>,<action>)  loop over classes list
##     to find generating sets or Iso/Automorphisms up to inner automorphisms
##
##  classes is a list of records like the ones returned from
##  MorMaxFusClasses.
##
##  params is a record containing optional components:
##  gens  generators that are to be mapped
##  from  preimage group (that contains gens)
##  to    image group (as it might be smaller than 'range')
##  free  free generators
##  rels  some relations that hold in from, given as list [word,order]
##  dom   a set of elements on which automorphisms act faithful
##  aut   Subgroup of already known automorphisms
##  condition function that must return `true' on the homomorphism.
##  setrun  If set to true approximate a run through sets by having the
##          class indices increasing.
##
##  action is a number whose bit-representation indicates the action to be
##  taken:
##  1     homomorphism
##  2     injective
##  4     surjective
##  8     find all (in contrast to one)
##
BindGlobal( "MorClassOrbs", function(G,C,R,D)
local i,cl,cls,rep,x,xp,p,b,g;
  i:=Index(G,C);
  if i>20000 or i<Size(D) then
    return List(DoubleCosetRepsAndSizes(G,C,D),j->j[1]);
  else
    if not IsBound(C!.conjclass) then
      cl:=[R];
      cls:=[R];
      rep:=[One(G)];
      i:=1;
      while i<=Length(cl) do
        for g in GeneratorsOfGroup(G) do
          x:=cl[i]^g;
          if not x in cls then
            Add(cl,x);
            AddSet(cls,x);
            Add(rep,rep[i]*g);
          fi;
        od;
        i:=i+1;
      od;
      SortParallel(cl,rep);
      C!.conjclass:=cl;
      C!.conjreps:=rep;
    fi;
    cl:=C!.conjclass;
    rep:=[];
    b:=BlistList([1..Length(cl)],[]);
    p:=1;
    repeat
      while p<=Length(cl) and b[p] do
        p:=p+1;
      od;
      if p<=Length(cl) then
        b[p]:=true;
        Add(rep,p);
        cls:=[cl[p]];
        for i in cls do
          for g in GeneratorsOfGroup(D) do
            x:=i^g;
            xp:=PositionSorted(cl,x);
            if not b[xp] then
              Add(cls,x);
              b[xp]:=true;
            fi;
          od;
        od;
      fi;
      p:=p+1;
    until p>Length(cl);

    return C!.conjreps{rep};
  fi;
end );

InstallGlobalFunction(MorClassLoop,function(range,clali,params,action)
local id,result,rig,dom,tall,tsur,tinj,thom,gens,free,rels,len,ind,cla,m,
      mp,cen,i,j,imgs,ok,size,l,hom,cenis,reps,repspows,sortrels,genums,wert,p,
      e,offset,pows,TestRels,pop,mfw,derhom,skip,cond,outerorder,setrun,
      finsrc;

  finsrc:=IsBound(params.from) and HasIsFinite(params.from)
          and IsFinite(params.from);
  len:=Length(clali);
  if ForAny(clali,i->Length(i)=0) then
    return []; # trivial case: no images for generator
  fi;

  id:=One(range);
  if IsBound(params.aut) then
    result:=params.aut;
    rig:=true;
    if IsBound(params.dom) then
      dom:=params.dom;
    else
      dom:=false;
    fi;
  else
    result:=[];
    rig:=false;
  fi;

  if IsBound(params.outerorder) then
    outerorder:=params.outerorder;
  else
    outerorder:=false;
  fi;

  if IsBound(params.setrun) then
    setrun:=params.setrun;
  else
    setrun:=false;
  fi;


  # extra condition?
  if IsBound(params.condition) then
    cond:=params.condition;
  else
    cond:=fail;
  fi;

  tall:=action>7; # try all
  if tall then
    action:=action-8;
  fi;
  derhom:=fail;
  tsur:=action>3; # test surjective
  if tsur then
    size:=Size(params.to);
    action:=action-4;
    if Index(range,DerivedSubgroup(range))>1 then
      derhom:=NaturalHomomorphismByNormalSubgroup(range,DerivedSubgroup(range));
    fi;
  fi;
  tinj:=action>1; # test injective
  if tinj then
    action:=action-2;
  fi;
  thom:=action>0; # test homomorphism

  if IsBound(params.gens) then
    gens:=params.gens;
  fi;

  if IsBound(params.rels) then
    free:=params.free;
    rels:=params.rels;
    if Length(rels)=0 then
      rels:=false;
    fi;
  elif thom then
    free:=GeneratorsOfGroup(FreeGroup(Length(gens)));
    mfw:=MorFroWords(free);
    # get some more
    if finsrc and Product(List(gens,Order))<2000 then
      for i in Cartesian(List(gens,i->[1..Order(i)])) do
        Add(mfw,Product(List([1..Length(gens)],z->free[z]^i[z])));
      od;
    fi;
    rels:=[];
    for i in mfw do
      p:=Order(MappedWord(i,free,gens));
      if p<>infinity then
        Add(rels,[i,p]);
      fi;
    od;
    if Length(rels)=0 then
      rels:=false;
    fi;
  else
    rels:=false;
  fi;

  pows:=[];
  if rels<>false then
    # sort the relators according to the generators they contain
    genums:=List(free,i->GeneratorSyllable(i,1));
    genums:=List([1..Length(genums)],i->Position(genums,i));
    sortrels:=List([1..len],i->[]);
    pows:=List([1..len],i->[]);
    for i in rels do
      l:=len;
      wert:=0;
      m:=[];
      for j in [1..NrSyllables(i[1])] do
        p:=genums[GeneratorSyllable(i[1],j)];
        e:=ExponentSyllable(i[1],j);
        Append(m,[p,e]); # modified extrep
        AddSet(pows[p],e);
        if p<len then
          wert:=wert+2; # conjugation: 2 extra images
          l:=Minimum(l,p);
        fi;
        wert:=wert+AbsInt(e);
      od;
      Add(sortrels[l],[m,i[2],i[2]*wert,[1,3..Length(m)-1],i[1]]);
    od;
    # now sort by the length of the relators
    for i in [1..len] do
      SortBy(sortrels[i],x->x[3]);
    od;
    if Length(pows)>0 then
      offset:=1-Minimum(List(Filtered(pows,i->Length(i)>0),
                            i->i[1])); # smallest occurring index
    fi;
    # test the relators at level tlev and set imgs
    TestRels:=function(tlev)
    local rel,k,j,p,start,gn,ex;

      if Length(sortrels[tlev])=0 then
        imgs:=List([tlev..len-1],i->reps[i]^(m[i][mp[i]]));
        imgs[Length(imgs)+1]:=reps[len];
        return true;
      fi;

      if IsPermGroup(range) then
        # test by tracing points
        for rel in sortrels[tlev] do
          start:=1;
          p:=start;
          k:=0;
          repeat
            for j in rel[4] do
              gn:=rel[1][j];
              ex:=rel[1][j+1];
              if gn=len then
                p:=p^repspows[gn][ex+offset];
              else
                p:=p/m[gn][mp[gn]];
                p:=p^repspows[gn][ex+offset];
                p:=p^m[gn][mp[gn]];
              fi;
            od;
            k:=k+1;
          # until we have the power or we detected a smaller potential order.
          until k>=rel[2] or (p=start and IsInt(rel[2]/k));
          if p<>start then
            return false;
          fi;
        od;
      fi;

      imgs:=List([tlev..len-1],i->reps[i]^(m[i][mp[i]]));
      imgs[Length(imgs)+1]:=reps[len];

      if tinj then
        return ForAll(sortrels[tlev],i->i[2]=Order(MappedWord(i[5],
                                      free{[tlev..len]}, imgs)));
      else
        return ForAll(sortrels[tlev],
                      i->IsInt(i[2]/Order(MappedWord(i[5],
                                          free{[tlev..len]}, imgs))));
      fi;

    end;
  else
    TestRels:=x->true; # to satisfy the code below.
  fi;

  pop:=false; # just to initialize

  # backtrack over all classes in clali
  l:=ListWithIdenticalEntries(len,1);
  ind:=len;
  while ind>0 do
    ind:=len;
    Info(InfoMorph,3,"step ",l);
    # test class combination indicated by l:
    cla:=List([1..len],i->clali[i][l[i]]);
    reps:=List(cla,Representative);
    skip:=false;
    if derhom<>fail then
      if not Size(Group(List(reps,i->Image(derhom,i))))=Size(Image(derhom)) then
#T The group `Image( derhom )' is abelian but initially does not know this;
#T shouldn't this be set?
#T Then computing the size on the l.h.s. may be sped up using `SubgroupNC'
#T w.r.t. the (abelian) group.
        skip:=true;
        Info(InfoMorph,3,"skipped");
      fi;
    fi;

    if pows<>fail and not skip then
      if rels<>false and IsPermGroup(range) then
        # and precompute the powers
        repspows:=List([1..len],i->[]);
        for i in [1..len] do
          for j in pows[i] do
            repspows[i][j+offset]:=reps[i]^j;
          od;
        od;
      fi;

      #cenis:=List(cla,i->Intersection(range,Centralizer(i)));
      # make sure we get new groups (we potentially add entries)
      cenis:=[];
      for i in cla do
        cen:=Intersection(range,Centralizer(i));
        if IsIdenticalObj(cen,Centralizer(i)) then
          m:=Size(cen);
          cen:=SubgroupNC(range,GeneratorsOfGroup(cen));
          SetSize(cen,m);
        fi;
        Add(cenis,cen);
      od;

      # test, whether a gen.sys. can be taken from the classes in <cla>
      # candidates.  This is another backtrack
      m:=[];
      m[len]:=[id];
      # positions
      mp:=[];
      mp[len]:=1;
      mp[len+1]:=-1;
      # centralizers
      cen:=[];
      cen[len]:=cenis[len];
      cen[len+1]:=range; # just for the recursion
      i:=len-1;

      # set up the lists
      while i>0 do
        #m[i]:=List(DoubleCosetRepsAndSizes(range,cenis[i],cen[i+1]),j->j[1]);
        m[i]:=MorClassOrbs(range,cenis[i],reps[i],cen[i+1]);
        mp[i]:=1;

        pop:=true;
        while pop and i<=len do
          pop:=false;
          while mp[i]<=Length(m[i]) and TestRels(i)=false do
            mp[i]:=mp[i]+1; #increment because of relations
            Info(InfoMorph,4,"early break ",i);
          od;
          if i<=len and mp[i]>Length(m[i]) then
            Info(InfoMorph,3,"early pop");
            pop:=true;
            i:=i+1;
            if i<=len then
              mp[i]:=mp[i]+1; #increment because of pop
            fi;
          fi;
        od;

        if pop then
          i:=-99; # to drop out of outer loop
        elif i>1 then
          cen[i]:=Centralizer(cen[i+1],reps[i]^(m[i][mp[i]]));
        fi;
        i:=i-1;
      od;

      if pop then
        Info(InfoMorph,3,"allpop");
        i:=len+2; # to avoid the following `while' loop
      else
        i:=1;
        Info(InfoMorph,3,"loop");
      fi;

      while i<len do
        if rels=false or TestRels(1) then
          if rels=false then
            # otherwise the images are set by `TestRels' as a side effect.
            imgs:=List([1..len-1],i->reps[i]^(m[i][mp[i]]));
            imgs[len]:=reps[len];
          fi;
          Info(InfoMorph,4,"orders: ",List(imgs,Order));

          # computing the size can be nasty. Thus try given relations first.
          ok:=true;

          if rels<>false then
            if tinj then
              ok:=ForAll(rels,i->i[2]=Order(MappedWord(i[1],free,imgs)));
            else
              ok:=ForAll(rels,i->IsInt(i[2]/Order(MappedWord(i[1],free,imgs))));
            fi;
          fi;

          # check surjectivity
          if tsur and ok then
            ok:= Size( SubgroupNC( range, imgs ) ) = size;
          fi;

          if ok and thom then
            Info(InfoMorph,3,"testing");
            imgs:=GroupGeneralMappingByImagesNC(params.from,range,gens,imgs);
            SetIsTotal(imgs,true);
            if tsur then
              SetIsSurjective(imgs,true);
            fi;
            ok:=IsSingleValued(imgs);
            if ok and tinj then
              ok:=IsInjective(imgs);
            fi;
          fi;

          if ok and cond<>fail then
            ok:=cond(imgs);
          fi;

          if ok then
            Info(InfoMorph,2,"found");
            # do we want one or all?
            if tall then
              if rig then
                if not imgs in result then
                  result:= GroupByGenerators( Concatenation(
                              GeneratorsOfGroup( result ), [ imgs ] ),
                              One( result ) );
                  # note its niceo
                  hom:=NiceMonomorphismAutomGroup(result,dom,gens);
                  SetNiceMonomorphism(result,hom);
                  SetIsHandledByNiceMonomorphism(result,true);

                  Size(result);
                  Info(InfoMorph,2,"new ",Size(result));
                fi;
              else
                Add(result,imgs);
                # can we deduce we got all?
                if outerorder<>false and Lcm(Concatenation([1],List(
                  Filtered(result,x->not IsInnerAutomorphism(x)),
                  x->First([2..outerorder],y->IsInnerAutomorphism(x^y)))))
                    =outerorder
                  then
                    Info(InfoMorph,1,"early break");
                    return result;
                fi;
              fi;
            else
              return imgs;
            fi;
          fi;
        fi;

        mp[i]:=mp[i]+1;
        while i<=len and mp[i]>Length(m[i]) do
          mp[i]:=1;
          i:=i+1;
          if i<=len then
            mp[i]:=mp[i]+1;
          fi;
        od;

        while i>1 and i<=len do
          while i<=len and TestRels(i)=false do
            Info(InfoMorph,4,"intermediate break ",i);
            mp[i]:=mp[i]+1;
            while i<=len and mp[i]>Length(m[i]) do
              Info(InfoMorph,3,"intermediate pop ",i);
              i:=i+1;
              if i<=len then
                mp[i]:=mp[i]+1;
              fi;
            od;
          od;

          if i<=len then # i>len means we completely popped. This will then
                        # also pop us out of both `while' loops.
            cen[i]:=Centralizer(cen[i+1],reps[i]^(m[i][mp[i]]));
            i:=i-1;
            #m[i]:=List(DoubleCosetRepsAndSizes(range,cenis[i],cen[i+1]),j->j[1]);
            m[i]:=MorClassOrbs(range,cenis[i],reps[i],cen[i+1]);
            mp[i]:=1;

          else
            Info(InfoMorph,3,"allpop2");
          fi;
        od;

      od;
    fi;

    # 'free for increment'
    l[ind]:=l[ind]+1;
    while ind>0 and l[ind]>Length(clali[ind]) do
      if setrun and ind>1 then
        # if we are running through sets, approximate by having the
        # l-indices increasing
        l[ind]:=Minimum(l[ind-1]+1,Length(clali[ind]));
      else
        l[ind]:=1;
      fi;
      ind:=ind-1;
      if ind>0 then
        l[ind]:=l[ind]+1;
      fi;
    od;
  od;

  return result;
end);


#############################################################################
##
#F  MorFindGeneratingSystem(<G>,<cl>) . .  find generating system with as few
##                      as possible generators from the first classes in <cl>
##
InstallGlobalFunction(MorFindGeneratingSystem,function(arg)
local G,cl,lcl,len,comb,combc,com,a,alltwo;
  G:=arg[1];
  cl:=arg[2];
  Info(InfoMorph,1,"FindGenerators");
  # throw out the 1-Class
  cl:=Filtered(cl,i->Length(i)>1 or Size(i[1].representative)>1);
  alltwo:=PrimeDivisors(Size(G))=[2];

  #create just a list of ordinary classes.
  lcl:=List(cl,i->Concatenation(List(i,j->j.classes)));
  len:=1;
  len:=Maximum(1,AbelianRank(G)-1);
  while true do
    len:=len+1;
    Info(InfoMorph,2,"Trying length ",len);
    # now search for <len>-generating systems
    comb:=UnorderedTuples([1..Length(lcl)],len);
    combc:=List(comb,i->List(i,j->lcl[j]));

    # test all <comb>inations
    com:=0;
    while com<Length(comb) do
      com:=com+1;
      # don't try only order 2 generators unless it's a 2-group
      if Set(Flat(combc[com]),i->Order(Representative(i)))<>[2] or
        alltwo then
        a:=MorClassLoop(G,combc[com],rec(to:=G),4);
        if Length(a)>0 then
          return a;
        fi;
      fi;
    od;
  od;
end);

#############################################################################
##
#F  Morphium(<G>,<H>,<DoAuto>) . . . . . . . .Find isomorphisms between G and H
##       modulo inner automorphisms. DoAuto indicates whether all
##       automorphism are to be found
##       This function thus does the main combinatoric work for creating
##       Iso- and Automorphisms.
##       It needs, that both groups are not cyclic.
##
InstallGlobalFunction(Morphium,function(G,H,DoAuto)
local combi,Gr,Gcl,Ggc,Hr,Hcl,bg,bpri,x,dat,
      gens,i,c,hom,elms,price,result,inns,bcl,vsu;

  if IsSolvableGroup(G) and CanEasilyComputePcgs(G) then
    gens:=MinimalGeneratingSet(G);
  else
    gens:=SmallGeneratingSet(G);
  fi;
  Gr:=MorRatClasses(G);
  Gcl:=MorMaxFusClasses(Gr);

  Ggc:=List(gens,i->First(Gcl,j->ForAny(j,j->ForAny(j.classes,k->i in k))));
  combi:=List(Ggc,i->Concatenation(List(i,i->i.classes)));
  price:=Product(combi,i->Sum(i,Size));
  Info(InfoMorph,1,"generating system ",Sum(Flat(combi),Size),
       " of price:",price,"");

  if ((not HasMinimalGeneratingSet(G) and price/Size(G)>10000)
     or Sum(Flat(combi),Size)>Size(G)/10 or IsSolvableGroup(G))
     and ValueOption("nogensyssearch")<>true then
    if IsSolvableGroup(G) then
      gens:=IsomorphismPcGroup(G);
      gens:=List(MinimalGeneratingSet(Image(gens)),
                 i->PreImagesRepresentative(gens,i));
      Ggc:=List(gens,i->First(Gcl,j->ForAny(j,j->ForAny(j.classes,k->i in k))));
      combi:=List(Ggc,i->Concatenation(List(i,i->i.classes)));
      bcl:=ShallowCopy(combi);
      SortBy(bcl,a->Sum(a,Size));
      bg:=gens;
      bpri:=Product(combi,i->Sum(i,Size));
      for i in [1..7*Length(gens)-12] do
        repeat
          for c in [1..Length(gens)] do
            if Random(1,3)<2 then
              gens[c]:=Random(G);
            else
              x:=bcl[Random(Filtered([1,1,1,1,2,2,2,3,3,4],k->k<=Length(bcl)))];
              gens[c]:=Random(Random(x));
            fi;
          od;
        until Index(G,SubgroupNC(G,gens))=1;
        Ggc:=List(gens,i->First(Gcl,
                  j->ForAny(j,j->ForAny(j.classes,k->i in k))));
        combi:=List(Ggc,i->Concatenation(List(i,i->i.classes)));
        Append(bcl,combi);
        SortBy(bcl,a->Sum(a,Size));
        price:=Product(combi,i->Sum(i,Size));
        Info(InfoMorph,3,"generating system of price:",price,"");
        if price<bpri then
          bpri:=price;
          bg:=gens;
        fi;
      od;

      gens:=bg;

    else
      gens:=MorFindGeneratingSystem(G,Gcl);
    fi;

    Ggc:=List(gens,i->First(Gcl,j->ForAny(j,j->ForAny(j.classes,k->i in k))));
    combi:=List(Ggc,i->Concatenation(List(i,i->i.classes)));
    price:=Product(combi,i->Sum(i,Size));
    Info(InfoMorph,1,"generating system of price:",price,"");
  fi;

  if not DoAuto then
    Hr:=MorRatClasses(H);
    Hcl:=MorMaxFusClasses(Hr);
  fi;

  vsu:=SomeVerbalSubgroups(G,H);
  if List(vsu[1],Size)<>List(vsu[2],Size) then
    # cannot be candidates
    return [];
  fi;

  # now test, whether it is worth, to compute a finer congruence
  # then ALSO COMPUTE NEW GEN SYST!
  # [...]

  if not DoAuto then
    combi:=[];
    for i in Ggc do
      c:=Filtered(Hcl,
           j->Set(j,k->k.size)=Set(i,k->k.size)
                and Length(j[1].classes)=Length(i[1].classes)
                and Size(j[1].class)=Size(i[1].class)
                and Size(j[1].representative)=Size(i[1].representative)
      # This test assumes maximal fusion among the rat.classes. If better
      # congruences are used, they MUST be checked here also!
        );
      if Length(c)<>1 then
        # Both groups cannot be isomorphic, since they lead to different
        # congruences!
        Info(InfoMorph,2,"different congruences");
        return fail;
      else
        Add(combi,c[1]);
      fi;
    od;
    combi:=List(combi,i->Concatenation(List(i,i->i.classes)));
  fi;

  # filter by verbal subgroups
  for i in [1..Length(gens)] do
    c:=Filtered([1..Length(vsu[1])],j->gens[i] in vsu[1][j]);
    c:=Filtered(combi[i],k->
         c=Filtered([1..Length(vsu[2])],j->Representative(k) in vsu[2][j]));
    if Length(c)<Length(combi[i]) then
      Info(InfoMorph,1,"images improved by verbal subgroup:",
      Sum(combi[i],Size)," -> ",Sum(c,Size));
      combi[i]:=c;
    fi;
  od;

  # combi contains the classes, from which the
  # generators are taken.

  #free:=GeneratorsOfGroup(FreeGroup(Length(gens)));
  #rels:=MorFroWords(free);
  #rels:=List(rels,i->[i,Order(MappedWord(i,free,gens))]);
  #result:=rec(gens:=gens,from:=G,to:=H,free:=free,rels:=rels);
  result:=rec(gens:=gens,from:=G,to:=H);

  if DoAuto then

    inns:=List(GeneratorsOfGroup(G),i->InnerAutomorphism(G,i));
    if Sum(Flat(combi),Size)<=MORPHEUSELMS then
      elms:=[];
      for i in Flat(combi) do
        if not ForAny(elms,j->Representative(i)=Representative(j)) then
          # avoid duplicate classes
          Add(elms,i);
        fi;
      od;
      elms:=Union(List(elms,AsList));
      Info(InfoMorph,1,"permrep on elements: ",Length(elms));

      Assert(2,ForAll(GeneratorsOfGroup(G),i->ForAll(elms,j->j^i in elms)));
      result.dom:=elms;
      inns:= GroupByGenerators( inns, IdentityMapping( G ) );

      hom:=NiceMonomorphismAutomGroup(inns,elms,gens);
      SetNiceMonomorphism(inns,hom);
      SetIsHandledByNiceMonomorphism(inns,true);

      result.aut:=inns;
    else
      elms:=false;
    fi;

    # catch case of simple groups to get outer automorphism orders
    # automorphism suffices.
    if IsNonabelianSimpleGroup(H) then
      dat:=DataAboutSimpleGroup(H);
      if IsBound(dat.fullAutGroup) then
        if dat.fullAutGroup[1]=1 then
          # all automs are inner.
          result:=rec(aut:=result.aut);
        else
          result.outerorder:=dat.fullAutGroup[1];
          result:=rec(aut:=MorClassLoop(H,combi,result,15));
        fi;
      else
        result:=rec(aut:=MorClassLoop(H,combi,result,15));
      fi;
    else
      result:=rec(aut:=MorClassLoop(H,combi,result,15));
    fi;

    if elms<>false then
      result.elms:=elms;
      result.elmsgens:=Filtered(gens,i->i<>One(G));
      inns:=SubgroupNC(result.aut,GeneratorsOfGroup(inns));
    fi;
    result.inner:=inns;
  else
    result:=MorClassLoop(H,combi,result,7);
  fi;

  return result;

end);

#############################################################################
##
#F  AutomorphismGroupAbelianGroup(<G>)
##
InstallGlobalFunction(AutomorphismGroupAbelianGroup,function(G)
local i,j,k,l,m,o,nl,nj,max,r,e,au,p,gens,offs;

  # trivial case
  if Size(G)=1 then
    au:= GroupByGenerators( [], IdentityMapping( G ) );
    i:=NiceMonomorphismAutomGroup(au,[One(G)],[One(G)]);
    SetNiceMonomorphism(au,i);
    SetIsHandledByNiceMonomorphism(au,true);
    SetIsAutomorphismGroup( au, true );
    SetIsFinite(au,true);
    return au;
  fi;

  # get standard generating system
  gens:=IndependentGeneratorsOfAbelianGroup(G);

  au:=[];
  # run by primes
  p:=PrimeDivisors(Size(G));
  for i in p do
    l:=Filtered(gens,j->IsInt(Order(j)/i));
    nl:=Filtered(gens,i->not i in l);

    #sort by exponents
    o:=List(l,j->LogInt(Order(j),i));
    e:=[];
    for j in Set(o) do
      Add(e,[j,l{Filtered([1..Length(o)],k->o[k]=j)}]);
    od;

    # construct automorphisms by components
    for j in e do
      nj:=Concatenation(List(Filtered(e,i->i[1]<>j[1]),i->i[2]));
      r:=Length(j[2]);

      # the permutations and addition
      if r>1 then
        Add(au,GroupHomomorphismByImagesNC(G,G,Concatenation(nl,nj,j[2]),
            #(1,2)
            Concatenation(nl,nj,j[2]{[2]},j[2]{[1]},j[2]{[3..Length(j[2])]})));
        Add(au,GroupHomomorphismByImagesNC(G,G,Concatenation(nl,nj,j[2]),
            #(1,..,n)
            Concatenation(nl,nj,j[2]{[2..Length(j[2])]},j[2]{[1]})));
        #for k in [0..j[1]-1] do
        k:=0;
          Add(au,GroupHomomorphismByImagesNC(G,G,Concatenation(nl,nj,j[2]),
              #1->1+i^k*2
              Concatenation(nl,nj,[j[2][1]*j[2][2]^(i^k)],
                                  j[2]{[2..Length(j[2])]})));
        #od;
      fi;

      # multiplications

      for k in List( Flat( GeneratorsPrimeResidues(i^j[1])!.generators ),
              Int )  do

        Add(au,GroupHomomorphismByImagesNC(G,G,Concatenation(nl,nj,j[2]),
            #1->1^k
            Concatenation(nl,nj,[j[2][1]^k],j[2]{[2..Length(j[2])]})));
      od;

    od;

    # the mixing ones
    for j in [1..Length(e)] do
      for k in [1..Length(e)] do
        if k<>j then
          nj:=Concatenation(List(e{Difference([1..Length(e)],[j,k])},i->i[2]));
          offs:=Maximum(0,e[k][1]-e[j][1]);
          if Length(e[j][2])=1 and Length(e[k][2])=1 then
            max:=Minimum(e[j][1],e[k][1])-1;
          else
            max:=0;
          fi;
          for m in [0..max] do
            Add(au,GroupHomomorphismByImagesNC(G,G,
               Concatenation(nl,nj,e[j][2],e[k][2]),
               Concatenation(nl,nj,[e[j][2][1]*e[k][2][1]^(i^(offs+m))],
                                    e[j][2]{[2..Length(e[j][2])]},e[k][2])));
          od;
        fi;
      od;
    od;
  od;

  if IsFpGroup(G) and IsWholeFamily(G) then
    # rewrite to standard generators
    gens:=GeneratorsOfGroup(G);
    au:=List(au,x->GroupHomomorphismByImagesNC(G,G,gens,List(gens,y->Image(x,y))));
  fi;

  for i in au do
    SetIsBijective(i,true);
    j:=MappingGeneratorsImages(i);
    if j[1]<>j[2] then
      SetIsInnerAutomorphism(i,false);
    fi;
    SetFilterObj(i,IsMultiplicativeElementWithInverse);
  od;

  au:= GroupByGenerators( au, IdentityMapping( G ) );
  SetIsAutomorphismGroup(au,true);
  SetIsFinite(au,true);

  SetInnerAutomorphismsAutomorphismGroup(au,TrivialSubgroup(au));

  if IsFinite(G) then
    SetIsFinite(au,true);
    SetIsGroupOfAutomorphismsFiniteGroup(au,true);
  fi;

  return au;
end);

#############################################################################
##
#F  IsomorphismAbelianGroups(<G>)
##
InstallGlobalFunction(IsomorphismAbelianGroups,function(G,H)
local o,p,gens,hens;

  # get standard generating system
  gens:=IndependentGeneratorsOfAbelianGroup(G);
  gens:=ShallowCopy(gens);

  # get standard generating system
  hens:=IndependentGeneratorsOfAbelianGroup(H);
  hens:=ShallowCopy(hens);

  o:=List(gens,Order);
  p:=List(hens,Order);

  SortParallel(o,gens);
  SortParallel(p,hens);

  if o<>p then
    return fail;
  fi;

  o:=GroupHomomorphismByImagesNC(G,H,gens,hens);
  SetIsBijective(o,true);

  return o;
end);

BindGlobal("OuterAutomorphismGeneratorsSimple",function(G)
local d,id,H,iso,aut,auts,i,all,hom,field,dim,P,diag,mats,gens,gal;
  if not IsNonabelianSimpleGroup(G) then
    return fail;
  fi;
  gens:=GeneratorsOfGroup(G);
  d:=DataAboutSimpleGroup(G);
  id:=d.idSimple;
  all:=false;
  if id.series="A" then
    if id.parameter=6 then
      # A6 is easy enough
      return fail;
    else
      H:=AlternatingGroup(id.parameter);
      if G=H then
        iso:=IdentityMapping(G);
      else
        iso:=IsomorphismGroups(G,H);
      fi;
      aut:=GroupGeneralMappingByImages(G,G,gens,
            List(gens,
              x->PreImagesRepresentative(iso,Image(iso,x)^(1,2))));
      auts:=[aut];
      all:=true;
    fi;

  elif id.series in ["L","2A"] or (id.series="C" and id.parameter[1]>2) then
    hom:=EpimorphismFromClassical(G);
    if hom=fail then return fail;fi;
    field:=FieldOfMatrixGroup(Source(hom));
    dim:=DimensionOfMatrixGroup(Source(hom));
    auts:=[];
    if Size(field)>2 then
      gal:=GaloisGroup(field);
      # Diagonal automorphisms
      if id.series="L" then
        P:=GL(dim,field);
      elif id.series="2A" then
        P:=GU(dim,id.parameter[2]);
        #gal:=GaloisGroup(GF(GF(id.parameter[2]),2));
      elif id.series="C" then
        if IsEvenInt(Size(field)) then
          P:=Source(hom);
        else
          P:=List(One(Source(hom)),ShallowCopy);
          for i in [1..Length(P)/2] do
            P[i][i]:=PrimitiveRoot(field);
          od;
          P:=ImmutableMatrix(field,P);
          if not ForAll(GeneratorsOfGroup(Source(hom)),
                     x->x^P in Source(hom)) then
            Error("changed shape!");
          elif P in Source(hom) then
            Error("P is in!");
          fi;

          P:=Group(Concatenation([P],GeneratorsOfGroup(Source(hom))));
          SetSize(P,Size(Source(hom))*2);
        fi;
      else
        Error("not yet done");
      fi;
      # Sufficiently many elements to get the mult. group
      aut:=Size(P)/Size(Source(hom));
      P:=GeneratorsOfGroup(P);
      if id.series="C" then
        # we know it's the first
        mats:=P{[1]};
        P:=P{[2..Length(P)]};
      else
        diag:=Group(One(field));
        mats:=[];
        while Size(diag)<aut do
          if not DeterminantMat(P[1]) in diag then
            diag:=ClosureGroup(diag,DeterminantMat(P[1]));
            Add(mats,P[1]);
          fi;
          P:=P{[2..Length(P)]};
        od;
      fi;
      auts:=Concatenation(auts,
        List(mats,s->GroupGeneralMappingByImages(G,G,gens,List(gens,x->
                  Image(hom,PreImagesRepresentative(hom,x)^s)))));

    else
      gal:=Group(()); # to force trivial
    fi;

    if Size(gal)>1 then
      # Galois
      auts:=Concatenation(auts,
        List(MinimalGeneratingSet(gal),
                s->GroupGeneralMappingByImages(G,G,gens,List(gens,x->
                  Image(hom,
--> --------------------

--> maximum size reached

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

[ Dauer der Verarbeitung: 0.36 Sekunden  (vorverarbeitet)  ]