Constructing Hamiltonian Cycles on a Square Grid

Background

A solution to New Scientist Enigma puzzle 1431 led me to consider how to generate all Hamiltonian Cycles for an N by N square grid.

The obvious way to generate Hamiltonian cycles for any graph is to find all Hamiltonian paths that end at a node adjacent to the start node. However, for a 6 by 6 grid, there are 458,696 Hamiltonian paths but only 1,072 Hamiltonian cycles. (Sloane A096969 and A003763)

The analysis in  R. Stoyan and V. Strehl: Enumeration of Hamiltonian Circuits in Rectangular Grids. Journal of Combinatorial Mathematics and Combin.Computing 21 (1996), 197-127, identifies an alternative way of identifying Hamiltonian cycles. For example, this figure shows a Hamiltonian Cycle for a 6 by 6 grid.

Stoyan and Strehl point out that each cell is either inside or outside the Hamiltonian Cycle, as shown here (yellow cells inside, white outside) Thus, Hamiltonian Cycles for an N by N graph can be associated with an "induced subgraph" of an (N-1) by (N-1) square graph.

Stoyan and Strehl show that any induced subgraph corresponds to a Hamiltonian Cycle if: A) It is a tree B) No 2 by 2 sub-grid of the extended grid matches one of the four patterns:

  Additionally, there are always $\frac{N^2-2}{2}$ inside cells and $\frac{(N-2)^2}{2}$ outside cells, and the four corner cells have to be inside.

These constrants can be used to find subgraphs that correspond to Hamiltonian cycles more efficiently than the search of Hamiltonian paths. For N=6, there are $\binom{(N-1)^2-4}{\frac{(N-2)^2}{2}} = \binom{21}{8} = 203,490$ combinations of outside cells that need to be tested.

However, the algorithm does not scale well. For N=8, there are $\binom{45}{18} =  1,715,884,494,940$ combinations of outside cells.

The performance could be improved by using the observation that there cannot be 2 contiguous outside cells on the boundary. However the algorithm as it stands runs in less than 1 second for N=6, and the suggested improvement would not significantly reduce the number of cases that would need to be tested for N=8.

Python Implementation

A Python implementation of an algorithm based on these constraints is presented below. The search for the four patterns of 2 by 2 sub-cells does not cover the extended grid, as it turns out that this is not necessary, at least for the cases N=4 and N=6 (I haven't worked out whether this extends to higher values of N)



The Python program is composed of a number of functions:

graph(m) constructs a rectangular m by m grid graph. Each node is represented by a number $0..m^2-1$. A dictionary defines the nodes that each node is connected to.

valid(gr,M,outside,oppmask,perimeter) determines whether a configuration of outside cells satisfies the conditions described above, by testing each 2 by 2 array of subcells, and by testing that the outside cells are all connected to an outside cell on the perimeter. oppmask is the pre-calculated set of 2 by 2 subcells

hamiltonian(outside,M) converts a valid configuration of outside cells on an M by M grid to the equivalent Hamiltonian cycle on an M+1 by M+1 grid.

Hamiltonians(N) yields all the Hamiltonian cycles for an N by N grid.

If the program is run stand-alone, it calculates the Hamiltonian cycles for N=4 and N=6, and prints randomly selected cycles and induced subgraphs. Here are some examples of induced subgraphs for N=6 : 

    

from random import randint
from itertools import combinations as comb

    
  
def graph(m): # define a rectangular grid graph m by m
  g = { n:set([n+1,n-m,n+m]) & set(xrange(m*m)) for n in xrange(0,m*m,m)}
  g.update( { n:set([n-1,n+1,n-m,n+m]) & set(xrange(m*m))
            for n in xrange(m*m) if 0 < n%m < m-1  })
  g.update({ n:set([n-1,n-m,n+m]) & set(xrange(m*m))
            for n in xrange(m-1,m*m,m)})
  return g

def valid(gr,M,outside,oppmasks,perimeter):
  
  for op in oppmasks:
    s = op.mask & outside
    if len(s) in (0,4) \
       or s == op.oppsA \
       or s == op.oppsB:
      return False

  def find_outside_connected(graph, start, pathset):
    pathset.add(start)
    for node in graph[start]:
      if node in outside and node not in pathset:
        find_outside_connected(graph, node, pathset)

  outside_connected=set()
  for s in perimeter & outside:
    find_outside_connected(gr, s, outside_connected)

  if len(outside_connected)<((M-1)*(M-1))/2:
    return False
  
  return True
  
def hamiltonian(outside,M):
  inside = set(range(M*M))-set(outside)
  # convert to an M+1 by M+1 array to trace the Hamiltonian
  converted = set([x+x/M for x in inside])
  N=M+1
  ham=[0,1]
  pos=1
  move=1
  while pos != 0:

    if move==1: # going right, 
      if pos-N in converted: move=-N
      elif pos in converted: move=+1
      else: move=+N  
        
    elif move==N: # going down, 
      if pos in converted: move=1
      elif pos-1 in converted: move=+N
      else: move=-1  

    elif move== -N: # going up, 
      if pos-N-1 in converted: move=-1
      elif pos-N in converted: move=-N
      else: move=1  

    elif move== -1: # going left, 
      if pos-1 in converted: move=N
      elif pos-N-1 in converted: move=-1
      else: move=-N  
        
    else:
      raise Exception("no more cases")
     
    pos+=move  
    ham.append(pos)

  if len(ham) != N*N+1 or len(set(ham)) != N*N:
    print ham
    raise Exception("Invalid Hamiltonian")

  return ham
  

def Hamiltonians(N):
  M=N-1
  gr=graph(M)
  noncorners = set(range(M*M))-set((0,M-1,M*(M-1),M*M-1))
  allnodes=set(range(M*M))
  centre = set([x for x in xrange(M*M) if x/M not in (0,M-1) and x%M not in (0,M-1)])
  perimeter = allnodes-centre

  class Twobytwo: # class to hold the masks for 2 by 2 tests
    pass
  oppmasks=[]
  for i in xrange(M-1):
    for j in xrange(M-1):
      m=Twobytwo()
      m.mask = set((M*i+j,M*i+j+1,M*(i+1)+j,M*(i+1)+j+1))
      m.oppsA = set((M*i+j,M*(i+1)+j+1))
      m.oppsB = set((M*i+j+1,M*(i+1)+j))
      oppmasks.append(m)

  for outside in comb(noncorners,((M-1)*(M-1))/2 ):
    if valid(gr, M, set(outside),oppmasks, perimeter):
      h =  hamiltonian(outside,M)

      if __name__=="__main__" and  randint(0,99)%100==0:
        s="\n"
        for i in xrange(M):
          for j in xrange(M):
            s+= "  " if M*i+j in outside else "* "
          s+="\n"
        print s
        print h

      yield h  

if __name__=="__main__":
  from time import clock
  start= clock()
  for N in (4,6):
    solutions = [h for h in Hamiltonians(N)]
    print "N=",N,":", len(solutions),"solutions"
    print "N=",N,":", len(set([tuple(x) for x in solutions]))  , "distinct solutions" 

  print clock()-start, "seconds"