Quelle  grape.g   Sprache: unbekannt

##############################################################################
##
##  grape.g (Version 4.9.3)    GRAPE Library     Leonard Soicher
##
##  Copyright (C) 1992-2025 Leonard Soicher, School of Mathematical Sciences, 
##                      Queen Mary University of London, London E1 4NS, U.K.
##
# This version includes code by Jerry James (debugged by LS) 
# which allows a user to use bliss instead of nauty for computing 
# automorphism groups and canonical labellings of graphs.
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, see https://www.gnu.org/licenses/gpl.html
#

GRAPE_RANDOM := false; # Determines if certain random methods are to be used
                       # in  GRAPE  functions.
                       # The default is that these random methods are not
                       # used (GRAPE_RANDOM=false).
                       # If these random methods are used (GRAPE_RANDOM=true),
                       # this does not affect the correctness of the results,
                       # as documented, returned by GRAPE functions, but may
                       # influence the time taken and the actual (correct)
                       # values returned.  Due to improvements
                       # in GRAPE and in permutation group methods in GAP4,
                       # the use of random methods is rarely necessary,
                       # and should only be employed by GRAPE experts.

GRAPE_NRANGENS := 18;  # The number of random generators taken for a subgroup
         # when  GRAPE_RANDOM=true.

GRAPE_NAUTY := true;   # Use nauty when true, else use bliss.

GRAPE_DREADNAUT_EXE := 
   ExternalFilename(DirectoriesPackagePrograms("grape"),"dreadnaut"); 
   # filename of dreadnaut or dreadnautB executable

GRAPE_BLISS_EXE := ExternalFilename(DirectoriesSystemPrograms(),"bliss"); 
   # filename of bliss executable

GRAPE_DREADNAUT_INPUT_USE_STRING := false;
   # If true then use a string for the stream used for input
   # to dreadnaut/nauty, if false use a file for this. 
   # Using a string is faster than using a file, but may use
   # too much storage.

# The following variant of GAP's Exec is more flexible, and does not require a
# shell. That makes it more reliable on Windows resp. with Cygwin. Moreover,
# it allows to redirect input and output.
BindGlobal("GRAPE_Exec", function(cmd, args, instream, outstream)
  local dir, status;

  if not IsString(cmd) then
    Error("<cmd> must be a file path");
  fi;

  if not IsInputStream(instream) then
    Error("<instream> must be an input stream");
  fi;

  if not IsOutputStream(outstream) then
    Error("<outstream> must be an output stream");
  fi;

  # execute in the current directory
  dir := DirectoryCurrent();

  # execute the command
  status := Process(dir, cmd, instream, outstream, args);

  return status;
end);

BindGlobal("GRAPE_OrbitNumbers",function(G,n)
#
# Returns the orbits of  G  on  [1..n]  in the form of a record 
# containing the orbit representatives and a length  n  list 
# orbitNumbers,  such that  orbitNumbers[j]=i  means that point 
# j  is in the orbit of the i-th representative.
#
local i,j,orbnum,reps,im,norb,g,orb;
if not IsPermGroup(G) or not IsInt(n) then 
   Error("usage: GRAPE_OrbitNumbers( <PermGroup>, <Int> )");
fi;
orbnum:=[];
for i in [1..n] do
   orbnum[i]:=0;
od;
reps:=[];
norb:=0;
for i in [1..n] do
   if orbnum[i]=0 then      # new orbit
      Add(reps,i);
      norb:=norb+1;
      orbnum[i]:=norb;
      orb:=[i];
      for j in orb do 
  for g in GeneratorsOfGroup(G) do
     im:=j^g;
     if orbnum[im]=0 then 
        orbnum[im]:=norb;
        Add(orb,im); 
     fi; 
  od; 
      od; 
   fi;
od;
return rec(representatives:=reps,orbitNumbers:=orbnum);
end);

BindGlobal("GRAPE_NumbersToSets",function(vec)
#
# Returns a list of sets as described by the numbers in vec, i.e.
# i is in the j-th set iff vec[i]=j>0.
#
local list,i,j;
if not IsList(vec) then 
   Error("usage: GRAPE_NumbersToSets( <List> )");
fi;
if Length(vec)=0 then
   return [];
fi;
list:=[];
for i in [1..Maximum(vec)] do
   list[i]:=[];
od;
for i in [1..Length(vec)] do
   j:=vec[i];
   if j>0 then 
      Add(list[j],i);
   fi;
od;
for i in [1..Length(list)] do
   IsSSortedList(list[i]);
od;
return list;
end);

BindGlobal("GRAPE_IntransitiveGroupGenerators",function(arg)
local conjperm,i,newgens,gens1,gens2,max1,max2;
gens1:=arg[1];
gens2:=arg[2];
if IsBound(arg[3]) then
   max1:=arg[3];
else
   max1:=LargestMovedPoint(gens1);
fi;
if IsBound(arg[4]) then
   max2:=arg[4];
else
   max2:=LargestMovedPoint(gens2);
fi;
if not (IsList(gens1) and IsList(gens2) and IsInt(max1) and IsInt(max2)) then 
   Error(
   "usage: GRAPE_IntransitiveGroupGenerators( <List>, <List> [,<Int> [,<Int> ]] )");
fi;
if Length(gens1)<>Length(gens2) then
   Error("Length(<gens1>) <> Length(<gens2>)");
fi;
conjperm:=PermList(Concatenation(List([1..max2],x->x+max1),[1..max1]));
newgens:=[];
for i in [1..Length(gens1)] do
   newgens[i]:=gens1[i]*(gens2[i]^conjperm);
od;
return newgens;
end);

BindGlobal("ProbablyStabilizer",function(G,pt)
#
# Returns a subgroup of  Stabilizer(G,pt),  which is very often 
# the full stabilizer. In fact, if GRAPE_RANDOM=false, then it
# is guaranteed to be the full stabilizer.
#
local sch,orb,x,y,im,k,gens,g,stabgens,i;
if not IsPermGroup(G) or not IsInt(pt) then
   Error("usage: ProbablyStabilizer( <PermGroup>, <Int> )");
fi;
if not GRAPE_RANDOM or HasSize(G) or HasStabChainMutable(G) or IsAbelian(G) then
   return Stabilizer(G,pt);
fi;
#
# At this point we know that  G  is non-abelian.  In particular,
# G  has at least two generators.
#
# First, we make a Schreier vector of permutations for the orbit  pt^G.
#
gens:=GeneratorsOfGroup(G);
sch:=[];
orb:=[pt];
sch[pt]:=();
for x in orb do
   for g in gens do
      im:=x^g;
      if not IsBound(sch[im]) then
  sch[im]:=g;
  Add(orb,im);
      fi;
   od;
od; 
# Now make  gens  into a new randomish generating sequence for  G.
gens:=ShallowCopy(GeneratorsOfGroup(G));
k:=Length(gens);
for i in [k+1..GRAPE_NRANGENS] do
   gens[i]:=gens[Random([1..k])];
od;
for i in [1..Length(gens)] do
   gens[i]:=Product(gens);
od;
# Now make a list  stabgens  of random elements of the stabilizer of  pt.
x:=();
stabgens:=[];
for i in [1..GRAPE_NRANGENS] do
   x:=x*Random(gens);
   im:=pt^x;
   while im<>pt do
      x:=x/sch[im];
      im:=im/sch[im];
   od;
   if x<>() then
      Add(stabgens,x);
   fi;
od;
return Group(stabgens,());
end);

BindGlobal("ProbablyStabilizerOrbitNumbers",function(G,pt,n)
#
# Returns the "orbit numbers" record for a subgroup of  Stabilizer(G,pt), 
# in its action on  [1..n].  
# This subgroup is very often the full stabilizer, and in fact, 
# if  GRAPE_RANDOM=false,  then it is guaranteed to be the full stabilizer. 
#
if not IsPermGroup(G) or not IsInt(pt) or not IsInt(n) then
   Error(
   "usage: ProbablyStabilizerOrbitNumbers( <PermGroup>, <Int>, <Int>  )");
fi;
return GRAPE_OrbitNumbers(ProbablyStabilizer(G,pt),n);
end);

BindGlobal("GRAPE_RepWord",function(gens,sch,r)
#
# Given a sequence  gens  of group generators, and a  (word type)
# schreier vector  sch  made using  gens,  this function returns a 
# record containing the orbit representative for  r  (wrt  sch),  and
# a word in  gens  taking this representative to  r. 
# (We assume  sch  includes the orbit of  r.)
#
local word,w;
word:=[]; 
w:=sch[r];
while w > 0 do
   Add(word,w); 
   r:=r/gens[w]; 
   w:=sch[r];
od;
return rec(word:=Reversed(word),representative:=r);
end);
   
BindGlobal("NullGraph",function(arg)
#
# Returns a null graph with  n  vertices and group  G=arg[1].
# If  arg[2]  is bound then  n=arg[2],  otherwise  n  is the maximum 
# largest moved point of the generators of  G.
# The  names,  autGroup,  maximumClique,  minimumVertexColouring,
# and  canonicalLabelling  components of the 
# returned null graph are left unbound; however, the  isSimple  
# component is set (to true).
#
local G,n,gamma,nadj,sch,orb,i,j,k,im,gens;
G:=arg[1];
if not IsPermGroup(G) or (IsBound(arg[2]) and not IsInt(arg[2])) then
   Error("usage: NullGraph( <PermGroup>, [, <Int> ] )");
fi;
n:=LargestMovedPoint(GeneratorsOfGroup(G)); 
if IsBound(arg[2]) then
   if arg[2] < n  then
      Error("<arg[2]> too small");
   fi;
   n:=arg[2];
fi;
gamma:=rec(isGraph:=true,order:=n,group:=G,schreierVector:=[],
    adjacencies:=[],representatives:=[],isSimple:=true);
#
# Calculate  gamma.representatives,  gamma.schreierVector,  and
# gamma.adjacencies.
#
sch:=gamma.schreierVector; 
gens:=GeneratorsOfGroup(gamma.group); 
nadj:=0;
for i in [1..n] do 
   sch[i]:=0; 
od;
for i in [1..n] do
   if sch[i]=0 then      # new orbit
      Add(gamma.representatives,i);
      nadj:=nadj+1;
      sch[i]:=-nadj;     # tells where to find the adjacency set.
      gamma.adjacencies[nadj]:=[];
      orb:=[i];
      for j in orb do 
         for k in [1..Length(gens)] do
            im:=j^gens[k];
            if sch[im]=0 then 
               sch[im]:=k; 
               Add(orb,im); 
     fi; 
  od; 
      od; 
   fi;
od;
gamma.representatives:=Immutable(gamma.representatives);
gamma.schreierVector:=Immutable(gamma.schreierVector);
return gamma;
end);

