Schinzel minimal

schinzel.py

# https://enigmaticcode.wordpress.com/2013/10/15/enigma-136-twelve-point-square/

from itertools import count, product, combinations
import gmpy2
 
F = gmpy2.mpq
 
# 1/2
h = F(1, 2)
 
# centre of the square and distance from it
(xt, yt, rt2) = (h, h, h)
 

# the smallest circle we haven't found yet
nmin = 3
# maximum radius to check
rmax = 20
# numbers which have already been printed on the console
printed = []
 
# canonical position of centres
def canonical(x, y):
  (x, y) = (x % 1, y % 1)
  if y > h: y = 1 - y
  if x > h: x = 1 - x
  if y > x: (x, y) = (y, x)
  return (x, y)
 
# record the smallest circle
r = dict()

 
for z in count(1):
  # min/max distance (squared)
  if z > rmax: break

  (m2, M2) = ((z - 1) ** 2, z ** 2)
 
  # collect lattice points
  lattice = list()
  for (x, y) in product(range(-z, z+1), repeat=2):
    # calculate distances
    (x0, y0, x1, y1) = (x ** 2, y ** 2, (x - 1) ** 2, (y - 1) ** 2)
    if any(m2 <= d2 < M2 for d2 in (x0 + y0, x0 + y1, x1 + y0, x1 + y1)):
      lattice.append((x, y))
 
  # we only need to check points in one sector of the circle
  while nmin in r: nmin += 1
  if nmin > 16:
    (k, check) = (8, list((x, y) for (x, y) in lattice if not(y < 0 or x < 0 or x + 1 < y)))
  elif nmin > 8:
    (k, check) = (4, list((x, y) for (x, y) in lattice if not(y < 0 or x < 0)))
  elif nmin > 4:
    (k, check) = (2, list((x, y) for (x, y) in lattice if not(y < 0)))
  else:
    (k, check) = (1, lattice)
 

  print("[checking radius <", z, "] [nmin=", nmin,"] [", len(check), " points in 1/", k, " sector] ...")
  for radius in range(rmax):
      if radius in r and radius not in printed:
          u = r[radius]
          print(radius, ": ", u[0]) 
          printed.append(radius)
  # choose two lattice points
  p = len(check)
  for (i, j) in combinations(range(0, p - 2), 2):
    ((x1, y1), (x2, y2)) = (check[i], check[j])
    # the parameters for the perpendicular bisector of p1 and p2 (Ax + By = C)
    a1 = 2 * (x2 - x1)
    b1 = 2 * (y2 - y1)
    c1 = (x1 ** 2 - x2 ** 2) + (y1 ** 2 - y2 ** 2)
    # does it pass through the origin square
    d2 = F(((a1 * xt) + (b1 * yt) + c1) ** 2, a1 ** 2 + b1 ** 2)
    if d2 > rt2: continue
 
    # choose the remaning lattice point
    for k in range(j + 1, p - 1):
      (x3, y3) = check[k]
      # the parameters for the perpendicular bisector of p1 and p3 (Ax + By = C)
      a2 = 2 * (x3 - x1)
      b2 = 2 * (y3 - y1)
      # check they are not parallel
      d = a1 * b2 - a2 * b1
      if d == 0: continue
      c2 = (x1 ** 2 - x3 ** 2) + (y1 ** 2 - y3 ** 2)      
      # do they intersect in the square
      y0 = F(a2 * c1 - a1 * c2, d)
      if not(0 <= y0 <= 1): continue
      x0 = F(b1 * c2 - b2 * c1, d)
      if not(0 <= x0 <= 1): continue
      # work out the radius of the circle
      r2 = (x0 - x1) ** 2 + (y0 - y1) ** 2
      if not(m2 <= r2 < M2): continue
 
      # count lattice points with the same distance
      n = sum(1 for (x, y) in lattice if (x0 - x) ** 2 + (y0 - y) ** 2 == r2)
      if n not in r or r2 < r[n][0]:
        # record (distance squared, canonical centre)
        (x0, y0) = canonical(x0, y0)
        r[n] = (r2, x0, y0)
        #print("-> n=", n, " r=", sqrt(r2), " (r^2=", r2, ") (", x0, ", ", y0, ")")

            

for radius in range(rmax):
    if radius in r:
        u = r[radius]
        print(radius, ": ", u[0])