#######################################################################
#
# Functions for computing fundamental groups, certain quotients 
# of fundamental groups, and covers of finite 2-dimensional simplicial 
# complexes.
#
# These functions are to be used within GAP>=4.4 with GRAPE>=4.3 installed.
#
# Copyright (C) Leonard H. Soicher 1999-2000,2004,2009,2012.
# 
# See the paper: Sarah Rees and Leonard H. Soicher,
# `An algorithmic approach to fundamental groups and covers of 
# combinatorial cell complexes', J. Symbolic Comp. 29 (2000), 59-77.
# This paper contains the theory behind the algorithms used here, as well
# as an example of the use of some GRAPE2.31/GAP3 implementations of these
# algorithms for clique complexes. 
# The current implementations of the algorithms are somewhat different.
# Note that we can now handle *general* 2-dimensional simplicial
# complexes as well as clique complexes -- see the documentation for
# the function `FundamentalRecordSimplicialComplex'.
#
# Updated July, 2004 to be compatible with GRAPE 4.2.
# Updated November, 2009 to be compatible with GAP 4.4.
# Updated September, 2012 to be compatible with GRAPE 4.6.1.
#

LoadPackage("grape");

SymmetricCentralizerTransitiveGroup := function(G)
#
# Let  G  be a transitive permutation group on  [1..n].
# Then this function returns the centralizer of  G  in the 
# symmetric group  S_n.
#
local n,i,j,transversal,im,g,orb,gens,orbs,p,pimages;
if not IsPermGroup(G) then 
   Error("usage: SymmetricCentralizerTransitiveGroup( <PermGroup> )");
fi;
n:=Maximum(LargestMovedPoint(GeneratorsOfGroup(G)),1);
transversal:=[()];
orb:=[1];
for i in orb do 
   for g in GeneratorsOfGroup(G) do
      im:=i^g;
      if not IsBound(transversal[im]) then
	 #
	 # new point in orbit.
	 #
	 Add(orb,im); 
	 transversal[im]:=transversal[i]*g;
      fi; 
   od; 
od; 
if Length(orb)<>n then
   Error("<G> not transitive");
fi;
#
# Now  transversal[i]  is an element of  G  taking  1  to  i
# (for i=1,...,n),  with transversal[1]=().
# 
gens:=[];                   # generators of the centralizer so far.
orbs:=List([1..n],x->[x]);  # orbits of the centralizer so far.  
			    #  orbs  is maintained as a list of sets.
for i in [2..n] do
   if ForAny(orbs,x->i=x[1]) then
      #
      #  i  is the least element in its orbit.
      # We thus need to determine if there is an element  p  in 
      # the centralizer of  G,  such that  1^p = i.  
      # This  p,  if it exists, is uniquely determined, and is
      # added to  gens.
      #
      p:=[];
      pimages:=BlistList([1..n],[]);
      for j in [1..n] do
	 p[j]:=i^transversal[j];
	 pimages[p[j]]:=true;
      od;
      if SizeBlist(pimages)=n then
	 #
	 #  p  defines a permutation of  [1..n].
	 #
	 p:=PermList(p);
	 if ForAll(GeneratorsOfGroup(G),g->g*p=p*g) then
	    #
	    #  p  centralizes  G.
	    #
	    Add(gens,p);
	    orbs:=List(Orbits(Group(gens,()),[1..n]),Set);
	 fi;
      fi;
   fi;
od;
return Group(gens,());
end;

