#############################################################################
##
#W  clashom.gi                  GAP library                  Alexander Hulpke
##
##
#Y  (C) 1999 School Math and Comp. Sci., University of St Andrews, Scotland
#Y  Copyright (C) 2002 The GAP Group
##
##  This file contains functions that compute the conjugacy classes of a
##  finite group by homomorphic images.
##  Literature: A.H: Conjugacy classes in finite permutation groups via
##  homomorphic images, MathComp, to appear.
##


# get classes from classical group if possible.
BindGlobal("ClassesFromClassical",function(G)
local H,hom,d,cl;
  if IsPermGroup(G) and (IsNaturalAlternatingGroup(G)  or
    IsNaturalSymmetricGroup(G)) then
    return ConjugacyClasses(G); # there is a method for this
  fi;
  if not IsSimpleGroup(PerfectResiduum(G)) then
    return fail;
  fi;
  d:=DataAboutSimpleGroup(PerfectResiduum(G));
  if d.idSimple.series<>"L" then
    return fail;
  fi;
  hom:=EpimorphismFromClassical(G);
  if hom=fail then
    return fail;
  fi;
  # so far matrix classes only implemented for GL/SL
  if not (IsNaturalGL(Source(hom)) or IsNaturalSL(Source(hom))) then
    return fail;
  fi;

  cl:=ConjugacyClasses(Source(hom));
  cl:=List(cl,x->Image(hom,Representative(x)));
  cl:=ConjugacyClassesByRandomSearch(G:seed:=cl);
  return cl;
end);

#############################################################################
##
#F  ClassRepsPermutedTuples(<g>,<ran>)
##
##  computes representatives of the colourbars with colours selected from
##  <ran>.
BindGlobal("ClassRepsPermutedTuples",function(g,ran)
local anz,erg,pat,pat2,sym,nrcomp,coldist,stab,dc,i,j,k,sum,schieb,lstab,
      stabs,p;
  anz:=NrMovedPoints(g);
  sym:=SymmetricGroup(anz);
  erg:=[];
  stabs:=[];
  for nrcomp in [1..anz] do
    # all sorted colour distributions
    coldist:=Combinations(ran,nrcomp);
    for pat in OrderedPartitions(anz,nrcomp) do
      Info(InfoHomClass,3,"Pattern: ",pat);
      # compute the partition stabilizer
      stab:=[];
      sum:=0;
      for i in pat do
	schieb:=MappingPermListList([1..i],[sum+1..sum+i]);
	sum:=sum+i;
	stab:=Concatenation(stab,
	        List(GeneratorsOfGroup(SymmetricGroup(i)),j->j^schieb));
      od;
      stab:=Subgroup(sym,stab);
      dc:=List(DoubleCosetRepsAndSizes(sym,stab,g),i->i[1]);

      # compute expanded pattern
      pat2:=[];
      for i in [1..nrcomp] do
        for j in [1..pat[i]] do
	  Add(pat2,i);
	od;
      od;

      for j in dc do
	# the new bar's stabilizer
	lstab:=Intersection(g,ConjugateSubgroup(stab,j));
	p:=Position(stabs,lstab);
	if p=fail then
	  Add(stabs,lstab);
	else
	  lstab:=stabs[p];
	fi;
	# the new bar
	j:=Permuted(pat2,j);
        for k in coldist do
	  Add(erg,[List(j,i->k[i]),lstab]);
	od;
      od;
    od;
  od;
  return erg;
end);

