# Connecting the Dots: Counting n²-gons in the n×n Grid

The recent volume of MAA's Mathematics Magazine had an article "Connecting the Dots: Maximal Polygons on a Square Grid" by Sam Chow, Ayla Gafni, and Paul Gafni about making (very convex) $$\n^2\$$-gons where each vertex is a different point of the $$\n \times n\$$ grid.

One is not allowed to have two consecutive segments be collinear, since that would merge two sides of the polygon.

For example, the following images are from Figures 2 and 10 in the Chow, Gafni, Gafni paper, which show a $$\49\$$-gon on the $$\7 \times 7\$$ grid of dots and a $$\100\$$-gon on the $$\10 \times 10\$$ grid of dots:  In the paper, the authors show that it is possible to form these $$\n^2\$$-gons for all integers $$\n \geq 2\$$ except for $$\n=3\$$ and $$\n=5\$$.

# Challenge

This challenge will have you count the number of "essentially different" $$\n^2\$$-gons on the $$\n \times n\$$ grid, where two polygons are considered the same if you can reflect or rotate one to match the other.

To score this, I will run the code on my machine—a 2017 MacBook with an Intel chip and 8GB of RAM—for five minutes with increasing sizes of $$\n\$$. Your score will be the maximum value of $$\n\$$ that your program can produce, and ties will be broken by the amount of time it took to get to that value.

# Small test cases

(From the Chow, Gafni, Gafni paper.)

• When $$\n = 2\$$ there is only one $$\4\$$-gon, the square.
• When $$\n = 3\$$ there are no $$\9\$$-gons.
• When $$\n = 4\$$ there is only one $$\16\$$-gon: • When $$\n = 5\$$ there are no $$\25\$$-gons.
• What do the dashed lines mean? Apr 10 at 20:52
• @Anush, the dashed lines are—I think—a way of showing the recursive structure of this construction. They indicate the boundary between one "layer" of squares and the next, but they're not relevant for this challenge. Apr 10 at 20:55
• Are there any test cases? (also, the 'if you can reflect or rotate one to match the other' will make this very annoying to solve) Apr 12 at 13:34
• I don't expect that the "reflect or rotate" criterion should be too annoying to implement, actually. One brute-force way to do it is to try out all eight symmetries of the square and take the one that is lexicographically earliest in some sense. Apr 12 at 16:37
• @thedefault. I added some small test cases. Apr 12 at 16:47

# Python 3 + Z3, N=5

Not 100% sure if this is correct. Basically a standard SAT/SMT coding of the Hamiltonian cycle problem, with additional clauses to block intersecting diagonals and collinear segments. Each time a solution is found all of its distinct rotations/flips are blocked. I've run it up to N=6 (2120 solutions found), which took 30 minutes.

Edit: I just realized that this only allows edges between adjacent grid points, but this restriction wasn't intended by the challenge. Oh well, I'll leave it up for now.

import collections
import heapq
import itertools
import math
import random
import time
import sys
import z3

show_paths = len(sys.argv) >= 2 and 'show' in sys.argv

if show_paths:
import numpy as np
import pylab as plt
from matplotlib import collections  as mc

def show_path(path):
lines = list(zip(path, path[1:] + path[:1]))
c = np.array([(1, 0, 0, 1), (0, 1, 0, 1), (0, 0, 1, 1)])
lc = mc.LineCollection(lines, colors=c, linewidths=2)
fig, ax = plt.subplots()
ax.autoscale()
ax.margins(0.1)
plt.show()

def compute_min_distances(C, distance_matrix, starting_points, banned):
min_distances = collections.defaultdict(lambda: float('inf'))
open_set = []
for i in starting_points:
state = (i, C[i])
min_distances[state] = 1
open_set.append((1, state))
while open_set:
d, state = open_set.pop()
if d > min_distances[state]:
continue
i, delta = state
x0, y0 = c0 = C[i]
for j, _ in distance_matrix[i]:
if j in banned:
continue
x1, y1 = c1 = C[j]
delta1 = (x1 - x0, y1 - y0)
if delta in (delta1, (x0 - x1, y0 - y1)):
continue
state = (j, delta1)
if d + 1 < min_distances[state]:
min_distances[state] = d + 1
heapq.heappush(open_set, (d + 1, state))
distance_array = [float('inf') for _ in distance_matrix]
distance_array = 0
for (i, _), d in min_distances.items():
distance_array[i] = min(distance_array[i], d)
return distance_array

