#############################################################################
##
#W  grppcfp.gi                  GAP library                      Bettina Eick
##
#Y  Copyright (C)  1997,  Lehrstuhl D fuer Mathematik,  RWTH Aachen, Germany
#Y  (C) 1998 School Math and Comp. Sci., University of St.  Andrews, Scotland
#Y  Copyright (C) 2002 The GAP Group
##
##  This file contains some functions to convert a pc group into an
##  fp group and vice versa.
##
Revision.grppcfp_gi :=
    "@(#)$Id: grppcfp.gi,v 4.44.2.3 2005/11/30 03:05:31 gap Exp $";

#############################################################################
##
#F  PcGroupFpGroup( F )
#F  PcGroupFpGroupNC( F )
##
InstallGlobalFunction( PcGroupFpGroup, function( F )
    return PolycyclicFactorGroup(
        FreeGroupOfFpGroup( F ),
        RelatorsOfFpGroup( F ) );
end );

InstallGlobalFunction( PcGroupFpGroupNC, function( F )
    return PolycyclicFactorGroupNC(
        FreeGroupOfFpGroup( F ),
        RelatorsOfFpGroup( F ) );
end );

#############################################################################
##
#F  IsomorphismFpGroupByPcgs( pcgs, str )
##
InstallGlobalFunction( IsomorphismFpGroupByPcgs, function( pcgs, str )
    local n, F, gens, rels, i, pis, exp, t, h, rel, comm, j, H, phi;

    n:=Length(pcgs);
    F    := FreeGroup( n, str );
    if n=0 then
      phi:=GroupHomomorphismByImagesNC(GroupOfPcgs(pcgs),F/[],[],[]);
      SetIsBijective( phi, true );
      return phi;
    fi;
    gens := GeneratorsOfGroup( F );
    pis  := RelativeOrders( pcgs );
    rels := [ ];
    for i in [1..n] do

        # the power
        exp := ExponentsOfRelativePower( pcgs, i ){[i+1..n]};
        t   := One( F );
        for h in [i+1..n] do
            t := t * gens[h]^exp[h-i];
        od;
        rel := gens[i]^pis[i] / t;
        Add( rels, rel );

        # the commutators
        for j in [i+1..n] do
            comm := Comm( pcgs[j], pcgs[i] );
            exp := ExponentsOfPcElement( pcgs, comm ){[i+1..n]};
            t   := One( F );
            for h in [i+1..n] do
                t := t * gens[h]^exp[h-i];
            od;
            rel := Comm( gens[j], gens[i] ) / t;
            Add( rels, rel );
        od;
    od;
    H := F / rels;
    SetSize(H,Product(RelativeOrders(pcgs)));
    phi := 
      GroupHomomorphismByImagesNC( GroupOfPcgs(pcgs), H, AsList( pcgs ),
                                        GeneratorsOfGroup( H ) );

    SetIsBijective( phi, true );
    return phi;
    
end );

#############################################################################
##
#M  IsomorphismFpGroupByCompositionSeries( G, str )
##
InstallOtherMethod( IsomorphismFpGroupByCompositionSeries, "pc groups",
               true, [IsGroup and CanEasilyComputePcgs,IsString], 0,
function( G,nam )
  return IsomorphismFpGroupByPcgs( Pcgs(G), nam );
end);

#############################################################################
##
#O  IsomorphismFpGroup( G )
##
InstallOtherMethod( IsomorphismFpGroup, "pc groups",
               true, [IsGroup and CanEasilyComputePcgs,IsString], 0,
function( G,nam )
  return IsomorphismFpGroupByPcgs( Pcgs( G ), nam);
end );

#############################################################################
##
#O  IsomorphismFpGroupByGeneratorsNC( G )
##
InstallMethod(IsomorphismFpGroupByGeneratorsNC,"pcgs",
  IsFamFamX,[IsGroup,IsPcgs,IsString],0,
function( G,p,nam )
  # this test now is obsolete but extremely cheap.
  if Product(RelativeOrders(p))<Size(G) then
    Error("pcgs does not generate the group");
  fi;
  return IsomorphismFpGroupByPcgs( p, nam);
end );