BindGlobal("CompleteGraph",function(arg)
#
# Returns a complete graph with  n  vertices and group  G=arg[1].
# If  arg[2]  is bound then  n=arg[2],  otherwise  n  is the maximum 
# largest moved point of the generators of the permutation group  G.
# If the boolean argument  arg[3]  is bound and has 
# value true then the complete graph will have all possible loops, 
# otherwise it will have no loops (the default).
#
# The  names,  autGroup,  maximumClique,  minimumVertexColouring,
# and  canonicalLabelling  components of the 
# returned complete graph are left unbound; however, the  isSimple  
# component is set (appropriately).
#
local G,n,gamma,i,mustloops;
G:=arg[1];
if not IsPermGroup(G) or (IsBound(arg[2]) and not IsInt(arg[2]))
        or (IsBound(arg[3]) and not IsBool(arg[3])) then
   Error("usage: CompleteGraph( <PermGroup>, [, <Int> [, <Bool> ]] )");
fi;
n:=LargestMovedPoint(GeneratorsOfGroup(G)); 
if IsBound(arg[2]) then
   if arg[2] < n  then
      Error("<arg[2]> too small");
   fi;
   n:=arg[2];
fi;
if IsBound(arg[3]) then
   mustloops:=arg[3];
else 
   mustloops:=false;
fi;
gamma:=NullGraph(G,n);
if gamma.order=0 then
   return gamma;
fi;
if mustloops then 
   gamma.isSimple:=false;
fi;
for i in [1..Length(gamma.adjacencies)] do
   gamma.adjacencies[i]:=[1..n];
   if not mustloops then
      RemoveSet(gamma.adjacencies[i],gamma.representatives[i]);
   fi;
od;
return gamma;
end);

BindGlobal("GRAPE_Graph",function(arg)
#
# First suppose that  arg[5]  is unbound or has value  false.
# Then  L=arg[2]  is a list of elements of a set  S  on which 
# G=arg[1]  acts with action  act=arg[3].  Also  rel=arg[4]  is a boolean
# function defining a  G-invariant relation on  S  (so that 
# for  g in G,  rel(x,y)  iff  rel(act(x,g),act(y,g)) ). 
# Then function  GRAPE_Graph  returns the graph  gamma  with vertex-names
# Immutable(Concatenation(Orbits(G,L,act))),  and  x  is joined to  y
# in  gamma  iff  rel(VertexName(gamma,x),VertexName(gamma,y)).
#
# If  arg[5]  has value  true  then it is assumed that  L=arg[2] 
# is invariant under  G=arg[1]  with action  act=arg[3]. Then
# the function  GRAPE_Graph  behaves as above, except that  gamma.names
# becomes an immutable copy of  L.
#
local G,L,act,rel,invt,gamma,vertexnames,i,reps,H,orb,x,y,adj;
G:=arg[1];
L:=arg[2];
act:=arg[3];
rel:=arg[4];
if IsBound(arg[5]) then
   invt:=arg[5];
else
   invt:=false;
fi;
if not (IsGroup(G) and IsList(L) and IsFunction(act) and IsFunction(rel) 
 and IsBool(invt)) then
   Error("usage: GRAPE_Graph( <Group>, <List>, <Function>, <Function> [, <Bool> ] )");
fi;
if invt then
   vertexnames:=Immutable(L);
else
   vertexnames:=Immutable(Concatenation(Orbits(G,L,act)));
fi;
gamma:=NullGraph(Action(G,vertexnames,act),Length(vertexnames));
Unbind(gamma.isSimple);
gamma.names:=vertexnames;
if not GRAPE_RANDOM then
   if (HasSize(G) and Size(G)<>infinity) or 
      (IsPermGroup(G) and HasStabChainMutable(G)) or
      (HasIsNaturalSymmetricGroup(G) and IsNaturalSymmetricGroup(G)) then
      StabChainOp(gamma.group,rec(limit:=Size(G)));
   fi;
fi;
reps:=gamma.representatives;
for i in [1..Length(reps)] do
   H:=ProbablyStabilizer(gamma.group,reps[i]);
   x:=vertexnames[reps[i]];
   if IsTrivial(H) then  
      gamma.adjacencies[i]:=Filtered([1..gamma.order],j->rel(x,vertexnames[j]));
   else
      adj:=[];
      for orb in OrbitsDomain(H,[1..gamma.order]) do
  y:=vertexnames[orb[1]];
  if rel(x,y) then
     Append(adj,orb);
  fi;
      od;
      Sort(adj);
      gamma.adjacencies[i]:=adj;
   fi;
   IsSSortedList(gamma.adjacencies[i]);
od;
return gamma;
end);

DeclareOperation("Graph",[IsGroup,IsList,IsFunction,IsFunction]);
InstallMethod(Graph,"for use in GRAPE with 4 parameters",
   [IsGroup,IsList,IsFunction,IsFunction],0,GRAPE_Graph);
DeclareOperation("Graph",[IsGroup,IsList,IsFunction,IsFunction,IsBool]);
InstallMethod(Graph,"for use in GRAPE with 5 parameters",
   [IsGroup,IsList,IsFunction,IsFunction,IsBool],0,GRAPE_Graph);

BindGlobal("JohnsonGraph",function(n,e)
#
# Returns the Johnson graph, whose vertices are the e-subsets
# of {1,...,n},  with x joined to y iff  Intersection(x,y)
# has size  e-1.
#
local rel,J;
if not IsInt(n) or not IsInt(e) then 
   Error("usage: JohnsonGraph( <Int>, <Int> )");
fi;
if e<0 or n<e then
   Error("must have 0 <= <e> <= <n>");
fi;
rel := function(x,y)
   return Length(Intersection(x,y))=e-1; 
end;
J:=Graph(SymmetricGroup(n),Combinations([1..n],e),OnSets,rel,true);
J.isSimple:=true;
return J;
end);

BindGlobal("HammingGraph",function(d,q)
#
# Where d and q are positive integers, this function returns the
# Hamming graph H(d,q), defined as follows. The set of vertices
# (actually vertex-names) is the set of all d-tuples of elements
# of [1..q], with vertices v and w joined by an edge iff their
# Hamming distance is 1. The group associated with the returned
# graph is  S_q wr S_d  in its product action on the vertices.
#
local W,projection,embedding,moved,act,rel;
if not IsPosInt(d) or not IsPosInt(q) then 
   Error("usage: HammingGraph( <PosInt>, <PosInt> )"); 
fi;
if q=1 then
   # special trivial case
   return Graph(Group(()),[ListWithIdenticalEntries(d,1)],
      function(x,g) return x; end,function(x,y) return false; end,true);
fi;
W:=WreathProductImprimitiveAction(SymmetricGroup([1..q]),SymmetricGroup([1..d]));
projection:=Projection(W);
embedding:=Embedding(W,d+1);
moved:=List([1..d],i->MovedPoints(Image(Embedding(W,i))));
act := function(x,g)
# Product action of  g  on d-tuple  x.
local bb,b,a,y,i;
bb:=g^projection;
b:=bb^embedding;
a:=g*b^(-1); # so g factorises as a*b in W
y:=[];
for i in [1..d] do
   y[i^bb]:=PositionSorted(moved[i],moved[i][x[i]]^a);
od;
return y;
end;
rel := function(x,y)
# boolean function returning true iff the d-tuples  x  and  y
# are at Hamming distance 1. 
local i,count;
count:=0; 
for i in [1..d] do
   if x[i]<>y[i] then
      count:=count+1;
      if count>1 then 
         return false;
      fi;
   fi;
od; 
return count=1;
end;
# Now construct and return the Hamming graph.
return Graph(W,Tuples([1..q],d),act,rel,true);
end); 

BindGlobal("IsGraph",function(obj)
#
# Returns  true  iff  obj  is a (GRAPE) graph.
#
return IsRecord(obj) and IsBound(obj.isGraph) and obj.isGraph=true
   and IsBound(obj.group) and IsBound(obj.schreierVector); 
end);

BindGlobal("CopyGraph",function(gamma)
#
# Returns a "structural" copy  delta  of the graph  gamma,  and
# also ensures that the appropriate components of  delta  are immutable. 
# 
local delta;
if not IsGraph(gamma) then
   Error("usage: CopyGraph( <Graph> )");
fi;
delta:=ShallowCopy(gamma);
delta.adjacencies:=StructuralCopy(delta.adjacencies);
delta.representatives:=Immutable(delta.representatives);
delta.schreierVector:=Immutable(delta.schreierVector);
if IsBound(delta.names) then
   delta.names:=Immutable(delta.names);
fi;
if IsBound(delta.maximumClique) then
   delta.maximumClique:=Immutable(delta.maximumClique);
fi;
if IsBound(delta.minimumVertexColouring) then
   delta.minimumVertexColouring:=Immutable(delta.minimumVertexColouring);
fi;
Unbind(delta.canonicalLabelling); # for safety
return delta;
end);

BindGlobal("OrderGraph",function(gamma)
#
# returns the order of  gamma.
#
if not IsGraph(gamma) then
   Error("usage: OrderGraph( <Graph> )");
fi;
return gamma.order;
end);

DeclareOperation("Vertices",[IsRecord]);
# to avoid the clash with `Vertices' defined in the xgap package
InstallMethod(Vertices,"for GRAPE graph",[IsRecord],0, 
function(gamma)
#
# Returns the vertex-set of graph  gamma.
#
if not IsGraph(gamma) then
   TryNextMethod();
fi;
return [1..gamma.order];
end);