def point_orientations(p, n):
for (x, y) in (p, p[::-1]):
for x1 in (x, n - 1 - x):
for y1 in (y, n - 1 - y):
yield (x1, y1)

def path_orientations(path, n):
orientations = list(zip(*map(lambda p: point_orientations(p, n), path)))
for o in orientations:
i = o.index((0, 0))
if o[i - 1] == (0, 1):
o = o[::-1]
i = (len(o) - 1) - i
if o[(i + 1) % len(o)] != (0, 1):
continue
yield o[i:] + o[:i]

start_time = time.time()
for N in itertools.count(2):
print("starting n = %d..." % N)
solver = z3.Solver()
C = [(x, y) for x in range(N) for y in range(N)]
C_inv = {c : i for i, c in enumerate(C)}
bit_width = math.ceil(math.log(N * N + 1, 2))
P = [z3.BitVec('P_%d' % i, bit_width) for i in range(N * N)]
V = [z3.BitVec('V_%d' % i, bit_width) for i in range(N * N)]
for v in V:
A = [[] for _ in range(N * N)]
for i, (x, y) in enumerate(C):
for dx, dy in itertools.product(range(-1, 2), repeat=2):
if not (dx or dy):
continue
try:
A[i].append(C_inv[(x + dx, y + dy)])
except KeyError:
pass
H = [[(j, z3.Bool('H_%d_%d' % (i, j))) for j in A[i]] for i in range(N * N)]
start = C_inv[(0, 1)]
min_distances = compute_min_distances(C, H, [start], {0})
end0, end1 = C_inv[(1, 0)], C_inv[(1, 1)]
min_distances_rev = compute_min_distances(C, H, [end0, end1], {0, start})
for i, (d0, d1) in enumerate(zip(min_distances, min_distances_rev)):
if i not in (0, start):
solver.add(z3.ULE(P[i], N * N - d1))
H_t = [[] for _ in range(N * N)]
for i, row in enumerate(H):
x0, y0 = C[i]
for j, v in row:
H_t[j].append(v)
if i == 0:
e = P[j] == 1
elif j == 0:
e = P[i] == (N * N - 1)
else:
e = P[j] == (P[i] + 1)
clauses = [V[i] == j, e]
x1, y1 = C[j]
dy = y1 - y0
dx = x1 - x0
blocked = []
if dy and dx:
x_m = min(x0, x1)
y_m = min(y0, y1)
if dy == dx:
blocked.append(((x_m + 1, y_m), (x_m, y_m + 1)))
else:
blocked.append(((x_m, y_m), (x_m + 1, y_m + 1)))
blocked.append(blocked[-1][::-1])
x2 = x1 + dx
y2 = y1 + dy
if all(0 <= c < N for c in (x2, y2)):
blocked.append(((x1, y1), (x2, y2)))
for (c0, c1) in blocked:
k = C_inv[c0]
l = C_inv[c1]
clauses.append(z3.Not(next(w for m, w in H[k] if m == l)))
for row in H_t:
solver.add(z3.PbEq([(v, 1) for v in row], 1))
solutions = set()
old_len = 0
while solver.check() == z3.sat:
model = solver.model()
successors = [int(str(model.eval(v))) for v in V]
path = 
while len(path) < N * N:
path.append(successors[path[-1]])
path = tuple(C[i] for i in path)

if show_paths:
show_path(path)

orientations = set(path_orientations(path, N))
assert not any(reoriented in solutions for reoriented in orientations)
for reoriented in orientations:
successor = {}
for p0, p1 in zip(reoriented, reoriented[1:] + reoriented[:1]):
successor[C_inv[p0]] = C_inv[p1]
solver.add(z3.Not(z3.And(*[v == successor[i] for i, v in enumerate(V)])))
s = "%d solutions found so far... (%.1f seconds elapsed)" % (len(solutions), time.time() - start_time)
sys.stdout.write(s + max(0, old_len - len(s)) * ' ' + '\r')
sys.stdout.flush()
old_len = len(s)
sys.stdout.write(' ' * old_len + '\r')
sys.stdout.flush()
end_time = time.time()
print("%d solution%s found for n = %d in %.2g seconds" % (len(solutions), 's' * (len(solutions) != 1), N, end_time - start_time))
start_time = end_time