#############################################################################
##
#F  ConjugacyClassesSubwreath(<F>,<M>,<n>,<autT>,<T>,<Lloc>,<comp>,<emb>,<proj>)
##
InstallGlobalFunction(ConjugacyClassesSubwreath,
function(F,M,n,autT,T,Lloc,components,embeddings,projections)
local clT,	# classes T
      lcl,	# Length(clT)
      clTR,	# classes under other group (autT,centralizer)
      fus,	# class fusion
      sci,	# |centralizer_i|
      oci,	# |reps_i|
      i,j,k,l,	# loop
      pfus,	# potential fusion
      op,	# operation of F on components
      ophom,	# F -> op
      clF,	# classes of F
      clop,	# classes of op
      bars,	# colour bars
      barsi,    # partial bars
      lallcolors,# |all colors|
      reps,Mproj,centralizers,centindex,emb,pi,varpi,newreps,newcent,
      newcentindex,centimages,centimgindex,C,p,P,selectcen,select,
      cen,eta,newcentlocal,newcentlocalindex,d,dc,s,t,elm,newcen,shift,
      cengen,b1,ore,
      # as in paper
      colourbar,newcolourbar,possiblecolours,potentialbars,bar,colofclass,
      clin,clout,
      etas,	# list of etas
      opfun,	# operation function
      r,rp,	# op-element complement in F
      cnt,
      brp,bcen,
      centralizers_r, # centralizers of r
      newcent_r,# new list to buid
      centrhom, # projection \rest{centralizer of r}
      localcent_r, # image
      cr,
      isdirprod,# is just M a direct product
      genpos,	# generator index
      genpos2,
      gen,	# generator
      stab,	# stabilizer
      stgen,	# local stabilizer generators
      trans,
      repres,
      img,
      limg,
      con,
      pf,
      orb,	# orbit
      orpo,	# orbit position
      minlen,	# minimum orbit length
      remainlen,#list of remaining lengths
      gcd,	# gcd of remaining orbit lengths
      stabtrue,
      diff,
      possible,
      combl,
      smacla,
      smare,
      ppos,
      maxdiff,
      again,	# run orbit again to get all
      trymap,	# operation to try
      skip,	# skip (if u=ug)
      ug,	# u\cap u^{gen^-1}
      scj,	# size(centralizers[j])
      dsz;	# Divisors(scj);


  Info(InfoHomClass,1,
       "ConjugacyClassesSubwreath called for almost simple group of size ",
        Size(T));
  isdirprod:=Size(M)=Size(autT)^n;

  # classes of T
  clT:=ClassesFromClassical(T);
  if clT=fail then
    clT:=ConjugacyClassesByRandomSearch(T);
  fi;

  clT:=List(clT,i->[Representative(i),Centralizer(i)]);
  lcl:=Length(clT);
  Info(InfoHomClass,1,"found ",lcl," classes in almost simple");
  clTR:=List(clT,i->ConjugacyClass(autT,i[1]));

  # possible fusion under autT
  fus:=List([1..lcl],i->[i]);
  for i in [1..lcl] do
    sci:=Size(clT[i][2]);
    # we have taken a permutation representation that  prolongates to autT!
    oci:=CycleStructurePerm(clT[i][1]); 

    # we have tested already the smaller-# classes
    pfus:=Filtered([i+1..lcl],j->CycleStructurePerm(clT[j][1])=oci and
      Size(clT[j][2])=sci);
    pfus:=Difference(pfus,fus[i]);
    if Length(pfus)>0 then
      Info(InfoHomClass,3,"possible fusion ",pfus);
      for j in pfus do
        if clT[j][1] in clTR[i] then
	  fus[i]:=Union(fus[i],fus[j]);
	  # fuse the entries
	  for k in fus[i] do
	    fus[k]:=fus[i];
	  od;
	fi;
      od;
    fi;
  od;
  fus:=Set(fus); # throw out duplicates
  colofclass:=List([1..lcl],i->PositionProperty(fus,j->i in j));
  Info(InfoHomClass,2,"fused to ",Length(fus)," colours");

  # get the allowed colour bars
  ophom:=ActionHomomorphism(F,components,OnSets,"surjective");
  op:=Image(ophom);
  lallcolors:=Length(fus);
  bars:=ClassRepsPermutedTuples(op,[1..lallcolors]);

  Info(InfoHomClass,1,"classes in normal subgroup");
  # inner classes
  reps:=[One(M)];
  centralizers:=[M];
  centindex:=[1];
  colourbar:=[[]];

  Mproj:=[];
  varpi:=[];

  for i in [1..n] do
    Info(InfoHomClass,1,"component ",i);
    barsi:=Set(Immutable(List(bars,j->j[1]{[1..i]})));
    emb:=embeddings[i];
    pi:=projections[i];
    Add(varpi,ActionHomomorphism(M,Union(components{[1..i]}),"surjective"));
    Add(Mproj,Image(varpi[i],M));
    newreps:=[];
    newcent:=[];
    newcentindex:=[];
    centimages:=[];
    centimgindex:=[];
    newcolourbar:=[];

    etas:=[]; # etas for the centralizers

    # fuse centralizers that become the same
    for j in [1..Length(centralizers)] do
      C:=Image(pi,centralizers[j]);
      p:=Position(centimages,C);
      if p=fail then
        Add(centimages,C);
	p:=Length(centimages);
      fi;
      Add(centimgindex,p);

      # #force 'centralizers[j]' to have its base appropriate to the component
      # # (this will speed up preimages)
      # cen:=centralizers[j];
      # d:=Size(cen);
      # cen:=Group(GeneratorsOfGroup(cen),());
      # StabChain(cen,rec(base:=components[i],size:=d));
      # centralizers[j]:=cen;
      # etas[j]:=ActionHomomorphism(cen,components[i],"surjective");

    od;
    Info(InfoHomClass,2,Length(centimages)," centralizer images");

    # consider previous centralizers
    for j in [1..Length(centimages)] do
      # determine all reps belonging to this centralizer
      C:=centimages[j];
      selectcen:=Filtered([1..Length(centimgindex)],k->centimgindex[k]=j);
      Info(InfoHomClass,2,"Number ",j,": ",Length(selectcen),
            " previous centralizers to consider");
      
      # 7'
      select:=Filtered([1..Length(centindex)],k->centindex[k] in selectcen);
      # Determine the addable colours
      if i=1 then 
	possiblecolours:=[1..Length(fus)];
      else
	possiblecolours:=[];
	#for k in select do
	#  bar:=colourbar[k];
	k:=1;
	while k<=Length(select) 
	  and Length(possiblecolours)<lallcolors do
	  bar:=colourbar[select[k]];
	  potentialbars:=Filtered(bars,j->j[1]{[1..i-1]}=bar);
	  UniteSet(possiblecolours,
	           potentialbars{[1..Length(potentialbars)]}[1][i]);
	  k:=k+1;
	od;

      fi;

      for k in Union(fus{possiblecolours}) do
	# double cosets
	if Size(C)=Size(T) then
	  dc:=[One(T)];
	else

	  Assert(1,IsSubgroup(T,clT[k][2]));
	  Assert(1,IsSubgroup(T,C));

	  dc:=List(DoubleCosetRepsAndSizes(T,clT[k][2],C),i->i[1]);
	fi;
	for t in selectcen do
	  # continue partial rep. 

#	  #force 'centralizers[j]' to have its base appropriate to the component
#	  # (this will speed up preimages)
#	  if not (HasStabChainMutable(cen) 
#	     and i<=Length(centralizers)
#	     and BaseStabChain(StabChainMutable(cen))[1] in centralizers[i])
#	    then
#	    d:=Size(cen);
#	    cen:= Group( GeneratorsOfGroup( cen ), One( cen ) );
#	    StabChain(cen,rec(base:=components[i],size:=d));
#	    #centralizers[t]:=cen;
#	  fi;

	  cen:=centralizers[t];

	  if not IsBound(etas[t]) then
	    if Number(etas,i->IsBound(i))>500 then
	      for d in
		Filtered([1..Length(etas)],i->IsBound(etas[i])){[1..500]} do
		Unbind(etas[d]);
	      od;
	    fi;
	    etas[t]:=ActionHomomorphism(cen,components[i],"surjective");
	  fi;
	  eta:=etas[t];

	  select:=Filtered([1..Length(centindex)],l->centindex[l]=t);
	  Info(InfoHomClass,3,"centralizer nr.",t,", ",
	       Length(select)," previous classes");
	  newcentlocal:=[];
	  newcentlocalindex:=[];

	  for d in dc do
	    for s in select do
	      # test whether colour may be added here
	      bar:=Concatenation(colourbar[s],[colofclass[k]]);
	      bar:=ShallowCopy(colourbar[s]);
	      Add(bar,colofclass[k]);
	      MakeImmutable(bar);
	      #if ForAny(bars,j->j[1]{[1..i]}=bar) then
	      if bar in barsi then
		# new representative
		elm:=reps[s]*Image(emb,clT[k][1]^d);
		if elm in Mproj[i] then
		  # store the new element
		  Add(newreps,elm);
		  Add(newcolourbar,bar);
		  if i<n then # we only need the centralizer for further 
		              # components
		    newcen:=ClosureGroup(Lloc,
		              List(GeneratorsOfGroup(clT[k][2]),g->g^d));
		    p:=Position(newcentlocal,newcen);
		    if p=fail then
		      Add(newcentlocal,newcen);
		      p:=Length(newcentlocal);
		    fi;
		    Add(newcentlocalindex,p);
		  else
		    Add(newcentlocalindex,1); # dummy, just for counting
		  fi;
		#else
		#  Info(InfoHomClass,5,"not in");
		fi;

	      #else
	      #	Info(InfoHomClass,5,bar,"not minimal");
	      fi;
	      # end the loops from step 9
	    od;
	  od;
	  Info(InfoHomClass,2,Length(newcentlocalindex),
	       " new representatives");

	  if i<n then # we only need the centralizer for further components

	    # Centralizer preimages
	    shift:=[];
	    for l in [1..Length(newcentlocal)] do
	      P:=PreImage(eta,Intersection(Image(eta),newcentlocal[l]));

	      p:=Position(newcent,P);
	      if p=fail then
		Add(newcent,P);
		p:=Length(newcent);
	      fi;
	      shift[l]:=p;
	    od;

	    # move centralizer indices to global
	    for l in newcentlocalindex do
	      Add(newcentindex,shift[l]);
	    od;

	  fi;

	# end the loops from step 6,7 and 8
	od;
      od;
    od;

    centralizers:=newcent;
    centindex:=newcentindex;
    reps:=newreps;
    colourbar:=newcolourbar;
    # end the loop of step 2.
  od;

  Info(InfoHomClass,1,Length(reps)," classreps constructed");

  # further fusion among bars
  newreps:=[];
  Info(InfoHomClass,2,"computing centralizers");
  for bar in bars do
    b1:=Immutable(bar[1]);
    select:=Filtered([1..Length(reps)],i->colourbar[i]=b1);
    if Length(select)>1 then
      Info(InfoHomClass,2,"test ",Length(select)," classes for fusion");
    fi;
    newcentlocal:=[];
    for i in [1..Length(select)] do
      if not ForAny(newcentlocal,j->reps[select[i]] in j) then
	#AH we could also compute the centralizer
        C:=Centralizer(F,reps[select[i]]);
	Add(newreps,[reps[select[i]],C]);
	if i<Length(select) and Size(bar[2])>1 then
	  # there are other reps with the same bar left and the bar
	  # stabilizer is bigger than M
	  if not IsBound(bar[2]!.colstabprimg) then
	    # identical stabilizers have the same link. Therefore store the
	    # preimage in them
	    bar[2]!.colstabprimg:=PreImage(ophom,bar[2]);
	  fi;
	  # any fusion would take place in the stabilizer preimage
	  # we know that C must fix the bar, so it is the centralizer there.
	  r:=ConjugacyClass(bar[2]!.colstabprimg,reps[select[i]],C);
	  Add(newcentlocal,r);
	fi;
      fi;
    od;
  od;

  Info(InfoHomClass,1,"fused to ",Length(newreps)," inner classes");
  clF:=newreps;
  clin:=ShallowCopy(clF);
  Assert(1,Sum(clin,i->Index(F,i[2]))=Size(M));
  clout:=[];

  # outer classes

  clop:=Filtered(ConjugacyClasses(op),i->Order(Representative(i))>1);

  for k in clop do
    Info(InfoHomClass,1,"lifting class ",Representative(k));

    r:=PreImagesRepresentative(ophom,Representative(k));
    # try to make r of small order
    rp:=r^Order(Representative(k));
    rp:=RepresentativeAction(M,Concatenation(components),
                                  Concatenation(OnTuples(components[1],rp^-1),
				  Concatenation(components{[2..n]})),OnTuples);
    if rp<>fail then
      r:=r*rp;
    else
      Info(InfoHomClass,2,
           "trying random modification to get large centralizer");
      cnt:=LogInt(Size(autT),2)*10;
      brp:=();
      bcen:=Size(Centralizer(F,r));
      repeat
        rp:=Random(M);
	cengen:=Size(Centralizer(M,r*rp));
        if cengen>bcen then
	  bcen:=cengen;
	  brp:=rp;
	  cnt:=LogInt(Size(autT),2)*10;
	else
	  cnt:=cnt-1;
	fi;
      until cnt<0;
      r:=r*brp;
      Info(InfoHomClass,2,"achieved centralizer size ",bcen);
    fi;
    Info(InfoHomClass,2,"representative ",r);
    cr:=Centralizer(M,r);

    # first look at M-action
    reps:=[One(M)];
    centralizers:=[M];
    centralizers_r:=[cr];
    for i in [1..n] do;
      newreps:=[];
      newcent:=[];
      newcent_r:=[];
      opfun:=function(a,m)
               return Comm(r,m)*a^m;
             end;

      for j in [1..Length(reps)] do
	scj:=Size(centralizers[j]);
	dsz:=0;
	centrhom:=ActionHomomorphism(centralizers_r[j],components[i],
	            "surjective");
	localcent_r:=Image(centrhom);
	Info(InfoHomClass,4,i,":",j);
	Info(InfoHomClass,3,"acting: ",Size(centralizers[j])," minimum ",
	      Int(Size(Image(projections[i]))/Size(centralizers[j])),
	      " orbits.");
	# compute C(r)-classes
	clTR:=[];
	for l in clT do
	  Info(InfoHomClass,4,"DC",Index(T,l[2])," ",Index(T,localcent_r));
	  dc:=DoubleCosetRepsAndSizes(T,l[2],localcent_r);
	  clTR:=Concatenation(clTR,List(dc,i->l[1]^i[1]));
	od;

	orb:=[];
	for p in [1..Length(clTR)] do

	  repres:=PreImagesRepresentative(projections[i],clTR[p]);
	  if i=1 or isdirprod
	     or reps[j]*RestrictedPermNC(repres,components[i]) 
	            in Mproj[i] then
	    stab:=Centralizer(localcent_r,clTR[p]);
	    if Index(localcent_r,stab)<Length(clTR)/10 then
	      img:=Orbit(localcent_r,clTR[p]);
	      #ensure Representative is in first position
	      if img[1]<>clTR[p] then
	        genpos:=Position(img,clTR[p]);
		img:=Permuted(img,(1,genpos));
	      fi;
	    else
	      img:=ConjugacyClass(localcent_r,clTR[p],stab);
	    fi;
	    Add(orb,[repres,PreImage(centrhom,stab),img,localcent_r]);
	  fi;
	od;
	clTR:=orb;

	#was:
	#clTR:=List(clTR,i->ConjugacyClass(localcent_r,i));
	#clTR:=List(clTR,j->[PreImagesRepresentative(projections[i],
	#                                            Representative(j)),
	#	         PreImage(centrhom,Centralizer(j)),
	#		 j]);

	# put small classes to the top (to be sure to hit them and make
	# large local stabilizers)
	Sort(clTR,function(a,b) return Size(a[3])<Size(b[3]);end);

	Info(InfoHomClass,3,Length(clTR)," local classes");

	cengen:=GeneratorsOfGroup(centralizers[j]);
	#cengen:=Filtered(cengen,i->not i in localcent_r);

	while Length(clTR)>0 do

	  # orbit algorithm on classes
	  stab:=clTR[1][2];
	  orb:=[clTR[1]];
	  #repres:=RestrictedPermNC(clTR[1][1],components[i]);
	  repres:=clTR[1][1];
	  trans:=[One(M)];
	  select:=[2..Length(clTR)];

	  orpo:=1;
	  minlen:=Size(orb[1][3]);
	  possible:=false;
	  stabtrue:=false;
	  pf:=infinity;
	  maxdiff:=Size(T);
	  again:=0;
	  trymap:=false;
	  ug:=[];
	  # test whether we have full orbit and full stabilizer
	  while Size(centralizers[j])>(Sum(orb,i->Size(i[3]))*Size(stab)) do
	    genpos:=1;
	    while genpos<=Length(cengen) and
	      Size(centralizers[j])>(Sum(orb,i->Size(i[3]))*Size(stab)) do
	      gen:=cengen[genpos];
	      skip:=false;
	      if trymap<>false then
	        orpo:=trymap[1];
	        gen:=trymap[2];
		trymap:=false;
	      elif again>0 then
		if not IsBound(ug[genpos]) then
		  ug[genpos]:=Intersection(centralizers_r[j],
				   ConjugateSubgroup(centralizers_r[j],gen^-1));
		fi;
		if again<500 and ForAll(GeneratorsOfGroup(centralizers_r[j]),
		          i->i in ug[genpos])
		 then
		  # the random elements will give us nothing new
		  skip:=true;
		else
		  # get an element not in ug[genpos]
		  repeat
		    img:=Random(centralizers_r[j]);
		  until not img in ug[genpos] or again>=500;
		  gen:=img*gen;
	        fi;
	      fi;

	      if not skip then

		img:=Image(projections[i],opfun(orb[orpo][1],gen));

		smacla:=select;

		if not stabtrue then
		  p:=PositionProperty(orb,i->img in i[3]);
		  ppos:=fail;
                else
		  # we have the stabilizer and thus are only interested in
		  # getting new elements.
		  ppos:=PositionProperty(select,
			   i->Size(clTR[i][3])<=maxdiff and img in clTR[i][3]);
		  if ppos=fail then
		    p:="ignore"; #to avoid the first case
		  else
		    ppos:=select[ppos]; # get the right value
		    p:=fail; # go to first case
		  fi;
		fi;

		if p=fail then
		  if ppos=fail then
		    p:=First(select,
			   i->Size(clTR[i][3])<=maxdiff and img in clTR[i][3]);
                  else
		    p:=ppos;
		  fi;

		  RemoveSet(select,p);
		  Add(orb,clTR[p]);

		  #change the transversal element to map to the representative
		  con:=trans[orpo]*gen;
		  limg:=opfun(repres,con);
		  con:=con*PreImagesRepresentative(centrhom,
			   RepresentativeAction(localcent_r,
						 Image(projections[i],limg),
						 Representative(clTR[p][3])));
		  Assert(1,Image(projections[i],opfun(repres,con))
			   =Representative(clTR[p][3]));
		  Add(trans,con);
		  for stgen in GeneratorsOfGroup(clTR[p][2]) do
		    Assert( 1, IsOne( Image( projections[i],
				   opfun(repres,con*stgen/con)/repres ) ) );
		    stab:=ClosureGroup(stab,con*stgen/con);
		  od;

		  # compute new minimum length

		  if Length(select)>0 then
		    remainlen:=List(clTR{select},i->Size(i[3]));
		    gcd:=Gcd(remainlen);
		    diff:=minlen-Sum(orb,i->Size(i[3]));

		    if diff<0 then
		      # only go through this if the orbit actually grew
		      # larger
		      minlen:=Sum(orb,i->Size(i[3]));
		      repeat
			if dsz=0 then
			  dsz:=DivisorsInt(scj);
			fi;
			while not minlen in dsz do
			  # minimum gcd multiple to get at least the
			  # smallest divisor
			  minlen:=minlen+
			            (QuoInt((First(dsz,i->i>=minlen)-minlen-1),
				            gcd)+1)*gcd;
			od;

			# now try whether we actually can add orbits to make up
			# that length
			diff:=minlen-Sum(orb,i->Size(i[3]));
			Assert(1,diff>=0);
			# filter those remaining classes small enough to make
			# up the length
			smacla:=Filtered(select,i->Size(clTR[i][3])<=diff);
			remainlen:=List(clTR{smacla},i->Size(i[3]));
			combl:=1;
			possible:=false;
			if diff=0 then
			  possible:=fail;
			fi;
			while gcd*combl<=diff 
			      and combl<=Length(remainlen) and possible=false do
			  if NrCombinations(remainlen,combl)<100 then
			    possible:=ForAny(Combinations(remainlen,combl),
					     i->Sum(i)=diff);
			  else
			    possible:=fail;
			  fi;
			  combl:=combl+1;
			od;
			if possible=false then
			  minlen:=minlen+gcd;
			fi;
		      until possible<>false;
		    fi; # if minimal orbit length grew

		    Info(InfoHomClass,5,"Minimum length of this orbit ",
		         minlen," (",diff," missing)");

                  fi;

		  if minlen*Size(stab)=Size(centralizers[j]) then
		    #Assert(1,Length(smacla)>0);
		    maxdiff:=diff;
		    stabtrue:=true;
                  fi;

		elif not stabtrue then
		  # we have an element that stabilizes the conjugacy class.
		  # correct this to an element that fixes the representative.
		  # (As we have taken already the centralizer in
		  # centralizers_r, it is sufficient to correct by
		  # centralizers_r-conjugation.)
		  con:=trans[orpo]*gen;
		  limg:=opfun(repres,con);
		  con:=con*PreImagesRepresentative(centrhom,
			   RepresentativeAction(localcent_r,
						 Image(projections[i],limg),
						 Representative(orb[p][3])));
		  stab:=ClosureGroup(stab,con/trans[p]);
		  if Size(stab)*2*minlen>Size(centralizers[j]) then
		    Info(InfoHomClass,3,
		         "true stabilizer found (cannot grow)");
		    minlen:=Size(centralizers[j])/Size(stab);
		    maxdiff:=minlen-Sum(orb,i->Size(i[3]));
		    stabtrue:=true;
		  fi;
		fi;

		if stabtrue then

		  smacla:=Filtered(select,i->Size(clTR[i][3])<=maxdiff);

		  if Length(smacla)<pf then
		    pf:=Length(smacla);
		    remainlen:=List(clTR{smacla},i->Size(i[3]));

		    Info(InfoHomClass,3,
			"This is the true orbit length (missing ",
			maxdiff,")");

		    if Size(stab)*Sum(orb,i->Size(i[3]))
		        =Size(centralizers[j]) then
                      maxdiff:=0;

		    elif Sum(remainlen)=maxdiff then
		      Info(InfoHomClass,2,
			  "Full possible remainder must fuse");
		      orb:=Concatenation(orb,clTR{smacla});
		      select:=Difference(select,smacla);

		    else
		      # test whether there is only one possibility to get
		      # this length
		      if Length(smacla)<20 and
		       Sum(List([1..Minimum(Length(smacla),
		                    Int(maxdiff/gcd+1))],
			   x-> NrCombinations(smacla,x)))<10000 then
			# get all reasonable combinations
			smare:=[1..Length(smacla)]; #range for smacla
			combl:=Concatenation(List([1..Int(maxdiff/gcd+1)],
					      i->Combinations(smare,i)));
			# pick those that have the correct length
			combl:=Filtered(combl,i->Sum(remainlen{i})=maxdiff);
			if Length(combl)>1 then
			  Info(InfoHomClass,3,"Addendum not unique (",
			  Length(combl)," possibilities)");
			  if (maxdiff<10 or again>0) 
			    and ForAll(combl,i->Length(i)<=5) then
			    # we have tried often enough, now try to pick the
			    # right ones 
			    possible:=false;
			    combl:=Union(combl);
			    combl:=smacla{combl};
			    genpos2:=1;
			    smacla:=[];
			    while possible=false and Length(combl)>0 do
			      img:=Image(projections[i],
				opfun(clTR[combl[1]][1],cengen[genpos2]));
			      p:=PositionProperty(orb,i->img in i[3]);
			      if p<>fail then
				# it is!
				Info(InfoHomClass,4,"got one!");

				# remember the element to try
				trymap:=[p,(cengen[genpos2]*
				  PreImagesRepresentative(
				    RestrictedMapping(projections[i],
				      centralizers[j]),
				    RepresentativeAction(
				    orb[p][4],
				    img,Representative(orb[p][3]))  ))^-1];

				Add(smacla,combl[1]);
				combl:=combl{[2..Length(combl)]};
				if Sum(clTR{smacla},i->Size(i[3]))=maxdiff then
				  # bingo!
				  possible:=true;
				fi;
			      fi;
			      genpos2:=genpos2+1;
			      if genpos2>Length(cengen) then
				genpos2:=1;
				combl:=combl{[2..Length(combl)]};
			      fi;
			    od;
			    if possible=false then
			      Info(InfoHomClass,4,"Even test failed!");
			    else
			      orb:=Concatenation(orb,clTR{smacla});
			      select:=Difference(select,smacla);
			      Info(InfoHomClass,3,"Completed orbit (hard)");
			    fi;
			  fi;
			else
			  combl:=combl[1];
			  orb:=Concatenation(orb,clTR{smacla{combl}});
			  select:=Difference(select,smacla{combl});
			  Info(InfoHomClass,3,"Completed orbit");
			fi;
		      fi;
		    fi;
		  fi;

	        fi;
	      else
	        Info(InfoHomClass,5,"skip");
	      fi; # if not skip

	      genpos:=genpos+1;
	    od;
	    orpo:=orpo+1;
	    if orpo>Length(orb) then
	      Info(InfoHomClass,3,"Size factor:",EvalF(
	      (Sum(orb,i->Size(i[3]))*Size(stab))/Size(centralizers[j])),
	      " orbit consists of ",Length(orb)," suborbits, iterating");

	      orpo:=1;
	      again:=again+1;
	    fi;
	  od;
	  Info(InfoHomClass,2,"Stabsize = ",Size(stab),
		", centstabsize = ",Size(orb[1][2]));
	  clTR:=clTR{select};

	  Info(InfoHomClass,2,"orbit consists of ",Length(orb)," suborbits,",
	       Length(clTR)," classes left.");

	  Info(InfoHomClass,3,List(orb,i->Size(i[2])));
	  Info(InfoHomClass,4,List(orb,i->Size(i[3])));

          # select the orbit element with the largest local centralizer
	  orpo:=1;
	  p:=2;
	  while p<=Length(orb) do
	    if IsBound(trans[p]) and Size(orb[p][2])>Size(orb[orpo][2]) then
	      orpo:=p;
	    fi;
	    p:=p+1;
	  od;
	  if orpo<>1 then
	    Info(InfoHomClass,3,"switching to orbit position ",orpo);
	    repres:=opfun(repres,trans[orpo]);
	    #repres:=RestrictedPermNC(clTR[1][1],repres);
	    stab:=stab^trans[orpo];
	  fi;


	  Assert(1,ForAll(GeneratorsOfGroup(stab),
                j -> IsOne( Image(projections[i],opfun(repres,j)/repres) )));

	  # correct stabilizer to element stabilizer
	  Add(newreps,reps[j]*RestrictedPermNC(repres,components[i]));
	  Add(newcent,stab);
	  Add(newcent_r,orb[orpo][2]);
	od;

      od;
      reps:=newreps;
      centralizers:=newcent;
      centralizers_r:=newcent_r;

      Info(InfoHomClass,2,Length(reps)," representatives");
    od;

    select:=Filtered([1..Length(reps)],i->reps[i] in M);
    reps:=reps{select};
    reps:=List(reps,i->r*i);
    centralizers:=centralizers{select};
    centralizers_r:=centralizers_r{select};
    Info(InfoHomClass,1,Length(reps)," in M");

    # fuse reps if necessary
    cen:=PreImage(ophom,Centralizer(k));
    newreps:=[];
    newcentlocal:=[];
    for i in [1..Length(reps)] do
      bar:=CycleStructurePerm(reps[i]);
      ore:=Order(reps[i]);
      newcentlocal:=Filtered(newreps,
		     i->Order(Representative(i))=ore and
		     i!.elmcyc=bar);
      if not ForAny(newcentlocal,j->reps[i] in j) then
        C:=Centralizer(cen,reps[i]);
	# AH can we use centralizers[i] here ? 
	Add(clF,[reps[i],C]);
	Add(clout,[reps[i],C]);
	bar:=ConjugacyClass(cen,reps[i],C);
	bar!.elmcyc:=CycleStructurePerm(reps[i]);
	Add(newreps,bar);
      fi;
    od;
    Info(InfoHomClass,1,"fused to ",Length(newreps)," classes");
  od;
  
  Assert(1,Sum(clout,i->Index(F,i[2]))=Size(F)-Size(M));

  Info(InfoHomClass,2,Length(clin)," inner classes, total size =",
        Sum(clin,i->Index(F,i[2])));
  Info(InfoHomClass,2,Length(clout)," outer classes, total size =",
        Sum(clout,i->Index(F,i[2])));
  Info(InfoHomClass,3," Minimal ration for outer classes =",
	EvalF(Minimum(List(clout,i->Index(F,i[2])/(Size(F)-Size(M)))),30));

  Info(InfoHomClass,1,"returning ",Length(clF)," classes");

  Assert(1,Sum(clF,i->Index(F,i[2]))=Size(F));
  return clF;

end);