DeclareOperation("IsVertex",[IsRecord,IsObject]);
InstallMethod(IsVertex,"for GRAPE graph",[IsRecord,IsObject],0, 
function(gamma,v)
#
# Returns  true  iff  v  is vertex of  gamma.
#
if not IsGraph(gamma) then
   TryNextMethod();
fi;
return IsInt(v) and v >= 1 and v <= gamma.order;
end);

BindGlobal("AssignVertexNames",function(gamma,names)
#
# Assign vertex names for  gamma,  so that the (external) name of 
# vertex  i  becomes  names[i],  by making  gamma.names  an immutable 
# copy of  names.
#
if not IsGraph(gamma) or not IsList(names) then
   Error("usage: AssignVertexNames( <Graph>, <List> )");
fi;
if Length(names)<>gamma.order then 
   Error("Length(<names>) <> gamma.order");
fi;
gamma.names:=Immutable(names);
end);

BindGlobal("VertexName",function(gamma,v)
#
# Returns the (external) name of the vertex  v  of  gamma.
#
if not IsGraph(gamma) or not IsInt(v) then
   Error("usage: VertexName( <Graph>, <Int> )");
fi;
if IsBound(gamma.names) then 
   return Immutable(gamma.names[v]);
else 
   return v;
fi;
end);

BindGlobal("VertexNames",function(gamma)
#
# Returns the list of (external) names of the vertices of  gamma. 
#
if not IsGraph(gamma) then
   Error("usage: VertexNames( <Graph> )");
fi;
if IsBound(gamma.names) then 
   return Immutable(gamma.names);
else 
   return Immutable([1..gamma.order]);
fi;
end);

BindGlobal("VertexDegree",function(gamma,v)
#
# Returns the vertex (out)degree of vertex  v  in the graph  gamma.
#
local rw,sch;
if not IsGraph(gamma) or not IsInt(v) then
   Error("usage: VertexDegree( <Graph>, <Int> )");
fi;
if v<1 or v>gamma.order then
   Error("<v> is not a vertex of <gamma>");
fi;
sch:=gamma.schreierVector;
rw:=GRAPE_RepWord(GeneratorsOfGroup(gamma.group),sch,v);
return Length(gamma.adjacencies[-sch[rw.representative]]); 
end);

BindGlobal("VertexDegrees",function(gamma)
#
# Returns the set of vertex (out)degrees for the graph  gamma.
#
local adj,degs;
if not IsGraph(gamma) then
   Error("usage: VertexDegrees( <Graph> )");
fi;
degs:=[];
for adj in gamma.adjacencies do
   AddSet(degs,Length(adj));
od;
return degs;
end);

BindGlobal("IsVertexPairEdge",function(gamma,x,y)
#
# Assuming that  x,y  are vertices of  gamma,  returns true
# iff  [x,y]  is an edge of  gamma.
#
local w,sch,gens;
sch:=gamma.schreierVector;
gens:=GeneratorsOfGroup(gamma.group);
w:=sch[x];
while w > 0 do
   x:=x/gens[w];
   y:=y/gens[w];
   w:=sch[x];
od;
return y in gamma.adjacencies[-w];
end);

DeclareOperation("IsEdge",[IsRecord,IsObject]);
InstallMethod(IsEdge,"for GRAPE graph",[IsRecord,IsObject],0, 
function(gamma,e)
#
# Returns  true  iff  e  is an edge of  gamma.
#
if not IsGraph(gamma) then
   TryNextMethod();
fi;
if not IsList(e) or Length(e)<>2 or not IsVertex(gamma,e[1])
   or not IsVertex(gamma,e[2]) then
   return false;
fi;
return IsVertexPairEdge(gamma,e[1],e[2]);
end);

BindGlobal("Adjacency",function(gamma,v)
#
# Returns (a copy of) the set of vertices of  gamma  adjacent to vertex  v.
#
local w,adj,rw,gens,sch;
sch:=gamma.schreierVector;
if sch[v] < 0 then 
   return ShallowCopy(gamma.adjacencies[-sch[v]]);
fi;
gens:=GeneratorsOfGroup(gamma.group);
rw:=GRAPE_RepWord(gens,sch,v);
adj:=gamma.adjacencies[-sch[rw.representative]]; 
for w in rw.word do 
   adj:=OnTuples(adj,gens[w]); 
od;
return SSortedList(adj);
end);

BindGlobal("IsSimpleGraph",function(gamma)
#
# Returns  true  iff graph  gamma  is simple (i.e. has no loops and 
# if [x,y] is an edge then so is [y,x]).  Also sets the isSimple 
# field of  gamma  if this field was not already bound.
#
local adj,i,x,H,orb;
if not IsGraph(gamma) then 
   Error("usage: IsSimpleGraph( <Graph> )");
fi;
if IsBound(gamma.isSimple) then
   return gamma.isSimple;
fi;
for i in [1..Length(gamma.adjacencies)] do
   adj:=gamma.adjacencies[i];
   x:=gamma.representatives[i];
   if x in adj then    # a loop exists
      gamma.isSimple:=false;
      return false;
   fi;
   H:=ProbablyStabilizer(gamma.group,x);
   for orb in OrbitsDomain(H,adj) do
      if not IsVertexPairEdge(gamma,orb[1],x) then
  gamma.isSimple:=false;
  return false;
      fi;
   od;
od;
gamma.isSimple:=true;
return true;
end);

BindGlobal("DirectedEdges",function(gamma)
#
# Returns the set of directed (ordered) edges of  gamma.
#
local i,j,edges;
if not IsGraph(gamma) then 
   Error("usage: DirectedEdges( <Graph> )");
fi;
edges:=[];
for i in [1..gamma.order] do
   for j in Adjacency(gamma,i) do
      Add(edges,[i,j]);
   od;
od;
IsSSortedList(edges);  # edges is a set.
return edges;
end);

BindGlobal("UndirectedEdges",function(gamma)
#
# Returns the set of undirected edges of  gamma,  which must be 
# a simple graph.
#
local i,j,edges;
if not IsGraph(gamma) then 
   Error("usage: UndirectedEdges( <Graph> )");
fi;
if not IsSimpleGraph(gamma) then
   Error("<gamma> must be a simple graph");
fi;
edges:=[];
for i in [1..gamma.order-1] do
   for j in Adjacency(gamma,i) do
      if i<j then 
  Add(edges,[i,j]);
      fi;
   od;
od;
IsSSortedList(edges);  # edges is a set.
return edges;
end);

BindGlobal("AddEdgeOrbit",function(arg) 
#
# Let  gamma=arg[1]  and  e=arg[2].
# If  arg[3]  is bound then it is assumed to be  Stabilizer(gamma.group,e[1]).
# This procedure adds edge orbit  e^gamma.group  to the edge-set of  gamma.
#
local w,word,sch,gens,gamma,e,x,y,orb,u,v;
gamma:=arg[1]; 
e:=arg[2];
if not IsGraph(gamma) or not IsList(e) 
    or (IsBound(arg[3]) and not IsPermGroup(arg[3])) then
   Error("usage: AddEdgeOrbit( <Graph>, <List>, [, <PermGroup> ] )");
fi;
if Length(e)<>2 or not IsVertex(gamma,e[1]) or not IsVertex(gamma,e[2]) then
   Error("invalid <e>");
fi;
sch:=gamma.schreierVector;
gens:=GeneratorsOfGroup(gamma.group);
x:=e[1];
y:=e[2];
w:=sch[x];
word:=[];
while w > 0 do
   Add(word,w);
   x:=x/gens[w];
   y:=y/gens[w];
   w:=sch[x];
od;
if not(y in gamma.adjacencies[-sch[x]]) then
   #  e  is not an edge of  gamma
   if not IsBound(arg[3]) then
      orb:=Orbit(Stabilizer(gamma.group,x),y);
   else
      if ForAny(GeneratorsOfGroup(arg[3]),x->e[1]^x<>e[1]) then
  Error("<arg[3]>  not equal to  Stabilizer(<gamma.group>,<e[1]>)");
      fi;
      orb:=[];
      for u in Orbit(arg[3],e[2]) do
  v:=u;
  for w in word do
     v:=v/gens[w];
  od;
  Add(orb,v);
      od;
   fi;
   UniteSet(gamma.adjacencies[-sch[x]],orb);
   if e[1]=e[2] then
      gamma.isSimple:=false;
   elif IsBound(gamma.isSimple) and gamma.isSimple then
      if not IsVertexPairEdge(gamma,e[2],e[1]) then 
  gamma.isSimple:=false; 
      fi;
   else 
      Unbind(gamma.isSimple);
   fi;
   Unbind(gamma.autGroup);
   Unbind(gamma.canonicalLabelling);
   Unbind(gamma.maximumClique);
   Unbind(gamma.minimumVertexColouring);
fi;   
end);

BindGlobal("RemoveEdgeOrbit",function(arg) 
#
# Let  gamma=arg[1]  and  e=arg[2].
# If  arg[3]  is bound then it is assumed to be  Stabilizer(gamma.group,e[1]).
# This procedure removes the edge orbit  e^gamma.group  from the edge-set 
# of  gamma, if this orbit exists, and otherwise does nothing.
#
local w,word,sch,gens,gamma,e,x,y,orb,u,v;
gamma:=arg[1]; 
e:=arg[2];
if not IsGraph(gamma) or not IsList(e) 
    or (IsBound(arg[3]) and not IsPermGroup(arg[3])) then
   Error("usage: RemoveEdgeOrbit( <Graph>, <List>, [, <PermGroup> ] )");
fi;
if Length(e)<>2 or not IsVertex(gamma,e[1]) or not IsVertex(gamma,e[2]) then
   Error("invalid <e>");
fi;
sch:=gamma.schreierVector;
gens:=GeneratorsOfGroup(gamma.group);
x:=e[1];
y:=e[2];
w:=sch[x];
word:=[];
while w > 0 do
   Add(word,w);
   x:=x/gens[w];
   y:=y/gens[w];
   w:=sch[x];
od;
if y in gamma.adjacencies[-sch[x]] then
   #  e  is an edge of  gamma
   if not IsBound(arg[3]) then
      orb:=Orbit(Stabilizer(gamma.group,x),y);
   else
      if ForAny(GeneratorsOfGroup(arg[3]),x->e[1]^x<>e[1]) then
  Error("<arg[3]>  not equal to  Stabilizer(<gamma.group>,<e[1]>)");
      fi;
      orb:=[];
      for u in Orbit(arg[3],e[2]) do
  v:=u;
  for w in word do
     v:=v/gens[w];
  od;
  Add(orb,v);
      od;
   fi;
   SubtractSet(gamma.adjacencies[-sch[x]],orb);
   if IsBound(gamma.isSimple) and gamma.isSimple then
      if IsVertexPairEdge(gamma,e[2],e[1]) then 
  gamma.isSimple:=false; 
      fi;
   else 
      Unbind(gamma.isSimple);
   fi;
   Unbind(gamma.autGroup);
   Unbind(gamma.canonicalLabelling);
   Unbind(gamma.maximumClique);
   Unbind(gamma.minimumVertexColouring);
fi;   
end);

