On 4-regular square-complementary graphs of large girth (full code)

# This file name: squcodg.g
#
# This code searches exhaustively for square-complementary (squco) graphs that are
# also d-regular and of girth >=5 (equivalently, d-regular and of order d^2+d+1).
# This cannot happen for odd d, since the graph would have an odd number of vertices 
# of odd degree. For d=2, it does find the 7-cycle (the only such graph) and for d=4 
# it does not find any, thus proving that there are no 4-regular squco graphs of 
# girth 5, which solves an open problem in Discrete Mathematics 327 (2014) 62-75. 
# For d = 6, it is still impractical to do an exhaustive search. 
#
# (C) 2020 Ariadna Juárez-Valencia and Miguel Pizaña.
# License: GPLv3 (https://www.gnu.org/licenses/gpl-3.0.html).
#
# Requisites: GAP 4.10.2+ (https://www.gap-system.org/), 
#             YAGS 0.0.5+ (http://xamanek.izt.uam.mx/yags/)
#
# Usage: Put this file in your GAP's working directory, and at the gap prompt type: 
#
# Read("squcodg.g"); RegularSqucoOfLargeGirth([]);

RequirePackage("yags");
SetInfoLevel(YAGSInfo.InfoClass,1);  # For progress reporting. 
#SetInfoLevel(YAGSInfo.InfoClass,3); # Use this for full progress reporting.

#################################
# Searching for d-regular graphs.
# Reload this file whenever you 
# change this value.
d:=4; 
#################################

# The square-complement of a graph
SquCo:=function(g) return ComplementGraph(g^2); end;

# Returns true when g is square-complementary.
IsSquCo:=function(g) 
  local g1;
  g1:=SquCo(g);
  return IsIsomorphicGraph(g,g1); 
end;

# This is the initial scaffolding, a graph on d^2+d+1 vertices 
# consisting of a tree on d^2+1 vertices with all internal 
# vertices of degree d, and d additional isolated vertices. 
# All the d-regular squco graphs of girth at least 5 must contain 
# this substructure.
H:=TreeGraph([d,d-1]);
H:=DisjointUnion(H,DiscreteGraph(d));
if d=4 then 
 SetCoordinates(H,
  [ [ 2, -279 ], [ -281, 77 ], [ -152, -149 ], [ 170, -161 ], [ 284, 69 ], 
  [ -164, 187 ], [ -194, 115 ], [ -195, 38 ], [ -166, -33 ], [ -111, -89 ], 
  [ -39, -119 ], [ 37, -120 ], [ 109, -91 ], [ 165, -36 ], [ 195, 35 ], 
  [ 196, 112 ], [ 166, 184 ], [ -109, 242 ], [ -36, 271 ], [ 40, 270 ], 
  [ 112, 240 ] ]);#this is for a nice drawing.
fi;

# Candidate edges for completing the sought graphs. These are essentially
# Combinations([d+2..d^2+d+1],2), that is, all the possible edges among the 
# leaves of the original tree and the d extra isolated vertices, except for 
# some edges that form triangles (thus girth<5) with the initial scaffolding. 
# The order was chosen heuristically for efficiency in the search.
U:=Concatenation(Combinations([d^2+2 .. d^2+d+1],2),
                    List(Cartesian([d^2+2..d^2+d+1],[d+2..d^2+1]),Set),
                    Combinations([d+2..d^2+1],2));
# These are the edges that form triangles with the scaffolding:
forbidden:= [];
for k in [1..d] do
  Append(forbidden,Combinations([(d+2)+(k-1)*(d-1)..2*d+(k-1)*(d-1)],2));
od;
# Triangle-forming edges are removed from the candidate edges. 
U:=Filtered(U, e-> not e in forbidden);