InstallGlobalFunction(ConjugacyClassesFittingFreeGroup,function(G)
local cs,	# chief series of G
      i,	# index cs
      cl,	# list [classrep,centralizer]
      hom,	# G->G/cs[i]
      M,	# cs[i-1]
      N,	# cs[i]
      subN,	# maximan normal in M over N
      csM,	# orbit of nt in M under G
      n,	# Length(csM)
      T,	# List of T_i
      Q,	# Action(G,T)
      Qhom,	# G->Q and F->Q
      S,	# PreImage(Qhom,Stab_Q(1))
      S1,	# Action of S on T[1]
      deg1,	# deg (s1)
      autos,	# automorphism for action
      arhom,	# autom permrep list
      Thom,	# S->S1
      T1,	# T[1] Thom
      w,	# S1\wrQ
      wbas,	# base of w
      emb,	# embeddings of w
      proj,	# projections of wbas
      components, # components of w
      reps,	# List reps in G for 1->i in Q
      F,	# action of G on M/N
      Fhom,	# G -> F
      FQhom,	# Fhom*Qhom
      genimages,# G.generators Fhom
      img,	# gQhom
      gimg,	# gFhom
      act,	# component permcation to 1
      j,k,	# loop
      C,	# Ker(Fhom)
      clF,	# classes of F
      ncl,	# new classes
      FM,	# normal subgroup in F, Fhom(M)
      FMhom,	# M->FM
      dc,	# double cosets
      jim,	# image of j
      Cim,
      CimCl,
      p,
      l,lj,
      l1,
      elm,
      zentr,
      onlysizes,
      good,bad,
      lastM;

  onlysizes:=ValueOption("onlysizes");
  # we assume the group has no solvable normal subgroup. Thus we get all
  # classes by lifts via nonabelian factors and can disregard all abelian
  # factors.

  # we will give classes always by their representatives in G and
  # centralizers by their full preimages in G.

  cs:= ChiefSeriesOfGroup( G );

  # the first step is always simple
  if HasAbelianFactorGroup(G,cs[2]) then
    # try to get the largest abelian factor
    i:=2;
    while i<Length(cs) and HasAbelianFactorGroup(G,cs[i+1]) do
      i:=i+1;
    od;
    cs:=Concatenation([G],cs{[i..Length(cs)]});
    # now cs[1]/cs[2] is the largest abelian factor

    cl:=List(RightTransversal(G,cs[2]),i->[i,G]);
  else
    # compute the classes of the simple nonabelian factor by random search
    hom:=NaturalHomomorphismByNormalSubgroupNC(G,cs[2]);
    cl:=ConjugacyClasses(Image(hom));
    cl:=List(cl,i->[PreImagesRepresentative(hom,Representative(i)),
		    PreImage(hom,StabilizerOfExternalSet(i))]);
  fi;
  lastM:=cs[2];

  for i in [3..Length(cs)] do
    # we assume that cl contains classreps/centralizers for G/cs[i-1]
    # we want to lift to G/cs[i]
    M:=cs[i-1];
    N:=cs[i];
    Info(InfoHomClass,1,i,":",Index(M,N),";  ",Size(N));
    if HasAbelianFactorGroup(M,N) then
      Info(InfoHomClass,2,"abelian factor ignored");
    else
      # nonabelian factor. Now it means real work.

      # 1) compute the action for the factor

      # first, we obtain the simple factors T_i/N.
      # we get these as intersections of the conjugates of the subnormal
      # subgroup

      csM:=CompositionSeries(M); # stored attribute
      if not IsSubset(csM[2],N) then
	# the composition series goes the wrong way. Now take closures of
	# its steps with N to get a composition series for M/N, take the
	# first proper factor for subN.
        n:=3;
	subN:=fail;
	while n<=Length(csM) and subN=fail do
	  subN:=ClosureGroup(N,csM[n]);
	  if Index(M,subN)=1 then
	    subN:=fail; # still wrong
	  fi;
	  n:=n+1;
	od;
      else
	subN:=csM[2];
      fi;
      
      if IsNormal(G,subN) then

	# only one -> Call standard process

	Fhom:=fail;
	# is this an almost top factor?
	if Index(G,M)<10 then
	  Thom:=NaturalHomomorphismByNormalSubgroupNC(G,subN);
	  T1:=Image(Thom,M);
	  S1:=Image(Thom);
	  if Size(Centralizer(S1,T1))=1 then
	    deg1:=NrMovedPoints(S1);
	    Info(InfoHomClass,2,
	      "top factor gives conjugating representation, deg ",deg1);

	    Fhom:=Thom;
	  fi;
	else
	  Thom:=NaturalHomomorphismByNormalSubgroupNC(M,subN);
	  T1:=Image(Thom,M);
	fi;

	if Fhom=fail then
	  autos:=List(GeneratorsOfGroup(G),
		    i->GroupHomomorphismByImagesNC(T1,T1,GeneratorsOfGroup(T1),
		      List(GeneratorsOfGroup(T1),
			    j->Image(Thom,PreImagesRepresentative(Thom,j)^i))));

	  # find (probably another) permutation rep for T1 for which all
	  # automorphisms can be represented by permutations
	  arhom:=AutomorphismRepresentingGroup(T1,autos);
	  S1:=arhom[1];
	  deg1:=NrMovedPoints(S1);
	  Fhom:=GroupHomomorphismByImagesNC(G,S1,GeneratorsOfGroup(G),arhom[3]);
	fi;


	C:=KernelOfMultiplicativeGeneralMapping(Fhom);
	F:=Image(Fhom,G);

	clF:=ClassesFromClassical(F);
	if clF=fail then
	  clF:=ConjugacyClassesByRandomSearch(F);
	fi;

	clF:=List(clF,j->[Representative(j),StabilizerOfExternalSet(j)]);

      else
	csM:=Orbit(G,subN); # all conjugates
	n:=Length(csM);

	if n=1 then
	  Error("this cannot happen");
	  T:=M;
	fi;

	T:=Intersection(csM{[2..Length(csM)]}); # one T_i
	if Length(GeneratorsOfGroup(T))>5 then
	  T:=Group(SmallGeneratingSet(T));
	fi;

	T:=Orbit(G,T); # get all the t's
	# now T[1] is a complement to csM[1] in G/N.
	
	# now compute the operation of G on M/N
	Qhom:=ActionHomomorphism(G,T,"surjective");
	Q:=Image(Qhom,G);
	S:=PreImage(Qhom,Stabilizer(Q,1)); 

	# find a permutation rep. for S-action on T[1]
	Thom:=NaturalHomomorphismByNormalSubgroupNC(T[1],N);
	T1:=Image(Thom);
	if not IsSubset([1..NrMovedPoints(T1)],
	                 MovedPoints(T1)) then
	  Thom:=Thom*ActionHomomorphism(T1,MovedPoints(T1),"surjective");
	fi;
	T1:=Image(Thom,T[1]);
	autos:=List(GeneratorsOfGroup(S),
		  i->GroupHomomorphismByImagesNC(T1,T1,GeneratorsOfGroup(T1),
		    List(GeneratorsOfGroup(T1),
			  j->Image(Thom,PreImagesRepresentative(Thom,j)^i))));

	# find (probably another) permutation rep for T1 for which all
	# automorphisms can be represented by permutations
	arhom:=AutomorphismRepresentingGroup(T1,autos);
	S1:=arhom[1];
	deg1:=NrMovedPoints(S1);
	Thom:=GroupHomomorphismByImagesNC(S,S1,GeneratorsOfGroup(S),arhom[3]);

	T1:=Image(Thom,T[1]);

	# now embed into wreath
	w:=WreathProduct(S1,Q);
	wbas:=DirectProduct(List([1..n],i->S1));
	emb:=List([1..n+1],i->Embedding(w,i));
	proj:=List([1..n],i->Projection(wbas,i));
	components:=WreathProductInfo(w).components;

	# define isomorphisms between the components
	reps:=List([1..n],i->
		PreImagesRepresentative(Qhom,RepresentativeAction(Q,1,i)));

	genimages:=[];
	for j in GeneratorsOfGroup(G) do
	  img:=Image(Qhom,j);
	  gimg:=Image(emb[n+1],img);
	  for k in [1..n] do
	    # look at part of j's action on the k-th factor.
	    # we get this by looking at the action of
	    #   reps[k] *   j    *   reps[k^img]^-1
	    # 1   ->    k  ->  k^img    ->           1
	    # on the first component. 
	    act:=reps[k]*j*(reps[k^img]^-1);
	    # this must be multiplied *before* permuting
	    gimg:=ImageElm(emb[k],ImageElm(Thom,act))*gimg;
	    gimg:=RestrictedPermNC(gimg,MovedPoints(w));
	  od;
	  Add(genimages,gimg);
	od;

	F:=Subgroup(w,genimages);
	if AssertionLevel()>0 then
	  Fhom:=GroupHomomorphismByImages(G,F,GeneratorsOfGroup(G),genimages);
	  Assert(1,fail<>Fhom);
	else
	  Fhom:=GroupHomomorphismByImagesNC(G,F,GeneratorsOfGroup(G),genimages);
	fi;
	C:=KernelOfMultiplicativeGeneralMapping(Fhom);

	Info(InfoHomClass,1,"constructed Fhom");

	# 2) compute the classes for F

	if n>1 then
          #if IsPermGroup(F) and NrMovedPoints(F)<18 then
	  #  # the old Butler/Theissen approach is still OK
	  #  clF:=[];
	  #  for j in 
	  #   Concatenation(List(RationalClasses(F),DecomposedRationalClass)) do
	  #    Add(clF,[Representative(j),StabilizerOfExternalSet(j)]);
	  #  od;
	  #else
	    FM:=F;
	    for j in components do
	      FM:=Stabilizer(FM,j,OnSets);
	    od;

	    clF:=ConjugacyClassesSubwreath(F,FM,n,S1,
		  Action(FM,components[1]),T1,components,emb,proj);
	  #fi;
	else
	  FM:=Image(Fhom,M);
	  Info(InfoHomClass,1,
	      "classes by random search in almost simple group");

	  clF:=ClassesFromClassical(F);
	  if clF=fail then
	    clF:=ConjugacyClassesByRandomSearch(F);
	  fi;

	  clF:=List(clF,j->[Representative(j),StabilizerOfExternalSet(j)]);
	fi;
      fi; # true orbit of T.

      Assert(1,Sum(clF,i->Index(F,i[2]))=Size(F));
      Assert(1,ForAll(clF,i->Centralizer(F,i[1])=i[2]));

      # 3) combine to form classes of sdp 

      # the length(cl)=1 gets rid of solvable stuff on the top we got ``too
      # early''.
      if IsSubgroup(N,KernelOfMultiplicativeGeneralMapping(Fhom)) then
	Info(InfoHomClass,1,
	    "homomorphism is faithful for relevant factor, take preimages");
	if Size(N)=1 and onlysizes=true then
	  cl:=List(clF,i->[PreImagesRepresentative(Fhom,i[1]),Size(i[2])]);
	else
	  cl:=List(clF,i->[PreImagesRepresentative(Fhom,i[1]),
			    PreImage(Fhom,i[2])]);
        fi;
      else
	Info(InfoHomClass,1,"forming subdirect products");

	FM:=Image(Fhom,lastM);
	FMhom:=RestrictedMapping(Fhom,lastM);
	if Index(F,FM)=1 then
	  Info(InfoHomClass,1,"degenerated to direct product");
	  ncl:=[];
	  for j in cl do
	    for k in clF do
	      # modify the representative with a kernel elm. to project
	      # correctly on the second component
	      elm:=j[1]*PreImagesRepresentative(FMhom,
			  LeftQuotient(Image(Fhom,j[1]),k[1]));
	      zentr:=Intersection(j[2],PreImage(Fhom,k[2]));
	      Assert(2,ForAll(GeneratorsOfGroup(zentr),
		      i->Comm(i,elm) in N));
	      Add(ncl,[elm,zentr]);
	    od;
	  od;

	  cl:=ncl;

	else

	  # first we add the centralizer closures and sort by them
	  # (this allows to reduce the number of double coset calculations)
	  ncl:=[];
	  for j in cl do
	    Cim:=Image(Fhom,j[2]);
	    CimCl:=Cim;
	    #CimCl:=ClosureGroup(FM,Cim); # should be unnecessary, as we took
	    # the full preimage
	    p:=PositionProperty(ncl,i->i[1]=CimCl);
	    if p=fail then
	      Add(ncl,[CimCl,[j]]);
	    else
	      Add(ncl[p][2],j);
	    fi;
	  od;

	  Qhom:=NaturalHomomorphismByNormalSubgroupNC(F,FM);
	  Q:=Image(Qhom);
	  FQhom:=Fhom*Qhom;

	  # now construct the sdp's
	  cl:=[];
	  for j in ncl do
	    lj:=List(j[2],i->Image(FQhom,i[1]));
	    for k in clF do
	      # test whether the classes are potential mates
	      elm:=Image(Qhom,k[1]);
	      if not ForAll(lj,i->RepresentativeAction(Q,i,elm)=fail) then

		#l:=Image(Fhom,j[1]);
		      
		if Index(F,j[1])=1 then
		  dc:=[()];
		else
		  dc:=List(DoubleCosetRepsAndSizes(F,k[2],j[1]),i->i[1]);
		fi;
		good:=0;
		bad:=0;
		for l in j[2] do
		  jim:=Image(FQhom,l[1]);
		  for l1 in dc do
		    elm:=k[1]^l1;
		    if Image(Qhom,elm)=jim then
		      # modify the representative with a kernel elm. to project
		      # correctly on the second component
		      elm:=l[1]*PreImagesRepresentative(FMhom,
				  LeftQuotient(Image(Fhom,l[1]),elm));
		      zentr:=PreImage(Fhom,k[2]^l1);
		      zentr:=Intersection(zentr,l[2]);

		      Assert(2,ForAll(GeneratorsOfGroup(zentr),
			      i->Comm(i,elm) in N));

		      Info(InfoHomClass,4,"new class, order ",Order(elm),
			  ", size=",Index(G,zentr));
		      Add(cl,[elm,zentr]);
		      good:=good+1;
		    else
		      Info(InfoHomClass,5,"not in");
		      bad:=bad+1;
		    fi;
		  od;
		od;
		Info(InfoHomClass,4,good," good, ",bad," bad of ",Length(dc));
	      fi;
	    od;
	  od;
        fi; # real subdirect product

      fi; # else Fhom not faithful on factor

      # uff. That was hard work. We're finally done with this layer.
      lastM:=N;
    fi; # else nonabelian
    Info(InfoHomClass,1,"so far ",Length(cl)," classes computed");
  od;

  if Length(cs)<3 then
    Info(InfoHomClass,1,"Fitting free factor returns ",Length(cl)," classes");
  fi;
  Assert( 1, Sum( List( cl, pair -> Size(G) / Size( pair[2] ) ) ) = Size(G) );
  return cl;
end);