• What is the output for N=6? Apr 12 at 19:30

Work in progress... Still far over the 5 minute deadline set by OP, but already sharing some code to show what I have done so far.

Output:

N=6
Found 2956807 solutions in 12953858177 iterations.
Running time 41 minutes 17 seconds.


(Note that the previous version reported 9358354 without filtering rotation/mirroring. It appears that version had too aggressive pruning of the search tree. Unfortunately, less pruning in this version increased the running time a lot)

The strategy is to simply brute-force all possible paths as fast as possible and apply some pruning rules to cut the search tree short. The most effective pruning rule is to enforce that edge points are visited in clock-wise order.

For efficiency, I use a lookup table for all possible lines starting from a each grid point. And another lookup table to check what other lines are blocked/crossed by a certain line. I use 1D arrays instead of 2D containers/arrays.

My first experiments were in Java, which I later ported to C++ in the hope that that would be faster. Guess what, this only saved a few percent. Modern Java is very fast and comparable to C++ (for this use-case at least).

Contrary the code of @user1502040, my code not only considers lines between adjacent points, but between any two points (as long as there are no other intermediate points on that line). That gives some nice arbitrary paths. For example: As suggested by @PeterKagey, for each solution found I generate all 8 permutations (mirroring/rotating) and use lexicographical sorting to see if we found the "lowest" solution of these eight. If not, we ignore the solution and trust we will find the lowest solution later during the search.

Below is the unpolished, work-in-progress code. While the program is running, it will update the output/progress.html file every 5 seconds to show a diagram with the current search path. If you open this file in a browser, the browser will automatically keep refreshing the page.

In addition, it will save an output/solution-xyz.html file every few solutions. This enables you to sample a few example solutions.

To be continued...

// To compile:
//    \$ g++ -Wall --std=c++11 -DNDEBUG -O3 ConnectingTheDots.cpp \
//          -o ConnectingTheDots
//
// Examples:
//    N=2
//    Found 1 solutions in 4 iterations.
//    Running time 0 minutes 0 seconds.
//
//    N=3
//    Found 0 solutions in 9 iterations.
//    Running time 0 minutes 0 seconds.
//
//    N=4
//    Found 1 solutions in 422 iterations.
//    Running time 0 minutes 0 seconds.
//
//    N=5
//    Found 0 solutions in 376217 iterations.
//    Running time 0 minutes 0 seconds.
//
//    N=6
//    Found 2956807 solutions in 12953858177 iterations.
//    Running time 41 minutes 17 seconds.
//
// Improvement opportunities:
//
// - Parallelize: a brute-force search is trivial to run massively parallel on
//   all available CPU cores. Use POSIX pthread to create worker threads.
//
// - Data locality: lookup tables currently use 32-bit integers and have large
//   ROW_SIZEs. This reduces data locality and thus memory caching. Ideas:
//       - Consider using 8-bit resp. 16-bit integers (for N <= 8)
//       - Store data consecutively in array and use row pointers.
//
// - Roll back: making a step and rolling it back involves selectively saving
//   and restoring parts of the state. The overhead of this bookkeeping is
//   large. Perhaps it is faster to blindly copy the entire state.
//
//   To make the copy fast, we could condense the state into a smaller memory
//   footprint (at the cost of some bit masking and shifting). For example, the
//   pointAvailable[] and lineAvailable[] boolean arrays could be bitmaps which
//   makes them 32x times smaller.
//
// - Island pruning rule: when the path waves from an edge point inward and
//   then back to the next edge point, it may not enclose an unvisited point.
//   These islands can never be visited. Beware that the island must be really
//   enclosed. Just checking neighbouring points is not enough, since sharp
//   "spikes" can still reach far away isolated points.

#include <stdlib.h>
#include <stdio.h>
#include <algorithm> // std::min
#include <time.h>
#include <sys/stat.h>   // mkdir

// User-specified dimensions
const int N = 6;