# In the case d=4, the candidate edges are these:
# U:=[ 
#   [ 18, 19 ], [ 18, 20 ], [ 18, 21 ], [ 19, 20 ], [ 19, 21 ], 
#   [ 20, 21 ], [ 6, 18 ],  [ 7, 18 ],  [ 8, 18 ],  [ 9, 18 ],  #10 
#   [ 10, 18 ], [ 11, 18 ], [ 12, 18 ], [ 13, 18 ], [ 14, 18 ], 
#   [ 15, 18 ], [ 16, 18 ], [ 17, 18 ], [ 6, 19 ],  [ 7, 19 ],  #20 
#   [ 8, 19 ],  [ 9, 19 ],  [ 10, 19 ], [ 11, 19 ], [ 12, 19 ], 
#   [ 13, 19 ], [ 14, 19 ], [ 15, 19 ], [ 16, 19 ], [ 17, 19 ], #30 
#   [ 6, 20 ],  [ 7, 20 ],  [ 8, 20 ],  [ 9, 20 ],  [ 10, 20 ], 
#   [ 11, 20 ], [ 12, 20 ], [ 13, 20 ], [ 14, 20 ], [ 15, 20 ], #40 
#   [ 16, 20 ], [ 17, 20 ], [ 6, 21 ],  [ 7, 21 ],  [ 8, 21 ], 
#   [ 9, 21 ],  [ 10, 21 ], [ 11, 21 ], [ 12, 21 ], [ 13, 21 ], #50 
#   [ 14, 21 ], [ 15, 21 ], [ 16, 21 ], [ 17, 21 ], [ 6, 9 ], 
#   [ 6, 10 ],  [ 6, 11 ],  [ 6, 12 ],  [ 6, 13 ],  [ 6, 14 ],  #60 
#   [ 6, 15 ],  [ 6, 16 ],  [ 6, 17 ],  [ 7, 9 ],   [ 7, 10 ], 
#   [ 7, 11 ],  [ 7, 12 ],  [ 7, 13 ],  [ 7, 14 ],  [ 7, 15 ],  #70 
#   [ 7, 16 ],  [ 7, 17 ],  [ 8, 9 ],   [ 8, 10 ],  [ 8, 11 ], 
#   [ 8, 12 ],  [ 8, 13 ],  [ 8, 14 ],  [ 8, 15 ],  [ 8, 16 ],  #80 
#   [ 8, 17 ],  [ 9, 12 ],  [ 9, 13 ],  [ 9, 14 ],  [ 9, 15 ], 
#   [ 9, 16 ],  [ 9, 17 ],  [ 10, 12 ], [ 10, 13 ], [ 10, 14 ], #90 
#   [ 10, 15 ], [ 10, 16 ], [ 10, 17 ], [ 11, 12 ], [ 11, 13 ], 
#   [ 11, 14 ], [ 11, 15 ], [ 11, 16 ], [ 11, 17 ], [ 12, 15 ], #100 
#   [ 12, 16 ], [ 12, 17 ], [ 13, 15 ], [ 13, 16 ], [ 13, 17 ], 
#   [ 14, 15 ], [ 14, 16 ], [ 14, 17 ] ];

# Auxiliary function that transforms lists into permutations 
# e.g.: [2,5,6,8] -> (2,5,6,8).
permcycle:=function(L) 
  local M,i;
  if Length(L)<=1 then return (); fi;
  M:=[1..Maximum(L)];
  for i in [1..Length(L)-1] do 
    M[L[i]]:=L[i+1];
  od;
  M[L[Length(L)]]:= L[1];
  return PermList(M);
end;

# We define a list of groups, 'grp', that acts on the d^2+d+1 vertices and 
# leave the scaffolding invariant, hence those are subgroups of the automorphism 
# group of the scaffolding. Actually, only the last d^2 vertices 
# (i.e. [d+2 .. d^2+d+1]) are considered in these permutations, since will only 
# use them to act on the candidate edges. 
# 'grp[1]' permutes the d additional isolated vertices freely.
# 'grp[k]' for k in [2 .. d+1] permutes the (k-1)-th bunch of sibling leaves 
#          of the tree.
# 'grp[d+2]' permutes the d bunches of leaves among themselves.
#
# The corresponding 'grpe[k]' groups, are the inherited permutation groups acting 
# on the indices (positions) of the candidate edges 'U'.
grp:=[];
grpe:=[];
grp[1]:=SymmetricGroup([d^2+2..d^2+d+1]);
grpe[1]:=Action(grp[1],U,OnSets);
for k in [1..d] do 
    grp[k+1]:=SymmetricGroup([(d+2)+(k-1)*(d-1)..2*d+(k-1)*(d-1)]);
    grpe[k+1]:=Action(grp[k+1],U,OnSets);