BindGlobal("EdgeOrbitsGraph",function(arg)
#
# Let  G=arg[1],  E=arg[2].
# Returns the (directed) graph with vertex-set {1,...,n} and edge-set 
# the union over  e in E  of  e^G,  where  n=arg[3]  if  arg[3]  is bound,
# and  n=LargestMovedPoint(GeneratorsOfGroup(G))  otherwise.
# (E can consist of just a singleton edge.)
#
local G,E,n,gamma,e;
G:=arg[1];
E:=arg[2];
if IsBound(E[1]) and IsInt(E[1]) then 
   # assume  E  consists of a single edge.
   E:=[E];
fi;
if IsBound(arg[3]) then 
   n:=arg[3];
else
   n:=LargestMovedPoint(GeneratorsOfGroup(G));
fi;
if not IsPermGroup(G) or not IsList(E) or not IsInt(n) then
   Error("usage: EdgeOrbitsGraph( <PermGroup>, <List> [, <Int> ] )");
fi;
gamma:=NullGraph(G,n);
for e in E do 
   AddEdgeOrbit(gamma,e); 
od;
return gamma;
end);

BindGlobal("GeneralizedOrbitalGraphs",function(arg)
#
# Let  G=arg[1]. Then  G  must be a non-trivial permutation group,
# acting transitively on  [1..n],  where n:=LargestMovedPoint(G).
#
# The optional second parameter  k=arg[2]  (default: false)  must
# be true or false or a non-negative integer. 
#
# Then this function returns a list of distinct generalized orbital
# graphs for  G  (where a *generalized orbital graph* for  G  is a 
# (simple) graph with vertex set [1..n] and edge-set a union of 
# zero or more G-orbits of 2-subsets of  [1..n]).
#
# If  k=true  then *all*  the generalized orbital graphs 
# for  G  are in the list,  if  k=false  (the default)  then all
# the non-null generalized orbital graphs for  G  are in the list,
# and if  k  is a non-negative integer then the list consists of
# all the generalized orbital graphs for  G  whose edge-set is the
# union of exactly  k  G-orbits of 2-subsets of  [1..n].
#
# The group associated with each returned graph in the list is  G. 
#
local G,k,comb,combinations,n,H,result,reps,i,L,M,mm;
if not (Length(arg) in [1,2]) then
   Error("must have 1 or 2 arguments");
fi;
G:=arg[1];
if IsBound(arg[2]) then
   k:=arg[2];
else
   k:=false;
fi;
if not (IsPermGroup(G) and (IsBool(k) or IsInt(k))) then
   Error("usage: GeneralizedOrbitalGraphs( <PermGroup> [, <Bool> or <Int> ] )");
fi;
n:=LargestMovedPoint(G);
if n=0 or not IsTransitive(G,[1..n]) then
   Error("<G> must be a non-trivial transitive group on [1..LargestMovedPoint( <G> )]");
fi;
if not ((k in [false,true]) or (IsInt(k) and k>=0)) then 
   Error("<k> must be in  [false,true]  or be a non-negative integer");
fi;
H:=Stabilizer(G,1);
reps:=Set(List(OrbitsDomain(H,[2..n]),Minimum));
#
# Now make a duplicate-free list  L  of the graphs
# with vertex-set  [1..n]  and edge-set the union
# of a nondiagonal G-orbital and its paired orbital.
# At the same time, make a list  M  whose i-th element is
# a list of edges of  L[i],  such that the union of the
# G-orbits of the edges in  M[i]  is the edge-set of  L[i].
#
L:=[];
M:=[];
for i in [1..Length(reps)] do
   if ForAll(L,x->not IsEdge(x,[1,reps[i]])) then
      mm:=[[1,reps[i]],[reps[i],1]];
      Add(L,EdgeOrbitsGraph(G,mm));
      Add(M,mm);
   fi;
od;
result:=[];
if k in [false,true] then
   combinations:=Combinations(M);
else  
   # k is a non-negative integer
   combinations:=Combinations(M,k);
fi;
for comb in combinations do
   if k<>false or comb<>[] then
      Add(result,EdgeOrbitsGraph(G,Concatenation(comb)));
   fi;
od;
return result;
end);

BindGlobal("CollapsedAdjacencyMat",function(arg)
#
# Returns the collapsed adjacency matrix  A  for  gamma=arg[2]  wrt  
# group  G=arg[1],  assuming  G <= Aut(gamma). 
# The rows and columns of  A  are indexed by the orbits 
# orbs[1],...,orbs[n], say, of  G  on the vertices of  
# gamma, and the entry  A[i][j]  of  A  is defined as follows:
#    Let  reps[i]  be a representative of the  i-th  G-orbit  orbs[i].
#    Then  A[i][j] equals the number of neighbours (in  gamma)
#    of  reps[i]  in  orbs[j]. 
# Note that this definition does not depend on the choice of 
# representative  reps[i].
#
# *** New for Grape 2.3: In the special case where this function 
# is given just one argument, then we must have  gamma=arg[1], 
# we must have  gamma.group  transitive on the vertices of  gamma,
# and then the returned collapsed adjacency matrix for  gamma  is 
# w.r.t. the stabilizer in  gamma.group  of  1.  Additionally 
# [1]=orbs[1].  This feature is to 
# conform with the definition of collapsed adjacency matrix in 
# Praeger and Soicher, "Low Rank Representations and Graphs for 
# Sporadic Groups", CUP, Cambridge, 1997.  (In GRAPE we allow a collapsed
# adjacency matrix to be more general, as we can collapse w.r.t. to
# an arbitrary subgroup of  Aut(gamma),  and  gamma  need not 
# even be vertex-transitive.)  
#
local G,gamma,orbs,i,j,n,A,orbnum,reps;
if Length(arg)=1 then
   gamma:=arg[1];
   if not IsGraph(gamma) then
      Error("usage: CollapsedAdjacencyMat( [<PermGroup>,] <Graph>)");
   fi;
   if gamma.order=0 then 
      return []; 
   fi; 
   if not IsTransitive(gamma.group,[1..gamma.order]) then
      Error(
       "<gamma.group> not transitive on vertices of single argument <gamma>"); 
   fi;
   G := Stabilizer(gamma.group,1);
else   
   G := arg[1];
   gamma := arg[2];
   if not IsPermGroup(G) or not IsGraph(gamma) then
      Error("usage: CollapsedAdjacencyMat( [<PermGroup>,] <Graph> )");
   fi;
   if gamma.order=0 then
      return [];
   fi;
fi;
orbs:=GRAPE_OrbitNumbers(G,gamma.order);
orbnum:=orbs.orbitNumbers;
reps:=orbs.representatives;
n:=Length(reps);
A:=NullMat(n,n);
for i in [1..n] do 
   for j in Adjacency(gamma,reps[i]) do
      A[i][orbnum[j]]:=A[i][orbnum[j]]+1;
   od;
od;
return A;
end);

BindGlobal("OrbitalDigraphColadjMats",function(arg)
#
# This function returns a sequence of collapsed adjacency 
# matrices for the the orbital digraphs of the (assumed) transitive  
# G=arg[1].  The matrices are collapsed w.r.t.  Stabilizer(G,1),  so
# that these are collapsed adjacency matrices in the sense of 
# Praeger and Soicher, "Low Rank Representations and Graphs for 
# Sporadic Groups", CUP, Cambridge, 1997. 
# The matrices are collapsed w.r.t. a fixed ordering of the G-suborbits,
# with the trivial suborbit  [1]  coming first.
#
# If  arg[2]  is bound, then it is assumed to be  Stabilizer(G,1).
#
local G,H,orbs,deg,i,j,k,n,coladjmats,A,orbnum,reps,gamma;
G:=arg[1];
if not IsPermGroup(G) or (IsBound(arg[2]) and not IsPermGroup(arg[2])) then
   Error("usage: OrbitalDigraphColadjMats( <PermGroup> [, <PermGroup> ] )");
fi;
if IsBound(arg[2]) then
   H:=arg[2];
   if ForAny(GeneratorsOfGroup(H),x->1^x<>1) then
      Error("<H> does not fix the point 1");
   fi;
else
   H:=Stabilizer(G,1);
fi;
deg:=Maximum(LargestMovedPoint(GeneratorsOfGroup(G)),1);
if not IsTransitive(G,[1..deg]) then
   Error("<G> not transitive");
fi;
gamma:=NullGraph(G,deg);
orbs:=GRAPE_OrbitNumbers(H,gamma.order);
orbnum:=orbs.orbitNumbers;
reps:=orbs.representatives;
if reps[1]<>1 then # this cannot happen!
   Error("internal error");
fi;
n:=Length(reps);
coladjmats:=[];
for i in [1..n] do
   AddEdgeOrbit(gamma,[1,reps[i]],H);
   A:=NullMat(n,n);
   for j in [1..n] do 
      for k in Adjacency(gamma,reps[j]) do
  A[j][orbnum[k]]:=A[j][orbnum[k]]+1;
      od;
   od;
   coladjmats[i]:=A;
   if i < n then
      RemoveEdgeOrbit(gamma,[1,reps[i]],H);
   fi;
od;
return coladjmats;
end);

BindGlobal("OrbitalGraphColadjMats",OrbitalDigraphColadjMats);
# for backward compatibility