// Calculated dimensions
const int MAX_POS = N*N;
const int MAX_LINE = MAX_POS*MAX_POS;

/** ---------------------------------------------------------------------------
* Class to check whether two line segments intersect.
* 2D integer coordinates only.
*
* Reference:
*  https://www.geeksforgeeks.org/check-if-two-given-line-segments-intersect/
*  http://www.dcs.gla.ac.uk/~pat/52233/slides/Geometry1x1.pdf
*/
class Intersect {
public:
/**
* Check whether lines p1q1 and p2q2 cross eachother.
*
* @return true if crossed
*/
static bool isCrossed(int p1, int q1, int p2, int q2) {
assert(p1 >= 0 && p1 < MAX_POS);
assert(q1 >= 0 && q1 < MAX_POS);
assert(p2 >= 0 && p2 < MAX_POS);
assert(q2 >= 0 && q2 < MAX_POS);
assert(p1 != q1);
assert(p2 != q2);

// endpoints may touch
if (p1 == p2 || p1 == q2 || q1 == p2 || q1 == q2) return false;

// convert point IDs to 2D coordinates
int p1x = p1 % N; int p1y = p1 / N;
int q1x = q1 % N; int q1y = q1 / N;
int p2x = p2 % N; int p2y = p2 / N;
int q2x = q2 % N; int q2y = q2 / N;

// Find the four orientations needed for general and
// special cases
int o1 = orientation(p1x, p1y, q1x, q1y, p2x, p2y);
int o2 = orientation(p1x, p1y, q1x, q1y, q2x, q2y);
int o3 = orientation(p2x, p2y, q2x, q2y, p1x, p1y);
int o4 = orientation(p2x, p2y, q2x, q2y, q1x, q1y);

// General case
if (o1 != o2 && o3 != o4) return true;

// Special Cases
// p1, q1 and p2 are colinear and p2 lies on segment p1q1
if (o1 == 0 && onSegment(p1x, p1y, p2x, p2y, q1x, q1y)) return true;

// p1, q1 and q2 are colinear and q2 lies on segment p1q1
if (o2 == 0 && onSegment(p1x, p1y, q2x, q2y, q1x, q1y)) return true;

// p2, q2 and p1 are colinear and p1 lies on segment p2q2
if (o3 == 0 && onSegment(p2x, p2y, p1x, p1y, q2x, q2y)) return true;

// p2, q2 and q1 are colinear and q1 lies on segment p2q2
if (o4 == 0 && onSegment(p2x, p2y, q1x, q1y, q2x, q2y)) return true;

return false; // Doesn't fall in any of the above cases
}

private:
/**
* To find orientation of ordered triplet (p, q, r).
*
* @return 0 --> p, q and r are colinear
*       1 --> Clockwise
*       2 --> Counterclockwise
*/
static int orientation(int px, int py, int qx, int qy, int rx, int ry) {
int val = (qy - py)*(rx - qx) - (qx - px)*(ry - qy);

if (val == 0) return 0; // colinear
return (val > 0) ? 1 : 2; // clock or counterclock wise
}

/**
* Given three colinear points p, q, r, the function checks if
* point q lies on line segment 'pr'.
*
* @return true if point q lies on line segment pr
*/
static bool onSegment(int px, int py, int qx, int qy, int rx, int ry) {
if (qx <= std::max(px, rx) && qx >= std::min(px, rx) &&
qy <= std::max(py, ry) && qy >= std::min(py, ry)) {
return true;
}
return false;
}
};

/** ---------------------------------------------------------------------------
* Lookup table with all possible connections.
*/
class LinesLookupTable {
public:
/**
* Build look up table.
*/
void build() {
assert(ROW_SIZE >= (N*N)-1);
lines = new int[MAX_POS*ROW_SIZE];

// iterate point i over grid
for (int i = 0; i < MAX_POS; i++) {
int x = i%N;
int y = i/N;

int nrLines = 0;

// iterate point j over grid
for (int j = 0; j < MAX_POS; j++) {
if (i == j) continue;

int dx = (j%N) - x;
int dy = (j/N) - y;

// not allowed to have another point on the line
bool isValid = true;
for (int div = 2; div < N; div++) {
if ((abs(dx) % div) == 0 &&
(abs(dy) % div) == 0) {
isValid = false;
break;
}
}
if (!isValid) continue;

// found a valid line
assert(nrLines < ROW_SIZE);
lines[i*ROW_SIZE + nrLines++] = j;
}

// write end of list marker
assert(nrLines < ROW_SIZE);
lines[i*ROW_SIZE + nrLines] = -1;
}
}

/**
* Get list of lines connecting specified position.
* @param position   start position
* @return array with connecting positions, terminated by -1
*/
const int *getLines(int position) const {
assert(position >= 0 && position < MAX_POS);
return &lines[position*ROW_SIZE];
}

private:
static const int  ROW_SIZE = 64;
int               *lines;
};