FundamentalRecordSimplicialComplex := function(arg)
#
# This function returns a so-called fundamental record for the 
# 2-dimensional simplicial complex  K  defined by  arg[1].
#
# If  gamma=arg[1]  is a graph, then it must be non-empty,
# simple and connected, and the simplicial complex  K  is
# taken to be the clique complex  K(gamma).  (The 0-, 1-, and
# 2-simplices of  K(gamma)  are respectively the vertices, edges, 
# and triangles of  gamma.)
# 
# Otherwise   arg[1]  must be a list  [gamma,triangles],  such that   
# gamma  is a non-empty simple, connected graph, and triangles 
# is a list of triangles of  gamma  (given as sets of vertices). 
# Then our simplicial complex  K  has 1-skeleton  gamma,  and the 
# set of 2-simplices of  K  is the union of the  gamma.group-orbits 
# of the elements of  triangles. 
#
# If  arg[2]  is bound then it must be the string  "abelianized"
# (or "abelianised" if you prefer).
#
# If  p=arg[3]  is bound then it must be a positive prime integer.
# 
# A fundamental record for  K  consists of the following components:  
#
#   `group' If  arg[2]  is unbound then this component is a finitely 
#           presented group isomorphic to the fundamental group of  K. 
#           If  arg[2]  is bound, but  arg[3]  is unbound, then this 
#           component is a finitely presented group isomorphic to the 
#           abelianized fundamental group of  K.  If  p=arg[3]  is bound
#           then `group' is an abelian group isomorphic to  H_1(K,GF(p)).
#
#  `edgeLabels'  is a generating list of elements of the above  group,
#                in 1-1 correspondence with the directed edges
#                of  gamma,  such that  edgeLabels[i]  corresponds to  
#                DirectedEdges(gamma)[i]. 
#                The edge labels are used to construct covers of  K.
#                There is a spanning tree  T  of  gamma,  such that
#                each directed edge of  T  corresponds to the identity. 
#                Furthermore, the product of the edge labels going 
#                around a directed 2-simplex is the identity.
#   
#  `abelianized' is equal to true if  arg[2]  is bound, otherwise false.
#
#  `prime'       is equal to  p=arg[3]  if  arg[3]  is bound,
#                otherwise `prime' is unbound.    
#               
local gamma,triangles,abelianized,p,identity,G,gens,delta,unusededges, 
      e,f,i,j,ii,jj,kk,labels,label,ngens,newgen,rel,rels,pres,presgens,
      mult,x,y,z,here,next,visited,edgestoaddQ,inedgestoaddQ,
      firstunused,unused,directededges,F,zero,one,edgelinks,
      edgelinksbound,two_simplices;
 
firstunused := function(i)
#
# This function returns the least  j>=i  such that
# directededges[j] is in the set (defined by) 
# unusededges.  If there is no such j, returns `fail'.
# Assumes that  i>=1. 
#
local j;
for j in [i..Length(directededges)] do
   if unusededges[directededges[j][1]][directededges[j][2]] then
      return j;
   fi;
od;
return fail;
end;

if not (IsGraph(arg[1]) or IsList(arg[1])) then 
   Error("usage: FundamentalRecordSimplicialComplex( <Graph> or <List>",
         " [, \"abelianized\" [, <Int> ]] )"); 
fi;
if IsGraph(arg[1]) then
   gamma:=arg[1];
   edgelinksbound:=false;
else
   gamma:=arg[1][1];
   triangles:=arg[1][2];
   if ForAny(triangles,x->not IsSSortedList(x) or Length(x)<>3) then
      Error("the elements of <triangles> must be sets of size 3");
   fi;
   edgelinks:=List([1..gamma.order],i->List([1..gamma.order],x->[]));
   edgelinksbound:=true;
   two_simplices:=Concatenation(Orbits(gamma.group,triangles,OnSets));
   for x in two_simplices do 
      for y in x do
         e:=Difference(x,[y]);
         Add(edgelinks[e[1]][e[2]],y);
      od;
   od;
   Unbind(two_simplices);
   for i in [1..gamma.order-1] do
      for j in [i+1..gamma.order] do
         edgelinks[i][j]:=Set(edgelinks[i][j]);
         edgelinks[j][i]:=edgelinks[i][j];
      od;
   od;
fi;
if IsBound(arg[2]) then 
   if arg[2]<>"abelianized" and arg[2]<>"abelianised" then
      Error("<arg[2]>, if bound, must be \"abelianized\" (or \"abelianised\")");
   fi;
   abelianized:=true;
   identity:=[];
