1
$\begingroup$

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.

hexagonal coincidence lattices

  • 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:

True/False output

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()
$\endgroup$

1 Answer 1

1
$\begingroup$

In divround, pint and qint must be the integer parts of pfrac and qfrac, However, as I understand, in Python, converting a float to an int always rounds it towards 0, i.e., positive floats are rounded down, but negative floats are rounded up, so int(-2.7) is -2, not -3. Instead of

pint, qint = [int(thing) for thing in (pfrac, qfrac)]

you should have

pint, qint = [math.floor(thing) for thing in (pfrac, qfrac)]

Now that I think of it, it should be enough for divround to simply return Point(round(pfrac), round(qfrac)). The result D in euclid_2021 would not necessarily be the closest Eisenstein integer to A/B, but it would be such that ∥D−A/B∥<1, which is enough for the algorithm. (Both the "real" and the "u-maginary" parts of D−A/B would be between −0.5 and 0.5, so its norm would be at most sqrt(0.5^2+0.5^2+0.5∗0.5) = sqrt(0.75).)

$\endgroup$
2
  • $\begingroup$ Yes, using math.floor() it now works for all the cases listed including (9, 2, 2, 6). I'll explore the alternative rounding scheme and compare both of my python implementations on a much larger set of test cases. $\endgroup$
    – uhoh
    Jan 3 at 20:52
  • $\begingroup$ I'd planned on finding a way to independently check my result to be 100% confident my script was providing absolutely correct results in all cases before accepting, but I have not yet found a suitable round tuit and it's certainly not a necessary condition for acceptance of this answer. :-) $\endgroup$
    – uhoh
    Mar 2 at 22:56

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.