BindGlobal("LocalInfo",function(arg)
#
# Calculates  "local info"  for  gamma=arg[1]  from point of view of vertex  
# set (or list or singleton vertex)  V=arg[2].
#
# Returns record containing the  "layer numbers"  for gamma w.r.t.
# V,  as well as the the local diameter and local girth of gamma w.r.t.  V.
# ( layerNumbers[i]=j>0 if vertex  i  is in layer[j]  (i.e. at distance
# j-1  from  V),  layerNumbers[i]=0  if vertex  i
# is not joined by a (directed) path from some element of  V.)
# Also, if a local parameter  ci[V], ai[V], or bi[V]  exists then
# this information is recorded in  localParameters[i+1][1], 
# localParameters[i+1][2], or localParameters[i+1][3], 
# respectively (otherwise a -1 is recorded). 
#
# *** If  gamma  is not simple then local girth and the local parameters
# may not be what you think. The local girth has no real meaning if 
# |V| > 1.
#
# *** But note: If arg[3] is bound and arg[3] > 0
# then the procedure stops after  layers[arg[3]]  has been determined.
#
# If  arg[4]  is bound (a set or list or singleton vertex), then the
# procedure stops when the first layer containing a vertex in  arg[4]  is 
# complete.  Moreover, if  arg[4]  is bound then the local info record 
# contains a  distance  field, whose value (if  arg[3]=0)  is  
# min{ d(v,w) | v in V, w in arg[4] }.
#
# If  arg[5]  is bound then it is assumed to be a subgroup of Aut(gamma)
# stabilising  V  setwise.
#
local gamma,V,layers,localDiameter,localGirth,localParameters,i,j,x,y,next,
      nprev,nhere,nnext,sum,orbs,orbnum,laynum,lnum,
      stoplayer,stopvertices,distance,loc,reps,layerNumbers;
gamma:=arg[1];
V:=arg[2];
if IsInt(V) then 
   V:=[V];
fi;
if not IsGraph(gamma) or not IsList(V) then 
   Error("usage: LocalInfo( <Graph>, <Int> or <List>, ... )");
fi;
if not IsSSortedList(V) then
   V:=SSortedList(V);
fi;
if V=[] or not IsSubset([1..gamma.order],V) then 
   Error("<V> must be non-empty set of vertices of <gamma>");
fi;
if IsBound(arg[3]) then 
   stoplayer:=arg[3]; 
   if not IsInt(stoplayer) or stoplayer < 0 then
      Error("<stoplayer> must be integer >= 0");
   fi;
else 
   stoplayer:=0; 
fi;
if IsBound(arg[4]) then 
   stopvertices:=arg[4]; 
   if IsInt(stopvertices) then
      stopvertices:=[stopvertices];
   fi;
   if not IsSSortedList(stopvertices) then
      stopvertices:=SSortedList(stopvertices);
   fi;
   if not IsSubset([1..gamma.order],stopvertices) 
      then
  Error("<stopvertices> must be a set of vertices of <gamma>");
   fi;
else 
   stopvertices:=[]; 
fi;
if IsBound(arg[5]) then 
   if not IsPermGroup(arg[5]) then 
      Error("<arg[5]> must be a permutation group (<= Stab(<V>)");
   fi;
   orbs:=GRAPE_OrbitNumbers(arg[5],gamma.order);
else
   if Length(V)=1 then
      if IsBound(gamma.autGroup) then
  orbs:=ProbablyStabilizerOrbitNumbers(gamma.autGroup,V[1],gamma.order);
      else
  orbs:=ProbablyStabilizerOrbitNumbers(gamma.group,V[1],gamma.order);
      fi;
   else
      orbs:=rec(representatives:=[1..gamma.order],
  orbitNumbers:=[1..gamma.order]); 
   fi;
fi;
orbnum:=orbs.orbitNumbers;
reps:=orbs.representatives;
laynum:=[];
for i in [1..Length(reps)] do 
   laynum[i]:=0;
od;
localGirth:=-1; 
distance:=-1;
localParameters:=[]; 
next:=[];
for i in V do
   AddSet(next,orbnum[i]);
od;
sum:=Length(next);
for i in next do 
   laynum[i]:=1;
od;
layers:=[]; 
layers[1]:=next; 
i:=1; 
if Length(Intersection(V,stopvertices)) > 0 then
   stoplayer:=1; 
   distance:=0;
fi;
while stoplayer<>i and Length(next)>0 do
   next:=[];
   for x in layers[i] do 
      nprev:=0; 
      nhere:=0; 
      nnext:=0;
      for y in Adjacency(gamma,reps[x]) do
  lnum:=laynum[orbnum[y]];
  if i>1 and lnum=i-1 then 
     nprev:=nprev+1;
  elif lnum=i then 
     nhere:=nhere+1;
  elif lnum=i+1 then
     nnext:=nnext+1;
  elif lnum=0 then
     AddSet(next,orbnum[y]); 
     nnext:=nnext+1;
     laynum[orbnum[y]]:=i+1;
  fi;
      od;
      if (localGirth=-1 or localGirth=2*i-1) and nprev>1 then 
  localGirth:=2*(i-1); 
      fi;
      if localGirth=-1 and nhere>0 then 
  localGirth:=2*i-1; 
      fi;
      if not IsBound(localParameters[i]) then 
  localParameters[i]:=[nprev,nhere,nnext];
      else
  if nprev<>localParameters[i][1] then 
     localParameters[i][1]:=-1; 
  fi;
  if nhere<>localParameters[i][2] then 
     localParameters[i][2]:=-1; 
  fi;
  if nnext<>localParameters[i][3] then 
     localParameters[i][3]:=-1; 
  fi;
      fi;
   od;
   if Length(next)>0 then 
      i:=i+1; 
      layers[i]:=next; 
      for j in stopvertices do 
  if laynum[orbnum[j]]=i then 
     stoplayer:=i; 
     distance:=i-1;
            break;
  fi;
      od;
      sum:=sum+Length(next); 
   fi;
od;
if sum=Length(reps) then 
   localDiameter:=Length(layers)-1; 
else 
   localDiameter:=-1; 
fi;
# now change  orbnum  to give the layer numbers instead of orbit numbers.
layerNumbers:=orbnum;
for i in [1..gamma.order] do
   layerNumbers[i]:=laynum[orbnum[i]];
od;
loc:=rec(layerNumbers:=layerNumbers,localDiameter:=localDiameter,
  localGirth:=localGirth,localParameters:=localParameters);
if Length(stopvertices) > 0 then 
   loc.distance:=distance;
fi;
return loc;
end);

BindGlobal("LocalInfoMat",function(A,rows)
#
# Calculates local info on a graph using a collapsed adjacency matrix 
#  A  for that graph.
# This local info is from the point of view of the set of vertices
# represented by the set  rows  of row indices of  A.
# The elements of  layers[i]  will be the row indices representing 
# the vertices of the i-th layer.  
# No  distance  field will be calculated.
#
# *** If  A  is not the collapsed adjacency matrix for a simple graph 
# then  localGirth  and localParameters may not be what you think.
# If  rows  does not represent a single vertex then  localGirth  has 
# no real meaning. 
#
local layers,localDiameter,localGirth,localParameters,i,j,x,y,next,
      nprev,nhere,nnext,sum,laynum,lnum,n;
if IsInt(rows) then 
   rows:=[rows];
fi;
if not IsMatrix(A) or not IsList(rows) then 
   Error("usage: LocalInfoMat( <Matrix>, <Int> or <List> )");
fi;
if not IsSSortedList(rows) then
   rows:=SSortedList(rows);
fi;
n:=Length(A);
if rows=[] or not IsSubset([1..n],rows) then 
   Error("<rows> must be non-empty set of row indices");
fi;
laynum:=ListWithIdenticalEntries(n,0);
localGirth:=-1; 
localParameters:=[]; 
next:=ShallowCopy(rows); 
for i in next do 
   laynum[i]:=1;
od;
layers:=[]; 
layers[1]:=next; 
i:=1; 
sum:=Length(rows);
while Length(next)>0 do
   next:=[];
   for x in layers[i] do 
      nprev:=0; 
      nhere:=0; 
      nnext:=0;
      for y in [1..n] do
  j:=A[x][y];
  if j>0 then
     lnum:=laynum[y];
     if i>1 and lnum=i-1 then 
        nprev:=nprev+j;
     elif lnum=i then 
        nhere:=nhere+j;
     elif lnum=i+1 then
        nnext:=nnext+j;
     elif lnum=0 then
        AddSet(next,y); 
        nnext:=nnext+j;
        laynum[y]:=i+1;
     fi;
  fi;
      od;
      if (localGirth=-1 or localGirth=2*i-1) and nprev>1 then 
  localGirth:=2*(i-1); 
      fi;
      if localGirth=-1 and nhere>0 then 
  localGirth:=2*i-1; 
      fi;
      if not IsBound(localParameters[i]) then 
  localParameters[i]:=[nprev,nhere,nnext];
      else
  if nprev<>localParameters[i][1] then 
     localParameters[i][1]:=-1; 
  fi;
  if nhere<>localParameters[i][2] then 
     localParameters[i][2]:=-1; 
  fi;
  if nnext<>localParameters[i][3] then 
     localParameters[i][3]:=-1; 
  fi;
      fi;
   od;
   if Length(next)>0 then 
      i:=i+1; 
      layers[i]:=next; 
      sum:=sum+Length(next); 
   fi;
od;
if sum=n then 
   localDiameter:=Length(layers)-1; 
else 
   localDiameter:=-1; 
fi;
return rec(layerNumbers:=laynum,localDiameter:=localDiameter,
    localGirth:=localGirth,localParameters:=localParameters);
end);