else
   abelianized:=false;
   pres:=PresentationFpGroup(FreeGroup(0));
   identity:=Identity(pres);
fi;
if IsBound(arg[3]) then 
   p:=arg[3];
   if not IsInt(p) then
      Error("usage: FundamentalRecordSimplicialComplex( <Graph> or <List>",
            " [, \"abelianized\" [, <Int> ]] )"); 
   fi;
   if p < 0 or (not IsPrime(p)) then
      Error("<arg[3]> must be a positive prime");
   fi;
   F:=GF(p);
   zero:=Zero(F);
   one:=One(F);
elif abelianized then
   p:=0;
   zero:=0;
   one:=1;
fi;
if gamma.order=0 then 
   Error("<gamma> must have at least one vertex");
fi;
if not IsSimpleGraph(gamma) then 
   Error("<gamma> must be a simple graph");
fi;
gamma:=NewGroupGraph(Group([],()),gamma);
directededges:=DirectedEdges(gamma);
delta:=NullGraph(Group([],()),gamma.order);
#
# From here on within this function, both  gamma.group  and  delta.group  
# are trivial. 
#
labels:=List([1..gamma.order],x->[]);
edgestoaddQ:=[];
inedgestoaddQ:=List([1..gamma.order],x->BlistList([1..gamma.order],[]));
#
# labels[x][y]  will store the edge-label of the (directed) edge  [x,y] 
# of  delta.
#
# inedgestoaddQ[x][y]=true  iff  [x,y]  is an edge of  gamma, 
# but not of delta,  such that  [x,y]  (together with  [y,x])  extends a  
# 2-claw  of  delta  to a triangle  {x,y,z}  which is a 2-simplex of 
# gamma.  The queue  edgestoaddQ  stores such edges  [x,y]. 
# Important note:  inedgestoaddQ[x][y]=true  iff  inedgestoaddQ[y][x]=true,  
# in which case just one of  [x,y]  and  [y,x]  is stored on  edgestoaddQ.
# 
unusededges:=List([1..gamma.order],x->BlistList([1..gamma.order],
                                      gamma.adjacencies[x]));
#
# The set of directed edges of  gamma  is the disjoint union of the 
# directed edges of  delta  and the edges  [x,y]  with  unusededges[x][y]=true.
# Note that  unusededges[x][y]=true  iff  unusededges[y][x]=true.
# 
# We now perform a BFS to make  delta  a spanning tree  T  for  gamma.
#
next:=[1];
visited:=BlistList([1..gamma.order],[1]);
while next<>[] do
   here:=next;
   next:=[];
   for x in here do 
      for y in gamma.adjacencies[x] do
	 if not visited[y] then
	    Add(next,y);
	    visited[y]:=true;
            labels[x][y]:=StructuralCopy(identity);
            labels[y][x]:=StructuralCopy(identity);
            unusededges[x][y]:=false;
            unusededges[y][x]:=false;
            AddSet(delta.adjacencies[x],y);
            AddSet(delta.adjacencies[y],x);
            #
            # now maintain  edgestoaddQ.
            #
            if edgelinksbound then 
               jj:=edgelinks[x][y];
            else
               jj:=gamma.adjacencies[y];
            fi;
            for z in Intersection(delta.adjacencies[x],jj) do
               inedgestoaddQ[y][z]:=true;
               inedgestoaddQ[z][y]:=true;
               Add(edgestoaddQ,[y,z]);
            od;
         fi;
      od;
   od;
od;
if ForAny(visited,x->x=false) then
   Error("<gamma> must be connected");
fi;
ngens:=0;
if not abelianized or p=0 then
   rels:=[];
