为什么爱森斯坦整数对欧几里得求最大公分母的方法的这种实现会因为这一点而卡住?

计算科学 Python 算法 麻木的 复数
2021-12-16 06:49:43

我的数学 SE 问题确定一对旋转六边形格子中的重合点是否最接近原点?解释了我遇到的问题。我不会在这里详细再现整个事情,但这张图片以图形方式说明了问题。

六边形重合格

  • 右侧:定义闭合点 (5, 6) 和开圆点 (5, 3) 的重合点(见红色箭头)是该重合点阵中距原点最近的六个点之一。
  • 左侧:由 (6, 5) 和 (5, 3) (红色箭头)类似定义的重合点远不是离原点最近的重合点。

唯一的区别是我们通过将闭合点阵旋转大约 6 度将 (6, 5) 替换为 (5, 6)。

我的问题的答案建议将众所周知的欧几里得算法的修改应用于该问题。

它现在应用于两个复数Eisenstein 整数,而不是两个实整数,这是复平面中点的六边形分布。有关更多信息,请参阅:

为了尝试从那个答案到 Python 中实现这个非常非常好的和清晰描述的分步方法,我定义了一个Point类并为该类提供了以下方法:.subtract(P), .multiply(P), .divround(P) , .prod(P) and.sub(P)` 尽我所能复制该答案中描述的内容。

请注意,一个Point对象P由两个整数 P.p定义,P.q并且所有方法都对这些整数进行操作并返回Point. 大多数操作是加法和乘法,唯一暂时进入浮点数的是.divround(P).

下面是一个包含 300 多个 (p, q), (k, l) 的测试程序,由函数测试euclid_2021(p, q, k, l)

到目前为止,不太严格的测试表明该算法给出了我认为正确的答案,至少它似乎与我之前编写的一个非常长、丑陋、看起来很疯狂的脚本非常吻合,但我没有已知的正确方法来检查两者。

但是在这里,我只询问一个明显的无限循环,该方法仅针对以下 300 多个案例中的一个陷入困境。

Group_1Group_3快速运行并返回表示它是否最接近的布尔值。

Group_2只要我省略第二种情况(9、2、2、6),一切都很好但是对于这种情况,循环会无限期地运行。

关于为什么这个案例失败的任何数学想法?

当仅测试“坏”情况(9、2、2、6)并插入以下打印行时:

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

该算法只是在两种状态之间来回循环:

(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)

第 1、2、3 组的真/假输出,但跳过“坏”情况:

真/假输出

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()
1个回答

divround,pintqint必须是 and 的整数部分pfracqfrac但是,据我所知,在 Python 中,将 a 转换float为 anint总是将其舍入为 0,即正浮点数向下舍入,但负浮点数向上舍入,int(-2.7)-2 也是如此,不是-3。代替

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

你应该有

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

divround现在想起来,简单的应该就够了return Point(round(pfrac), round(qfrac))结果 D ineuclid_2021不一定是最接近 A/B 的 Eisenstein 整数,但它会使得 ∥D-A/B∥<1,这对于算法来说已经足够了。(D−A/B 的“实数”和“虚数”部分都在 -0.5 和 0.5 之间,因此其范数最多为 sqrt(0.5^2+0.5^2+0.5∗0.5) =平方(0.75)。)