## Lifting code, using new format and compatible with matrix groups

#############################################################################
##
#F  FFClassesVectorSpaceComplement( <N>, <p>, <gens>, <howmuch> )
##
##  This function creates a record  containing information about a complement
##  in <N> to the span of <gens>.
##
# field, dimension, subgenerators (as vectors),howmuch
BindGlobal("FFClassesVectorSpaceComplement",function(field,r, Q,howmuch )
local   zero,  one,  ran,  n,  nan,  cg,  pos,  i,  j,  v;
    
    one:=One( field);  zero:=Zero(field);
    ran:=[ 1 .. r ];
    n:=Length( Q );    nan:=[ 1 .. n ];
    
    cg:=rec( matrix        :=[  ],
	       one           :=one,
               baseComplement:=ShallowCopy( ran ),
               commutator    :=0,
               centralizer   :=0,
               dimensionN    :=r,
               dimensionC    :=n );

    if n = 0  or  r = 0  then
        cg.inverse:=NullMapMatrix;
	cg.projection    :=IdentityMat( r, one );
	cg.needed    :=[];
        return cg;
    fi;
    
    for i  in nan  do
        cg.matrix[ i ]:=Concatenation( Q[ i ], zero * nan );
        cg.matrix[ i ][ r + i ]:=one;
    od;
    TriangulizeMat( cg.matrix );
    pos:=1;
    for v  in cg.matrix  do
        while v[ pos ] = zero  do
            pos:=pos + 1;
        od;
        RemoveSet( cg.baseComplement, pos );
        if pos <= r  then  cg.commutator :=cg.commutator  + 1;
                     else  cg.centralizer:=cg.centralizer + 1;  fi;
    od;

    if howmuch=1 then
      return Immutable(cg);
    fi;

    cg.needed        :=[  ];
    cg.projection    :=IdentityMat( r, one );

    # Find a right pseudo inverse for <Q>.
    Append( Q, cg.projection );
    Q:=MutableTransposedMat( Q );
    TriangulizeMat( Q );
    Q:=TransposedMat( Q );
    i:=1;
    j:=1;
    while i <= r  do
        while j <= n and Q[ j ][ i ] = zero  do
            j:=j + 1;
        od;
        if j <= n and Q[ j ][ i ] <> zero  then
            cg.needed[ i ]:=j;
        else

            # If <Q> does  not  have full rank, terminate when the bottom row
            # is reached.
            i:=r;

        fi;
        i:=i + 1; 
    od;
    
    if IsEmpty( cg.needed )  then
        cg.inverse:=NullMapMatrix;
    else
        cg.inverse:=Q{ n + ran }
                       { [ 1 .. Length( cg.needed ) ] };
	cg.inverse:=ImmutableMatrix(field,cg.inverse,true);
    fi;
    if IsEmpty( cg.baseComplement )  then
        cg.projection:=NullMapMatrix;
    else

        # Find a base change matrix for the projection onto the complement.
        for i  in [ 1 .. cg.commutator ]  do
            cg.projection[ i ][ i ]:=zero;
        od;
        Q:=[  ];
        for i  in [ 1 .. cg.commutator ]  do
            Q[ i ]:=cg.matrix[ i ]{ ran };
        od;
        for i  in [ cg.commutator + 1 .. r ]  do
            Q[ i ]:=ListWithIdenticalEntries( r, zero );
            Q[ i ][ cg.baseComplement[ i-r+Length(cg.baseComplement) ] ]
             :=one;
        od;
        cg.projection:=cg.projection ^ Q;
        cg.projection:=cg.projection{ ran }{ cg.baseComplement };
	cg.projection:=ImmutableMatrix(field,cg.projection,true);
        
    fi;

    return cg;
end);

#############################################################################
##
#F  VSDecompCentAction( <N>, <h>, <C>, <howmuch> )
##
##  Given a homomorphism C -> N, c |-> [h,c],  this function determines (a) a
##  vector space decomposition N =  [h,C] + K with  projection onto K and (b)
##  the  ``kernel'' S <  C which plays   the role of  C_G(h)  in lemma 3.1 of
##  [Mecky, Neub\"user, Bull. Aust. Math. Soc. 40].
##
BindGlobal("VSDecompCentAction",function( pcgs, h, C, field,howmuch )
local   i,  tmp,  v,x,cg;

  i:=One(field);
  x:=List( C, c -> ExponentsOfPcElement(pcgs,Comm( h, c ) )*i);
#Print(Number(x,IsZero)," from ",Length(x),"\n");

  cg:=FFClassesVectorSpaceComplement(field,Length(pcgs),x,howmuch);
  tmp:=[  ];
  for i  in [ cg.commutator + 1 ..
	      cg.commutator + cg.centralizer ]  do
    v:=cg.matrix[ i ];
    tmp[ i - cg.commutator ]:=PcElementByExponentsNC( C,
	      v{ [ cg.dimensionN + 1 ..
		  cg.dimensionN + cg.dimensionC ] } );
  od;
  Unbind(cg.matrix);
  cg.cNh:=tmp;
  return cg;
end);

#############################################################################
##
#F  LiftClassesEANonsolvGeneral( <H>,<N>,<NT>,<cl> )
##
BindGlobal("LiftClassesEANonsolvGeneral",
  function( H, Npcgs, cl, hom, pcisom,solvtriv,fran)
    local  classes,    # classes to be constructed, the result
           correctingelement,
           field,      # field over which <N> is a vector space
	   one,
           h,          # preimage `cl.representative' under <hom>
	   cg,
           cNh,        # centralizer of <h> in <N>
           C,  gens,   # preimage `Centralizer( cl )' under <hom>
           r,          # dimension of <N>
           ran,        # constant range `[ 1 .. r ]'
           aff,        # <N> as affine space
           imgs,  M,   # generating matrices for affine operation
           orb,        # orbit of affine operation
           rep,# set of classes with canonical representatives
           c,  i, # loop variables
	   PPcgs,denomdepths,
	   reduce,
	   corr,
	   correctionfactor,
	   stabfac,
	   stabfacgens,
	   stabfacimg,
	   stabrad,
	   sz,gpsz,subsz,solvsz,
	   orblock,
	   stage,
	   b,j,
	   p,vp,genum,
	   st,gpe,
	   fe,epi,
	   repword,repwords,radidx,img,
	   radrange,farange,
	   ratio,
	   reps,
	   deno,docorrect,
	   k,failcnt,orpo,
	   stabilizergen,stabstack,
	   comm,s,stab;# for class correction

  correctingelement:=function(h,rep,fe)
  local comm;
    comm:=Comm( h, fe ) * Comm( rep, fe );
    comm:= ExponentsOfPcElement(Npcgs,comm)*one;
    ConvertToVectorRep(comm,field);
    comm := List(comm * cg.inverse,Int);
    comm:=PcElementByExponentsNC
      ( Npcgs, Npcgs{ cg.needed }, -comm );
    fe:=fe*comm;
    return fe;
  end;

  h := cl[1];

  field := GF( RelativeOrders( Npcgs )[ 1 ] );
  one:=One(field);
  PPcgs:=ParentPcgs(NumeratorOfModuloPcgs(Npcgs));
  denomdepths:=ShallowCopy(DenominatorOfModuloPcgs(Npcgs)!.depthsInParent);
  Add(denomdepths,Length(PPcgs)+1); # one
  
  # Determine the subspace $[h,N]$ and calculate the centralizer of <h>.
  #cNh := ExtendedPcgs( DenominatorOfModuloPcgs( N!.capH ),
  #               VSDecompCentAction( N, h, N!.capH ) );

  #oldcg:=KernelHcommaC(Npcgs,h,NumeratorOfModuloPcgs(Npcgs),2);

  #cg:=VSDecompCentAction( Npcgs, h, NumeratorOfModuloPcgs(Npcgs),field,2 );
  cg:=VSDecompCentAction( Npcgs, h, Npcgs,field,2 );
#Print("complen =",Length(cg.baseComplement)," of ",cg.dimensionN,"\n");
#if Length(Npcgs)>5 then Error("cb"); fi;

  cNh:=cg.cNh;

  r := Length( cg.baseComplement );
  ran := [ 1 .. r ];
  
  # Construct matrices for the affine operation on $N/[h,N]$.
  Info(InfoHomClass,4,"space=",Size(field),"^",r);
  if Size(field)^r>3*10^8 then Error("too large");fi;
  aff := ExtendedVectors( field ^ r );

  gens:=Concatenation(cl[2],Npcgs,cl[3]); # all generators
#if ForAny(gens,x->Order(x)=1) then Error("HUH2"); fi;
  gpsz:=cl[5];

  solvsz:=cl[6];

  radidx:=Length(Npcgs)+Length(cl[2]);

  imgs := [  ];
  for c  in gens  do
    M := [  ];
    for i  in [ 1 .. r ]  do
	M[ i ] := Concatenation( ExponentsConjugateLayer( Npcgs,
	      Npcgs[ cg.baseComplement[ i ] ] , c )
	      * cg.projection, [ Zero( field ) ] );
    od;
    M[ r + 1 ] := Concatenation( ExponentsOfPcElement
			  ( Npcgs, Comm( h, c ) ) * cg.projection,
			  [ One( field ) ] );

    M:=ImmutableMatrix(field,M);
    Add( imgs, M );
  od;

#if Size(field)^r>10^7 then Error("BIG");fi;

  # now compute orbits, being careful to get stabilizers in steps
  #orbreps:=[];
  #stabs:=[];

  orb:=OrbitsRepsAndStabsVectorsMultistage(gens{[1..radidx]},
	imgs{[1..radidx]},pcisom,solvsz,solvtriv,
	gens{[radidx+1..Length(gens)]},
	imgs{[radidx+1..Length(gens)]},cl[4],hom,gpsz,OnRight,aff);

  classes:=[];
  deno:=DenominatorOfModuloPcgs(Npcgs);
  for b in orb do
    rep := PcElementByExponentsNC( Npcgs, Npcgs{ cg.baseComplement },
		    b.rep{ ran } );
    #Assert(2,ForAll(GeneratorsOfGroup(stabsub),i->Comm(i,h*rep) in NT));
    stabrad:=ShallowCopy(b.stabradgens);
#Print("startdep=",List(stabrad,x->DepthOfPcElement(PPcgs,x)),"\n");
#if ForAny(stabrad,x->Order(x)=1) then Error("HUH3"); fi;
    stabfacgens:=b.stabfacgens;
    stabfacimg:=b.stabfacimgs;

    # correct generators. Partially in Pc Image
    if Length(cg.needed)>0 then

      stabrad:=List(stabrad,x->correctingelement(h,rep,x));
      # must make proper pcgs -- correction does not preserve that
      stabrad:=TFMakeInducedPcgsModulo(PPcgs,stabrad,denomdepths);

      # we change by radical elements, so the images in the factor don't
      # change
      stabfacgens:=List(stabfacgens,x->correctingelement(h,rep,x));

    fi;

    correctionfactor:=Characteristic(field)^Length(cg.needed);
    subsz:=b.subsz/correctionfactor;
    c := [h * rep,stabrad,stabfacgens,stabfacimg,subsz,
	   b.stabrsubsz/correctionfactor];
    Assert(3,Size(Group(Concatenation(DenominatorOfModuloPcgs(Npcgs),
       stabrad,stabfacgens)))=subsz);

    Add(classes,c);

  od;

  return classes;

end);

#############################################################################
##
#F  LiftClassesEANonsolvCentral( <H>, <N>, <cl> )
##
# the version for pc groups implicitly uses a pc-group orbit-stabilizer
# algorithm. We can't  do this but have to use a more simple-minded
# orbit/stabilizer approach.
BindGlobal("LiftClassesEANonsolvCentral",
  function( H, Npcgs, cl,hom,pcisom,solvtriv,fran )
local  classes,    	# classes to be constructed, the result
	field,         	# field over which <Npcgs> is a vector space
	o,
	n,r,		# dimensions
	space,
	com,
	comms,
	mats,
	decomp,
	reduce,
	gens,
	radidx,
	stabrad,stabfacgens,stabfacimg,stabrsubsz,relo,orblock,fe,st,
	orb,rep,reps,repword,repwords,p,stabfac,img,vp,genum,gpsz,
	subsz,solvsz,i,j,
	v,
	h,          	# preimage `cl.representative' under <hom>
	C,		# preimage `Centralizer( cl )' under <hom>
	w,    		# coefficient vectors for projection along $[h,N]$
	c;          	# loop variable

  field := GF( RelativeOrders( Npcgs )[ 1 ] );
  h := cl[1];
  #reduce:=ReducedPermdegree(C);
  #if reduce<>fail then
  #  C:=Image(reduce,C);
  #  Info(InfoHomClass,4,"reduced to deg:",NrMovedPoints(C));
  #  h:=Image(reduce,h);
  #  N:=ModuloPcgs(SubgroupNC(C,Image(reduce,NumeratorOfModuloPcgs(N))),
#	          SubgroupNC(C,Image(reduce,DenominatorOfModuloPcgs(N))));
#  fi;

  # centrality still means that conjugation by c is multiplication with
  # [h,c] and that the complement space is generated by commutators [h,c]
  # for a generating set {c|...} of C.

  field:=GF(RelativeOrders(Npcgs)[1]);
  n:=Length(Npcgs);
  o:=One(field);
  stabrad:=Concatenation(cl[2],Npcgs);
  radidx:=Length(stabrad);
  stabfacgens:=cl[3];
  stabfacimg:=cl[4];
  gpsz:=cl[5];
  subsz:=gpsz;
  solvsz:=cl[6];
  stabfac:=TrivialSubgroup(Image(hom));

  gens:=Concatenation(stabrad,stabfacgens); # all generators
  # commutator space basis

  comms:=List(gens,c->o*ExponentsOfPcElement(Npcgs,Comm(h,c)));
  List(comms,x->ConvertToVectorRep(x,field));
  space:=List(comms,ShallowCopy);
  TriangulizeMat(space);
  space:=Filtered(space,i->i<>Zero(i)); # remove spurious columns

  com:=BaseSteinitzVectors(IdentityMat(n,field),space);

  # decomposition of vectors into the subspace basis
  r:=Length(com.subspace);
  if r>0 then
    # if the subspace is trivial, everything stabilizes

    decomp:=Concatenation(com.subspace,com.factorspace)^-1;
    decomp:=decomp{[1..Length(decomp)]}{[1..r]};
    decomp:=ImmutableMatrix(field,decomp);

    # build matrices for the affine action
    mats:=[];
    for w in comms do
      c:=IdentityMat(r+1,o);
      c[r+1]{[1..r]}:=w*decomp; # translation bit
      c:=ImmutableMatrix(field,c);
      Add(mats,c);
    od;

    #subspace affine enumerator
    v:=ExtendedVectors(field^r);

    # orbit-stabilizer algorithm solv/nonsolv version
    #C := Stabilizer( C, v, v[1],GeneratorsOfGroup(C), mats,OnPoints );

    # assume small domain -- so no bother with bitlist
    orb:= [v[1]];
    reps:=[One(gens[1])];
    repwords:=[[]];
    stabrad:=[];
    stabrsubsz:=Size(solvtriv);

    vp:=1;

    for genum in [radidx,radidx-1..1] do
      relo:=RelativeOrders(pcisom!.sourcePcgs)[
	      DepthOfPcElement(pcisom!.sourcePcgs,gens[genum])];
      img:=orb[1]*mats[genum];
      repword:=repwords[vp];
      p:=Position(orb,img);
      if p=fail then
	for j in [1..Length(orb)*(relo-1)] do
	  img:=orb[j]*mats[genum];
	  Add(orb,img);
	  Add(reps,reps[j]*gens[genum]);
	  Add(repwords,repword);
	od;
      else
	rep:=gens[genum]/reps[p];
	Add(stabrad,rep);
	stabrsubsz:=stabrsubsz*relo;
      fi;

    od;
    stabrad:=Reversed(stabrad);

    Assert(1,solvsz=stabrsubsz*Length(orb));

    #nosolvable part
    orblock:=Length(orb);
    vp:=1;
    stabfacgens:=[];
    stabfacimg:=[];
    while vp<=Length(orb) do
      for genum in [radidx+1..Length(gens)] do
	img:=orb[vp]*mats[genum];
	rep:=reps[vp]*gens[genum];
	repword:=Concatenation(repwords[vp],[genum-radidx]);
	p:=Position(orb,img);
	if p=fail then
	  Add(orb,img);
	  Add(reps,rep);
	  Add(repwords,repword);
	  for j in [1..orblock-1] do
	    img:=orb[vp+j]*mats[genum];
    #if img in orb then Error("HUH");fi;
	    Add(orb,img);
	    Add(reps,reps[vp+j]*mats[genum]);
	    # repword stays the same!
	    Add(repwords,repword);
	  od;
	else
	  st:=rep/reps[p];
	  if Length(repword)>0 then
	    # build the factor group element
	    fe:=One(Image(hom));
	    for i in repword do
	      fe:=fe*cl[4][i];
	    od;
	    for i in Reversed(repwords[p]) do
	      fe:=fe/cl[4][i];
	    od;
	    if not fe in stabfac then
	      # not known -- add to generators
	      Add(stabfacgens,st);
	      Add(stabfacimg,fe);
	      stabfac:=ClosureGroup(stabfac,fe);
	    fi;
	  fi;
	fi;
      od;
      vp:=vp+orblock;
    od;

    subsz:=stabrsubsz*Size(stabfac);
  else
    stabrsubsz:=solvsz;
  fi;

  if Length(com.factorspace)=0 then
    classes :=[[h,stabrad,stabfacgens,stabfacimg,subsz,stabrsubsz]];
  else
    classes:=[];
    # enumerator of complement
    v:=field^Length(com.factorspace);
    for w in v do
      c := [h * PcElementByExponentsNC( Npcgs,w*com.factorspace),
	    stabrad,stabfacgens,stabfacimg,subsz,stabrsubsz];
      #if reduce<>fail then
  #	Add(classes,[PreImagesRepresentative(reduce,c[1]),
  #	  PreImage(reduce,c[2])]);
  #      else

  Assert(3,c[6]
    =Size(Group(Concatenation(c[2],DenominatorOfModuloPcgs(Npcgs)))));

      Add(classes,c);
  #      fi;
    od;
  fi;

#  Assert(1,ForAll(classes,i->i[1] in H and IsSubset(H,i[2])));
  return classes;
end);