fi;
unused:=firstunused(1);
while unused <> fail do
   #
   # Loop invariant:  delta  is a simple connected subgraph
   #                  of  gamma  on the same vertex-set,  
   #                  and the set of directed edges of  gamma   
   #                  is the disjoint union of the directed edges of 
   #                  delta  and the set (defined by)  unusededges. 
   #                  If not abelianized then  pres  is a presentation
   #                  for the fundamental group of the subcomplex  J  of  K
   #                  induced on the graph  delta.
   #                  If abelianized and p=0 then 
   #                  (the free abelian group of rank ngens)/rels 
   #                  is isomorphic to the abelianized fundamental 
   #                  group of  J.  If abelianized and p>0 then
   #                  GF(p)^ngens is isomorphic to  H_1(J,GF(p)). 
   #
   for e in edgestoaddQ do
      x:=e[1];
      y:=e[2];
      unusededges[x][y]:=false;
      unusededges[y][x]:=false;
      inedgestoaddQ[x][y]:=false;
      inedgestoaddQ[y][x]:=false;
      #
      #  [x,y]  completes some 2-claw in  delta  to a 2-simplex.
      #
      AddSet(delta.adjacencies[x],y);
      AddSet(delta.adjacencies[y],x);
      if ngens = 0 then
         labels[x][y]:=StructuralCopy(identity);
         labels[y][x]:=StructuralCopy(identity);
      else
         #
         #  label  [x,y]  and  [y,x],  and determine possible new relators.
         #
         ii:=Intersection(delta.adjacencies[x],delta.adjacencies[y]);
         if edgelinksbound then  
            IntersectSet(ii,edgelinks[x][y]);
         fi;
         if abelianized then
            labels[x][y]:=labels[x][ii[1]]+labels[ii[1]][y]; 
            labels[y][x]:=-labels[x][y]; 
         else
            labels[x][y]:=labels[x][ii[1]]*labels[ii[1]][y]; 
            labels[y][x]:=labels[x][y]^(-1);
         fi;
         i:=1;
         while i<Length(ii) and ngens>0 do 
            i:=i+1;
            if abelianized then
               rel:=labels[x][y]+labels[y][ii[i]]+labels[ii[i]][x]; 
               if rel<>zero*[1..ngens] then 
                  if p=0 then 
                     if not (-rel in rels) then
                        AddSet(rels,rel);
                     fi;
                  else  #  p>0  and non-trivial relator  rel  to process
                     if ngens=1 then
                        ngens:=0; 
                        for jj in [1..delta.order] do
                           for kk in delta.adjacencies[jj] do
                              labels[jj][kk]:=[];
                           od;
                        od;
                     else 
                        j:=1;
                        while rel[j]=zero do
                           j:=j+1;
                        od;
                        MultRowVector(rel,-rel[j]^(-1));
                        # 
                        # rel  now has  zero  in the first j-1 positions,
                        # and   -one  in the jth position.
                        #
                        ngens:=ngens-1;
                        #
                        # rewrite the edge-labels to quotient out by  rel
                        #
                        for jj in [1..delta.order] do
                           for kk in delta.adjacencies[jj] do
                              if kk>=jj then
                                 break;
                              fi;
                              label:=labels[jj][kk];
                              AddRowVector(label,rel,label[j]);
                              # now  label  has a zero in the j-th position.
                              if j=1 then
                                 LeftShiftRowVector(label,1);
                              else
                                 for e in [j..ngens] do
                                    label[e]:=label[e+1];
                                 od;
                                 Unbind(label[ngens+1]);
                              fi;
                              labels[kk][jj] := -label;
                           od;
                        od;
                     fi;
                  fi;
               fi;
            else  # not abelianized 
               rel:=labels[x][y]*labels[y][ii[i]]*labels[ii[i]][x]; 
               if rel<>identity and (not rel in rels) 
                                and (not rel^-1 in rels) then
	          AddRelator(pres,rel);
                  AddSet(rels,rel);
               fi;
            fi;
         od;
      fi;
      #
      # Determine the edges to add to  edgestoaddQ.
      #
      for f in [[x,y],[y,x]] do
         if edgelinksbound then
            jj:=edgelinks[f[1]][f[2]];
         else
            jj:=gamma.adjacencies[f[2]];
         fi;                                                  
         jj:=Intersection(jj,delta.adjacencies[f[1]]);
         SubtractSet(jj,delta.adjacencies[f[2]]);
         for z in jj do
            if not inedgestoaddQ[f[2]][z] then 
               inedgestoaddQ[f[2]][z]:=true;
               inedgestoaddQ[z][f[2]]:=true;
               Add(edgestoaddQ,[f[2],z]);
            fi;
         od;
      od;
   od;
   edgestoaddQ:=[];
   unused:=firstunused(unused);
   if unused<>fail then
      #  
      # We must increase the number of generators.
      #
      ngens:=ngens+1;
      x:=directededges[unused][1];
      y:=directededges[unused][2];
      #
      # We shall label  [x,y]  with a new generator.
      #
      if abelianized then  
         newgen:=zero*[1..ngens];
         newgen[ngens]:=one;
         if ngens=1 and p>0 then
            ConvertToVectorRep(newgen);
         fi;
         labels[x][y]:=newgen;
         labels[y][x]:=-newgen;
         # now maintain  labels  and  rels.
         for jj in [1..delta.order] do
            for kk in delta.adjacencies[jj] do
               labels[jj][kk][ngens]:=zero;
               if ngens=1 and p>0 then
                  ConvertToVectorRep(labels[jj][kk]);
               fi;
            od;
         od;                           
         if p=0 then 
            for rel in rels do 
               rel[ngens]:=zero;
            od;
         fi;   
      else
         AddGenerator(pres);
         newgen:=pres!.(ngens);
         labels[x][y]:=newgen;
         labels[y][x]:=newgen^-1;
      fi;
      unusededges[x][y]:=false;
      unusededges[y][x]:=false;
      AddSet(delta.adjacencies[x],y);
      AddSet(delta.adjacencies[y],x);
      #
      # Determine the edges to add to  edgestoaddQ.
      #
      for f in [[x,y],[y,x]] do
         if edgelinksbound then
            jj:=edgelinks[f[1]][f[2]];
         else
            jj:=gamma.adjacencies[f[2]];
         fi;
         for z in Intersection(jj,delta.adjacencies[f[1]]) do
            inedgestoaddQ[f[2]][z]:=true;
            inedgestoaddQ[z][f[2]]:=true;
            Add(edgestoaddQ,[f[2],z]);    
         od;
      od;
      unused:=firstunused(unused);
   fi;