/** ---------------------------------------------------------------------------
* Create line ID based on two points.
* We choose a consistent order of begin and end point to prevent
* double counting.
*
* @param p  begin point
* @param q  end point
*/
inline int makeLineId(int p, int q) {
assert(p >= 0 && p < MAX_POS);
assert(q >= 0 && q < MAX_POS);
return std::min(p, q)*MAX_POS + std::max(p, q);
}

/** ---------------------------------------------------------------------------
* Lookup table with all other lines that are blocked by a certain line.
*/
class BlockedLookupTable {
public:
BlockedLookupTable(LinesLookupTable *linesLookupTable_) :
linesLookupTable(linesLookupTable_) {
}

/**
* Build look up table.
*/
void build() {
blocked = new int[MAX_LINE*ROW_SIZE];

// iterate point i over grid
for (int i = 0; i < MAX_POS; i++) {
// iterate point j over all lines connecting i
const int *jp = linesLookupTable->getLines(i);
for (int j = *jp++; j != -1; j = *jp++) {
assert(j >= 0 && j < MAX_POS);
int lineId = makeLineId(i, j);
int nrBlocked = 0;

// iterate point p over grid
for (int p = 0; p < MAX_POS; p++) {
// iterate point q over all lines connecting p
const int *qp = linesLookupTable->getLines(p);
for (int q = *qp++; q != -1; q = *qp++) {
assert(q >= 0 && q < MAX_POS);
int lineId2 = makeLineId(p, q);
if (lineId == lineId2) continue;

// check if line i-j crosses line p-q
if (Intersect::isCrossed(i, j, p, q)) {
assert(nrBlocked < ROW_SIZE);
blocked[lineId*ROW_SIZE + nrBlocked++] = lineId2;
}
}
}

// block lines colinear with line i-j
int dx = (j%N) - (i%N);
int dy = (j/N) - (i/N);

// extend out from i
int x = (i%N) - dx;
int y = (i/N) - dy;
if (x >= 0 && x < N && y >= 0 && y < N) {
int k = x + y*N;
int lineId2 = makeLineId(i, k);
assert(nrBlocked < ROW_SIZE);
blocked[lineId*ROW_SIZE + nrBlocked++] = lineId2;
}

// extend out from j
x = (j%N) + dx;
y = (j/N) + dy;
if (x >= 0 && x < N && y >= 0 && y < N) {
int k = x + y*N;
int lineId2 = makeLineId(j, k);
assert(nrBlocked < ROW_SIZE);
blocked[lineId*ROW_SIZE + nrBlocked++] = lineId2;
}

// write end of list marker
assert(nrBlocked < ROW_SIZE);
blocked[lineId*ROW_SIZE + nrBlocked] = -1;
}
}
}

/**
* Get list of lines blocked by the specified line
* @param lineId  line
* @return array with blocked lines, terminated by -1
*/
const int *getBlockedLines(int lineId) {
assert(lineId >= 0 && lineId < MAX_LINE);
assert(lineId/MAX_POS < lineId%MAX_POS);
return &blocked[lineId*ROW_SIZE];
}

private:
static const int  ROW_SIZE = 512;
LinesLookupTable   *linesLookupTable;
int               *blocked;
};

