My Math SE question determining if a coincident point in a pair of rotated hexagonal lattices is closest to the origin? explains the problem I have. I won't reproduce the whole thing in detail here, but this image illustrates the problem graphically.
- On the right: the coincident point defined the closed dot point (5, 6) and the open circle point (5, 3) (see red arrow) is (one of the) six closest points to the origin in this coincidence lattice.
- On the left: the coincident similarly defined by (6, 5) and (5, 3) (red arrow) is far from being the closest coincident point to the origin.
The only difference is that we've swapped (6, 5) for (5, 6) by rotating the closed dot lattice by about six degrees.
The answer to my question proposes applying a modification of the well known Euclid algorithm to the problem.
Rather than two real integers, it's now applied to two complex Eisenstein integers which is a hexagonal distribution of points in the complex plane. For more see:
In order to try to implement this very very nicely and clearly described step-by-step method from that answer into Python, I've defined a Point class and given the class the following methods: .subtract(P), .multiply(P), .divround(P), .prod(P)and.sub(P)` reproducing to the best of my abilities what's described in that answer.
Note that a Point object P is defined by two integers P.p and P.q and all of the methods operate on these and return fresh instances of Point. Mostly the operations are addition and multiplication, the only one going off into floats temporarily is .divround(P).
Below is a test program with over 300 cases for (p, q), (k, l) to be tested by the function euclid_2021(p, q, k, l).
So far less-than-rigorous testing suggests that the algorithm is giving what I think is the right answer, at least it seems to agree very nicely with a very long, ugly, crazy looking script that I'd put together earlier, but I don't have a known-correct method to check either against.
But here I am only asking about an apparently endless loop that the method gets stuck in for only one of the 300+ cases below.
Group_1 and Group_3 run quickly and return booleans representing if it's closest or not.
For Group_2 as long as I omit the 2nd case (9, 2, 2, 6) everything is okay. But for that case the loop keeps running indefinitely.
Any thoughts mathematically on why this one case fails?
When testing only the "bad" case (9, 2, 2, 6) with the following print line inserted:
def euclid_2021_with_print(p, q, k, l):
A = Point(p, q)
B = Point(k, l)
if A.abs < B.abs:
A, B = B, A
while B.nonzero:
D = A.divround(B)
B, A = A.subtract(D.multiply(B)), B
print((D.p, D.q), (A.p, A.q), (B.p, B.q)) ### Print inside while loop
return A.isone
the algorithm simply cycles back and forth between two states:
(0, 0) (4, -3) (-2, 3)
(0, 0) (-2, 3) (4, -3)
(0, 0) (4, -3) (-2, 3)
(0, 0) (-2, 3) (4, -3)
(0, 0) (4, -3) (-2, 3)
(0, 0) (-2, 3) (4, -3)
True/False output for Groups 1, 2, 3 but skipping the "bad" case:
import numpy as np
import matplotlib.pyplot as plt
class Point():
ones = ((1, 0), (-1, 0), (0, 1), (0, -1), (-1, 1), (1, -1))
u = 0.5 * (1 + np.sqrt(3) * 1j)
def __init__(self, p, q):
self.p = int(p)
self.q = int(q)
self.pq = (self.p, self.q)
self.nonzero = not ((self.p == 0) and (self.q == 0))
self.isone = self.pq in self.ones
self.abs = np.sqrt(self.p**2 + self.q**2)
self.xy = self.p + self.u * self.q
self.x, self.y = np.real(self.xy), np.imag(self.xy)
def subtract(self, P):
return Point(self.p - P.p, self.q - P.q)
def multiply(self, P):
return Point(self.p * P.p - self.q * P.q,
self.p * P.q + self.q * P.p + self.q * P.q)
def divround(self, P):
(x, y), (z, t) = self.pq, P.pq
ptop, qtop = x * (z + t) + y*t, (y*z - x*t)
bot = z**2 + z*t + t**2
pfrac, qfrac = ptop/bot, qtop/bot
pint, qint = [int(thing) for thing in (pfrac, qfrac)]
four = ((pint, qint), (pint, qint+1), (pint+1, qint), (pint+1, qint+1))
ndsqs = []
for (p, q) in four:
dp, dq = pfrac - p, qfrac - q
ndsqs.append( -1 * (dp**2 + dp * dq + dq**2))
return Point(*four[np.argmax(ndsqs)])
def prod(self, B):
a, b, c, d = self.x, self.y, B.x, B.y
return Point(a*c - b*d, a*d + b*c + b*d)
def sub(self, B):
a, b, c, d = self.x, self.y, B.x, B.y
return Point(a-c, b-d)
def euclid_2021(p, q, k, l):
A = Point(p, q)
B = Point(k, l)
if A.abs < B.abs:
A, B = B, A
while B.nonzero:
D = A.divround(B)
B, A = A.subtract(D.multiply(B)), B
return A.isone
group_1 = [[4, 0, 3, 0], [4, 0, 0, 3], [0, 4, 3, 0], [0, 4, 0, 3],
[5, 0, 3, 1], [5, 0, 1, 3], [0, 5, 3, 1], [0, 5, 1, 3],
[6, 0, 4, 1], [6, 0, 1, 4], [0, 6, 4, 1], [0, 6, 1, 4],
[6, 0, 3, 2], [6, 0, 2, 3], [0, 6, 3, 2], [0, 6, 2, 3],
[7, 0, 5, 0], [7, 0, 0, 5],
[5, 3, 5, 0], [5, 3, 0, 5], [3, 5, 5, 0], [3, 5, 0, 5],
[0, 7, 5, 0], [0, 7, 0, 5], [7, 0, 4, 2], [7, 0, 2, 4],
[5, 3, 4, 2], [5, 3, 2, 4], [3, 5, 4, 2], [3, 5, 2, 4],
[0, 7, 4, 2], [0, 7, 2, 4],
[7, 0, 3, 3],
[5, 3, 3, 3], [3, 5, 3, 3],
[0, 7, 3, 3],
[8, 0, 6, 0], [8, 0, 0, 6], [0, 8, 6, 0], [0, 8, 0, 6],
[8, 0, 4, 3], [8, 0, 3, 4], [0, 8, 4, 3], [0, 8, 3, 4],
[9, 0, 6, 1], [9, 0, 1, 6], [0, 9, 6, 1], [0, 9, 1, 6],
[10, 0, 7, 1], [10, 0, 1, 7], [0, 10, 7, 1], [0, 10, 1, 7],
[10, 0, 6, 2], [10, 0, 2, 6], [0, 10, 6, 2], [0, 10, 2, 6],
[11, 0, 8, 0], [11, 0, 0, 8], [0, 11, 8, 0], [0, 11, 0, 8],
[11, 0, 7, 2], [11, 0, 2, 7], [0, 11, 7, 2], [0, 11, 2, 7],
[11, 0, 6, 3], [11, 0, 3, 6], [0, 11, 6, 3], [0, 11, 3, 6],
[12, 0, 9, 0], [12, 0, 0, 9], [0, 12, 9, 0], [0, 12, 0, 9],
[12, 0, 8, 1], [12, 0, 1, 8], [0, 12, 8, 1], [0, 12, 1, 8],
[12, 0, 8, 2], [12, 0, 2, 8], [0, 12, 8, 2], [0, 12, 2, 8],
[12, 0, 7, 3], [12, 0, 3, 7], [0, 12, 7, 3], [0, 12, 3, 7],
[12, 0, 6, 4], [12, 0, 4, 6], [0, 12, 6, 4], [0, 12, 4, 6],
[12, 0, 5, 5], [0, 12, 5, 5],
[2, 1, 2, 0], [2, 1, 0, 2], [1, 2, 2, 0], [1, 2, 0, 2],
[3, 1, 2, 1], [3, 1, 1, 2], [1, 3, 2, 1], [1, 3, 1, 2],
[4, 1, 2, 2], [1, 4, 2, 2],
[5, 1, 4, 0], [5, 1, 0, 4], [1, 5, 4, 0], [1, 5, 0, 4],
[6, 1, 5, 0], [6, 1, 0, 5], [1, 6, 5, 0], [1, 6, 0, 5],
[7, 1, 5, 1], [7, 1, 1, 5], [1, 7, 5, 1], [1, 7, 1, 5],
[8, 1, 6, 1], [8, 1, 1, 6], [1, 8, 6, 1], [1, 8, 1, 6],
[8, 1, 5, 2], [8, 1, 2, 5], [1, 8, 5, 2], [1, 8, 2, 5],
[8, 1, 4, 3], [8, 1, 3, 4], [1, 8, 4, 3], [1, 8, 3, 4],
[9, 1, 7, 0],
[9, 1, 5, 3], [9, 1, 3, 5],
[9, 1, 0, 7],
[6, 5, 7, 0],
[6, 5, 5, 3], [6, 5, 3, 5],
[6, 5, 0, 7], [5, 6, 7, 0],
[5, 6, 5, 3], [5, 6, 3, 5],
[5, 6, 0, 7],
[1, 9, 7, 0],
[1, 9, 5, 3], [1, 9, 3, 5],
[1, 9, 0, 7],
[9, 1, 6, 2], [9, 1, 2, 6],
[6, 5, 6, 2], [6, 5, 2, 6],
[5, 6, 6, 2], [5, 6, 2, 6],
[1, 9, 6, 2], [1, 9, 2, 6],
[9, 1, 4, 4],
[6, 5, 4, 4], [5, 6, 4, 4],
[1, 9, 4, 4],
[10, 1, 8, 0], [10, 1, 0, 8], [1, 10, 8, 0], [1, 10, 0, 8],
[10, 1, 7, 1], [10, 1, 1, 7], [1, 10, 7, 1], [1, 10, 1, 7],
[10, 1, 6, 3], [10, 1, 3, 6], [1, 10, 6, 3], [1, 10, 3, 6],
[10, 1, 5, 4], [10, 1, 4, 5], [1, 10, 5, 4], [1, 10, 4, 5],
[11, 1, 8, 1], [11, 1, 1, 8], [1, 11, 8, 1], [1, 11, 1, 8],
[11, 1, 6, 4], [11, 1, 4, 6], [1, 11, 6, 4], [1, 11, 4, 6],
[11, 1, 5, 5], [1, 11, 5, 5],
[2, 2, 2, 1], [2, 2, 1, 2],
[4, 2, 4, 0], [4, 2, 0, 4], [2, 4, 4, 0], [2, 4, 0, 4],
[5, 2, 4, 1], [5, 2, 1, 4], [2, 5, 4, 1], [2, 5, 1, 4],
[6, 2, 4, 2], [6, 2, 2, 4], [2, 6, 4, 2], [2, 6, 2, 4],
[6, 2, 3, 3], [2, 6, 3, 3],
[7, 2, 6, 0], [7, 2, 0, 6], [2, 7, 6, 0], [2, 7, 0, 6],
[7, 2, 5, 2], [7, 2, 2, 5], [2, 7, 5, 2], [2, 7, 2, 5],
[7, 2, 4, 3], [7, 2, 3, 4], [2, 7, 4, 3], [2, 7, 3, 4],
[8, 2, 7, 0], [8, 2, 5, 3], [8, 2, 3, 5], [8, 2, 0, 7],
[2, 8, 7, 0],
[2, 8, 5, 3], [2, 8, 3, 5],
[2, 8, 0, 7],
[8, 2, 6, 1], [8, 2, 1, 6], [2, 8, 6, 1], [2, 8, 1, 6],
[8, 2, 4, 4], [2, 8, 4, 4],
[9, 2, 7, 1], [9, 2, 1, 7], [2, 9, 7, 1], [2, 9, 1, 7]]
group_2 = [[9, 2, 6, 2], [9, 2, 2, 6], [2, 9, 6, 2], [2, 9, 2, 6]]
group_3 = [[9, 2, 5, 4], [9, 2, 4, 5], [2, 9, 5, 4], [2, 9, 4, 5],
[10, 2, 8, 0], [10, 2, 0, 8], [2, 10, 8, 0], [2, 10, 0, 8],
[10, 2, 8, 1], [10, 2, 1, 8], [2, 10, 8, 1], [2, 10, 1, 8],
[10, 2, 7, 2], [10, 2, 2, 7], [2, 10, 7, 2], [2, 10, 2, 7],
[10, 2, 6, 3], [10, 2, 3, 6], [2, 10, 6, 3], [2, 10, 3, 6],
[4, 3, 4, 1], [4, 3, 1, 4], [3, 4, 4, 1], [3, 4, 1, 4],
[4, 3, 3, 2], [4, 3, 2, 3], [3, 4, 3, 2], [3, 4, 2, 3],
[6, 3, 6, 0], [6, 3, 0, 6], [3, 6, 6, 0], [3, 6, 0, 6],
[6, 3, 4, 3], [6, 3, 3, 4], [3, 6, 4, 3], [3, 6, 3, 4],
[7, 3, 6, 1], [7, 3, 1, 6], [3, 7, 6, 1], [3, 7, 1, 6],
[8, 3, 7, 0], [8, 3, 5, 3], [8, 3, 3, 5], [8, 3, 0, 7],
[3, 8, 7, 0], [3, 8, 5, 3], [3, 8, 3, 5], [3, 8, 0, 7],
[8, 3, 7, 1], [8, 3, 1, 7], [3, 8, 7, 1], [3, 8, 1, 7],
[8, 3, 6, 2], [8, 3, 2, 6], [3, 8, 6, 2], [3, 8, 2, 6],
[9, 3, 8, 0], [9, 3, 0, 8], [3, 9, 8, 0], [3, 9, 0, 8],
[9, 3, 7, 2], [9, 3, 2, 7], [3, 9, 7, 2], [3, 9, 2, 7],
[9, 3, 6, 3], [9, 3, 3, 6], [3, 9, 6, 3], [3, 9, 3, 6],
[9, 3, 5, 4], [9, 3, 4, 5], [3, 9, 5, 4], [3, 9, 4, 5],
[4, 4, 5, 0], [4, 4, 0, 5],
[4, 4, 4, 2], [4, 4, 2, 4],
[4, 4, 3, 3],
[5, 4, 6, 0], [5, 4, 0, 6], [4, 5, 6, 0], [4, 5, 0, 6],
[5, 4, 5, 1], [5, 4, 1, 5], [4, 5, 5, 1], [4, 5, 1, 5],
[6, 4, 6, 1], [6, 4, 1, 6], [4, 6, 6, 1], [4, 6, 1, 6],
[6, 4, 5, 2], [6, 4, 2, 5], [4, 6, 5, 2], [4, 6, 2, 5],
[7, 4, 7, 0],
[7, 4, 5, 3], [7, 4, 3, 5],
[7, 4, 0, 7],
[4, 7, 7, 0],
[4, 7, 5, 3], [4, 7, 3, 5],
[4, 7, 0, 7],
[7, 4, 6, 2], [7, 4, 2, 6], [4, 7, 6, 2], [4, 7, 2, 6],
[7, 4, 4, 4], [4, 7, 4, 4],
[8, 4, 8, 0], [8, 4, 0, 8], [4, 8, 8, 0], [4, 8, 0, 8],
[8, 4, 7, 1], [8, 4, 1, 7], [4, 8, 7, 1], [4, 8, 1, 7],
[8, 4, 6, 3], [8, 4, 3, 6], [4, 8, 6, 3], [4, 8, 3, 6],
[8, 4, 5, 4], [8, 4, 4, 5], [4, 8, 5, 4], [4, 8, 4, 5],
[5, 5, 6, 1], [5, 5, 1, 6],
[5, 5, 5, 2], [5, 5, 2, 5],
[7, 5, 8, 0], [7, 5, 0, 8], [5, 7, 8, 0], [5, 7, 0, 8],
[7, 5, 7, 1], [7, 5, 1, 7], [5, 7, 7, 1], [5, 7, 1, 7],
[7, 5, 6, 3], [7, 5, 3, 6], [5, 7, 6, 3], [5, 7, 3, 6],
[7, 5, 5, 4], [7, 5, 4, 5], [5, 7, 5, 4], [5, 7, 4, 5],
[6, 6, 7, 1], [6, 6, 1, 7],
[6, 6, 6, 3], [6, 6, 3, 6],
[6, 6, 5, 4], [6, 6, 4, 5]]
G1 = [euclid_2021(*pqkl) for pqkl in group_1]
G2 = [euclid_2021(*pqkl) for pqkl in [group_2[i] for i in (0, 2, 3)]]
G3 = [euclid_2021(*pqkl) for pqkl in group_3]
plt.plot(G1)
plt.plot(G2, '-r', linewidth=6)
plt.plot(G3)
plt.ylim(-0.1, 1.1)
plt.show()