od; 
if not abelianized then
   G:=FpGroupPresentation(pres);
   if Length(GeneratorsOfGroup(G))<>ngens then
      Error("Internal problem: GeneratorsOfGroup(<G>) inconsistent");
   fi;   
   if ngens=0 then
      labels:=ListWithIdenticalEntries(Length(directededges),Identity(G));
   else
      presgens:=GeneratorsOfPresentation(pres);
      gens:=GeneratorsOfGroup(G);
      labels:=List(directededges,
                   e->MappedWord(labels[e[1]][e[2]],presgens,gens));
   fi;
   return rec(group:=G, edgeLabels:=labels, abelianized:=false);
fi;
#
# now must have abelianized=true
#
if p=0 then
   if ngens=0 then
      G:=FreeGroup(0);
      labels:=ListWithIdenticalEntries(Length(directededges),Identity(G));
   else 
      G:=FreeGroup(ngens);
      gens:=GeneratorsOfGroup(G);
      rels:=Concatenation(List(Combinations(gens,2),x->Comm(x[1],x[2])),
            List(rels,x->Product(List([1..ngens],i->gens[i]^x[i]))));
      G:=G/rels;
      if Length(GeneratorsOfGroup(G))<>ngens then
         Error("Internal problem: GeneratorsOfGroup(<G>) inconsistent");
      fi;   
      gens:=GeneratorsOfGroup(G);
      labels:=List(directededges,e->labels[e[1]][e[2]]);
      Apply(labels,x->Product(List([1..ngens],i->gens[i]^x[i])));
   fi; 
   return rec(group:=G, edgeLabels:=labels, abelianized:=true);