od;
perm1:=();
perm2:=();
for k in [1..d-1] do 
  perm1:=perm1*permcycle([d+1+k, d+d+k .. d+1+(d-1)*(d-1)+k]);
  perm2:=perm2*permcycle([d+1+k,d+d+k]);
od;
grp[d+2]:=Group([perm1,perm2]);
grpe[d+2]:=Action(grp[d+2],U,OnSets);

# These are the groups 'grp' in the case d=4.
# grp[1]:=Group([(18,19,20,21),(18,19)]);
# grp[2]:=Group([(6,7,8),(6,7)]);
# grp[3]:=Group([(9,10,11),(9,10)]);
# grp[4]:=Group([(12,13,14),(12,13)]);
# grp[5]:=Group([(15,16,17),(15,16)]);
# grp[6]:=Group([(6,9)(7,10)(8,11),(6,9,12,15)(7,10,13,16)(8,11,14,17)]);

# Number of candidate edges.
NU:=Length(U);

# Checking function for the backtracking. It must return 'true' if the 
# current L is a feasible solution (still has a chance to be completed 
# to the sought example). Otherwise, it returns 'false'. L is the list 
# of (indices of) the candidate edges currently being considered to extend 
# the initial scaffolding 'H'. 
chk:=function(L,H) 
  local x,i,len,G,last,remaining,Lmin;
  len:=Length(L);
  if len=0 then return true; fi; # Empty list is OK.
  last:=L[len];
  # L is a set of (indices of) edges, so we can assume 
  # it is given in increasing order:
  if len>=2 and last <= L[len-1] then return false; fi; 

  G:=AddEdges(H,U{L});# Scaffolding + edges described by L.

  # Degree d+1 is not allowed:
  if MaxDegree(G)>d then return false; fi;
  # Make sure that there are still enough candidate edges for the vertices 
  # to achieve the desired degree d:
  for x in [d+2..d^2+d+1] do 
    remaining:=Length(Filtered(U{[last+1..NU]},z->x in z));
    if remaining+VertexDegree(G,x)<d then return false; fi;
  od; 
  # Cycles of length <5 are not allowed.
  if Girth(G)<5 then return false; fi;
  # Check whether the current configuration L has already been considered in a 
  # previous isomorphic configuration. Here, 'previous' means in lexicographic 
  # order which is also the searching order.
  for i in [1..Length(grp)] do 
    Lmin:=Minimum(Orbit(grpe[i],L,OnSets));
    if Lmin<L then 
       # Reports when a previously considered isomorphic configuration is found:
      Info(YAGSInfo.InfoClass,1,L,"~",Lmin,"\n"); 
      return false;
    fi; 
  od;
  # Otherwise it is OK to keep trying.
  return true;
end;

# Tests whether the current configuration is already a final solution.
done:=function(L,H) 
    local G;
    G:=AddEdges(H,U{L});
    if MinDegree(G)<d then return false; fi;
    return IsSquCo(G);
end;

# This function does all the job with L=[] (initial configuration). This returns 
# the set of (indices of) edges that added to the scaffolding, produces the 
# sought graph.  Or it returns 'fail' if no such graph exists.
# With a different L, it can start searching at the given configuration L.
# Useful when we need to restart the computation or when we want to 
# parallelize the search using several GAP sessions.
RegularSqucoOfLargeGirth:=function(L) 
   return Backtrack(L,[1..NU],chk,done,H);
end;