Given a square 5x5 grid of points, find five circles that together pass through all 25 points.

fivecircles.py

#!/usr/bin/env/python
#   Author: David Radcliffe (dave@gotmath.com)
#   Date: 20 December 2013
#   This code is in the public domain.

# Problem: Given a square 5x5 grid of points, find five circles that together pass through all 25 points.
# There are 21 essentially different solutions. Two solutions are considered the same if they differ only
# by a rotation or reflection. This script finds all of these solutions.
# The solutions are rendered graphically at http://www.gotmath.com/doc/fivecircles.html


from math import sqrt
from time import time
start_time = time()

N = 6
num_circles = 6
filename = "circle_covers_%d_%d.csv" % (N, num_circles)
N2 = N*N
eps = 1e-6
circles = set()
arrangements = [set() for k in range(num_circles + 1)]


# Are the points (x1, y1), (x2, y2), and (x3, y3) collinear?
def collinear(x1, y1, x2, y2, x3, y3):
    return (x2 - x1) * (y3 - y1) == (x3 - x1) * (y2 - y1)


# Euclidean distance between (x1, y1) and (x2, y2)
def dist(x1, y1, x2, y2):
    return sqrt((x1 - x2)**2 + (y1 - y2)**2)


# Use Ptolemy's theorem to determine whether the points (x1, y1), (x2, y2), (x3, y3), and (x4, y4)
# are concyclic (lie on a common circle).
def concyclic(x1, y1, x2, y2, x3, y3, x4, y4):
    a = dist(x1, y1, x2, y2) * dist(x3, y3, x4, y4)
    b = dist(x1, y1, x3, y3) * dist(x2, y2, x4, y4)
    c = dist(x1, y1, x4, y4) * dist(x2, y2, x3, y3)
    return abs(a + b + c - 2*max(a,b,c)) < eps


# Get the (x,y) coordinates of the i-th grid point.
def coordinates(i):
    return i % N, i / N


# Find all circles that pass through 3 or more points of the grid.
# Each circle is represented as a string of 0s and 1s. The 1s indicate the grid points that lie on the circle.
def get_circles():
    grid = ['0'] * N2
    for i1 in range(N2-2):
        x1, y1 = coordinates(i1)
        grid[i1] = '1'
        for i2 in range(i1+1, N2-1):
            x2, y2 = coordinates(i2)
            grid[i2] = '1'
            for i3 in range(i2+1, N2):
                x3, y3 = coordinates(i3)
                grid[i3] = '1'
                if not collinear(x1, y1, x2, y2, x3, y3):
                    for i4 in range(N2):
                        x4, y4 = coordinates(i4)
                        if concyclic(x1, y1, x2, y2, x3, y3, x4, y4):
                            grid[i4] = '1'
                    g = ''.join(grid)
                    if not(g in circles):
                        circles.add(g)
                grid = ['0']*N2
                grid[i1] = '1'
                grid[i2] = '1'
            grid[i2] = '0'
        grid[i1] = '0'
VerticalFlip = {}
Transpose = {}


def vertical_flip(g):
    if g in VerticalFlip:
        return VerticalFlip[g]
    flipped = ''.join([g[k*N:(k+1)*N] for k in xrange(N-1,-1,-1)])
    VerticalFlip[g] = flipped
    return flipped


def transpose(g):
    if g in Transpose:
        return Transpose[g]
    transposed = ''.join(g[y*N+x] for x in range(N) for y in range(N))
    Transpose[g] = transposed
    return transposed


def symmetrize(arrangement):
    arr1 = sorted(arrangement)
    arr2 = sorted(map(vertical_flip, arr1))
    arr3 = sorted(map(transpose, arr2))
    arr4 = sorted(map(vertical_flip, arr3))
    arr5 = sorted(map(transpose, arr4))
    arr6 = sorted(map(vertical_flip, arr5))
    arr7 = sorted(map(transpose, arr6))
    arr8 = sorted(map(vertical_flip, arr7))
    return tuple(min(arr1, arr2, arr3, arr4, arr5, arr6, arr7, arr8))


def weight(arrangement):
    n = len(arrangement)
    return sum(any(arrangement[i][j] == '1' for i in xrange(n)) for j in xrange(N2))


def get_arrangements(k):
    if k == 0:
        arrangements[0].add(())
    else:
        min_weight = 1 + (k*N2-1)/num_circles
        for arrangement in arrangements[k-1]:
            for circle in circles:
                new_arrangement = arrangement + (circle,)
                if weight(new_arrangement) >= min_weight:
                    new_arrangement = symmetrize(new_arrangement)
                    arrangements[k].add(new_arrangement)


# Find the center and radius of a circle given three points on the circle.
# Formula from https://en.wikipedia.org/wiki/Circumscribed_circle#Cartesian_coordinates
def center_and_radius(circle):
    i = circle.find("1")
    ax, ay = coordinates(i)
    i = circle.find("1", i+1)
    bx, by = coordinates(i)
    i = circle.find("1", i+1)
    cx, cy = coordinates(i)
    a = ax * ax + ay * ay
    b = bx * bx + by * by
    c = cx * cx + cy * cy
    d = 2.0 * (ax * (by - cy) + bx * (cy - ay) + cx * (ay - by))
    x = (a * (by - cy) + b * (cy - ay) + c * (ay - by)) / d
    y = (a * (cx - bx) + b * (ax - cx) + c * (bx - ax)) / d
    r = dist(x, y, ax, ay)
    return x, y, r

get_circles()

# for k in range(6):
#     get_arrangements(k)
#     print k, len(arrangements[k])
# print "Time: %d seconds" % (time() - start_time)

for k in range(num_circles+1):
    get_arrangements(k)
    print "%d arrangements with %d circles" % (len(arrangements[k]), k)

outfile = open(filename, "wt")
outlines = [','.join("x%d,y%d,r%d" % (k, k, k) for k in xrange(1,1+num_circles))]

for arr in arrangements[num_circles]:
    line = ','.join(list(str(x) for circle in arr
                         for x in center_and_radius(circle)))
    outlines.append(line)

outfile.writelines("\n".join(outlines))
print "\n".join(outlines)
print 'Time: %s seconds' % (time() - start_time)