Connecting the Dots: Counting n²-gons in the n×n Grid
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