#############################################################################
##
#F  InitEpimorphismSQ( F )
##
InstallGlobalFunction( InitEpimorphismSQ, function( F )
local g, gens, r, rels, ng, nr, pf, pn, pp, D, P, M, Q, I, A, G, min,
gensA, relsA, gensG, imgs, prei, i, j, k, l, norm, index, diag, n,genu;

  if IsFpGroup(F) then
    gens := GeneratorsOfGroup( FreeGroupOfFpGroup( F ) );
    ng   := Length( gens );
    genu:=List(gens,i->GeneratorSyllable(i,1));
    genu:=List([1..Maximum(genu)],i->Position(genu,i));
    rels := RelatorsOfFpGroup( F );
    nr   := Length( rels );

	# build the relation matrix for the commutator  quotient  group
	M := [];
	for i in [ 1..Maximum( nr, ng ) ] do
	    M[i] := List( [ 1..ng ], i->0 );
	    if i <= nr then
	      r := rels[i];
	      for j in [1..NrSyllables(r)] do
		g := GeneratorSyllable(r,j);
		k:=genu[g];
		M[i][k] := M[i][k] + ExponentSyllable(r,j);
	      od;
	    fi;
	od;

    # compute normal form
    norm := NormalFormIntMat( M,15 );
    D := norm.normal; 
    P := norm.rowtrans; 
    Q := norm.coltrans; 
    I := Q^-1;
    min := Minimum( Length(D), Length(D[1]) );
    diag := List( [1..min], x -> D[x][x] );
    if ForAny( diag, x -> x = 0 ) then
        Info(InfoSQ,1,"solvable quotient is infinite");
        return false;
    fi;

    # compute pc presentation for the finite quotient
    n := Filtered( diag, x -> x <> 1 );
    n := Length( Flat( List( n, x -> FactorsInt( x ) ) ) );
    A := FreeGroup(IsSyllableWordsFamily, n );
	gensA := GeneratorsOfGroup( A );

	index := [];
    relsA := []; 
	g := 1;	
    pf := [];
	for i in [ 1..ng ] do
	    if D[i][i] <> 1 then
	        index[i] := g;
	        pf[i] := TransposedMat( Collected( FactorsInt( D[i][i] ) ) );
            pf[i] := rec( factors := pf[i][1],
                          powers  := pf[i][2] );
	        for j in [ 1..Length( pf[i].factors ) ] do
		        pn := pf[i].factors[j];
          	    pp := pf[i].powers [j];
		        for k in [ 1..pp ] do
		            relsA[g] := []; 
                    relsA[g][g] := gensA[g]^pn;
		            for l in [ 1..g-1 ] do
		                relsA[g][l] := gensA[g]^gensA[l]/gensA[g];
		            od;
		            if j <> 1 or k <> 1 then
		                relsA[g-1][g-1] := relsA[g-1][g-1]/gensA[g];
		            fi;
		            g := g + 1;
		        od;
	        od;
	    fi;
	od;

    relsA := Flat( relsA );
    A     := A / relsA;

    # compute corresponding pc group
	G := PcGroupFpGroup( A );
    gensG := Pcgs( G );
  
    # set up epimorphism F -> A -> G
	imgs  := [];
	for i in [ 1..ng ] do
	  imgs[i] := One( G );
	  for j in [ 1..ng ] do
	    if Q[i][j] <> 0 and D[j][j] <> 1 then
	      imgs[i] := imgs[i] * gensG[index[j]]^( Q[i][j] mod D[j][j] );
	    fi;
	  od;
	od;

    # compute preimages
	prei := [];
	for i in [ 1..ng ] do
		if D[i][i] <> 1 then
		    r := One( FreeGroupOfFpGroup( F ) );
		    for j in [ 1..ng ] do
		    	if imgs[j] <> One( G ) then
		    	    r := r * gens[j] ^ ( I[i][j] mod Order( imgs[j] ) );
		    	fi;
		    od;
		    g := index[i];
		    for j in [ 1..Length( pf[i].factors ) ] do
		    	pn := pf[i].factors[j];
		    	pp := pf[i].powers [j];
		    	for k in [ 1..pp ] do
		    	    prei[g] := r; 
                    g := g + 1; 
                    r := r ^ pn;
		    	od;
		    od;
		fi;
	od;

    return rec( source := F,
                image  := G,
                imgs   := imgs,
                prei   := prei );
  elif IsMapping(F) then
    if IsWholeFamily(Image(F)) then
      return rec(source:=Source(F),
		image:=Parent(Image(F)), # parent will replace full group
		                         # with other gens.
		imgs:=List(GeneratorsOfGroup(Source(F)),
			    i->Image(F,i)));
    else
      # ensure the image group is the whole family
      gensG:=Pcgs(Image(F));
      G:=GroupByPcgs(gensG);
      return rec(source:=Source(F),
		image:=G,
		imgs:=List(GeneratorsOfGroup(Source(F)),
			    i->PcElementByExponentsNC(FamilyPcgs(G),
			         ExponentsOfPcElement(gensG,Image(F,i)))));
    fi;
  fi;
  Error("Syntax!");
end );