fi;
#
# now must have abelianized=true and p>0.
#
if ngens=0 then  
   G:=TrivialGroup(IsPcGroup);
   labels:=ListWithIdenticalEntries(Length(directededges),Identity(G));
else
   G:=AbelianGroup(ListWithIdenticalEntries(ngens,p));
   gens:=GeneratorsOfGroup(G);
   labels:=List(directededges,e->labels[e[1]][e[2]]);
   Apply(labels,x->List(x,IntFFE));
   Apply(labels,x->Product(List([1..ngens],i->gens[i]^x[i])));
fi;
return rec(group:=G, edgeLabels:=labels, abelianized:=true, prime:=p);
end;

LiftedAutomorphism := function(arg)
#
# Given an automorphism  aut=arg[3]  of the connected simple graph  
# gamma=arg[1]  (which must have at least one vertex),  
# this function attempts to lift  aut  to an automorphism 
# (permuting the fibres) of the (assumed) connected m-fold
# cover  arg[2]=delta  of  gamma.  If  D=arg[4]  is bound, then it
# is assumed to be a subgroup of the deck group of  delta.
#
# This function returns `fail' if  aut  does not lift to 
# an automorphism of  delta;  otherwise a lift of  aut  is returned. 
# 
# We assume that the  i-th  vertex of  delta  projects
# to vertex  ((i-1) div m)+1  of  gamma  (as it would
# if produced by the function  CoveringGraph  (applied to  gamma)). 
#
# *Warning*  This function does only partial checking of the assumptions 
# on  gamma,  delta,  and  aut  (and  D  if  arg[4]  is bound),
# which must have the properties described above.
#
local gamma,delta,aut,D,liftedaut,m,rep,reps,x,y,here,next,visited,
      adjliftedautx,c,int,canlift,liftedautimages,projection,leastlift;
gamma:=arg[1];
delta:=arg[2];
aut:=arg[3];
if IsBound(arg[4]) then 
   D:=arg[4];
else
   D:=Group([],());
fi;
if not (IsGraph(gamma) and IsGraph(delta) 
	and IsPerm(aut) and IsPermGroup(D)) then 
   Error(
   "usage: LiftedAutomorphism( <Graph>, <Graph>, <Perm> [, <PermGroup> ] )");
fi;
if gamma.order=0 then
   Error("<gamma> must have at least one vertex");
fi;
m:=delta.order/gamma.order;  #  delta  is an m-fold cover of  gamma.
if not IsInt(m) or m=0 then
   Error("<delta> is not a connected cover of <gamma>");
fi;
if not IsSimpleGraph(gamma) then
   Error("<gamma> not a simple graph");
fi;
if not IsSimpleGraph(delta) then
   Error("<delta> not a simple graph");
fi;

projection := function(v)
#
# Returns the projection in  gamma  of the vertex  v  of delta.
#
return QuoInt(v-1,m)+1;
end;

leastlift := function(v)
#
# Returns the least element in  [1..delta.order]  that projects to 
# vertex  v  of  gamma. 
#
return m*(v-1)+1;
end;