#############################################################################
##
#F  LiftClassesEATrivRep
##
BindGlobal("LiftClassesEATrivRep",
  function( H, Npcgs, cl, fants,hom, pcisom,solvtriv)
    local  classes,    # classes to be constructed, the result
           h,field,one,solvsz,radidx,gens,imgs,M,bas,
	   gpsz,c,i,npcgsact,usent,dim,found,nsgens,ntgens,nsimgs,mo,
	   basis, ssidx,cpidx,compl,pcgsimgs,pcgssel,
	   sel,pcgs,fasize,nsfgens,fgens,a,norb,fstab,rep,reps,frep,freps,
	   orb,p,rsgens,el,img,j,basinv,newo,orbslev,ssd,result,o,subs,orbsub,
	   sgens,sfgens,z,minvecs,orpo,norpo,maxorb,
	   IteratedMinimizer,OrbitMinimizer,Minimizer,miss;

  npcgsact:=function(c)
    local M,i;
    M := [  ];
    for i  in [ 1 .. dim ]  do
	M[ i ] := ExponentsConjugateLayer( Npcgs,
	Npcgs[ i ] , c )*one;
    od;
    M:=ImmutableMatrix(field,M);
    return M;
  end;

  pcgs:=MappingGeneratorsImages(pcisom)[1];
  field:=GF(RelativeOrders(Npcgs)[1]);
  one:=One(field);
  dim:=Length(Npcgs);

  # action of group
  h := cl[1];
  gens:=Concatenation(cl[2],Npcgs,cl[3]); # all generators
  fgens:=Concatenation(ListWithIdenticalEntries(
	    Length(Npcgs)+Length(cl[2]),One(Range(hom))),cl[4]);
  gpsz:=cl[5];

  solvsz:=cl[6];

  radidx:=Length(Npcgs)+Length(cl[2]);
  imgs := [  ];
  for c  in gens  do
    Add( imgs, npcgsact(c));
  od;
  sel:=Filtered([1..Length(imgs)],x->Order(imgs[x])>1);

  usent:=0;
  found:=0;
  while usent<Length(fants) do
    usent:=usent+1;
    nsfgens:=NormalIntersection(fants[usent],Group(cl[4]));
    fasize:=Size(nsfgens);
    nsfgens:=SmallGeneratingSet(nsfgens);
    nsgens:=List(nsfgens,x->PreImagesRepresentative(hom,x));
    nsimgs:=List(Concatenation(pcgs,nsgens),npcgsact);
    mo:=GModuleByMats(nsimgs,field);
    if not MTX.IsIrreducible(mo) then
      # split space as direct sum under normal sub -- clifford Theory
      o:=MTX.BasesMinimalSubmodules(mo);
      if Length(o)>50 then
	o:=o{Set(List([1..50],x->Random([1..Length(o)])))};
      fi;

      for i in Filtered([1..Length(o)],
	  x->(mo.dimension mod Length(o[x])=0) and Length(o[x])>found) do
	# subspace and images as orbit
	bas:=o[i];
	ssd:=Length(bas);
	if found<ssd and Size(field)^ssd<3*10^7 then
	  Info(InfoHomClass,2,"Trying subspace ",ssd," in ",mo.dimension);
	  orbsub:=Orbit(Group(imgs{sel}),bas,OnSubspacesByCanonicalBasis);
	  if Length(orbsub)*Length(bas)<>Length(bas[1]) then
	    subs:=MTX.InducedActionSubmodule(mo,bas);
	    subs:=MTX.Homomorphisms(subs,mo);
	    orbsub:=Filtered(subs,x->x in orbsub);
	  fi;
	  if Length(orbsub)*Length(bas)=Length(bas[1]) then
	    found:=ssd;
	    el:=[orbsub,bas,fasize,nsgens,nsimgs,nsfgens,mo];

	    subs:=List([1..Length(orbsub)],x->[(x-1)*ssd+1..x*ssd]);
	    bas:=ImmutableMatrix(field,Concatenation(orbsub)); # this is the new basis
	    basinv:=bas^-1;
	  else
	    Info(InfoHomClass,3,"failed ",Length(orbsub));
	  fi;
	fi;
      od;
    fi;

  od;

  if found=0 then
    Info(InfoHomClass,3,"basic case");
    #Error("BASIC");
    return fail;
  else
    ssd:=found;
    #el is [orbsub,bas,fasize,nsgens,nsimgs,nsfgens,mo];
    orbsub:=el[1];
    bas:=el[2];
    fasize:=el[3];
    nsgens:=el[4];
    nsimgs:=el[5];
    nsfgens:=el[6];
    mo:=el[7];
    Info(InfoHomClass,2,"Using subspace ",ssd," in ",mo.dimension);

    subs:=List([1..Length(orbsub)],x->[(x-1)*ssd+1..x*ssd]);
    bas:=ImmutableMatrix(field,Concatenation(orbsub)); # this is the new basis
    basinv:=bas^-1;

  fi;

  imgs:=List(imgs,x->bas*x*basinv); # write wrt new basis

  # now determine N-orbits, stepwise

  solvtriv:=Subgroup(Range(pcisom),
      List(DenominatorOfModuloPcgs(Npcgs),x->ImagesRepresentative(pcisom,x)));

  orb:=[rec(len:=1,rep:=Zero(bas[1]),
	stabfacgens:=nsgens,
	stabfacimgs:=nsfgens,
	# only generators in factor
	stabradgens:=Filtered(pcgs,x->not x in DenominatorOfModuloPcgs(Npcgs)),
        stabrsubsz:=Size(Image(pcisom)),
	subsz:=fasize*Product(RelativeOrders(pcgs))
		   )];

  orbslev:=[];
  maxorb:=1;
  for i in [1..Length(subs)] do
    norb:=[];
    el:=Elements(VectorSpace(field,IdentityMat(Length(bas),field){subs[i]}));
    for o in orb do
      newo:= OrbitsRepsAndStabsVectorsMultistage(
	     o.stabradgens,List(o.stabradgens,x->bas*npcgsact(x)*basinv),
	     pcisom,o.stabrsubsz,solvtriv,
	     o.stabfacgens,List(o.stabfacgens,x->bas*npcgsact(x)*basinv),
	     o.stabfacimgs,hom,o.subsz,OnRight,
	     el);
      for j in newo do
	if j.len>maxorb then maxorb:=j.len;fi;
	if i>1 then
	  j.len:=j.len*o.len;
	  j.rep:=Concatenation(o.rep{[1..(i-1)*ssd]},j.rep{[(i-1)*ssd+1..Length(j.rep)]});
	  MakeImmutable(j.rep);
	fi;
	Add(norb,j);
      od;
    od;
    Info(InfoHomClass,3,"Level ",i," , ",Length(norb)," orbits");
    orb:=norb;
    Add(orbslev,ShallowCopy(orb));
  od;

  IteratedMinimizer:=function(vec,allcands)
  local i,stops,a,cands,mapper,fmapper,stabfacgens,stabradgens,stabfacimgs,
	range,lcands,lvec;
    cands:=allcands;
    mapper:=One(Source(hom));
    fmapper:=One(Range(hom));
    stabfacgens:=nsgens;
    stabfacimgs:=nsfgens;
    stabradgens:=pcgs;
    for i in [1..Length(subs)] do
      range:=[1..i*ssd];
      lcands:=Filtered(orbslev[i],
	x->ForAny(cands,y->y.rep{range}=x.rep{range}));
      lvec:=Concatenation(vec{range},Zero(vec{[i*ssd+1..Length(vec)]}));
      result:=OrbitMinimumMultistage(stabradgens,
	   List(stabradgens,x->bas*npcgsact(x)*basinv),
	   pcisom,0,0, # solvable size not needed
	   stabfacgens,
	   List(stabfacgens,x->bas*npcgsact(x)*basinv),
	   stabfacimgs,
	   hom,0, #gpsz not needed
	   OnRight,lvec,maxorb,#Maximum(List(lcands,x->x.len)),
	   Set(List(lcands,x->x.rep)));
      a:=First(lcands,x->x.rep{range}=result.min{range});
      mapper:=mapper*result.elm;
      fmapper:=fmapper*result.felm;
      vec:=vec*bas*npcgsact(result.elm)*basinv; # map vector to so far canonical
      # not all classes are feasible
      Assert(1,ForAny(cands,x->x.rep{range}=vec{range}));
      cands:=Filtered(cands,x->x.rep{range}=vec{range}); 
      stabradgens:=a.stabradgens;
      stabfacgens:=a.stabfacgens;
      stabfacimgs:=a.stabfacimgs;
    od;
    if Length(cands)<>1 then Error("nonunique");fi;
    return rec(elm:=mapper,felm:=fmapper,min:=vec,nclass:=cands[1]);
  end;

  pcgsimgs:=List(pcgs,x->bas*npcgsact(x)*basinv);
  pcgssel:=Filtered([1..Length(pcgs)],x->not IsOne(pcgsimgs[x]));
  nsimgs:=List(nsgens,x->bas*npcgsact(x)*basinv);


  OrbitMinimizer:=function(vec,allcands)
  local a;
    
  if false and allcands[1].len>1 then
    Error();
  fi;
    a:=OrbitMinimumMultistage(pcgs,pcgsimgs,pcisom,0,0,
	nsgens,nsimgs,nsfgens,hom,0,
	OnRight,vec,allcands[1].len,minvecs);
    a.nclass:=First(allcands,x->x.rep=a.min);
    return a;

  end;

  orpo:=NewDictionary(orb[Length(orb)].rep,true,field^Length(orb[1].rep));
  for p in [1..Length(orb)] do
    AddDictionary(orpo,orb[p].rep,p);
  od;

  # now do an orbit algorithm on orb. As the orbit is short no need for
  # two-step. 

  newo:=[];
  while Length(orb)>0 do
    # pick new one
    p:=First([1..Length(orb)],x->IsBound(orb[x]));
    norb:=[orb[p]];
    norpo:=[];
    norpo[p]:=1;
    el:=Filtered(orb,x->x.len=orb[p].len);
    minvecs:=Set(List(el,x->x.rep));
#el:=orbslev[3];
    if orb[p].len>30000 then
      Minimizer:=IteratedMinimizer;
    else
      Minimizer:=OrbitMinimizer;
    fi;

    # as Rad <=N we can assume that the radical part of the stabilizer
    # is known
    rsgens:=ShallowCopy(orb[p].stabradgens);
    a:=Difference([1..Length(gens)],sel);
    sgens:=Concatenation(orb[p].stabfacgens,gens{a});
    sfgens:=Concatenation(orb[p].stabfacimgs,fgens{a});
    fstab:=Group(sfgens);
    reps:=[One(Source(hom))];
    freps:=[One(Range(hom))];
    Unbind(orb[p]);

    # factor missing from stop
    miss:=cl[5]/(norb[1].len*Size(fstab)*norb[1].stabrsubsz);

    i:=1;
    while i<=Length(norb) and miss>1 do
      for j in sel do
	img:=OnRight(norb[i].rep,imgs[j]);
	img:=Minimizer(img,el);

	rep:=reps[i]*gens[j]*img.elm;
	frep:=freps[i]*fgens[j]*img.felm;
	p:=LookupDictionary(orpo,img.min);
	#p:=PositionProperty(norb,x->x.rep=img.min);
	if p=fail then
	  Error("unknown minimum");
	elif IsBound(norpo[p]) then
	  # old point
	  p:=norpo[p];
	  if miss>=2 then
	Assert(1,norb[i].rep*imgs[j]*bas*npcgsact(img.elm)*basinv=norb[p].rep);
    #Print("A",i," ",j," ",Length(el),"\n");
	    # old point -- stabilize
	    a:=frep/freps[p];
	    if not a in sfgens then
	      Add(sgens,rep/reps[p]);
	      Add(sfgens,a);
	      miss:=miss*Size(fstab);
	      fstab:=ClosureGroup(fstab,a);
	      miss:=miss/Size(fstab);
#Print("miss1:",EvalF(miss)," ",i," of ",Length(norb),"\n");

	    fi;
	  fi;
	else
  #Print("B",i," ",j," ",Length(el),"\n");
	  # new point
	  #p:=PositionProperty(orb,x->x.rep=img.min);
	Assert(1,norb[i].rep*imgs[j]*bas*npcgsact(img.elm)*basinv=orb[p].rep);
	  Add(norb,orb[p]);
	  norpo[p]:=Length(norb);
	  Add(reps,rep);
	  Add(freps,frep);
	  miss:=miss*(Length(norb)-1)/Length(norb);
#Print("miss3:",EvalF(miss)," ",i," of ",Length(norb),"\n");
	  # add conjugate stabilizer
	  #Append(rsgens,List(orb[p].stabradgens,x->rep*x/rep));
	  for z in [1..Length(orb[p].stabfacgens)] do
	    a:=frep*orb[p].stabfacimgs[z]/frep;
	    if not a in fstab then
	      Add(sgens,rep*orb[p].stabfacgens[z]/rep);
	      Add(sfgens,a);
	      miss:=miss*Size(fstab);
	      fstab:=ClosureGroup(fstab,a);
	      miss:=miss/Size(fstab);
#Print("miss2:",miss,"\n");
	    fi;
	  od;
	  Unbind(orb[p]);
	fi;
      od;
      i:=i+1;
    od;
if miss<>1 then 
  # something is dodgy -- fallback to the default algorithm
  return fail;Error("HEH?");fi;
    Info(InfoHomClass,3,"Fused ",Length(norb),"*",norb[1].len," ",
      Number(orb,x->IsBound(x))," left");
    Assert(0,ForAll(rsgens,x->norb[1].rep*bas*npcgsact(x)*basinv=norb[1].rep));
    Assert(0,ForAll(sgens,x->norb[1].rep*bas*npcgsact(x)*basinv=norb[1].rep));
#if ForAny(rsgens,x->Order(x)=1) then Error("HUH5"); fi;

    a:=[h*PcElementByExponents(Npcgs,norb[1].rep*bas),rsgens,sgens,sfgens,
	cl[5]/Length(norb)/norb[1].len, norb[1].stabrsubsz];

#rsgens:=List(rsgens,x->ImageElm(pcisom,x));
#if rsgens<>InducedPcgsByGenerators(FamilyPcgs(Range(pcisom)),rsgens) then
#  Error("nonpcgs!");
#fi;

    Add(newo,a);

  od;
  return newo;
end);