#############################################################################
##
#F  LiftEpimorphismSQ( epi, M, c )
##
InstallGlobalFunction( LiftEpimorphismSQ, function( epi, M, c )
    local F, G, pcgsG, n, H, pcgsH, d, gensf, pcgsN, htil, gtil, mtil,mtilinv,
          w, e, g, m, i, A, V, rel, l, v, mats, j, t, mat, k, elms, imgs,
          lift, null, vec, new, U, sol, sub, elm, r,tval,tvalp,
	  ex,pos,i1,genid,rels,reln,stopi;

    F := epi.source;
    gensf := GeneratorsOfGroup( FreeGroupOfFpGroup( F ) );
    r := Length( gensf );

    genid:=[];
    for i in [1..r] do
      genid[GeneratorSyllable(gensf[i],1)]:=i;
    od;

    d := M.dimension;

    G := epi.image;
    pcgsG := Pcgs( G );
    n := Length( pcgsG );

    H := ExtensionNC( G, M, c );
    pcgsH := Pcgs( H );
    pcgsN := InducedPcgsByPcSequence( pcgsH, pcgsH{[n+1..n+d]} );


    htil := pcgsH{[1..n]};
    gtil := [];
    mtil := [];
    mtilinv:=[];
    for w in epi.imgs do
      e := ExponentsOfPcElement( pcgsG, w );
      g := PcElementByExponentsNC( pcgsH, htil, e );
      Add( gtil, g );
      m := ImmutableMatrix(M.field, IdentityMat( d, M.field ) );
      for i in [1..n] do
	  m := m * M.generators[i]^e[i];
      od;
      Add( mtil,m);
      #Add( mtilinv, ImmutableMatrix(M.field,m^-1 ));
    od;
    mtilinv:=List(mtil,i->i^-1);

    # set up inhom eq
    A := List( [1..r*d], x -> [] );
    V := [];

    # for each relator of G add 

    rels:=RelatorsOfFpGroup(F);

    stopi:=[4,8,15,30,200];
    AddSet(stopi,Length(rels));
    for reln in [1..Length(rels)] do
      if IsInt(reln/100) then
	Info(InfoSQ,2,reln);
      fi;

      rel:=rels[reln];
      l := NrSyllables( rel );

      # right hand side
      # was: v := MappedWord( rel, gensf, gtil );
      v:=One(gtil[1]);
      for i in [1..l] do
	j := genid[GeneratorSyllable(rel,i)];
	ex:=ExponentSyllable(rel,i);
	if ex<0 then
	  v:=v/gtil[j]^(-ex);
	else
	  v:=v*gtil[j]^ex;
	fi;
      od;

      v := ExponentsOfPcElement( pcgsN, v ) * One( M.field );
      Append( V, v );
  
      # left hand side
      mats := ListWithIdenticalEntries( r,
		  Immutable( NullMat( d, d, M.field ) ) );

      # ahulpke, 28-feb-00: it seems to be much more clever, to run
      # through this loop backwards. Then `MappedWord' can be replaced by
      # a multiplication
      # Similarly the iterated calls to `Subword' are very expensive - better
      # use the internal syllable indexing

      # tval is the product from position i on, tvalp the product from
      # position i+1 on (the tval of the last round)
      tval:=One(mats[1]);

      for i in [l,l-1..1] do
	j := genid[GeneratorSyllable(rel,i)];
	ex:=ExponentSyllable(rel,i);
	if ex<0 then
	  pos:=false;
	  ex:=-ex;
	else
	  pos:=true;
	fi;
	for i1 in [1..ex] do
	  tvalp:=tval;
	  if pos then
	    tval:=mtil[j]*tval;
	    mat:=tvalp;
	    mats[j] := mats[j] + mat;
	  else
	    tval:=mtilinv[j]*tval;
	    mat := tval;
	    mats[j] := mats[j] - mat;
	  fi;
	od;
      od;

      for i in [1..r] do
	  for j in [1..d] do
	      k := d * (i-1) + j;
	      Append( A[k], mats[i][j] );
	  od;
      od; 

      # do these tests several times earlier to speed up
      if reln in stopi then
	sol := SolutionMat( A, V );
	# if there is no solution, then there is no lift
	if sol=fail then
#T return value should be fail?
	  if reln<Length(rels) then
	    Info(InfoSQ,3,"early break:",reln);
	  fi;
	  return false;
	fi;
      fi;
    od;

    # create lift
    elms := [];
    for i in [1..r] do
        sub := - sol{[d*(i-1)+1..d*i]}; 
        elm := PcElementByExponentsNC( pcgsN, sub );
        Add( elms, elm );
    od;
    imgs := List( [1..r], x -> gtil[x] * elms[x] ) ;
    lift := rec( source := F,
                 image  := H, 
                 imgs   := imgs ); 

    # in non-split case this is it
    if IsRowVector( c ) then return lift; fi;
  
    # otherwise check
    U    := Subgroup( H, imgs );
    if Size( U ) = Size( H ) 
     and c=0 then # c=0 is the ordinary case
      return lift;
    else
      lift:=false; # indicate the lift is no good
    fi;

    # this is not optimal - see Plesken
    null := NullspaceMat( A );
    Info(InfoSQ,2,"nullspace dimension:",Length(null));
    for vec in null do
        new  := vec + sol;
        elms := [];
        for i in [1..r] do
            sub := new{[d*(i-1)+1..d*i]}; 
            elm := PcElementByExponentsNC( pcgsN, sub );
            Add( elms, elm );
        od;
        imgs := List( [1..r], x -> gtil[x] * elms[x] );
        U    := Subgroup( H, imgs );
        if Size( U ) = Size( H ) then
	  if lift<>false then
	    Info(InfoSQ,2,"found one");
	    lift:=SubdirProdPcGroups(H,imgs,
				     lift.image,lift.imgs);
	    H:=lift[1];
	    imgs:=lift[2];
	  fi;
	  lift := rec( source := F,
		      image  := H, 
		      imgs   := imgs );
	  if c=0 then
	    return lift;
	  fi;
        fi;
    od;

    # give up 
    return lift; # if c=0 this is automatically false
end );

