Connecting the Dots: Counting n²-gons in the n×n Grid The recent volume of MAA's Mathematics Magazine had an article &quot;Connecting the Dots: Maximal Polygons on a Square Grid&quot; 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 fastest-code challenge will have you count the number of &quot;essentially different&quot; \$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.
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) &gt;= 2 and 'show' in sys.argv[1]

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.add_collection(lc)
        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 &gt; 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 &lt; 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] = 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(&quot;starting n = %d...&quot; % 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:
        solver.add(z3.ULT(v, N * N))
    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.UGE(P[i], d0))
            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 &lt;= c &lt; 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)))
            solver.add(z3.Implies(v, z3.And(*clauses)))
    solver.add(P[0] == 0)
    solver.add(P[start] == 1)
    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 = [0]
        while len(path) &lt; 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)
        solutions.add(path)
        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 = &quot;%d solutions found so far... (%.1f seconds elapsed)&quot; % (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(&quot;%d solution%s found for n = %d in %.2g seconds&quot; % (len(solutions), 's' * (len(solutions) != 1), N, end_time - start_time))
    start_time = end_time