InstallGlobalFunction(ConjugacyClassesViaRadical,function (G)
local r,	#radical
      f,	# G/r
      hom,	# G->f
      pcgs,mpcgs, #(modulo) pcgs
      pcisom,
      gens,
      ser,	# series
      radsize,len,ntrihom,
      mran,nran,fran,
      central,
      fants,
      d,
      solvtriv,
      M,N,	# normal subgrops
      ind,	# indices
      i,	#loop
      new,	# new classes
      cl,ncl;	# classes

  # it seems to be cleaner (and avoids deferring abelian factors) if we
  # factor out the radical first. (Note: The radical method for perm groups
  # stores the nat hom.!)
  ser:=FittingFreeLiftSetup(G);
  radsize:=Product(RelativeOrders(ser.pcgs));
  len:=Length(ser.pcgs);

  if radsize=1 then
    cl:=ConjugacyClassesFittingFreeGroup(G:onlysizes:=false);
    ncl:=[];
    for i in cl do
      r:=ConjugacyClass(G,i[1],i[2]);
      Add(ncl,r);
    od;
    return ncl;
  fi;

  pcgs:=ser.pcgs;
  pcisom:=ser.pcisom;
  fants:=[];

  # store centralizers as rep, centralizer generators in radical,
  # centralizer generators with nontrivial
  # radfactor image, corresponding radfactor images
  # the generators in the radical do not list the generators of the
  # current layer after immediate lifting.

  if radsize=Size(G) then
    # solvable case
    hom:=ser.factorhom;
    d:=MappingGeneratorsImages(hom);
    mran:=Filtered([1..Length(d[2])],x->not IsOne(d[2][x]));
    cl:=[[One(G),[],d[1]{mran},d[2]{mran},Size(G),Size(G)]];
  else
    # nonsolvable
    if radsize>1 then
      hom:=ser.factorhom;
      ntrihom:=true;
      f:=Image(hom);
      # we need centralizers
      cl:=ConjugacyClassesFittingFreeGroup(f:onlysizes:=false);
      fants:=Filtered(NormalSubgroups(f),x->Size(x)>1 and Size(x)<Size(f));
    else
      if IsPermGroup(G) then
	hom:=SmallerDegreePermutationRepresentation(G);
	ntrihom:=not IsOne(hom);;
      else
        hom:=IdentityMapping(G);
	ntrihom:=false;
      fi;
      f:=Image(hom);
      cl:=ConjugacyClassesFittingFreeGroup(f);
    fi;

    if ntrihom then
      ncl:=[];
      for i in cl do
	new:=[PreImagesRepresentative(hom,i[1])];
	if not IsInt(i[2]) then
	  Add(new,[]); # no generators in radical yet
	  gens:=SmallGeneratingSet(i[2]);
	  Add(new,
	    List(gens,x->PreImagesRepresentative(hom,x)));
          Add(new,gens);
	  #TODO: PreImage groups?
	  #Add(new,PreImage(hom,i[2]));
	  Add(new,radsize*Size(i[2]));
	  Add(new,radsize);
	fi;
        Add(ncl,new);
      od;
      cl:=ncl;

    fi;
  fi;

  Assert(3,ForAll(cl,x->x[6]=Size(Group(Concatenation(x[2],pcgs)))));

  for d in [2..Length(ser.depths)] do
    #M:=ser[i-1];
    #N:=ser[i];
    mran:=[ser.depths[d-1]..len];
    nran:=[ser.depths[d]..len];
    fran:=[mran,nran];

    mpcgs:=InducedPcgsByPcSequenceNC(pcgs,pcgs{mran}) mod
           InducedPcgsByPcSequenceNC(pcgs,pcgs{nran});

    central:= ForAll(GeneratorsOfGroup(G),
		i->ForAll(mpcgs,
		  j->DepthOfPcElement(pcgs,Comm(i,j))>=ser.depths[d]));
    
    # abelian factor, use affine methods
    Info(InfoHomClass,1,"abelian factor ",d,": ",
      Product(RelativeOrders(ser.pcgs){mran}), "->",
      Product(RelativeOrders(ser.pcgs){nran})," central:",central);

    ncl:=[];
    solvtriv:=Subgroup(Range(pcisom),
	List(DenominatorOfModuloPcgs(mpcgs),x->ImagesRepresentative(pcisom,x)));
    for i in cl do
      #Assert(2,ForAll(GeneratorsOfGroup(i[2]),j->Comm(i[1],j) in M));
      if (central or ForAll(Concatenation(i[2],i[3]),
		i->ForAll(mpcgs,
		  j->DepthOfPcElement(pcgs,Comm(i,j))>=ser.depths[d])) ) then
	Info(InfoHomClass,3,"central step");
	new:=LiftClassesEANonsolvCentral(G,mpcgs,i,hom,pcisom,solvtriv,fran);
      elif Length(fants)>0 and Order(i[1])=1 then
	# special case for trivial representetive
	new:=LiftClassesEATrivRep(G,mpcgs,i,fants,hom,pcisom,solvtriv);
	if new=fail then
	  new:=LiftClassesEANonsolvGeneral(G,mpcgs,i,hom,pcisom,solvtriv,fran);
	fi;
      else
	new:=LiftClassesEANonsolvGeneral(G,mpcgs,i,hom,pcisom,solvtriv,fran);
      fi;
      Assert(3,ForAll(new,x->x[6]
        =Size(Group(Concatenation(x[2],DenominatorOfModuloPcgs(mpcgs))))));

#if ForAny(new,x->x[2]<>TFMakeInducedPcgsModulo(pcgs,x[2],nran)) then Error("HUH6");fi;
#Print(List(new,x->List(x[2],y->DepthOfPcElement(pcgs,y))),"\n");

      #Assert(1,ForAll(new, i->ForAll(GeneratorsOfGroup(i[2]),j->Comm(j,i[1]) in N)));
      ncl:=Concatenation(ncl,new);
      Info(InfoHomClass,2,Length(new)," new classes (",Length(ncl)," total)");
    od;
    cl:=ncl;
    Info(InfoHomClass,1,"Now: ",Length(cl)," classes (",Length(ncl)," total)");
  od;

  if Order(cl[1][1])>1 then
    # the identity is not in first position
    Info(InfoHomClass,2,"identity not first, sorting");
    Sort(cl,function(a,b) return Order(a[1])<Order(b[1]);end);
  fi;

  Info(InfoHomClass,1,"forming classes");
  ncl:=[];
  for i in cl do
    if IsInt(i[2]) then
      r:=ConjugacyClass(G,i[1]);
      SetSize(r,Size(G)/i[2]);
    else
      d:=SubgroupByFittingFreeData(G,i[3],i[4],i[2]);
      Assert(1,Size(d)=i[5]);
      Assert(2,Centralizer(G,i[1])=d);
      SetSize(d,i[5]);
      r:=ConjugacyClass(G,i[1],d);
      SetSize(r,Size(G)/i[5]);
    fi;
    Add(ncl,r);
  od;

  # as this test is cheap, do it always
  if Sum(ncl,Size)<>Size(G) then
    Error("wrong classes");
  fi;

  cl:=ncl;

  return cl;
end);

#############################################################################
##
#F  LiftConCandCenNonsolvGeneral( <H>,<N>,<NT>,<cl> )
##
BindGlobal("LiftConCandCenNonsolvGeneral",
  function( H, Npcgs, reps, hom, pcisom,solvtriv,fran)
    local  nreps,      # new reps to be constructed, the result
           correctingelement,
	   minvec,
	   cano,       # element that will be canonical
	   cl,
           field,      # field over which <N> is a vector space
	   one,
           h,          # preimage `cl.representative' under <hom>
	   cg,
           cNh,        # centralizer of <h> in <N>
           C,  gens,   # preimage `Centralizer( cl )' under <hom>
           r,          # dimension of <N>
           ran,        # constant range `[ 1 .. r ]'
           aff,        # <N> as affine space
           imgs,  M,   # generating matrices for affine operation
           orb,        # orbit of affine operation
           rep,# set of classes with canonical representatives
           c,  i, # loop variables
	   PPcgs,denomdepths,
	   reduce,
	   corr,
	   correctionfactor,
	   censize,cenradsize,
	   stabfac,
	   stabfacgens,
	   stabfacimgs,
	   stabrad,
	   sz,gpsz,subsz,solvsz,
	   orblock,
	   stage,
	   b,j,x,
	   minimal,minstab,mappingelm,
	   p,vp,genum,
	   st,gpe,
	   fe,epi,
	   repword,repwords,radidx,img,
	   radrange,farange,
	   ratio,
	   deno,docorrect,
	   k,failcnt,orpo,
	   stabilizergen,stabstack,
	   sel,
	   comm,s,stab;# for class correction

  correctingelement:=function(h,rep,fe)
  local comm;
    comm:=Comm( h, fe ) * Comm( rep, fe );
    comm:= ExponentsOfPcElement(Npcgs,comm)*one;
    ConvertToVectorRep(comm,field);
    comm := List(comm * cg.inverse,Int);
    comm:=PcElementByExponentsNC
      ( Npcgs, Npcgs{ cg.needed }, -comm );
    fe:=fe*comm;
    return fe;
  end;

  mappingelm:=function(orb,pos)
  local mc,mcf,i;
    mc:=One(Source(hom));
    mcf:=One(Range(hom));
    for i in orb.repwords[pos] do
      mc:=mc*orb.gens[i];
      mcf:=mcf*orb.fgens[i];
    od;
    i:=pos mod orb.orblock;
    if i=0 then i:=orb.orblock;fi;
    mc:=orb.reps[i]*mc;
    return [mc,mcf];
  end;

  # all reps given must have the same canonical representative on the level
  # above. So they also all have the same centralizer and we can use this.
  cl:=reps[1];
  h := cl[3];

  field := GF( RelativeOrders( Npcgs )[ 1 ] );
  one:=One(field);
  PPcgs:=ParentPcgs(NumeratorOfModuloPcgs(Npcgs));
  denomdepths:=ShallowCopy(DenominatorOfModuloPcgs(Npcgs)!.depthsInParent);
  Add(denomdepths,Length(PPcgs)+1); # one
  
  # Determine the subspace $[h,N]$ and calculate the centralizer of <h>.
  #cNh := ExtendedPcgs( DenominatorOfModuloPcgs( N!.capH ),
  #               VSDecompCentAction( N, h, N!.capH ) );

  #oldcg:=KernelHcommaC(Npcgs,h,NumeratorOfModuloPcgs(Npcgs),2);

  #cg:=VSDecompCentAction( Npcgs, h, NumeratorOfModuloPcgs(Npcgs),field,2 );
  cg:=VSDecompCentAction( Npcgs, h, Npcgs,field,2 );
#Print("complen =",Length(cg.baseComplement)," of ",cg.dimensionN,"\n");
#if Length(Npcgs)>5 then Error("cb"); fi;

  cNh:=cg.cNh;

  r := Length( cg.baseComplement );
  ran := [ 1 .. r ];
  
  # Construct matrices for the affine operation on $N/[h,N]$.
  Info(InfoHomClass,4,"space=",Size(field),"^",r);
  if Size(field)^r>3*10^8 then Error("too large");fi;
  aff := ExtendedVectors( field ^ r );

  # Format for cl is:
  # 1:Element, 2: Conjugate element, 3: Element that will be canonical
  # in factor, 4:conjugator, 5:cenpcgs,
  # 6:cenfac, 7:cenfacimgs, 8:censize, 9:cenfacsize

  gens:=Concatenation(cl[5],Npcgs,cl[6]); # all generators
#if ForAny(gens,x->Order(x)=1) then Error("HUH2"); fi;
  gpsz:=cl[8];

  solvsz:=cl[8]/cl[9];

  radidx:=Length(Npcgs)+Length(cl[5]);

  imgs := [  ];
  for c  in gens  do
    M := [  ];
    for i  in [ 1 .. r ]  do
	M[ i ] := Concatenation( ExponentsConjugateLayer( Npcgs,
	      Npcgs[ cg.baseComplement[ i ] ] , c )
	      * cg.projection, [ Zero( field ) ] );
    od;
    M[ r + 1 ] := Concatenation( ExponentsOfPcElement
			  ( Npcgs, Comm( h, c ) ) * cg.projection,
			  [ One( field ) ] );

    M:=ImmutableMatrix(field,M);
    Add( imgs, M );
  od;

#if Size(field)^r>10^7 then Error("BIG");fi;

  # now compute orbits, being careful to get stabilizers in steps
  #orbreps:=[];
  #stabs:=[];

  # change reps to list of more general format
  nreps:=[];
  for x in reps do
    p:=rec(list:=x);
    p.vector:=ExponentsOfPcElement(Npcgs,LeftQuotient(h,x[2]))*One(field);
    p.exponents:=ShallowCopy(p.vector);
    ConvertToVectorRep(p.vector,field);
    p.vector:=p.vector*cg.projection;
    Add(p.vector,One(field));
    Add(nreps,p);
  od;
  reps:=nreps; 

  nreps:=[];
  sel:=[1..Length(reps)];
  while Length(sel)>0 do
    p:=sel[1]; # the one to do
    RemoveSet(sel,p);
    orb:=OrbitsRepsAndStabsVectorsMultistage(gens{[1..radidx]},
	  imgs{[1..radidx]},
	  pcisom,solvsz,solvtriv,
	  gens{[radidx+1..Length(gens)]},
	  imgs{[radidx+1..Length(gens)]},reps[p].list[7],hom,gpsz,OnRight,
	  aff:orbitseed:=reps[p].vector);
    orb:=orb[1];
    # find minimal element, mapper, stabilizer of minimal element
    minvec:=Minimum(orb.orbit);
    minimal:=mappingelm(orb,Position(orb.orbit,minvec));

    # get the real minimum, including N-Orbit
    corr:=ExponentsOfPcElement(Npcgs,
	                       LeftQuotient(h,reps[p].list[2]^minimal[1]));
    corr:=PcElementByExponentsNC(Npcgs,Npcgs{cg.needed},-corr*cg.inverse);
    minimal[1]:=minimal[1]*corr; # real minimizer

    # element that will be the canonical representative in the factor (tail
    # zeroed out)
    corr:=ExponentsOfPcElement(Npcgs,
	   LeftQuotient(h,reps[p].list[2]^minimal[1]));
    cano:=h*PcElementByExponents(Npcgs,corr);

    # this will not be a pcgs, but we induce later anyhow
    stabrad:=List(orb.stabradgens,x->x^minimal[1]); 
    stabfacgens:=List(orb.stabfacgens,x->x^minimal[1]);
    stabfacimgs:=List(orb.stabfacimgs,x->x^minimal[2]);

    censize:=orb.subsz;
    cenradsize:=orb.stabrsubsz;

    # correct generators. Partially in Pc Image
    if Length(cg.needed)>0 then

      rep:=LeftQuotient(h,reps[p].list[2]^minimal[1]);

      stabrad:=List(stabrad,x->correctingelement(h,rep,x));
      # must make proper pcgs -- correction does not preserve that
      stabrad:=TFMakeInducedPcgsModulo(PPcgs,stabrad,denomdepths);

      # we change by radical elements, so the images in the factor don't
      # change
      stabfacgens:=List(stabfacgens,x->correctingelement(h,rep,x));

      correctionfactor:=Characteristic(field)^Length(cg.needed);
      censize:=censize/correctionfactor;
      cenradsize:=cenradsize/correctionfactor;
    fi;

    # Format for cl is:
    # 1:Element, 2: Conjugate element, 3: Element that will be canonical
    # in factor, 4:conjugator, 5:cenpcgs,
    # 6:cenfac, 7:cenfacimgs, 8:censize, 9:cenfacsize
 
    Add(nreps,[reps[p].list[1],
               reps[p].list[2]^minimal[1],
	       cano,
	       reps[p].list[4]*minimal[1],
	       stabrad,stabfacgens,stabfacimgs,
	       censize,censize/cenradsize]);

    for cl in sel do
      b:=Position(orb.orbit,reps[cl].vector);
      if b<>fail then
	RemoveSet(sel,cl);
	# now find rep mapping 1 here
	b:=mappingelm(orb,b);
	b:=[LeftQuotient(b[1],minimal[1]),
			      LeftQuotient(b[2],minimal[2])];

	# get the real minimum, including N-Orbit
	corr:=ExponentsOfPcElement(Npcgs,
				  LeftQuotient(h,reps[cl].list[2]^b[1]));
	corr:=PcElementByExponentsNC(Npcgs,Npcgs{cg.needed},corr*cg.inverse);
	b[1]:=b[1]/corr; # real minimizer

	# Format for cl is:
	# 1:Element, 2: Conjugate element, 3: Element that will be canonical
	# in factor, 4:conjugator, 5:cenpcgs,
	# 6:cenfac, 7:cenfacimgs, 8:censize, 9:cenfacsiz
	Add(nreps,[reps[cl].list[1],
		  reps[cl].list[2]^b[1],
		  cano,
		  reps[cl].list[4]*b[1],
		  stabrad,stabfacgens,stabfacimgs,
		  censize,censize/cenradsize]);

      elif ValueOption("conjugacytest")=true then
	# in conj test this would mean fail
	return fail;
      fi;
    od;
  od;

  return nreps;

end);