#############################################################################
##
#F  BlowUpCocycleSQ( v, K, F )
##
InstallGlobalFunction( BlowUpCocycleSQ, function( v, K, F )
    local Q, B, vectors, hlp, i, k;

    if F = K then return v; fi;

    Q := AsField( K, F );
    B := Basis( Q );
    vectors:= BasisVectors( B );
    hlp := [];
    for i in [ 1..Length( v ) ] do
        for k in [ 1..Length( vectors ) ] do
        	Add( hlp, Coefficients( B, v[i] * vectors[k] )[1] );
        od;
    od;
    return hlp;
end );

#############################################################################
##
#F  TryModuleSQ( epi, M )
##
InstallGlobalFunction( TryModuleSQ, function( epi, M )
    local  C, lift, co, cb, cc, r, q, ccpos, ccnum, l, v, qi, c;

    # first try a split extension
    lift := LiftEpimorphismSQ( epi, M, 0 );
    if not IsBool( lift ) then return lift; fi;

    # get collector
    C := CollectorSQ( epi.image, M.absolutelyIrreducible, true );

    # compute the two cocycles
    co := TwoCocyclesSQ( C, epi.image, M.absolutelyIrreducible );

    # if there is one non split extension,  try all mod coboundaries
    if 0 < Length(co) then
        cb := TwoCoboundariesSQ( C, epi.image, M.absolutelyIrreducible );

        # use only those coboundaries which lie in <co>
        if 0 < Length(C.avoid)  then
            cb := SumIntersectionMat( co, cb )[2];
        fi;

        # convert them into row spaces
        if 0 < Length(cb)  then
            cc  := BaseSteinitzVectors( co, cb ).factorspace;
        else
            cc := co;
        fi;

        # try all non split extensions
        if 0 < Length(cc)  then
            r  := PrimitiveRoot( M.absolutelyIrreducible.field );
            q  := Size( M.absolutelyIrreducible.field );

            # loop over all vectors of <cc>
            for ccpos in [ 1 .. Length(cc) ]  do
                for ccnum in [ 0 .. q^(Length(cc)-ccpos)-1 ]  do
                    v := cc[Length(cc)-ccpos+1];
                    for l in [ 1 .. Length(cc)-ccpos ]  do
                        qi := QuoInt( ccnum, q^(l-1) );
                        if qi mod q <> q-1  then
                            v := v + r^(qi mod q) * cc[l];
                        fi;
                    od;

                    # blow cocycle up
                    c := BlowUpCocycleSQ( v, M.field, 
                         M.absolutelyIrreducible.field );

                    # try to lift epimorphism

		    lift := LiftEpimorphismSQ( epi, M, c);

		    # return if we have found a lift
		    if not IsBool( lift ) then return lift; fi;

                od;
            od;
        fi;
    fi;

    # give up
    return false;
end );