c:=leastlift(1^aut);
reps:=List(Orbits(D,[c..c+m-1]),x->x[1]); 
for rep in reps do 
   #
   # Determine if  aut  can lift to an element mapping  1  to  rep.
   # We build this lift in the list  liftedaut,  such that, if  liftedaut[i]
   # is bound, then it is the image of  i  under the lifted automorphism
   # determined so far. 
   #
   liftedaut:=[rep];
   liftedautimages:=BlistList([1..delta.order],liftedaut);
   canlift:=true;
   #
   # Now we do a BFS on  delta  to determine the images under liftedaut 
   # of the vertices > 1  of  delta.
   #
   next:=[1];
   visited:=BlistList([1..delta.order],[1]);
   while next<>[] and canlift do
      here:=next;
      next:=[];
      for x in here do 
	 adjliftedautx:=Adjacency(delta,liftedaut[x]);
	 for y in Adjacency(delta,x) do
	    if visited[y] then
	       if x<y and (not liftedaut[y] in adjliftedautx) then
		  #
		  # liftedaut  cannot be extended to an automorphism of  delta,
		  # since it would not map the edge  [x,y]  to an edge.
		  # (Note that since  delta  is simple we need only make
		  # this check when  x<y.)
		  #
		  canlift:=false;
	       fi;
	    else
	       Add(next,y);
	       visited[y]:=true;
	       #
	       # now determine  liftedaut[y]
	       #
	       c:=leastlift(projection(y)^aut);
	       int:=Intersection(adjliftedautx,[c..c+m-1]);
	       if Length(int)<>1 then
		  Error("<delta> not a cover of <gamma>,\n",
			"or <aut> not an automorphism of <gamma>");
	       fi;
	       liftedaut[y]:=int[1];
	       if liftedautimages[liftedaut[y]] then
		  #
		  # liftedaut  cannot define a permutation, since it would 
		  # map two different vertices to  liftedaut[y]. 
		  #
		  canlift:=false;
	       else
		  liftedautimages[liftedaut[y]]:=true;
	       fi;
	    fi;
	 od;
      od;
   od;
   if canlift then
      if SizeBlist(visited) < delta.order then
	 #
	 # Not all vertices of  delta  were visited in the BFS.
	 #
	 Error("<delta> not connected");
      fi;
      return PermList(liftedaut);
   fi;
od;
return fail;  
end;

CoveringGraph := function(arg) 
#
# This is a function which constructs a finite (topological) cover 
# of a graph  gamma,  when this cover is defined by a finite
# permutation representation of the fundamental group of  gamma,    
# and certain edge-labels of  gamma,  as described below.
#
# Let  gamma=arg[1]  be a connected simple graph having at least one 
# vertex,  G=arg[2]  a group,  and  edgelabels=arg[3]  a sequence of 
# elements of  G,  such that  Length(edgelabels)  is equal to
# the number of directed edges of  gamma.  We define a mapping 
#  label  from the set  D  of directed edges of  gamma  into  G  by   
#
#      label( DirectedEdges(gamma)[i] ) = edgelabels[i] 
#
# (DirectedEdges(gamma) is a GAP set), and we assume that the following
# properties hold:
#
#    (1)  G is generated by  label(D),  the set of edge-labels.
#    (2)  For each directed edge  [v,w]  in  D,  we have
#         label([v,w]) = (label([w,v])^-1. 
#    (3)  There is a spanning tree  T  of   gamma,  such that
#         label(d) = Identity(G),  for each directed edge  d   of  T.
#
# (These properties are assumed, but not checked, by this function.)
#
# Then  G  is a quotient of the fundamental group  
# pi_1(gamma) = pi_1(gamma,T)  of  gamma  w.r.t.  T  
# (where we regard  gamma  as a 1-dimensional simplicial complex).
#
# If the optional parameter  H=arg[4]  is not given, then  G
# must be a permutation group on  [1..m],  where  m  is defined 
# to be  Maximum(LargestMovedPoint(GeneratorsOfGroup(G)),1).
# In this case, the function  CoveringGraph  returns the m-fold cover  
# delta  of  gamma,  such that the vertices of  delta  are the 
# ordered pairs  {(v,i) | v in V(gamma), i in [1..m]},  and
# (v,i)  is joined by a (directed) edge to  (w,j)  in  delta  if and
# only if  [v,w]  is a (directed) edge of  gamma  and  i^label([v,w])=j. 
# (Note that this  delta  will be connected if and only if  G  is 
# transitive on  [1..m].)  Moreover, the vertices of the returned 
# graph  delta  are ordered so that the  i-th  vertex of  delta  
# projects to vertex  ((i-1) div m)+1 of  gamma.
#
# If the optional parameter  H=arg[4]  is given,  then it must
# be a subgroup of finite index in  G,  in which case 
# G  is replaced by the image of  G  acting on the right cosets 
# of  H,  and the  edge-labels  are replaced by their images under this 
# action.  Then the function  CoveringGraph  returns the cover  delta  
# as defined above for this new permutation group  G  and the 
# new edge-labels. (Note that this cover will be connected.)
#
local gamma,G,edgelabels,H,action,isconnecteddelta,D,DD,ddreps,
      delta,names,m,mm,x,directededges,i,j,b1,b2,lab,L;