/** ---------------------------------------------------------------------------
* Class for printing paths to console and rendering path diagrams in HTML file
*/
class Output {
public:
Output() :
updateProgressTime(0) {
}

/**
* Show current path.
*/
void printPath(const int *path, int pathLen) {
if (pathLen > 0) {
printf("%d", path);
for (int i = 1; i < pathLen; i++) {
printf(", %d", path[i]);
}
}
if (pathLen < MAX_POS+1) {
printf(", ...");
}
printf("\n");
}

/**
* Periodically dump path to console, and redraw diagram in HTML file.
*/
void updateProgress(const int *path, int pathLen) {
if (updateProgressInterval > 0 &&
time(NULL) >= updateProgressTime + updateProgressInterval) {

printPath(path, pathLen);

char filePath;
snprintf(filePath, sizeof(filePath), "%s/%s", outputDir,
progressFile);
drawPath(path, pathLen, filePath, true);

updateProgressTime = time(NULL);
}
}

/**
* Draw solution diagram in HTML file.
*/
void drawSolution(const int *path, int pathLen, int solutionId) {
char fileName;
snprintf(fileName, sizeof(fileName), solutionFile, solutionId);

char filePath;
snprintf(filePath, sizeof(filePath), "%s/%s", outputDir, fileName);

drawPath(path, pathLen, filePath, false);
printf("Written \"%s\"\n", filePath);
}

private:
static const int updateProgressInterval   = 5;  // in seconds (0=disable)
static constexpr const char *outputDir    = "output";
static constexpr const char *progressFile = "progress.html";
static constexpr const char *solutionFile = "solution-%04d.html";

time_t   updateProgressTime;

/**
* Ensure output directory exists.
*/
void makeOutputDir() {
struct stat st = {0};
if (stat(outputDir, &st) == -1) {
if (mkdir(outputDir, 0775)) {
perror("Failed to create output directory");
exit(-1);
}
}
}

/**
* Draw current (partial) path in SVG diagram and embed in HTML file.
* @param filePath   directory and file name of HTML file
* @param refresh    if true the HTML file will periodically reload in the
*                   browser
*/
void drawPath(const int *path, int pathLen, const char *filePath,
bool refresh) {
makeOutputDir();

char tempFilePath;
snprintf(tempFilePath, sizeof(tempFilePath), "%s.tmp", filePath);
FILE *fp = fopen(tempFilePath, "w");
if (fp == NULL) {
perror("Failed to create file");
return;
}
fprintf(fp, "<!DOCTYPE html>\n");
fprintf(fp, "<html>\n");
if (refresh) {
fprintf(fp, "<meta http-equiv=\"refresh\" content=\"%d\" />\n",
updateProgressInterval+1);
}

fprintf(fp, "<body>\n");
fprintf(fp, "<svg xmlns=\"http://www.w3.org/2000/svg\" version=\"1.1\""
" width=\"%.1f\" height=\"%.1f\">\n", N*100.0, N*100.0);
for (int y = 0; y < N; y++) {
for (int x = 0; x < N; x++) {
fprintf(fp, "<circle cx=\"%.1f\" cy=\"%.1f\" r=\"5.0\" "
"fill=\"black\" />\n", 50.0 + x*100.0, 50.0 + y*100.0);
}
}
for (int i = 0; i+1 < pathLen; i++) {
int x1 = path[i]%N;
int y1 = path[i]/N;
int x2 = path[i+1]%N;
int y2 = path[i+1]/N;
fprintf(fp, "<polyline points=\"%.1f,%1.f %.1f,%.1f\" "
"style=\"fill:none;stroke:#000000;stroke-width:2.0;\" "
"/>\n",
50.0 + x1*100.0, 50.0 + y1*100.0,
50.0 + x2*100.0, 50.0 + y2*100.0);
}
fprintf(fp, "</svg>\n");
fprintf(fp, "</body>\n");
fprintf(fp, "</html>\n");
fclose(fp);
rename(tempFilePath, filePath);
}
};