#############################################################################
##
#F  AllModulesSQ( epi, M )
##
InstallGlobalFunction( AllModulesSQ, function( epi, M,onlyact )
local  C, lift, co, cb, cc, r, q, ccpos, ccnum, l, v, qi,
       c,all,cnt,total,i,j,iter,sel,dim;

    iter:=onlyact<Length(Pcgs(epi.image)); # are we running in iteration?

    all:=epi;

    if not iter then
      # first try a split extension
      # the -1 indicates we want *all* sdps
      lift := LiftEpimorphismSQ( epi, M, -1 );
      if not IsBool( lift ) then 
	all:=lift;
	Info(InfoSQ,2,"semidirect ",Size(all.image)/Size(epi.image)," found");
      fi;
    fi;

    # get collector
    dim:=M.absolutelyIrreducible.dimension;
    C := CollectorSQ( epi.image, M.absolutelyIrreducible, true );

    # compute the two cocycles
    co := TwoCocyclesSQ( C, epi.image, M.absolutelyIrreducible );

    # if there is one non split extension,  try all mod coboundaries
    if 0 < Length(co) then
        cb := TwoCoboundariesSQ( C, epi.image, M.absolutelyIrreducible );

	q:=false;
	if iter and Length(cb)>0 then
	  # we only want those cocycles, which are trivial for the extra
	  # generators
	  # find those indices which can have nontrivial cocycles
	  r:=Length(Pcgs(epi.image));
	  v:=[1..dim];
	  sel:=[];
	  for i in [1..r] do
	    for j in [1..Minimum(i,onlyact)] do
	      UniteSet(sel,((i^2-i)/2+j-1)*dim+v);
	    od;
	  od;
	  v:=IdentityMat(Length(co[1]),M.absolutelyIrreducible.field){sel};
	  v:=ImmutableMatrix(M.absolutelyIrreducible.field,v);

	  r:=SumIntersectionMat(v,co)[2];
	  if Length(r)<Length(co) then
	    Info(InfoSQ,1,"don't need all cocycles/reduced cohomology");
	    co:=r;
	    q:=true; # use as flag whether it got changed
	  fi;
	fi;

        # use only those coboundaries which lie in <co>
        if 0 < Length(C.avoid) or q then
            cb := SumIntersectionMat( co, cb )[2];
        fi;

        # representatives for basis for the 2-cohomology
        if 0 < Length(cb)  then
            cc  := BaseSteinitzVectors( co, cb ).factorspace;
        else
            cc := co;
        fi;

        # try all non split extensions
        if 0 < Length(cc)  then

            r  := PrimitiveRoot( M.absolutelyIrreducible.field );
            q  := Size( M.absolutelyIrreducible.field );

	    total:=Int(q^Length(cc)/(q-1)); # approximately
	    cnt:=0;
            # loop over all vectors of <cc>
            for ccpos in [ 1 .. Length(cc) ]  do
                for ccnum in [ 0 .. q^(Length(cc)-ccpos)-1 ]  do
		  cnt:=cnt+1;
		  if cnt mod 10 =0 then
		    CompletionBar(InfoSQ,2,"cocycle loop: ",cnt/total);
		  fi;
		  v := cc[Length(cc)-ccpos+1];
		  for l in [ 1 .. Length(cc)-ccpos ]  do
		    qi := QuoInt( ccnum, q^(l-1) );
		    if qi mod q <> q-1  then
		      v := v + r^(qi mod q) * cc[l];
		    fi;
		  od;

		  # blow cocycle up
		  c := BlowUpCocycleSQ( v, M.field, 
			M.absolutelyIrreducible.field );

		  # try to lift epimorphism

		  lift := LiftEpimorphismSQ( epi, M, c);

		  # return if we have found a lift
		  if not IsBool( lift ) then 
		    lift:=SubdirProdPcGroups(all.image,all.imgs,
					      lift.image,lift.imgs);
		    all:=rec(source:=epi.source,
			      image:=lift[1],
			      imgs:=lift[2]);
		    Info(InfoSQ,2,"locally ",Size(all.image)/Size(epi.image),
			  " found");
		  fi;

                od;
            od;
        fi;
	CompletionBar(InfoSQ,2,"cocycle loop: ",false);
    fi;

    # return all lifts
    return all;
end );