gamma:=arg[1];
G:=arg[2];
edgelabels:=arg[3];
if IsBound(arg[4]) then 
   H:=arg[4];
   if not (IsGraph(gamma) and IsGroup(G) 
	   and IsList(edgelabels) and IsGroup(H)) then
      Error("usage: CoveringGraph( <Graph>, <PermGroup>, <List> ), or\n",
	    "CoveringGraph( <Graph>, <Group>, <List>, <Group> )"); 
   fi;
   if Length(GeneratorsOfGroup(G))=0 then
      # Convert  G  to the trivial permutation group to
      # avoid possible problems with a zero-generator FpGroup.
      # Also, convert  H  and  edgelabels.
      G:=Group([],());
      H:=Subgroup(G,[]);
      edgelabels:=ListWithIdenticalEntries(Length(edgelabels),Identity(G));
   fi;
   if IsFpGroup(G) then
      action:=FactorCosetAction(G,H);
   else
      action:=ActionHomomorphism(G,RightCosets(G,H),OnRight);
   fi;
   G:=Image(action);
   isconnecteddelta:=true;
else
   if not (IsGraph(gamma) and IsPermGroup(G) and IsList(edgelabels)) then
      Error("usage: CoveringGraph( <Graph>, <PermGroup>, <List> ), or\n",
	    "CoveringGraph( <Graph>, <Group>, <List>, <Group> )"); 
   fi;
   action:=IdentityMapping(G);
fi;
if gamma.order=0 then 
   Error("<gamma> must have at least one vertex");
fi;
if not IsSimpleGraph(gamma) then 
   Error("<gamma> not a simple graph");
fi;
if not IsConnectedGraph(gamma) then 
   Error("<gamma> not a connected graph");
fi;
m:=Maximum(LargestMovedPoint(GeneratorsOfGroup(G)),1);
names:=Cartesian(VertexNames(gamma),[1..m]);
if not IsBound(isconnecteddelta) then
   isconnecteddelta:=IsTransitive(G,[1..m]);
fi;
if isconnecteddelta then
   #
   # compute the deck group of  delta
   #
   DD:=SymmetricCentralizerTransitiveGroup(G);
   D:=Action(DD,[1..m*gamma.order],
		function(x,g) 
		   return x+(names[x][2]^g-names[x][2]); 
		end); 
else
   DD:=Group([],());
   D:=Group([],());
fi;
#
# From here,  D  is a fixed subgroup of the deck group of  delta,
# with  DD  the image of  D  acting on the fibre  [1..m].
#
delta:=NullGraph(D,Length(names));
AssignVertexNames(delta,names);
directededges:=DirectedEdges(gamma);
if Length(edgelabels)<>Length(directededges) then
   Error("<edgelabels> and DirectedEdges(<gamma>) must have same length"); 
fi;
ddreps:=Set(List(OrbitsDomain(DD,[1..m]),Minimum)); 
mm:=Length(ddreps);
for i in [1..Length(directededges)] do
   lab:=Image(action,edgelabels[i]);
   b1:=(directededges[i][1]-1)*mm;
   b2:=(directededges[i][2]-1)*m;
   for j in [1..mm] do
      AddSet(delta.adjacencies[b1+j],b2+ddreps[j]^lab);
   od;
od;
if isconnecteddelta then
   #
   #  delta  is connected, and so we try to lift generators of  
   #  gamma.group  to automorphisms of  delta. 
   #
   L:=Filtered(List(GeneratorsOfGroup(gamma.group),
	       g->LiftedAutomorphism(gamma,delta,g,D)),x->x<>fail);
   if L<>[] then
      delta:=NewGroupGraph(Group(Concatenation(GeneratorsOfGroup(D),L)),delta);
   fi;
fi;
return delta;
end;