/** ---------------------------------------------------------------------------
* Search all possible paths.
*/
class Solver {
public:
Solver(
LinesLookupTable *linesLookupTable_,
BlockedLookupTable *blockedLookupTable_,
Output *output_):

linesLookupTable(linesLookupTable_),
blockedLookupTable(blockedLookupTable_),
output(output_) {
}

/**
* Search all possible paths
*/
void solve() {
buildNextEdgePointTable();

pointAvailable = new bool[MAX_POS];
for (int i = 0; i < MAX_POS; i++) {
pointAvailable[i] = (nextEdgePoint[i] == -1);
}
pointAvailable[N] = false; // end of edge
pointAvailable[nextEdgePoint] = true;

lineAvailable = new bool[MAX_LINE];
for (int i = 0; i < MAX_LINE; i++) {
lineAvailable[i] = true;
}

path = new int[MAX_POS+1];
pathLen = 0;

nrIterations = 0;
nrSolutions = 0;
path[pathLen++] = 0;
pointAvailable = false;
recurseNextStep(0);

printf("Found %llu solutions in %llu iterations.\n",
nrSolutions, nrIterations);
}

private:
LinesLookupTable    *linesLookupTable;
BlockedLookupTable  *blockedLookupTable;
Output              *output;

int                 *nextEdgePoint;
bool                *pointAvailable;
bool                *lineAvailable;
unsigned long long   nrSolutions;
int                 *path;
int                  pathLen;
unsigned long long   nrIterations;

/**
* Build nextEdgePoint[] lookup table
*/
void buildNextEdgePointTable() {
nextEdgePoint = new int[MAX_POS];
for (int i = 0; i < MAX_POS; i++) {
nextEdgePoint[i] = -1;
}
int prev = 0;
for (int x = 1; x < N; x++) {
int y = 0;
int point = y*N + x;
nextEdgePoint[prev] = point;
prev = point;
}
for (int y = 1; y < N; y++) {
int x = N-1;
int point = y*N + x;
nextEdgePoint[prev] = point;
prev = point;
}
for (int x = N-2; x >= 0; x--) {
int y = N-1;
int point = y*N + x;
nextEdgePoint[prev] = point;
prev = point;
}
for (int y = N-2; y >= 1; y--) {
int x = 0;
int point = y*N + x;
nextEdgePoint[prev] = point;
prev = point;
}
}

/**
* Iterate over all possible next steps starting at current position.
* @param i  current 1D position
*/
void recurseNextStep(int i) {
nrIterations++;
if (pathLen == MAX_POS) {
// We are ready to make the last hop. Check if we can close the loop.
closeLoop(i);
} else {
int blockedLines;
int nrBlockedLines = 0;

// Iterate over all possible next steps 'j' starting at
// current position 'i'.
const int *jp = linesLookupTable->getLines(i);
for (int j = *jp++; j >= 0; j = *jp++) {
assert(j >= 0 && j < MAX_POS);

// target point still available?
if (pointAvailable[j]) {

// line to target point still available?
int lineId = makeLineId(i, j);
assert(lineId >= 0 && lineId < MAX_LINE);
if (lineAvailable[lineId]) {

// found a valid next step
path[pathLen++] = j;

// mark points and lines as used, but keep track of
// changes so we can roll them back later
lineAvailable[lineId] = false;
const int *bp = blockedLookupTable->getBlockedLines(lineId);
for (int lineId2 = *bp++; lineId2 != -1; lineId2 = *bp++) {
assert(lineId2 >= 0 && lineId2 < MAX_LINE);
if (lineAvailable[lineId2]) {
assert(nrBlockedLines < 256);
blockedLines[nrBlockedLines++] = lineId2;
lineAvailable[lineId2] = false;
}
}
pointAvailable[j] = false;

// make sure edge points are visited in order
if (nextEdgePoint[j] != -1) {
assert(!pointAvailable[nextEdgePoint[j]]);
pointAvailable[nextEdgePoint[j]] = true;
}

// search next step
recurseNextStep(j);

// roll back changes
if (nextEdgePoint[j] != -1) {
pointAvailable[nextEdgePoint[j]] = false;
}
pointAvailable[j] = true;
pathLen--;
while (nrBlockedLines > 0) {
nrBlockedLines--;
int lineId2 = blockedLines[nrBlockedLines];
assert(lineId2 >= 0 && lineId2 < MAX_LINE);
lineAvailable[lineId2] = true;
}
assert(lineId >= 0 && lineId < MAX_LINE);
lineAvailable[lineId] = true;
}
}
}

output->updateProgress(path, pathLen);
}
}
}

