圆的整数解的算法？

Question

我试图寻找方程的整数解：

y^2 + x^2 = 2n^2

如果我在wolfram alpha中搜索，它们几乎都会立即被发现，即使是非常大的n。当我实施暴力方法时，速度非常慢：

def psearch(n, count):
    for x in range(0, n):
        for y in range(0, n):
            if x*x + y*y == 2*n**2:
                print(x,y)
                count += 1
    return count

所以我假设有一种更快的方法来得到上面方程的所有整数解。我怎样才能在python中做到这一点，从而使它的运行时间大大降低？你知道吗

注：我见过this question但是它是关于在圆中寻找格点的，而不是圆方程的整数解。我还想找到具体的解决方案，而不仅仅是解决方案的数量。你知道吗

编辑：我仍然在寻找一个数量级更快的东西。下面是一个例子：n=5应该有12个整数解来找到那些应该在Wolfram alpha上搜索这个方程的解。你知道吗

编辑二：@victor zen对这个问题给出了惊人的答案。有人能想出一种方法来进一步优化他的解决方案吗？你知道吗

Answer 1

在第一个八分体y>0，x<y内搜索就足够了（四个解(±n, ±n)很明显，通过对称，一个解(x, y)产生8个拷贝(±x, ±y)，(±y, ±x)）。你知道吗

根据单调性，对于一个给定的y，至多有一个x。您可以通过循序渐进地跟随圆弧，减少y，然后调整x来找到它。如果您尽可能严格地维护条件x²+y²≤2n²，那么下面的代码将被优化为只使用基本整数算术（为了效率，使用2x而不是x）。你知道吗

    x, y, d= 2 * n, 2 * n, 0
    while y > 0:
        y, d= y - 2, d - y + 1
        if d < 0:
            x, d= x + 2, d + x + 1
        if d == 0:
            print(x >> 1, '² + ', y >> 1, '² = 2.', n, '²', sep='')

以下是1和100之间n的所有解决方案：

7² + 1² = 2.5²
14² + 2² = 2.10²
17² + 7² = 2.13²
21² + 3² = 2.15²
23² + 7² = 2.17²
28² + 4² = 2.20²
31² + 17² = 2.25²
35² + 5² = 2.25²
34² + 14² = 2.26²
41² + 1² = 2.29²
42² + 6² = 2.30²
46² + 14² = 2.34²
49² + 7² = 2.35²
47² + 23² = 2.37²
51² + 21² = 2.39²
56² + 8² = 2.40²
49² + 31² = 2.41²
63² + 9² = 2.45²
62² + 34² = 2.50²
70² + 10² = 2.50²
69² + 21² = 2.51²
68² + 28² = 2.52²
73² + 17² = 2.53²
77² + 11² = 2.55²
82² + 2² = 2.58²
84² + 12² = 2.60²
71² + 49² = 2.61²
79² + 47² = 2.65²
85² + 35² = 2.65²
89² + 23² = 2.65²
91² + 13² = 2.65²
92² + 28² = 2.68²
98² + 14² = 2.70²
103² + 7² = 2.73²
94² + 46² = 2.74²
93² + 51² = 2.75²
105² + 15² = 2.75²
102² + 42² = 2.78²
112² + 16² = 2.80²
98² + 62² = 2.82²
97² + 71² = 2.85²
113² + 41² = 2.85²
115² + 35² = 2.85²
119² + 17² = 2.85²
123² + 3² = 2.87²
119² + 41² = 2.89²
126² + 18² = 2.90²
119² + 49² = 2.91²
133² + 19² = 2.95²
137² + 7² = 2.97²
124² + 68² = 2.100²
140² + 20² = 2.100²

Answer 2

你可以通过只考虑一个象限和乘以4来优化这个算法。你知道吗

import math
def psearch(n, count):
  for x in range( 0 , 2*n  + 1):
    ysquare = 2*(n**2) - x * x
    if (ysquare <0):
      break
    y = int(math.sqrt(ysquare))

    if ysquare ==  y * y :
      print(x,y)
      count+=1

  return count


print(psearch(13241324,0) * 4)

输出

(1269716, 18682964)
(1643084, 18653836)
(11027596, 15134644)
(12973876, 13503476)
(13241324, 13241324)
(13503476, 12973876)
(15134644, 11027596)
(18653836, 1643084)
(18682964, 1269716)
36

Answer 3

在算法中，搜索所有可能的y值。这是不必要的。这里的诀窍是认识到

y^2 + x^2 = 2n^2

直接暗示

y^2 = 2n^2-x^2

这意味着你只需要检查2n^2-x^2是一个完美的正方形。你可以通过

y2 = 2*n*n - x*x 
#check for perfect square
y = math.sqrt(y2)
if int(y + 0.5) ** 2 == y2:
    #We have a perfect square.

另外，在您的算法中，您只检查x值到n。这是不正确的。由于y^2始终为正或零，我们可以通过将y^2设置为其最低值（即0）来确定需要检查的最高x值。因此，我们需要检查所有满足

x^2 <= 2n^2

减少到

abs(x) <= sqrt(2)*n.

将此与只检查顶部象限的优化结合起来，您就有了一个优化的psearch

def psearch(n):
    count = 0
    top = math.ceil(math.sqrt(2*n*n))
    for x in range(1, top):
        y2 = 2*n*n - x*x 
        #check for perfect square
        y = math.sqrt(y2)
        if int(y + 0.5) ** 2 == y2:
            count+=4
    return count