圆整数解的算法？

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 中做到这一点，以便它的运行时间要低得多？

注意：我见过这个问题，但是它是关于在圆内找到格点，而不是圆方程的整数解。我也有兴趣找到具体的解决方案，而不仅仅是解决方案的数量。

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

编辑 2：@victor zen 对这个问题给出了惊人的答案。谁能想到进一步优化他的解决方案的方法？

Answer 1

在您的算法中，您正在搜索所有可能的 y 值。这是不必要的。这里的诀窍是要意识到y^2 + x^2 = 2n^2直接意味着y^2 = 2n^2-x^2所以这意味着你只需要检查 2n^2-x^2 是一个完美的正方形。你可以这样做y2 = 2*n*n - x*x&nbsp;#check for perfect squarey = math.sqrt(y2)if int(y + 0.5) ** 2 == y2:&nbsp; &nbsp; #We have a perfect square.此外，在您的算法中，您只检查不超过 n 的 x 值。这是不正确的。由于 y^2 将始终为正数或零，我们可以通过将 y^2 设置为其最小值（即 0）来确定我们需要检查的最高 x 值。因此，我们需要检查所有满足的整数 x 值x^2 <= 2n^2这减少到abs(x) <= sqrt(2)*n.将此与仅检查上象限的优化相结合，您将获得优化的 psearchdef psearch(n):&nbsp; &nbsp; count = 0&nbsp; &nbsp; top = math.ceil(math.sqrt(2*n*n))&nbsp; &nbsp; for x in range(1, top):&nbsp; &nbsp; &nbsp; &nbsp; y2 = 2*n*n - x*x&nbsp;&nbsp; &nbsp; &nbsp; &nbsp; #check for perfect square&nbsp; &nbsp; &nbsp; &nbsp; y = math.sqrt(y2)&nbsp; &nbsp; &nbsp; &nbsp; if int(y + 0.5) ** 2 == y2:&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; count+=4&nbsp; &nbsp; return count

Answer 2

y>0在第一个八分圆内搜索就足够了x<y（四个解(±n, ±n)是显而易见的，并且通过对称性，一个解会(x, y)产生 8 个副本(±x, ±y), (±y, ±x)）。通过单调性，对于一个给定y的，至多有一个x。您可以通过逐步跟随圆弧，减少y然后调整来找到它x。如果您x²+y²≤2n²尽可能严格地保持条件，您会得到下面的代码，该代码已优化为仅使用初等整数算术（为了提高效率，2x使用 代替x）。&nbsp; &nbsp; x, y, d= 2 * n, 2 * n, 0&nbsp; &nbsp; while y > 0:&nbsp; &nbsp; &nbsp; &nbsp; y, d= y - 2, d - y + 1&nbsp; &nbsp; &nbsp; &nbsp; if d < 0:&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; x, d= x + 2, d + x + 1&nbsp; &nbsp; &nbsp; &nbsp; if d == 0:&nbsp; &nbsp; &nbsp; &nbsp; &nbsp; &nbsp; print(x >> 1, '² + ', y >> 1, '² = 2.', n, '²', sep='')以下是nbetween1和的所有解决方案100：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 3

您可以通过仅考虑一个象限并乘以 4 来优化此算法。import mathdef psearch(n, count):&nbsp; for x in range( 0 , 2*n&nbsp; + 1):&nbsp; &nbsp; ysquare = 2*(n**2) - x * x&nbsp; &nbsp; if (ysquare <0):&nbsp; &nbsp; &nbsp; break&nbsp; &nbsp; y = int(math.sqrt(ysquare))&nbsp; &nbsp; if ysquare ==&nbsp; y * y :&nbsp; &nbsp; &nbsp; print(x,y)&nbsp; &nbsp; &nbsp; count+=1&nbsp; return countprint(psearch(13241324,0) * 4)输出(1269716, 18682964)(1643084, 18653836)(11027596, 15134644)(12973876, 13503476)(13241324, 13241324)(13503476, 12973876)(15134644, 11027596)(18653836, 1643084)(18682964, 1269716)36