/**
* Try to close the loop back to position 0.
* @param position   current 1D position
*/
void closeLoop(int i) {
// Iterate over all possible next steps 'j' starting at
// current position 'i'.
const int *jp = linesLookupTable->getLines(i);
for (int j = *jp++; j >= 0; j = *jp++) {
assert(j >= 0 && j < MAX_POS);
if (j != 0) continue;

// line to target point still available?
int lineId = makeLineId(i, j);
if (lineAvailable[lineId]) {
path[pathLen++] = j;

// found a valid loop! Now check if this is the "lowest" of all
// equivalent paths.
if (isLowestSolution()) {
// Sample some random paths by drawing every few solutions.
// for N=6, we will find 2956807 solutions. Aim for 100 examples.
if ((nrSolutions % 29568) == 0) {
output->drawSolution(path, pathLen, nrSolutions);
}

nrSolutions++;
}

pathLen--;
}
}
}

/**
* Check whether the solution found is the "lowest" equivalent solution
* when sorted lexicographically.
*
* By rotating and mirroring each path can result in up to 8 equivalent
* paths. We choose the lexicographically lowest solution as the "main"
* solution.
*
* @returns true if the current path equals lowest equivalent path
*/
bool isLowestSolution() {
assert(pathLen >= MAX_POS);

int rotatedPath[MAX_POS+1];
memcpy(rotatedPath, path, pathLen*sizeof(rotatedPath));

for (int permutation = 1; permutation < 8; permutation++) {
// rotate clockwise
for (int i = 0; i < pathLen; i++) {
int x = rotatedPath[i]%N;
int y = rotatedPath[i]/N;

int rotX = N-1-y;
int rotY = x;

rotatedPath[i] = rotY*N + rotX;
}

// after 4 rotations we are upright again. Now mirror horizontally.
// note after mirroring the path runs counter-clockwise. We therefore
// also need to reverse the order of the points visited
if (permutation == 4) {
#ifndef NDEBUG
// assert we are upright again
for (int i = 0; i < pathLen; i++) {
assert(path[i] == rotatedPath[i]);
}
#endif
int mirroredPath[MAX_POS+1];
for (int i = 0; i < pathLen; i++) {
int x = rotatedPath[i]%N;
int y = rotatedPath[i]/N;

int rotX = N-1-x;
int rotY = y;

mirroredPath[pathLen-1-i] = rotY*N + rotX;
}
memcpy(rotatedPath, mirroredPath, pathLen*sizeof(rotatedPath));
}

if (comparePaths(path, rotatedPath) > 0) return false;
}

// of all 8 permutations, the current path is the "lowest" (or equal
// to some of the permutations in case of symmetric, but since the path
// generator will only generate that path once, it is not a problem).
return true;
}

/**
* Compare two paths.
* Assumes path A starts at point 0, path B can start anywhere.
*
* @return -1 if A < B, +1 if A > B
*/
int comparePaths(const int *pathA, const int *pathB) {
assert(pathA == 0);

// find point 0 in path B
int offset;
for (offset = 0; offset < MAX_POS && pathB[offset] != 0; offset++);
assert(pathB[offset] == 0);

// compare path A with path B
for (int i = 0, j = offset; i < MAX_POS; i++, j = (j+1) % MAX_POS) {
if (pathA[i] < pathB[j]) return -1;
if (pathA[i] > pathB[j]) return +1;
}
return 0;
}
};

/** ---------------------------------------------------------------------------
* Main program
*/
int main() {
LinesLookupTable    *linesLookupTable;
BlockedLookupTable  *blockedLookupTable;
Output              *output;
Solver              *solver;

time_t startTime = time(NULL);

linesLookupTable   = new LinesLookupTable();
blockedLookupTable = new BlockedLookupTable(linesLookupTable);
output             = new Output();
solver             = new Solver(linesLookupTable, blockedLookupTable,
output);
printf("N=%d\n", N);

linesLookupTable->build();
blockedLookupTable->build();
solver->solve();

int seconds = time(NULL) - startTime;
printf("Running time %d minutes %d seconds.\n", seconds/60, seconds%60);

return 0;
}


• I think C++ can be slow if you are not used to optimizing it. Could you try C? Or maybe someone else can once you post your code. May 20 at 7:28
• @Anush: I don't expect a lot of difference between C and C++ since I stayed away from containers and other STL abstractions. The code is basically "C-with-classes". I added some improvement suggestions at the top of the file. May 20 at 19:48