#############################################################################
##
#F  TryLayerSQ( epi, layer )
##
InstallGlobalFunction( TryLayerSQ, function( epi, layer )
    local field, dim, reps, rep, lift;

    # compute modules for prime
    field := GF(layer[1]);
    dim   := layer[2];
    reps  := IrreducibleModules( epi.image, field, dim );
    reps:=reps[2]; # the actual modules
        
    # loop over the representations
    for rep in reps do
        lift := TryModuleSQ( epi, rep );
        if not IsBool( lift ) then
           if not layer[3] or rep.dimension = dim then
               return lift;
           fi; 
        fi;
    od;
    
    # give up
    return false;
end );

#############################################################################
##
#F  EAPrimeLayerSQ( epi, prime )
##
InstallGlobalFunction( EAPrimeLayerSQ, function( epi, prime )
local field, dim, rep, lift,all,dims,allmo,mo,start,found,genum,genepi;

  # compute modules for prime
  field := GF(prime);
  start:=epi;
  dims:=List(CharacterDegrees(epi.image,prime),i->i[1]);

  genum:=Length(Pcgs(epi.image)); # number of generators of the starting
                                  # group. (We need to consider nontrivial
				  # cocycles only for those elements, as we
				  # only want to get one layer.)
  # build all modules
  allmo:=[];
  for dim in dims do
    rep  := IrreducibleModules( epi.image, field, dim );
    rep:=rep[2]; # the actual modules
    rep:=Filtered(rep,i->i.dimension=dim);
    Info(InfoSQ,1,"Dimension ",dim,", ",Length(rep)," modules");
    allmo[dim]:=rep;
  od;

  repeat # extend as long as possible
    all:=epi;
    genepi:=Length(Pcgs(epi.image));
    found:=false;
    for dim in dims do
      # loop over the representations
      for rep in [1..Length(allmo[dim])] do
	Info(InfoSQ,2,"Module representative ",dim," #",rep);
	mo:=allmo[dim][rep];

	# inflate to extra generators
	if genum<genepi then
	  mo:=GModuleByMats(Concatenation(mo.generators,
	     List([1..genepi-genum],
	          i->One(mo.generators[1]))),field);
	  if allmo[dim][rep].absolutelyIrreducible=allmo[dim][rep] then
	    mo.absolutelyIrreducible:=mo;
	  else
	    mo.absolutelyIrreducible:=GModuleByMats(
	      Concatenation(allmo[dim][rep].absolutelyIrreducible.generators,
	      List([1..genepi-genum],
		  i->One(allmo[dim][rep].absolutelyIrreducible.generators[1]))),
		  allmo[dim][rep].absolutelyIrreducible.field);
	  fi;
	fi;

	lift := AllModulesSQ( epi, mo,genum);
	if Size(lift.image)>Size(epi.image) then
	  found:=true;
	  lift:=SubdirProdPcGroups(all.image,all.imgs,
				    lift.image,lift.imgs);
	  all:=rec(source:=epi.source,
		    image:=lift[1],
		    imgs:=lift[2]);
	  Info(InfoSQ,1,"globally ",Size(all.image)/Size(start.image)," found");
        fi;
      od;
      
    od;
    epi:=all;
  until not found;

  return all;
end );