BindGlobal("InducedSubgraph",function(arg) 
#
# Returns the subgraph of  gamma=arg[1]  induced on the list  V=arg[2]  of
# distinct vertices of  gamma. 
# If  arg[3]  is unbound, then the trivial group is the group associated 
# with the returned induced subgraph. 
# If  arg[3]  is bound, this function assumes that  G=arg[3]   fixes  
# V  setwise, and is a group of automorphisms of the induced subgraph 
# when restricted to  V.  In this case, the image of  G  acting on  V  is 
# the group associated with the returned induced subgraph. 
#
# The i-th vertex of the induced subgraph corresponds to vertex V[i] of
# gamma,  with the i-th vertex-name of the induced subgraph being the 
# vertex-name in  gamma  of V[i].
#
local gamma,V,G,Ggens,gens,indu,i,j,W,VV,X;
gamma:=arg[1];
V:=arg[2];
if not IsGraph(gamma) or not IsList(V) 
    or (IsBound(arg[3]) and not IsPermGroup(arg[3])) then
   Error("usage: InducedSubgraph( <Graph>, <List>, [, <PermGroup> ] )");
fi;
VV:=SSortedList(V);
if Length(V)<>Length(VV) then
   Error("<V> must not contain repeated elements");
fi;
if not IsSubset([1..gamma.order],VV) then
   Error("<V> must be a list of vertices of <gamma>");
fi;
if IsBound(arg[3]) then 
   G:=arg[3];
else
   G:=Group([],());
fi;
W:=[];
for i in [1..Length(V)] do 
   W[V[i]]:=i; 
od;
Ggens:=GeneratorsOfGroup(G);
gens:=[]; 
for i in [1..Length(Ggens)] do 
   gens[i]:=[];
   for j in V do 
      gens[i][W[j]]:=W[j^Ggens[i]]; 
   od;
   gens[i]:=PermList(gens[i]);
od;
indu:=NullGraph(Group(gens,()),Length(V));
if IsBound(gamma.isSimple) and gamma.isSimple then 
   indu.isSimple:=true;
else 
   Unbind(indu.isSimple); 
fi;
for i in [1..Length(indu.representatives)] do
   X:=W{Intersection(VV,Adjacency(gamma,V[indu.representatives[i]]))};
   Sort(X);
   indu.adjacencies[i]:=X;
od;
if not IsBound(gamma.names) then
   indu.names:=Immutable(V);
else
   indu.names:=Immutable(gamma.names{V});
fi;
return indu;
end);

BindGlobal("Distance",function(arg)
#
# Let  gamma=arg[1],  X=arg[2],  Y=arg[3].
# Returns the distance  d(X,Y)  in the graph  gamma, where  X,Y 
# are singleton vertices or lists of vertices.
# (Returns  -1  if no (directed) path joins  X  to  Y  in  gamma.)
# If  arg[4]  is bound, then it is assumed to be a subgroup
# of  Aut(gamma)  stabilizing  X  setwise.
#
local gamma,X,Y;
gamma:=arg[1];
X:=arg[2];
if IsInt(X) then 
   X:=[X];
fi;
Y:=arg[3];
if IsInt(Y) then
   Y:=[Y];
fi;
if not (IsGraph(gamma) and IsList(X) and IsList(Y)) then 
   Error("usage: Distance( <Graph>, <Int> or <List>, ",
        "<Int> or <List> [, <PermGroup> ] )");
fi;
if IsBound(arg[4]) then 
   return LocalInfo(gamma,X,0,Y,arg[4]).distance;
else
   return LocalInfo(gamma,X,0,Y).distance;
fi;
end);

DeclareOperation("Diameter",[IsRecord]);
# to avoid the clash with `Diameter' defined in gap4r5
InstallMethod(Diameter,"for GRAPE graph",[IsRecord],0, 
function(gamma)
#
# Returns the diameter of  gamma. 
# A diameter of  -1  means that gamma is not (strongly) connected.  
#
local r,d,loc,reps;
if not IsGraph(gamma) then
   TryNextMethod();
fi;
if gamma.order=0 then
   Error("<gamma> has no vertices");
fi;
d:=-1;
if IsBound(gamma.autGroup) then
   reps:=GRAPE_OrbitNumbers(gamma.autGroup,gamma.order).representatives;
else 
   reps:=gamma.representatives;
fi;
for r in reps do 
   loc:=LocalInfo(gamma,r);
   if loc.localDiameter=-1 then
      return -1; 
   fi;
   if loc.localDiameter > d then
      d:=loc.localDiameter;
   fi;
od;
return d;
end);

DeclareOperation("Girth",[IsRecord]);
InstallMethod(Girth,"for GRAPE graph",[IsRecord],0, 
function(gamma)
#
# Returns the girth of  gamma,  which must be a simple graph. 
# A girth of  -1  means that gamma is a forest.  
#
local r,g,locgirth,stoplayer,adj,reps;
if not IsGraph(gamma) then
   TryNextMethod();
fi;
if gamma.order=0 then
   return -1;
fi;
if not IsSimpleGraph(gamma) then
   Error("<gamma> not a simple graph");
fi;
adj:=gamma.adjacencies[1];
if adj<>[] and Intersection(adj,Adjacency(gamma,adj[1]))<>[] then
   return 3;
fi;
g:=-1;
stoplayer:=0;
if IsBound(gamma.autGroup) then
   reps:=GRAPE_OrbitNumbers(gamma.autGroup,gamma.order).representatives;
else 
   reps:=gamma.representatives;
fi;
for r in reps do 
   locgirth:=LocalInfo(gamma,r,stoplayer).localGirth;
   if locgirth=3 then 
      return 3;
   fi;
   if locgirth<>-1 then
      if g=-1 or locgirth<g then
   g:=locgirth;
   stoplayer:=Int((g+1)/2)+1; # now no need for  LocalInfo  to create a 
                                     # layer beyond the  stoplayer-th  one.
      fi;
   fi;
od;
return g;
end);

BindGlobal("IsRegularGraph",function(gamma)
#
# Returns  true  iff the graph  gamma  is (out)regular.
#
local deg,i;
if not IsGraph(gamma) then 
   Error("usage: IsRegularGraph( <Graph> )");
fi;
if gamma.order=0 then 
   return true;
fi;
deg:=Length(gamma.adjacencies[1]);
for i in [2..Length(gamma.adjacencies)] do
   if deg <> Length(gamma.adjacencies[i]) then 
      return false;
   fi;
od;
return true;
end);

BindGlobal("IsNullGraph",function(gamma)
#
# Returns  true  iff the graph  gamma  has no edges.
#
local i;
if not IsGraph(gamma) then 
   Error("usage: IsNullGraph( <Graph> )");
fi;
for i in [1..Length(gamma.adjacencies)] do
   if Length(gamma.adjacencies[i])<>0 then 
      return false;
   fi;
od;
return true;
end);

BindGlobal("IsCompleteGraph",function(arg)
#
# Returns  true  iff the graph  gamma=arg[1]  is a complete graph.
# The optional boolean parameter  arg[2]  
# is true iff all loops must exist for  gamma  to be considered
# a complete graph (default: false); otherwise loops are ignored 
# (except to possibly set  gamma.isSimple). 
#
local deg,i,notnecsimple,gamma,mustloops;
gamma:=arg[1];
if IsBound(arg[2]) then 
   mustloops := arg[2];
else
   mustloops := false;
fi;
if not IsGraph(gamma) or not IsBool(mustloops) then 
   Error("usage: IsCompleteGraph( <Graph> [, <Bool> ] )");
fi;
notnecsimple := not IsBound(gamma.isSimple) or not gamma.isSimple;
for i in [1..Length(gamma.adjacencies)] do
   deg := Length(gamma.adjacencies[i]);
   if deg < gamma.order-1 then 
      return false;
   fi;
   if deg=gamma.order-1 then
      if mustloops then
  return false;
      fi;
      if notnecsimple and 
       (gamma.representatives[i] in gamma.adjacencies[i]) then
  gamma.isSimple := false;
  return false;
      fi;
   fi;
od;
return true;
end);

BindGlobal("IsLoopy",function(gamma)
#
# Returns  true  iff graph  gamma  has a loop.
#
local i;
if not IsGraph(gamma) then 
   Error("usage: IsLoopy( <Graph> )");
fi;
for i in [1..Length(gamma.adjacencies)] do
   if gamma.representatives[i] in gamma.adjacencies[i] then
      gamma.isSimple := false;
      return true;
   fi;
od;
return false;
end);

DeclareOperation("IsConnectedGraph",[IsRecord]);
InstallMethod(IsConnectedGraph,"for GRAPE graph",[IsRecord],0, 
function(gamma)
#
# Returns true iff  gamma  is (strongly) connected.
#
if not IsGraph(gamma) then
   TryNextMethod();
fi;
if gamma.order=0 then
   return true;
fi;
if IsSimpleGraph(gamma) then
   return LocalInfo(gamma,1).localDiameter > -1;
else
   return Diameter(gamma) > -1;
fi;
end);

DeclareOperation("ConnectedComponent",[IsRecord,IsPosInt]);
InstallMethod(ConnectedComponent,"for GRAPE graph",[IsRecord,IsPosInt],0, 
function(gamma,v)
#
# Returns the set of all vertices in  gamma  which can be reached by 
# a path starting at vertex  v.  The graph  gamma  must be simple.
#
local comp,laynum;
if not IsGraph(gamma) then
   TryNextMethod();
fi;
if not IsSimpleGraph(gamma) then
   Error("<gamma> not a simple graph");
fi;
if not IsVertex(gamma,v) then
   Error("<v> is not a vertex of <gamma>");
fi;
laynum:=LocalInfo(gamma,v).layerNumbers;
comp:=Filtered([1..gamma.order],j->laynum[j]>0);
IsSSortedList(comp);
return comp;
end);

DeclareOperation("ConnectedComponents",[IsRecord]);
InstallMethod(ConnectedComponents,"for GRAPE graph",[IsRecord],0, 
function(gamma)
#
# Returns the set of the vertex-sets of the connected components
# of  gamma,  which must be a simple graph.
#
local comp,used,i,j,x,cmp,laynum;
if not IsGraph(gamma) then
   TryNextMethod();
fi;
if not IsSimpleGraph(gamma) then
   Error("<gamma> not a simple graph");
fi;
comp:=[]; 
used:=BlistList([1..gamma.order],[]);
for i in [1..gamma.order] do 
   # Loop invariant: used[j]=true for all j<i. 
   if not used[i] then   # new component
      laynum:=LocalInfo(gamma,i).layerNumbers;
      cmp:=Filtered([i..gamma.order],j->laynum[j]>0);
      IsSSortedList(cmp);
      for x in Orbit(gamma.group,cmp,OnSets) do
  Add(comp,x);
  for j in x do 
     used[j]:=true;
  od;
      od;
   fi; 
od;
return SSortedList(comp);
end);