# canonical rep/centralizer
BindGlobal("TFCanonicalClassRepresentative",function (G,candidates)
local r,	#radical
      f,	# G/r
      hom,	# G->f
      prereps,  # fixed factor class reps preimages
      pcgs,mpcgs, #(modulo) pcgs
      pcisom,
      gens,
      ser,	# series
      radsize,len,ntrihom,
      mran,nran,fran,
      central,
      #fants,
      reps,
      nreps,
      fr,
      conj,
      d,
      solvtriv,
      select,sel,pos,
      M,N,	# normal subgrops
      ind,	# indices
      i,j,	#loop
      new,	# new classes
      classrange,
      cl,ncl;	# classes

  # it seems to be cleaner (and avoids deferring abelian factors) if we
  # factor out the radical first. (Note: The radical method for perm groups
  # stores the nat hom.!)
  ser:=FittingFreeLiftSetup(G);
  radsize:=Product(RelativeOrders(ser.pcgs));
  len:=Length(ser.pcgs);

  pcgs:=ser.pcgs;
  pcisom:=ser.pcisom;
  #fants:=fail;

  # store centralizers as rep, centralizer generators in radical,
  # centralizer generators with nontrivial
  # radfactor image, corresponding radfactor images
  # the generators in the radical do not list the generators of the
  # current layer after immediate lifting.

  if radsize=Size(G) then
    # solvable case
    hom:=ser.factorhom;
    d:=MappingGeneratorsImages(hom);
    mran:=Filtered([1..Length(d[2])],x->not IsOne(d[2][x]));
    # elm, elmconj, canonical factor, conjugator, cenpcgs,cenfac,cenfacimgs,censize,cenfacsiz
    reps:=List(candidates,x->[x,x,One(G),One(G),[],d[1]{mran},d[2]{mran},Size(G),Size(G)]);
  else
    # nonsolvable
    if radsize>1 then
      hom:=ser.factorhom;
      ntrihom:=true;
      f:=Image(hom);
      # we need centralizers
      #fants:=Filtered(NormalSubgroups(f),x->Size(x)>1 and Size(x)<Size(f));
    else
      if IsPermGroup(G) then
	hom:=SmallerDegreePermutationRepresentation(G);
	ntrihom:=not IsOne(hom);;
      else
        hom:=IdentityMapping(G);
	ntrihom:=false;
      fi;
      f:=Image(hom);
    fi;
    if not IsBound(f!.someClassReps) then
      f!.someClassReps:=[ConjugacyClass(f,One(f))]; # identity first
    fi;
    if HasConjugacyClasses(f) and
      Length(f!.someClassReps)<Length(ConjugacyClasses(f)) then
      # expand the list of stored factor classes for once.
      cl:=Filtered(ConjugacyClasses(f),x-> not ForAny(f!.someClassReps,
	   y->Order(Representative(x))=Order(Representative(y))
	   and Representative(y) in x));
      Append(f!.someClassReps,cl);
    fi;
    cl:=f!.someClassReps;

    classrange:=[1..Length(cl)];
    r:=ValueOption("candidatenums");
    if r<>fail and HasConjugacyClasses(G) then
      # candidatenums gives the numbers of some classes in G that should be
      # tried first (as they likely contain the element). Us this to reduce
      # conjugacy test in factor.
      if not IsBound(G!.radicalfactorclassmap) then
	G!.radicalfactorclassmap:=[];
      fi;
      fr:=G!.radicalfactorclassmap;
      for i in Filtered(r,x->not IsBound(fr[x])) do
	# compute images that are not yet known
	d:=ImagesRepresentative(hom,Representative(ConjugacyClasses(G)[i]));
	j:=First([1..Length(cl)],x->IsBound(cl[x]) and d in cl[x]);
	if j<>fail then
	  fr[i]:=j;
	fi;
      od;
      r:=Filtered(r,x->IsBound(fr[x]));
      r:=Set(fr{r}); # class numbers in factor
      classrange:=Concatenation(r,
		   Filtered(classrange,x->not x in r));

    fi;

    nreps:=[];
    for i in candidates do
      fr:=ImagesRepresentative(hom,i);
      conj:=fail;
      j:=0;
      while conj=fail and j<Length(classrange) do
	j:=j+1;
        if Order(fr)=Order(Representative(cl[classrange[j]])) then
	  conj:=RepresentativeAction(f,fr,Representative(cl[classrange[j]]));
	fi;
      od;

      if conj=fail then
	# not yet found, and classes of f were not known -- store this rep
	# image as canonical one for future use.
	j:=j+1;
	Add(cl,ConjugacyClass(f,fr));
	conj:=One(f);
      else
	j:=classrange[j];
      fi;
      # store fixed preimages of reps to avoid any impact of homomorphism.
      if not IsBound(f!.classpreimgs) then
	f!.classpreimgs:=[];
      fi;
      prereps:=f!.classpreimgs;
      if not IsBound(prereps[j]) then
	prereps[j]:=PreImagesRepresentative(hom,Representative(cl[j]));
      fi;

      r:=PreImagesRepresentative(hom,conj);
      d:=GeneratorsOfGroup(Centralizer(cl[j]));
      # Format for cl is:
      # 1:Element, 2: Conjugate element, 3: Element that will be canonical
      # in factor, 4:conjugator, 5:cenpcgs,
      # 6:cenfac, 7:cenfacimgs, 8:censize, 9:cenfacsize
      Add(nreps,[i,i^r,prereps[j],r,[],
        List(d,x->PreImagesRepresentative(hom,x)),d,
	radsize*Size(Centralizer(cl[j])), Size(Centralizer(cl[j]))]);
    od;
    reps:=nreps;

  fi;

  for d in [2..Length(ser.depths)] do
    #M:=ser[i-1];
    #N:=ser[i];
    mran:=[ser.depths[d-1]..len];
    nran:=[ser.depths[d]..len];
    fran:=[mran,nran];

    mpcgs:=InducedPcgsByPcSequenceNC(pcgs,pcgs{mran}) mod
           InducedPcgsByPcSequenceNC(pcgs,pcgs{nran});

    central:= ForAll(GeneratorsOfGroup(G),
		i->ForAll(mpcgs,
		  j->DepthOfPcElement(pcgs,Comm(i,j))>=ser.depths[d]));
    
    # abelian factor, use affine methods
    Info(InfoHomClass,1,"abelian factor ",d,": ",
      Product(RelativeOrders(ser.pcgs){mran}), "->",
      Product(RelativeOrders(ser.pcgs){nran})," central:",central);

    nreps:=[];
    solvtriv:=Subgroup(Range(pcisom),
	List(DenominatorOfModuloPcgs(mpcgs),x->ImagesRepresentative(pcisom,x)));
    select:=[1..Length(reps)];
    while Length(select)>0 do
      pos:=select[1];
      sel:=Filtered(select,x->reps[x][3]=reps[pos][3]);
      if ValueOption("conjugacytest")=true and Length(sel)<>2 then
	return fail;
      fi;
      Info(InfoHomClass,2,Length(sel),"in candidate group");
      select:=Difference(select,sel);
      new:=LiftConCandCenNonsolvGeneral(G,mpcgs,reps{sel},hom,pcisom,solvtriv,fran);
      # conj test 
      if new=fail then
	return new;
      fi;
      Append(nreps,new);
    od;
    reps:=nreps;
  od;

  # arrange back to same ordering as before.
  nreps:=[];
  for i in candidates do
    Add(nreps,First(reps,x->IsIdenticalObj(x[1],i)));
  od;

  return nreps;
end);

#############################################################################
##
#M  Centralizer( <G>, <e> ) . . . . . . . . . . . . . . using TF method
##
InstallMethod( CentralizerOp, "TF method:elm",IsCollsElms,
  [ IsGroup and IsFinite and HasFittingFreeLiftSetup,
  IsMultiplicativeElementWithInverse ], OVERRIDENICE,
function( G, e )
  if IsPermGroup(G) or IsPcGroup(G) then TryNextMethod();fi;
  if not e in G then 
    # TODO: form larger group containing e.
    TryNextMethod();
  fi;
  e:=TFCanonicalClassRepresentative(G,[e])[1];
  if e=fail then TryNextMethod();fi;
  return SubgroupByFittingFreeData(G,e[6],e[7],e[5]);
end );

#############################################################################
##
#M  Centralizer( <G>, <e> ) . . . . . . . . . . . . . . using TF method
##
InstallMethod( CentralizerOp, "TF method:subgroup",IsIdenticalObj,
  [ IsGroup and IsFinite and HasFittingFreeLiftSetup,
  IsGroup and IsFinite and HasGeneratorsOfGroup], 
  2*OVERRIDENICE,
function( G, S )
local c,e;
  if IsPermGroup(G) or IsPcGroup(G) then TryNextMethod();fi;
  c:=G;
  for e in GeneratorsOfGroup(S) do
    c:=Centralizer(c,e);
  od;
  return c;
end );

#############################################################################
##
#M  RepresentativeAction( <G>, <d>, <e>, <act> ) . . . . . using TF method
##
InstallOtherMethod( RepresentativeActionOp, "TF Method on elements",
  IsCollsElmsElmsX, 
  [ IsGroup and IsFinite and HasFittingFreeLiftSetup,
        IsMultiplicativeElementWithInverse,
	IsMultiplicativeElementWithInverse, IsFunction ], 
  OVERRIDENICE,
function ( G, d, e, act )
local c;
  if IsPermGroup(G) or IsPcGroup(G) then TryNextMethod();fi;
  if not (d in G and e in G) then 
    # TODO: form larger group containing e.
    TryNextMethod();
  fi;

  if act=OnPoints then #and d in G and e in G then
    c:=TFCanonicalClassRepresentative(G,[d,e]:conjugacytest);
    if c=fail then
      return fail;
    else
      if c[1][2]=c[2][2] then
	return c[1][4]/c[2][4]; # map via canonicals
      else
	return fail; # not conjugate
      fi;
    fi;
  fi;
  TryNextMethod();
end);

#############################################################################
##
#M  ConjugacyClasses( <G> ) . . . . . . . . . . . . . . . . . . of perm group
##
InstallMethod( ConjugacyClasses, "perm group", true, 
  [ IsPermGroup and IsFinite],OVERRIDENICE,
function( G )
local cl;
  if IsNaturalSymmetricGroup(G) or IsNaturalAlternatingGroup(G) then
    # there are better methods for Sn/An
    TryNextMethod();
  fi;

  cl:=ConjugacyClassesForSmallGroup(G);
  if cl<>fail then
    return cl;
  elif IsSimpleGroup( G ) then
    cl:=ClassesFromClassical(G);
    if cl=fail then
      cl:=ConjugacyClassesByRandomSearch( G );
    fi;
    return cl;
  else
    return ConjugacyClassesViaRadical(G);
  fi;
end );

#############################################################################
##
#M  ConjugacyClasses( <G> ) . . . . . . . . . . . . . . . . . . of perm group
##
InstallMethod( ConjugacyClasses, "TF Method", true,
  [ IsGroup and IsFinite and HasFittingFreeLiftSetup],OVERRIDENICE,
function(G)
  if IsPermGroup(G) or IsPcGroup(G) then TryNextMethod();fi;
  return ConjugacyClassesViaRadical(G);
end);

#############################################################################
#
# old code, kept for comparison

#############################################################################
##
#F  GeneralStepClEANSNonsolv( <H>,<N>,<NT>,<cl> )
##
BindGlobal("GeneralStepClEANSNonsolv",function( H, N,NT, cl )
    local  classes,    # classes to be constructed, the result
           field,      # field over which <N> is a vector space
	   one,
           h,          # preimage `cl.representative' under <hom>
           cNh,        # centralizer of <h> in <N>
           C,  gens,   # preimage `Centralizer( cl )' under <hom>
           r,          # dimension of <N>
           ran,        # constant range `[ 1 .. r ]'
           aff,        # <N> as affine space
           xset,       # affine operation of <C> on <aff>
	   eo,	       # xorbits/stabilizers
           imgs,  M,   # generating matrices for affine operation
           orb,        # orbit of affine operation
           rep,# set of classes with canonical representatives
           c,  i, # loop variables
	   reduce,
	   stabsub,
	   comm,s,stab;# for class correction
	  
  #NT:=AsSubgroup(H,NT);
  C := cl[2];
  field := GF( RelativeOrders( N )[ 1 ] );
  one:=One(field);
  h := cl[1];
  reduce:=ReducedPermdegree(C);
  if reduce<>fail then
    C:=Image(reduce,C);
    Info(InfoHomClass,4,"reduced to deg:",NrMovedPoints(C));
    h:=Image(reduce,h);
    NT:=Image(reduce,NT);
    N:=ModuloPcgs(SubgroupNC(C,Image(reduce,NumeratorOfModuloPcgs(N))),NT);
  fi;
  
  
  # Determine the subspace $[h,N]$ and calculate the centralizer of <h>.
  #cNh := ExtendedPcgs( DenominatorOfModuloPcgs( N!.capH ),
  #               KernelHcommaC( N, h, N!.capH ) );
  N!.capH:=N;
  cNh:=KernelHcommaC( N, h, N!.capH,2 );

  r := Length( N!.subspace.baseComplement );
  ran := [ 1 .. r ];
  
  # Construct matrices for the affine operation on $N/[h,N]$.
  aff := ExtendedVectors( field ^ r );

  if IsSolvableGroup(C) then
    gens:=Pcgs(C);
  else
    gens:=GeneratorsOfGroup(C);
  fi;
  imgs := [  ];
  for c  in gens  do
    M := [  ];
    for i  in [ 1 .. r ]  do
	M[ i ] := Concatenation( ExponentsConjugateLayer( N,
	      N[ N!.subspace.baseComplement[ i ] ] , c )
	      * N!.subspace.projection, [ Zero( field ) ] );
    od;
    M[ r + 1 ] := Concatenation( ExponentsOfPcElement
			  ( N, Comm( h, c ) ) * N!.subspace.projection,
			  [ One( field ) ] );

    M:=ImmutableMatrix(field,M);
    Add( imgs, M );
  od;
  xset := ExternalSet( C, aff, gens, imgs );

  classes := [  ];

  # NC is safe
  stabsub:=ClosureSubgroupNC(NT,cNh);
  SetActionKernelExternalSet(xset,stabsub);
  eo:=ExternalOrbitsStabilizers( xset );

  for orb in eo do
    rep := PcElementByExponentsNC( N, N{ N!.subspace.baseComplement },
		    Representative( orb ){ ran } );
    Assert(2,ForAll(GeneratorsOfGroup(stabsub),i->Comm(i,h*rep) in NT));

    # filter those we don't get anyhow.
    stab:=Filtered(GeneratorsOfGroup(StabilizerOfExternalSet(orb)),
                   i->not i in stabsub);

    comm := [  ];
    for s  in [ 1 .. Length( stab ) ]  do
	comm[ s ] := ExponentsOfPcElement( N,
	    Comm( rep, stab[ s ] ) * Comm( h, stab[ s ] ) )*one;
    od;
    comm:=ImmutableMatrix(field,comm);
    comm := comm * N!.subspace.inverse;
    for s  in [ 1 .. Length( comm ) ]  do
      stab[ s ] := stab[ s ] / PcElementByExponentsNC
	  ( N, N{ N!.subspace.needed }, comm[ s ] );
	  #( N!.capH, N!.capH{ N!.subspace.needed }, comm[ s ] );
      Assert(2,Comm(h*rep,stab[s]) in NT);
    od;

  # NC is safe
    stab:=ClosureSubgroupNC(stabsub,stab);

    if IsSolvableGroup(C) then
      SetIsSolvableGroup(stab,true);
    fi;
    c := [h * rep,stab];
    Assert(2,ForAll(GeneratorsOfGroup(stab),i->Comm(i,c[1]) in NT));

    if reduce<>fail then
      Add(classes,[PreImagesRepresentative(reduce,c[1]),
	PreImage(reduce,c[2])]);
    else
      Add(classes,c);
    fi;
  od;

  Assert(1,ForAll(classes,i->i[1] in H and IsSubset(H,i[2])));
  return classes;

end);