#############################################################################
##
#F  SQ( <F>, <...> ) / SolvableQuotient( <F>, <...> )
##
InstallGlobalFunction( SolvableQuotient, function ( F, primes )
    local G, epi, tup, lift, i, found, fac, j, p, iso;

    # initialise epimorphism
    epi := InitEpimorphismSQ(F);
    iso := IsomorphismSpecialPcGroup( epi.image );
    epi.image := Image( iso );
    epi.imgs := List( epi.imgs, x -> Image( iso, x ) );
    G   := epi.image;
    Info(InfoSQ,1,"init done, quotient has size ",Size(G));

    # if the commutator factor group is trivial return
    if Size( G ) = 1 then return epi; fi;

    # if <primes> is a list of tuples, it denotes a chief series
    if IsList( primes ) and IsList( primes[1] ) then
        Info(InfoSQ,2,"have chief series given");
        for tup in primes{[2..Length(primes)]} do
            Info(InfoSQ,1,"trying ", tup);
            tup[3] := true;
            lift := TryLayerSQ( epi, tup );
            if IsBool( lift ) then 
                return epi;
            else
                epi := ShallowCopy( lift );
                iso := IsomorphismSpecialPcGroup( epi.image );
                epi.image := Image( iso );
                epi.imgs := List( epi.imgs, x -> Image( iso, x ) );
                G   := epi.image;
            fi;
            Info(InfoSQ,1,"found quotient of size ", Size(G));
        od;
    fi;

    # if <primes> is a list of primes, we have to use try and error
    if IsList( primes ) and IsInt( primes[1] ) then
        found := true;
        i     := 1;
        while found and i <= Length( primes ) do
            p := primes[i];
            tup := [p, 0, false];
            Info(InfoSQ,1,"trying ", tup);
            lift := TryLayerSQ( epi, tup );
            if not IsBool( lift ) then
                epi := ShallowCopy( lift );
                iso := IsomorphismSpecialPcGroup( epi.image );
                epi.image := Image( iso );
                epi.imgs := List( epi.imgs, x -> Image( iso, x ) );
                G := epi.image;
                found := true;
                i := 1; 
            else
                i := i + 1;
            fi;
            Info(InfoSQ,1,"found quotient of size ", Size(G));
        od;
    fi;
                
    # if <primes> is an integer it is size we want
    if IsInt(primes)  then
        i := primes / Size( G );
        found := true;
        while i > 1 and found do
            fac := Collected( FactorsInt( i ) );
            found := false;
            j := 1;
            while not found and j <= Length( fac ) do
                fac[j][3] := false;
		Info(InfoSQ,1,"trying ", fac[j]);
                lift := TryLayerSQ( epi, fac[j] );
                if not IsBool( lift ) then
                    epi := ShallowCopy( lift );
                    iso := IsomorphismSpecialPcGroup( epi.image );
                    epi.image := Image( iso );
                    epi.imgs := List( epi.imgs, x -> Image( iso, x ) );
                    G := epi.image;
                    found := true;
                    i := primes / Size( G );
                else
                    j := j + 1;
                fi;
		Info(InfoSQ,1,"found quotient of size ", Size(G));
            od;
        od;
    fi;

    # this is the result - should be G only with setted epimorphism
    return epi;
end );

InstallGlobalFunction(EpimorphismSolvableQuotient,function(arg)
local g, sq, hom;
  g:=arg[1];
  sq:=CallFuncList(SQ,arg);
  hom:=GroupHomomorphismByImages(g,sq.image,GeneratorsOfGroup(g),sq.imgs);
  return hom;
end);