BindGlobal("ComponentLocalInfos",function(gamma)
#
# Returns a sequence of localinfos for the connected components of  
# gamma  (w.r.t. some vertex in each component).
# The graph  gamma  must be simple.
#
local comp,used,i,j,k,laynum;
if not IsGraph(gamma) then 
   Error("usage: ComponentLocalInfos( <Graph> )");
fi;
if not IsSimpleGraph(gamma) then
   Error("<gamma> not a simple graph");
fi;
comp:=[]; 
used:=BlistList([1..gamma.order],[]);
k:=0;
for i in [1..gamma.order] do 
   if not used[i] then   # new component
      k:=k+1; 
      comp[k]:=LocalInfo(gamma,i);
      laynum:=comp[k].layerNumbers;
      for j in [1..gamma.order] do
  if laynum[j] > 0 then
     used[j]:=true;
  fi;
      od;
   fi; 
od;
return comp;
end);

BindGlobal("Bicomponents",function(gamma)
#
# If  gamma  is bipartite, returns a length 2 list of
# bicomponents, or parts, of  gamma,  else returns the empty list.
# *** This function is for simple  gamma  only.
#
# Note: if gamma.order=0 this function returns [[],[]], and if 
# gamma.order=1 this function returns [[],[1]] (unlike GRAPE 2.2
# which returned [], which was inconsistent with considering 
# a zero vertex graph to be bipartite).
#
local bicomps,i,lnum,loc,locs;
if not IsGraph(gamma) then 
   Error("usage: Bicomponents( <Graph> )");
fi;
if not IsSimpleGraph(gamma) then
   Error("<gamma> not a simple graph");
fi;
bicomps:=[[],[]]; 
if gamma.order=0 then 
   return bicomps;
fi;
if IsNullGraph(gamma) then 
   return [[1..gamma.order-1],[gamma.order]];
fi;
locs:=ComponentLocalInfos(gamma);
for loc in locs do
   for i in [2..Length(loc.localParameters)] do
      if loc.localParameters[i][2]<>0 then
         #  gamma  not bipartite.
  return [];
      fi;
   od;
   for i in [1..Length(loc.layerNumbers)] do 
      lnum:=loc.layerNumbers[i];
      if lnum>0 then
  if lnum mod 2 = 1 then
     AddSet(bicomps[1],i);
  else
     AddSet(bicomps[2],i);
  fi;
      fi; 
   od;
od;
return bicomps;
end);

DeclareOperation("IsBipartite",[IsRecord]);
InstallMethod(IsBipartite,"for GRAPE graph",[IsRecord],0, 
function(gamma)
#
# Returns  true  iff  gamma  is bipartite. 
# *** This function is only for simple  gamma.
#
# Note: Now the one vertex graph is considered to be bipartite 
# (as well as the zero vertex graph). This is a change from the inconsistent 
# GRAPE 2.2 view that a zero vertex graph is bipartite, but not a one 
# vertex graph.
#
if not IsGraph(gamma) then
   TryNextMethod();
fi;
if not IsSimpleGraph(gamma) then
   Error("<gamma> not a simple graph");
fi;
return Length(Bicomponents(gamma))=2;
end);

BindGlobal("Layers",function(arg)
#
# Returns the list of vertex layers of  gamma=arg[1],  
# starting from  V=arg[2],  which may be a vertex list or singleton vertex. 
# Layers[i]  is the set of vertices at distance  i-1  from  V.
# If  arg[3]  is bound then it is assumed to be a subgroup 
# of  Aut(gamma)  stabilizing  V  setwise.
#
local gamma,V;
gamma:=arg[1];
V:=arg[2];
if IsInt(V) then
   V:=[V];
fi;
if not IsGraph(gamma) or not IsList(V) 
                      or (IsBound(arg[3]) and not IsPermGroup(arg[3])) then 
   Error("usage: Layers( <Graph>, <Int> or <List>, [, <PermGroup>] )");
fi;
if IsBound(arg[3]) then
   return GRAPE_NumbersToSets(LocalInfo(gamma,V,0,[],arg[3]).layerNumbers);
else
   return GRAPE_NumbersToSets(LocalInfo(gamma,V).layerNumbers);
fi;
end);

BindGlobal("LocalParameters",function(arg)
#
# Returns the local parameters of simple, connected  gamma=arg[1],  
# w.r.t to vertex list (or singleton vertex)  V=arg[2].
# The nonexistence of a local parameter is denoted by  -1.
# If  arg[3]  is bound then it is assumed to be a subgroup 
# of  Aut(gamma)  stabilizing  V  setwise.
#
local gamma,V,loc;
gamma:=arg[1];
V:=arg[2];
if IsInt(V) then
   V:=[V];
fi;
if not IsGraph(gamma) or not IsList(V) 
                      or (IsBound(arg[3]) and not IsPermGroup(arg[3])) then 
   Error("usage: LocalParameters( <Graph>, <Int> or <List>, [, <PermGroup>] )");
fi;
if not IsSimpleGraph(gamma) then
   Error("<gamma> not a simple graph");
fi;
if Length(V)>1 and not IsConnectedGraph(gamma) then
   Error("<gamma> not a connected graph");
fi;
if IsBound(arg[3]) then
   loc:=LocalInfo(gamma,V,0,[],arg[3]);
else
   loc:=LocalInfo(gamma,V);
fi;
if loc.localDiameter=-1 then
   Error("<gamma> not a connected graph");
fi;
return loc.localParameters;
end);

BindGlobal("GlobalParameters",function(gamma)
#
# Determines the global parameters of connected, simple graph  gamma.
# The nonexistence of a global parameter is denoted by  -1.
#
local i,j,k,reps,pars,lp,loc;
if not IsGraph(gamma) then 
   Error("usage: GlobalParameters( <Graph> )");
fi;
if gamma.order=0 then
   return [];
fi;
if not IsSimpleGraph(gamma) then
   Error("<gamma> not a simple graph");
fi;
if IsBound(gamma.autGroup) then
   reps:=GRAPE_OrbitNumbers(gamma.autGroup,gamma.order).representatives;
else 
   reps:=gamma.representatives;
fi;
loc:=LocalInfo(gamma,reps[1]);
if loc.localDiameter=-1 then
   Error("<gamma> not a connected graph");
fi;
pars:=loc.localParameters;
for i in [2..Length(reps)] do
   lp:=LocalInfo(gamma,reps[i]).localParameters;
   for j in [1..Maximum(Length(lp),Length(pars))] do
      if not IsBound(lp[j]) or not IsBound(pars[j]) then
  pars[j]:=[-1,-1,-1];
      else
  for k in [1..3] do
     if pars[j][k]<>lp[j][k] then
        pars[j][k]:=-1;
     fi;
  od;
      fi;
   od;
od;
return pars;
end);

DeclareOperation("IsDistanceRegular",[IsRecord]);
InstallMethod(IsDistanceRegular,"for GRAPE graph",[IsRecord],0, 
function(gamma)
#
# Returns  true  iff  gamma  is distance-regular 
# (a graph must be simple to be distance-regular).
#
local i,reps,pars,lp,loc,d;
if not IsGraph(gamma) then
   TryNextMethod();
fi;
if gamma.order=0 then
   return true;
fi;
if not IsSimpleGraph(gamma) then
   return false;
fi;
if IsBound(gamma.autGroup) then
   reps:=GRAPE_OrbitNumbers(gamma.autGroup,gamma.order).representatives;
else 
   reps:=gamma.representatives;
fi;
loc:=LocalInfo(gamma,reps[1]);
pars:=loc.localParameters;
d:=loc.localDiameter;
if d=-1 then  # gamma not connected
   return false;
fi;
if -1 in Flat(pars) then 
   return false;
fi;
for i in [2..Length(reps)] do
   loc:=LocalInfo(gamma,reps[i]);
   if loc.localDiameter<>d then 
      return false;
   fi;
   if pars <> loc.localParameters then
      return false;
   fi;
od;
return true;
end);

BindGlobal("DistanceSet",function(arg)
#
# Let  gamma=arg[1],  distances=arg[2],  V=arg[3].
# Returns the set of vertices  w  of  gamma,  such that  d(V,w)  is in
# distances (a list or singleton distance). 
# If  arg[4]  is bound, then it is assumed to be a subgroup
# of  Aut(gamma)  stabilizing  V  setwise.
#
local gamma,distances,V,maxlayer,distset,laynum,x,i;
gamma:=arg[1];
distances:=arg[2];
V:=arg[3];
if IsInt(distances) then   # assume  distances  consists of a single distance.
   distances:=[distances];
fi;
if not (IsGraph(gamma) and IsList(distances) and (IsList(V) or IsInt(V))) then
   Error("usage: DistanceSet( <Graph>, <Int> or <List>, ",
        "<Int> or <List> [, <PermGroup> ] )");
fi;
if not IsSSortedList(distances) then 
   distances:=SSortedList(distances);
fi;
distset:=[];
if Length(distances)=0 then
   return distset;
fi;
maxlayer:=Maximum(distances)+1;
if IsBound(arg[4]) then
   laynum:=LocalInfo(gamma,V,maxlayer,[],arg[4]).layerNumbers;
else
   laynum:=LocalInfo(gamma,V,maxlayer).layerNumbers;
fi;
for i in [1..gamma.order] do
   if laynum[i]-1 in distances then
      Add(distset,i);
   fi;
od;
IsSSortedList(distset);
return distset;
end);

BindGlobal("DistanceSetInduced",function(arg)
#
# Let  gamma=arg[1],  distances=arg[2],  V=arg[3].
# Returns the graph induced on the set of vertices  w  of  gamma,  
# such that  d(V,w)  is in distances (a list or singleton distance). 
# If  arg[4]  is bound, then it is assumed to be a subgroup
# of  Aut(gamma)  stabilizing  V  setwise.
#
local gamma,distances,V,distset,H;
gamma:=arg[1];
distances:=arg[2];
V:=arg[3];
if IsInt(distances) then   # assume  distances  consists of a single distance.
   distances:=[distances];
fi;
if IsInt(V) then
   V:=[V];
fi;
if IsBound(arg[4]) then
   H:=arg[4];
elif Length(V)=1 then
   H:=ProbablyStabilizer(gamma.group,V[1]);
else
   H:=Group([],());
fi;
if not (IsGraph(gamma) and IsList(distances) and IsList(V) and IsPermGroup(H))
  then
   Error("usage: DistanceSetInduced( <Graph>, ",
  "<Int> or <List>, <Int> or <List> [, <PermGroup> ] )");
fi;
distset:=DistanceSet(gamma,distances,V,H);
return InducedSubgraph(gamma,distset,H);
end);