# new version, no subspace
#############################################################################
##
#F  GeneralStepCanEANSNonsolv( <H>,<N>,<NT>,<C>,<reps> )
##
## canonical rep
BindGlobal("GeneralStepCanEANSNonsolv",function( H, N,NT, C,h,reps,repo,nostab )
  local SchreierVectorProduct, field, one, r, ran, gens, imgs, M, invimgs, repvec, repgps, newreps, aff, sel, i, repsofi, orb, rep, dict, q, stab, sti, stabgens, p, img, mi, os, a, mipo, mimap, map, ngrp, c, j;

  SchreierVectorProduct:=function(n)
  local w,q,a;
    w:=One(C);
    while n<>1 do
      q:=rep[n];
      w:=gens[q]*w;
      n:=LookupDictionary(dict,orb[n]*invimgs[q]);
    od;
    return w;
  end;

  #NT:=AsSubgroup(H,NT);
  field := GF( RelativeOrders( N )[ 1 ] );
  one:=One(field);
  #reduce:=ReducedPermdegree(C);
  #if reduce<>fail then
  #  C:=Image(reduce,C);
  #  Info(InfoHomClass,4,"reduced to deg:",NrMovedPoints(C));
  #  h:=Image(reduce,h);
  #  NT:=Image(reduce,NT);
  #  N:=ModuloPcgs(SubgroupNC(C,Image(reduce,NumeratorOfModuloPcgs(N))),NT);
  #fi;

  r := Length(N);
  ran := [ 1 .. r ];
  
  # Construct matrices for the affine operation on $N/[h,N]$.
  gens:=Filtered(GeneratorsOfGroup(C),i->not i in NT);
  if Length(gens)>20 then
    gens:=Filtered(SmallGeneratingSet(C),i->not i in NT);
  fi;
  imgs := [  ];
  for c  in gens  do
    M := [  ];
    for i  in ran  do
      #M[i]:=Concatenation(ExponentsConjugateLayer(N,N[i],c)*one,[Zero(field)]);
      M[i]:=Concatenation(ExponentsOfPcElement(N,N[i]^c)*one,[Zero(field)]);
    od;
    M[r+1]:=Concatenation(ExponentsOfPcElement(N,Comm(h,c))*one,[One(field)]);

    M:=ImmutableMatrix(field,M);
    Add( imgs, M );
  od;
  invimgs:=List(imgs,Inverse);

  # get vectors for reps
  repvec:=List(repo,i->Concatenation(
	   ExponentsOfPcElement(N,LeftQuotient(h,reps[i][1]))*one,[one]));
  for i in repvec do
    ConvertToVectorRep(i,field);
  od;
  repgps:=[];
  newreps:=[];
  aff:=field^(r+1);
  sel:=[1..Length(repo)];
  while Length(sel)>0 do
    i:=sel[1];
    repsofi:=reps[repo[i]];
    RemoveSet(sel,i);
    # since we want representatives as well, recode the orbit algorithm.
    orb:=[repvec[i]];
    rep:=[0];
    dict:=NewDictionary(repvec[i],true,aff);
    AddDictionary(dict,repvec[i],1);

    # get stabilizing generators
    q:=gens{Filtered([1..Length(gens)],i->orb[1]*imgs[i]=orb[1])};
    if q=gens or nostab then 
      stab:=C;
    else
      stab:=ClosureGroup(NT,q);
    fi;
    sti:=5;
    if nostab then sti:=-1;fi;
    stabgens:=[];
    p:=1;
    while p<=Length(orb) do
      for j in [1..Length(gens)] do
	img:=orb[p]*imgs[j];
	q:=LookupDictionary(dict,img);
	if q=fail then
	  Add(orb,img);
	  AddDictionary(dict,img,Length(orb));
	  Add(rep,j);
	elif Size(C)/Size(stab)>Length(orb) then
	  if sti=0 then
	    Add(stabgens,[p,j,q]);
	    if Random([1..QuoInt(Length(orb),5)])=1 then
	      os:=Random([1..Length(stabgens)]);
	      mi:=stabgens[os];
	      stabgens[os]:=stabgens[Length(stabgens)];
	      Unbind(stabgens[Length(stabgens)]);
	      os:=Size(stab);
	      stab:=ClosureGroup(stab,SchreierVectorProduct(mi[1])*gens[mi[2]]
					/ SchreierVectorProduct(mi[3]));
              if Size(stab)>os then
		sti:=1;
	      fi;
	    fi;
	  else
	    os:=Size(stab);
	    stab:=ClosureGroup(stab,SchreierVectorProduct(p)*gens[j]
				      / SchreierVectorProduct(q));
            if Size(stab)=os then
	      sti:=sti-1;
	    fi;
	  fi;
	fi;
      od;
      p:=p+1;
    od;
    # add missing schreier gens
    a:=Size(C)/Length(orb);
    while Size(stab)<a and not nostab do
      os:=Random([1..Length(stabgens)]);
      mi:=stabgens[os];
      stabgens[os]:=stabgens[Length(stabgens)];
      Unbind(stabgens[Length(stabgens)]);
      stab:=ClosureGroup(stab,SchreierVectorProduct(mi[1])*gens[mi[2]]
				/ SchreierVectorProduct(mi[3]));
    od;

    Info(InfoHomClass,3,"Orbit length ",Length(orb),
        " with ",Length(gens)," generators");
    mi:=Minimum(orb); # the ``canonical'' rep.
    mipo:=LookupDictionary(dict,mi);
    mimap:=SchreierVectorProduct(mipo); # element moving starter to minimal
    map:=mimap;
    stab:=stab^map; # stabilize minimal element
    mi:=PcElementByExponentsNC(N,mi{ran});

    Assert(1,ForAll(GeneratorsOfGroup(stab),x->Comm(x,h*mi) in NT));

    ngrp:=[[repo[i]],h*mi,stab];
    Add(repgps,ngrp);
    newreps[repo[i]]:=[repsofi[1]^map,repsofi[2]*map,Length(repgps)];
    Assert(1,LeftQuotient(h*mi*One(NT),repsofi[1]^map) in NT);

    for i in ShallowCopy(sel) do
      q:=LookupDictionary(dict,repvec[i]);
      if q<>fail then
	RemoveSet(sel,i);
	repsofi:=reps[repo[i]];
	Add(ngrp[1],repo[i]);
	map:=LeftQuotient(SchreierVectorProduct(q),mimap);
	newreps[repo[i]]:=[repsofi[1]^map,repsofi[2]*map,Length(repgps)];
	Assert(1,LeftQuotient(h*mi,repsofi[1]^map) in NT);
      fi;
    od;

  od;

  return [repgps,newreps];

end);

  
#############################################################################
##
#F  CentralStepClEANSNonsolv( <H>, <N>, <cl> )
##
# the version for pc groups implicitly uses a pc-group orbit-stabilizer
# algorithm. We can't  do this but have to use a more simple-minded
# orbit/stabilizer approach.
BindGlobal("CentralStepClEANSNonsolv",function( H, N, cl )
local  classes,    	# classes to be constructed, the result
	f,          	# field over which <N> is a vector space
	o,
	n,r,		# dimensions
	space,
	com,
	comms,
	mats,
	decomp,
	reduce,
	v,
	h,          	# preimage `cl.representative' under <hom>
	C,		# preimage `Centralizer( cl )' under <hom>
	w,    		# coefficient vectors for projection along $[h,N]$
	c;          	# loop variable

  C:=cl[2];
  h := cl[1];
  reduce:=ReducedPermdegree(C);
  if reduce<>fail then
    C:=Image(reduce,C);
    Info(InfoHomClass,4,"reduced to deg:",NrMovedPoints(C));
    h:=Image(reduce,h);
    N:=ModuloPcgs(SubgroupNC(C,Image(reduce,NumeratorOfModuloPcgs(N))),
	          SubgroupNC(C,Image(reduce,DenominatorOfModuloPcgs(N))));
  fi;

  # centrality still means that conjugation by c is multiplication with
  # [h,c] and that the complement space is generated by commutators [h,c]
  # for a generating set {c|...} of C.

  f:=GF(RelativeOrders(N)[1]);
  n:=Length(N);
  o:=One(f);
  # commutator space basis
  comms:=List(GeneratorsOfGroup(C),c->o*ExponentsOfPcElement(N,Comm(h,c)));
  List(comms,x->ConvertToVectorRep(x,f));
  space:=List(comms,ShallowCopy);
  TriangulizeMat(space);
  space:=Filtered(space,i->i<>Zero(i)); # remove spurious columns

  com:=BaseSteinitzVectors(IdentityMat(n,f),space);

  # decomposition of vectors into the subspace basis
  r:=Length(com.subspace);
  if r>0 then
    # if the subspace is trivial, everything stabilizes

    decomp:=Concatenation(com.subspace,com.factorspace)^-1;
    decomp:=decomp{[1..Length(decomp)]}{[1..r]};
    decomp:=ImmutableMatrix(f,decomp);

    # build matrices for the affine action
    mats:=[];
    for w in comms do
      c:=IdentityMat(r+1,o);
      c[r+1]{[1..r]}:=w*decomp; # translation bit
      c:=ImmutableMatrix(f,c);
      Add(mats,c);
    od;

    #subspace affine enumerator
    v:=ExtendedVectors(f^r);

    C := Stabilizer( C, v, v[1],GeneratorsOfGroup(C), mats,OnPoints );
  fi;
  Assert(1,Size(cl[2])/Size(C)=Size(f)^r);

  if Length(com.factorspace)=0 then
    if reduce<>fail then
      classes:=[[PreImagesRepresentative(reduce,h),PreImage(reduce,C)]];
    else
      classes:=[[h,C]];
    fi;
  else
    classes:=[];
    # enumerator of complement
    v:=f^Length(com.factorspace);
    for w in v do
      c := [h * PcElementByExponentsNC( N,w*com.factorspace),C ];
      if reduce<>fail then
	Add(classes,[PreImagesRepresentative(reduce,c[1]),
	  PreImage(reduce,c[2])]);
      else
	Add(classes,c);
      fi;
    od;
  fi;

  Assert(1,ForAll(classes,i->i[1] in H and IsSubset(H,i[2])));
  return classes;
end);

BindGlobal("ConjugacyClassesViaRadical_Old",function (G)
local r,	#radical
      f,	# G/r
      hom,	# G->f
      pcgs,mpcgs, #(modulo) pcgs
      ser,	# series
      M,N,	# normal subgrops
      ind,	# indices
      i,	#loop
      new,	# new classes
      cl,ncl;	# classes

  # it seems to be cleaner (and avoids deferring abelian factors) if we
  # factor out the radical first. (Note: The radical method for perm groups
  # stores the nat hom.!)
  ser:=PermliftSeries(G);
  pcgs:=ser[2];
  ser:=ser[1];
  r:=ser[1];

  if Size(r)<Size(G) then
    if Size(r)>1 then
      hom:=NaturalHomomorphismByNormalSubgroupNC(G,r);
      f:=Image(hom);
      # we need centralizers
      cl:=ConjugacyClassesFittingFreeGroup(f:onlysizes:=false);
    else
      hom:=SmallerDegreePermutationRepresentation(G); # old code
      f:=Image(hom);
      cl:=ConjugacyClassesFittingFreeGroup(f);
    fi;

    if not IsOne(hom) then
      ncl:=[];
      for i in cl do
	new:=[PreImagesRepresentative(hom,i[1])];
	if not IsInt(i[2]) then
	  Add(new,PreImage(hom,i[2]));
	fi;
        Add(ncl,new);
      od;
      cl:=ncl;
    fi;
  else
    cl:=[[One(G),G]];
  fi;

  for i in [2..Length(ser)] do
    M:=ser[i-1];
    N:=ser[i];
    
    # abelian factor, use affine methods
    Info(InfoHomClass,1,"abelian factor: ",Size(M),"->",Size(N));
    if pcgs=false then
      mpcgs:=ModuloPcgs(M,N);
    else
      mpcgs:=pcgs[i-1] mod pcgs[i];
    fi;

    ncl:=[];
    for i in cl do
      Assert(2,ForAll(GeneratorsOfGroup(i[2]),j->Comm(i[1],j) in M));
      if ForAll(GeneratorsOfGroup(i[2]),
		i->ForAll(mpcgs,j->Comm(i,j) in N)) then
	Info(InfoHomClass,3,"central step");
	new:=CentralStepClEANSNonsolv(G,mpcgs,i);
      else
	new:=GeneralStepClEANSNonsolv(G,mpcgs,AsSubgroup(G,N),i);
      fi;
      Assert(1,ForAll(new,
                  i->ForAll(GeneratorsOfGroup(i[2]),j->Comm(j,i[1]) in N)));
      Info(InfoHomClass,2,Length(new)," new classes");
      ncl:=Concatenation(ncl,new);
    od;
    cl:=ncl;
    Info(InfoHomClass,1,"Now: ",Length(cl)," classes");
  od;

  if Order(cl[1][1])>1 then
    # the idenity is not in first position
    Info(InfoHomClass,2,"identity not first, sorting");
    Sort(cl,function(a,b) return Order(a[1])<Order(b[1]);end);
  fi;

  Info(InfoHomClass,1,"forming classes");
  ncl:=[];
  for i in cl do
    if IsInt(i[2]) then
      r:=ConjugacyClass(G,i[1]);
      SetSize(r,Size(G)/i[2]);
    else
      Assert(1,Centralizer(G,i[1])=i[2]);
      r:=ConjugacyClass(G,i[1],i[2]);
    fi;
    Add(ncl,r);
  od;

  cl:=ncl;

  # temporary fix for wrong centralizers -- this code will go away anyhow
  # in next release
  if Sum(ncl,Size)<>Size(G) then
    ncl:=List(ncl,x->ConjugacyClass(G,Representative(x)));
    if Sum(ncl,Size)<>Size(G) then
      Error("wrong classes");
    fi;
  fi;

  return ncl;
end);





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