BindGlobal("DistanceGraph",function(gamma,distances)
#
# Returns graph  delta  with the same vertex-set, names, and group as 
# gamma,  and  [x,y]  is an edge of  delta  iff  d(x,y)  (in gamma)
# is in  distances. 
#
local r,delta,d,i;
if IsInt(distances) then
   distances:=[distances];
fi;
if not IsGraph(gamma) or not IsList(distances) then 
   Error("usage: DistanceGraph( <Graph>, <Int> or <List> )");
fi;
delta:=rec(isGraph:=true,order:=gamma.order,group:=gamma.group,
    schreierVector:=Immutable(gamma.schreierVector),adjacencies:=[],
    representatives:=Immutable(gamma.representatives));
if IsBound(gamma.names) then
   delta.names:=Immutable(gamma.names);
fi;
for i in [1..Length(delta.representatives)] do 
   delta.adjacencies[i]:=DistanceSet(gamma,distances,delta.representatives[i]);
od;
if not (0 in distances) and IsBound(gamma.isSimple) and gamma.isSimple then
   delta.isSimple:=true;
fi;
return delta;
end);

BindGlobal("ComplementGraph",function(arg)
#
# Returns the complement of the graph  gamma=arg[1]. 
# arg[2] is true iff loops/nonloops are to be complemented (default:false).
#
local gamma,comploops,i,delta,notnecsimple;
gamma:=arg[1];
if IsBound(arg[2]) then
   comploops:=arg[2];
else
   comploops:=false;
fi;
if not IsGraph(gamma) or not IsBool(comploops) then 
   Error("usage: ComplementGraph( <Graph> [, <Bool> ] )");
fi;
notnecsimple:=not IsBound(gamma.isSimple) or not gamma.isSimple;
delta:=rec(isGraph:=true,order:=gamma.order,group:=gamma.group,
    schreierVector:=Immutable(gamma.schreierVector),adjacencies:=[],
    representatives:=Immutable(gamma.representatives));
if IsBound(gamma.names) then
   delta.names:=Immutable(gamma.names);
fi;
if IsBound(gamma.autGroup) then
   delta.autGroup:=gamma.autGroup;
fi;
if IsBound(gamma.isSimple) then
   if gamma.isSimple and not comploops then
      delta.isSimple:=true;
   fi;
fi;
for i in [1..Length(delta.representatives)] do 
   delta.adjacencies[i]:=Difference([1..gamma.order],gamma.adjacencies[i]);
   if not comploops then
      RemoveSet(delta.adjacencies[i],delta.representatives[i]);
      if notnecsimple and (gamma.representatives[i] in gamma.adjacencies[i])
       then
  AddSet(delta.adjacencies[i],delta.representatives[i]);
      fi;
   fi;
od;
return delta;
end);

BindGlobal("PointGraph",function(arg)
#
# Assuming that  gamma=arg[1]  is simple, connected, and bipartite, 
# this function returns the connected component containing  
# v=arg[2]  of the distance-2  graph of  gamma=arg[1]  
# (default:  arg[2]=1,  unless  gamma has zero
# vertices, in which case a zero vertex graph is returned). 
# Thus, if  gamma  is the incidence graph of a (connected) geometry, and 
# v  represents a point, then the point graph of the geometry is returned.
#
local gamma,delta,bicomps,comp,v,gens,hgens,i,g,j,outer;
gamma:=arg[1];
if IsBound(arg[2]) then 
   v:=arg[2];
else
   v:=1;
fi;
if not IsGraph(gamma) or not IsInt(v) then
   Error("usage: PointGraph( <Graph> [, <Int> ])");
fi;
if gamma.order=0 then
   return CopyGraph(gamma);
fi;
bicomps:=Bicomponents(gamma);
if Length(bicomps)=0 or not IsSimpleGraph(gamma) 
       or not IsConnectedGraph(gamma) then
   Error("<gamma> not  simple,connected,bipartite");
fi;
if v in bicomps[1] then 
   comp:=bicomps[1];
else
   comp:=bicomps[2];
fi;
delta:=DistanceGraph(gamma,2);
# construct Schreier generators for the subgroup of  gamma.group 
# fixing  comp.
gens:=GeneratorsOfGroup(gamma.group);
hgens:=[];
for i in [1..Length(gens)] do
   g:=gens[i];
   if v^g in comp then
      AddSet(hgens,g);
      if IsBound(outer) then
  AddSet(hgens,outer*g/outer);
      fi;
   else    # g is an "outer" element
      if IsBound(outer) then
  AddSet(hgens,g/outer);
  AddSet(hgens,outer*g);
      else 
  outer:=g;
  for j in [1..i-1] do
     AddSet(hgens,outer*gens[j]/outer);
  od;
  g:=g^2;
  if g <> () then
     AddSet(hgens,g);
  fi;
      fi;
   fi;
od;
return InducedSubgraph(delta,comp,Group(hgens,()));
end);

BindGlobal("EdgeGraph",function(gamma)
#
# Returns the edge graph, also called the line graph, of the 
# (assumed) simple graph  gamma.
# This edge graph  delta  has the unordered edges of  gamma  as 
# vertices, and  e  is joined to  f  in delta precisely when 
# e<>f,  and  e,f  have a common vertex in  gamma.
#
local delta,i,j,k,edgeset,adj,r,e,f;
if not IsGraph(gamma) then 
   Error("usage: EdgeGraph( <Graph> )");
fi;
if not IsSimpleGraph(gamma) then
   Error("<gamma> not a simple graph");
fi;
edgeset:=UndirectedEdges(gamma);
delta:=NullGraph(Action(gamma.group,edgeset,OnSets),Length(edgeset));
delta.names:=Immutable(List(edgeset,e->List(e,i->VertexName(gamma,i))));
for i in [1..Length(delta.representatives)] do
   r:=delta.representatives[i];
   e:=edgeset[r];
   adj:=delta.adjacencies[i];
   for k in [1,2] do
      for j in Adjacency(gamma,e[k]) do
  f:=SSortedList([e[k],j]);
  if e<>f then
     AddSet(adj,PositionSet(edgeset,f));
  fi;
      od;
   od;
od;
return delta;
end);

BindGlobal("QuotientGraph",function(gamma,R)
#
# Returns the quotient graph  delta  of  gamma  defined by the 
# smallest  gamma.group  invariant equivalence relation  S 
# (on the vertices of  gamma)  containing the relation  R  (given
# as a list of ordered pairs of vertices of  gamma).
# The vertices of this quotient  delta  are the equivalence 
# classes of  S,  and  [X,Y]  is an edge of  delta  iff  
# [x,y]  is an edge of  gamma  for some  x in X,  y in Y.
#
local root,Q,F,V,W,i,j,r,q,x,y,names,gens,delta,g,h,m,pos;

root := function(x)
#
# Returns the root of the tree containing  x  in the forest represented 
# by  F,  and compresses the path travelled in this tree to find the root.
# F[x]=-x  if  x  is a root, else  F[x]  is the parent of  x.
#
local t,u;
t:=F[x];
while t>0 do 
   t:=F[t];
od;
# compress path
u:=F[x];
while u<>t do
   F[x]:=-t;
   x:=u;
   u:=F[x];
od;  
return -t;
end;

if not IsGraph(gamma) or not IsList(R) then
   Error("usage: QuotientGraph( <Graph>, <List> )");
fi;
if gamma.order<=1 or Length(R)=0 then
   delta:=CopyGraph(gamma);
   delta.names:=Immutable(List([1..gamma.order],i->[VertexName(gamma,i)]));
   return delta;
fi;
if IsInt(R[1]) then   # assume  R  consists of a single pair.
   R:=[R];
fi;
F:=[];
for i in [1..gamma.order] do
   F[i]:=-i;
od;
Q:=[];
for r in R do
   x:=root(r[1]);
   y:=root(r[2]);
   if x<>y then
      if x>y then
  Add(Q,x);
  F[x]:=y;
      else
  Add(Q,y);
  F[y]:=x;
      fi;
   fi;
od;
for q in Q do 
   for g in GeneratorsOfGroup(gamma.group) do
      x:=root(F[q]^g);
      y:=root(q^g);
      if x<>y then
  if x>y then
     Add(Q,x);
     F[x]:=y;
  else
     Add(Q,y);
     F[y]:=x;
  fi;
      fi;
   od;
od;
for i in Reversed([1..gamma.order]) do
   if F[i] < 0 then
      F[i]:=-F[i];
   else
      F[i]:=root(F[i]);
   fi;
od;
V:=SSortedList(F);
W:=[];
names:=[];
m:=Length(V);
for i in [1..m] do
   W[V[i]]:=i;
   names[i]:=[];
od;
for i in [1..gamma.order] do
   Add(names[W[F[i]]],VertexName(gamma,i));
od; 
gens:=[];
for g in GeneratorsOfGroup(gamma.group) do
   h:=[];
   for i in [1..m] do
      h[i]:=W[F[V[i]^g]];
   od;
   Add(gens,PermList(h));
od;
delta:=NullGraph(Group(gens,()),m);
delta.names:=Immutable(names);
IsSSortedList(delta.representatives);
for i in [1..gamma.order] do
   pos:=PositionSet(delta.representatives,W[F[i]]);
   if pos<>fail then
      for j in Adjacency(gamma,i) do
  AddSet(delta.adjacencies[pos],W[F[j]]);
      od;
   fi;
od;
if IsLoopy(delta) then
   delta.isSimple:=false;
elif IsBound(gamma.isSimple) and gamma.isSimple then
   delta.isSimple:=true;
else
   Unbind(delta.isSimple);
fi;
return delta;
end);
 
--> --------------------

--> maximum size reached

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