如何枚举满足附加约束条件的 x^2 + y^2 = z^2 - 1？

Question

假设N是一个数字 (10<=N<=10^5)。

我需要将它分成三个数字(x,y,z)，使得满足以下条件。

1. x<=y<=z 
2. x^2+y^2=z^2-1;
3. x+y+z<=N

我需要找出在一种方法下可以从给定数字中获得多少个组合。

我已经尝试了以下方法，但对于较大的数字来说需要很长时间，并导致超时。

int N= Int32.Parse(Console.ReadLine());
List<String> res = new List<string>();

//x<=y<=z
int mxSqrt = N - 2;
int a = 0, b = 0;
for (int z = 1; z <= mxSqrt; z++)
{
    a = z * z;
    for (int y = 1; y <= z; y++)
    {
        b = y * y;
        for (int x = 1; x <= y; x++)
        {
            int x1 = b + x * x;
            int y1 = a - 1;
            if (x1 == y1 && ((x + y + z) <= N))
            {
                res.Add(x + "," + y + "," + z);
            }
        }
    }
}
Console.WriteLine(res.Count());

我的问题：

我的解决方案需要更长的时间用于更大的数字（我认为是由于循环嵌套），如何改进它？

有没有更好的方法来解决同样的问题？

- A_Sk

---

关于一些优化：在最外层循环中执行j*j和i*i，而不是在最内层的k循环中执行。 - Peter B

如果你想要继续浪费 CPU 周期，那么请务必为每个 k 重新计算 j*j 和 i*i（这些与 k 无关)。 - Peter B

你只需要两个循环而不是三个，z = sqrt(....)。 - MrD at KookerellaLtd

一个提示：当要求以“x，y，z”给出时，请不要使用“i，j，k”作为循环变量，“x，y”作为其他值。 - Hans Kesting

个人而言，我会取一个变量，比如说y，然后计算“y^2 + 1”的质因数...这基本上给出了所有可能的x和z值...然后过滤掉不符合其他条件的值...(正如我所说，你可能可以预先计算所有这些值，它们与N无关)。 - MrD at KookerellaLtd

Answer 1

这里有一种方法可以枚举三元组，而不是穷举测试它们，使用数论，如此描述：https://mathoverflow.net/questions/29644/enumerating-ways-to-decompose-an-integer-into-the-sum-of-two-squares 由于这个数学方法花了我一些时间去理解和实现（收集了一些在上面提到的代码），而且我对这个主题并没有太多的权威感，所以我将其留给读者自己研究。这基于将数字表示为高斯整数共轭。 （a + bi）*（a-bi）= a ^ 2 + b ^ 2 。我们首先将数字z ^ 2-1分解成质数，将质数分解为高斯共轭，并找到不同的表达式，我们扩展和简化以获得a + bi，然后可以升级为a ^ 2 + b ^ 2 。

阅读平方和函数的好处是发现我们可以通过使用一个检查来排除任何包含形式为4k + 3的质数的奇次幂的候选者z^2 - 1。仅仅使用这个检查，我就能够使用下面的Rosetta质因数分解代码将Prune在10^5上的循环从214秒减少到19秒（在repl.it上）。

这里的实现只是一个演示。它没有处理或优化限制x和y。相反，它只是在枚举过程中进行。在这里尝试一下。

Python代码：

# https://math.stackexchange.com/questions/5877/efficiently-finding-two-squares-which-sum-to-a-prime
def mods(a, n):
    if n <= 0:
        return "negative modulus"
    a = a % n
    if (2 * a > n):
        a -= n
    return a

def powmods(a, r, n):
    out = 1
    while r > 0:
        if (r % 2) == 1:
            r -= 1
            out = mods(out * a, n)
        r /= 2
        a = mods(a * a, n)
    return out

def quos(a, n):
    if n <= 0:
        return "negative modulus"
    return (a - mods(a, n))/n

def grem(w, z):
    # remainder in Gaussian integers when dividing w by z
    (w0, w1) = w
    (z0, z1) = z
    n = z0 * z0 + z1 * z1
    if n == 0:
        return "division by zero"
    u0 = quos(w0 * z0 + w1 * z1, n)
    u1 = quos(w1 * z0 - w0 * z1, n)
    return(w0 - z0 * u0 + z1 * u1,
           w1 - z0 * u1 - z1 * u0)

def ggcd(w, z):
    while z != (0,0):
        w, z = z, grem(w, z)
    return w

def root4(p):
    # 4th root of 1 modulo p
    if p <= 1:
        return "too small"
    if (p % 4) != 1:
        return "not congruent to 1"
    k = p/4
    j = 2
    while True:
        a = powmods(j, k, p)
        b = mods(a * a, p)
        if b == -1:
            return a
        if b != 1:
            return "not prime"
        j += 1

def sq2(p):
    if p % 4 != 1:
      return "not congruent to 1 modulo 4"
    a = root4(p)
    return ggcd((p,0),(a,1))

# https://rosettacode.org/wiki/Prime_decomposition#Python:_Using_floating_point
from math import floor, sqrt

def fac(n):
    step = lambda x: 1 + (x<<2) - ((x>>1)<<1)
    maxq = long(floor(sqrt(n)))
    d = 1
    q = n % 2 == 0 and 2 or 3 
    while q <= maxq and n % q != 0:
        q = step(d)
        d += 1
    return q <= maxq and [q] + fac(n//q) or [n]

# My code...
# An answer for  https://dev59.com/glQJ5IYBdhLWcg3wT0LJ

from collections import Counter
from itertools import product
from sympy import I, expand, Add

def valid(ps):
  for (p, e) in ps.items():
    if (p % 4 == 3) and (e & 1):
      return False
  return True

def get_sq2(p, e):
  if p == 2:
    if e & 1:
      return [2**(e / 2), 2**(e / 2)]
    else:
      return [2**(e / 2), 0]
  elif p % 4 == 3:
    return [p, 0]
  else:
    a,b = sq2(p)
    return [abs(a), abs(b)]

def get_terms(cs, e):
  if e == 1:
    return [Add(cs[0], cs[1] * I)]
  res = [Add(cs[0], cs[1] * I)**e]
  for t in xrange(1, e / 2 + 1):
    res.append(
      Add(cs[0] + cs[1]*I)**(e-t) * Add(cs[0] - cs[1]*I)**t)
  return res

def get_lists(ps):
  items = ps.items()
  lists = []
  for (p, e) in items:
    if p == 2:
      a,b = get_sq2(2, e)
      lists.append([Add(a, b*I)])
    elif p % 4 == 3:
      a,b = get_sq2(p, e)
      lists.append([Add(a, b*I)**(e / 2)])
    else:
      lists.append(get_terms(get_sq2(p, e), e))
  return lists


def f(n):
  for z in xrange(2, n / 2):
    zz = (z + 1) * (z - 1)
    ps = Counter(fac(zz))
    is_valid = valid(ps)
    if is_valid:
      print "valid (does not contain a prime of form\n4k + 3 with an odd power)"
      print "z: %s, primes: %s" % (z, dict(ps))
      lists = get_lists(ps)
      cartesian = product(*lists)
      for element in cartesian:
        print "prime square decomposition: %s" % list(element)
        p = 1
        for item in element:
          p *= item
        print "complex conjugates: %s" % p
        vals = p.expand(complex=True, evaluate=True).as_coefficients_dict().values()
        x, y = vals[0], vals[1] if len(vals) > 1 else 0
        print "x, y, z: %s, %s, %s" % (x, y, z)
        print "x^2 + y^2, z^2-1: %s, %s" % (x**2 + y**2, z**2 - 1)
      print ''

if __name__ == "__main__":
  print f(100)

输出：

valid (does not contain a prime of form
4k + 3 with an odd power)
z: 3, primes: {2: 3}
prime square decomposition: [2 + 2*I]
complex conjugates: 2 + 2*I
x, y, z: 2, 2, 3
x^2 + y^2, z^2-1: 8, 8

valid (does not contain a prime of form
4k + 3 with an odd power)
z: 9, primes: {2: 4, 5: 1}
prime square decomposition: [4, 2 + I]
complex conjugates: 8 + 4*I
x, y, z: 8, 4, 9
x^2 + y^2, z^2-1: 80, 80

valid (does not contain a prime of form
4k + 3 with an odd power)
z: 17, primes: {2: 5, 3: 2}
prime square decomposition: [4 + 4*I, 3]
complex conjugates: 12 + 12*I
x, y, z: 12, 12, 17
x^2 + y^2, z^2-1: 288, 288

valid (does not contain a prime of form
4k + 3 with an odd power)
z: 19, primes: {2: 3, 3: 2, 5: 1}
prime square decomposition: [2 + 2*I, 3, 2 + I]
complex conjugates: (2 + I)*(6 + 6*I)
x, y, z: 6, 18, 19
x^2 + y^2, z^2-1: 360, 360

valid (does not contain a prime of form
4k + 3 with an odd power)
z: 33, primes: {17: 1, 2: 6}
prime square decomposition: [4 + I, 8]
complex conjugates: 32 + 8*I
x, y, z: 32, 8, 33
x^2 + y^2, z^2-1: 1088, 1088

valid (does not contain a prime of form
4k + 3 with an odd power)
z: 35, primes: {17: 1, 2: 3, 3: 2}
prime square decomposition: [4 + I, 2 + 2*I, 3]
complex conjugates: 3*(2 + 2*I)*(4 + I)
x, y, z: 18, 30, 35
x^2 + y^2, z^2-1: 1224, 1224

- גלעד ברקן

---

虽然对于 N = 10 和 N = 100 找到的组合数量与此处发布的其他解决方案相匹配，但是当涉及到 N = 1000 及更高时，计数远远不足。这是所有其他解决方案的缺陷还是这个解决方案的缺陷？ - cdlane

对于 f(1000)，上述代码会产生结果 x, y, z: 112, 476, 489，总和为1077。同样地，x, y, z: 242, 418, 483，总和为1143。我是否误解了如何运行它？ - cdlane

@cdlane 正如我在答案中所述，“这里的实现只是一个演示。它没有处理或优化限制x和y的方法。相反，它只是边枚举边进行。” - גלעד ברקן

@cdlane 我并不是想完全回答OP的要求。我的回答更多是为了建立其背后的思想。根据math.stackexchange的答案，遵循这个理论应该可以列举出所有可能的三元组。 - גלעד ברקן

明白了，我很欣赏它的意图以及如果完成的话它的潜在效率。我看到它获得了OP发布的赏金，但我想知道为什么它既没有满足赏金的标准。 - cdlane

Answer 2

我没有时间进行充分的测试，但似乎产生了与你的代码相同的结果（在100 -> 6个结果和1000 -> 55个结果）。

当 N=1000 时，时间为 2ms，而不使用 List 的情况下比你的 144ms 更快。

当 N=10000 时，时间为 28ms。

var N = 1000;
var c = 0;

for (int x = 2; x < N; x+=2)
{
    for (int y = x; y < (N - x); y+=2)
    {
        long z2 = x * x + y * y + 1;
        int z = (int) Math.Sqrt(z2);
        if (x + y + z > N)
            break;
        if (z * z == z2)
            c++;
    }
}

Console.WriteLine(c);

- NotFound

---

@Dukeling 是的，我也是这么想的，但因为谷歌给了我一个不正确的结果，所以我还是选择了它。我会更新答案。 - NotFound

你介意我更新我的答案并列出你的最新优化并给你信用吗？我使用你的解决方案在N = 100000上获得了9秒，如果使用并行，则需2.1秒，尽管我在自己的解决方案中获得了1.3秒。我不确定你的边界是否比我的更有效，但我认为你的解决方案可能是N范围的前50％中最好的。 - Mat

@Mat 确认无误。我也根据 @Dukeling 的反馈更改了计算方法，检查完美平方根。对于一些非常大的数字，似乎会得到不正确的结果（例如 3999680306388005621 来源）。 - NotFound

你应该注意的另一个重要细节是使用 long。如果没有错误，由于 int 溢出到负数，对于更高的 N 值，你将得到错误的结果。 - Mat

Answer 3

< p > x和y的范围是问题的重要部分。我个人使用了这个Wolfram Alpha查询，并检查了变量的确切形式。

感谢@Bleep-Bloop和评论中提供的非常优雅的边界优化，即x < n和x <= y < n - x。结果相同，时间几乎相同。

此外，由于x和y的可能值仅为正偶数，因此我们可以将循环迭代次数减半。

为了进一步优化，由于我们计算了x的上限，我们构建了一个所有可能的x值的列表，并进行并行计算。这在较高的N值上节省了大量时间，但对于较小的值速度较慢，因为需要并行处理的开销。

以下是最终代码：

非并行版本，使用int值：

List<string> res = new List<string>();
int n2 = n * n;

double maxX = 0.5 * (2.0 * n - Math.Sqrt(2) * Math.Sqrt(n2 + 1));

for (int x = 2; x < maxX; x += 2)
{
    int maxY = (int)Math.Floor((n2 - 2.0 * n * x - 1.0) / (2.0 * n - 2.0 * x));

    for (int y = x; y <= maxY; y += 2)
    {
        int z2 = x * x + y * y + 1;
        int z = (int)Math.Sqrt(z2);

        if (z * z == z2 && x + y + z <= n)
            res.Add(x + "," + y + "," + z);
    }
}

使用long值的并行版本：

using System.Linq;

...

// Use ConcurrentBag for thread safety
ConcurrentBag<string> res = new ConcurrentBag<string>();
long n2 = n * n;

double maxX = 0.5 * (2.0 * n - Math.Sqrt(2) * Math.Sqrt(n2 + 1L));

// Build list to parallelize
int nbX = Convert.ToInt32(maxX);
List<int> xList = new List<int>();
for (int x = 2; x < maxX; x += 2)
    xList.Add(x);

Parallel.ForEach(xList, x =>
{
    int maxY = (int)Math.Floor((n2 - 2.0 * n * x - 1.0) / (2.0 * n - 2.0 * x));

    for (long y = x; y <= maxY; y += 2)
    {
        long z2 = x * x + y * y + 1L;
        long z = (long)Math.Sqrt(z2);

        if (z * z == z2 && x + y + z <= n)
            res.Add(x + "," + y + "," + z);
    }
});

在 i5-8400 CPU 上单独运行时，我得到了以下结果：

N: 10； 解决方案：1； 经过时间：0.03 毫秒 （非并行，int）

N: 100； 解决方案：6； 经过时间：0.05 毫秒 （非并行，int）

N: 1000； 解决方案：55； 经过时间：0.3 毫秒 （非并行，int）

N: 10000； 解决方案：543； 经过时间：13.1 毫秒 （非并行，int）

N: 100000； 解决方案：5512； 经过时间：849.4 毫秒 （并行，long）

---

当 N 大于 36340 时，必须使用 long，因为当它被平方时，它会超过 int 的最大值。最后，当使用 int 时，当 N 约为 23000 时，并行版本开始比简单版本更好。

- Mat

---

我不确定 Wolphram Alpha 的展开和因式分解是否百分之百可靠。我曾经遇到过一个错误。 - גלעד ברקן

@גלעדברקן 并不总是100%正确的。但对于这个问题，我已经得到了正确的结果和更短的时间，所以我认为它是可以的。如果您发现更好的界限和/或解决方案，请随时告知我（我们）！ - Mat

Answer 4

在Python中有一个简单的改进（将其转换为基于C的代码更快的等效方法留给读者自己练习）。为了得到准确的计算时间，我在之前的运行中移除了打印解决方案本身的部分。

• 使用外循环来处理一个自由变量（我选择了z），它只受其与N的关系所限制。
• 使用内循环（我选择了y）受外循环索引的约束。
• 第三个变量可直接根据需求2进行计算。

计时结果:

-------------------- 10 
 1 solutions found in 2.3365020751953125e-05  sec.
-------------------- 100 
 6 solutions found in 0.00040078163146972656  sec.
-------------------- 1000 
 55 solutions found in 0.030081748962402344  sec.
-------------------- 10000 
 543 solutions found in 2.2078349590301514  sec.
-------------------- 100000 
 5512 solutions found in 214.93411707878113  sec.

对于大规模的情况，需要花费3:35的时间，再加上您整理和/或打印结果的时间。

如果您需要更快的代码（这仍然是一种相当暴力的方法），请研究丢番图方程和参数化，以生成给定z^2 - 1目标值的(y, x)对。

import math
import time

def break3(N):
    """
    10 <= N <= 10^5
    return x, y, z triples such that:
        x <= y <= z
        x^2 + y^2 = z^2 - 1        
        x + y + z <= N
    """

    """
    Observations:
    z <= x + y
    z < N/2
    """

    count = 0
    z_limit = N // 2
    for z in range(3, z_limit):

        # Since y >= x, there's a lower bound on y
        target = z*z - 1
        ymin = int(math.sqrt(target/2))
        for y in range(ymin, z):
            # Given y and z, compute x.
            # That's a solution iff x is integer.
            x_target = target - y*y
            x = int(math.sqrt(x_target))
            if x*x == x_target and x+y+z <= N:
                # print("solution", x, y, z)
                count += 1

    return count


test = [10, 100, 1000, 10**4, 10**5]
border = "-"*20

for case in test: 
    print(border, case)
    start = time.time()
    print(break3(case), "solutions found in", time.time() - start, "sec.")

- Prune

---

是的，我看到了无效的输出！在63个解决方案中，只有55个在总和范围内。一种简单的方法是单行检查，现在已经在答案中了。等我有半个小时的时间，我会解决根本问题（内部循环的上限）。 - Prune

我将保留这个解决方案。是的，我们可以适当限制上限，但在当前形式下，程序更易于阅读和维护。 - Prune

我不确定是否漏掉了什么，但是对于 N = 10 不是有两个解吗？(x, y, z) = (0, 0, 1) 或者 (2, 2, 3)。除非假设1不在解的范围内（OP从1开始循环，所以可能...）。 - Mat

@Mat 排除退化解是原始问题发布的一部分。 - Prune

我成功将你的循环从10^5的214秒缩短到了19秒（在repl.it上，可以查看我的答案：）。 - גלעד ברקן

我不需要运行它；我在欧拉计划中使用了那个算法的兄弟。我点赞了，干得好！ - Prune

Answer 5

#include<iostream>
#include<math.h>
int main()
{
    int N = 10000;
    int c = 0;
    for (int x = 2; x < N; x+=2)
    {
        for (int y = x; y < (N - x); y+=2)
        {
            auto z = sqrt(x * x + y * y + 1);
            if(x+y+z>N){
                break;
            }
            if (z - (int) z == 0)
            {
                c++;
            }
        }
    }
    std::cout<<c;
}

这是我的解决方案。在测试此问题的先前解决方案时，我发现x，y始终为偶数，而z为奇数。我不知道背后的数学本质，目前正在尝试弄清楚。"Original Answer"翻译成"最初的答案"。

- Justice_Lords

---

奇偶性是平方数模4的性质直接导致的结果。 - Prune

Answer 6

我希望用C#完成并覆盖所有测试用例，这些用例基于问题中提供的条件。

基本代码已转换为long以处理N<=100000的上限，并使用了我能想到的所有优化。我使用了@Mat的Wolfram Alpha查询的替代形式尽可能地进行预计算。此外，我还进行了最小完美平方测试，以避免在上限时调用数百万个sqrt()函数。

public static void Main()
{
    int c = 0;

    long N = long.Parse(Console.ReadLine());
    long N_squared = N * N;

    double half_N_squared = N_squared / 2.0 - 0.5;
    double x_limit = N - Math.Sqrt(2) / 2.0 * Math.Sqrt(N_squared + 1);

    for (long x = 2; x < x_limit; x += 2)
    {
        long x_squared = x * x + 1;

        double y_limit = (half_N_squared - N * x) / (N - x);

        for (long y = x; y < y_limit; y += 2)
        {
            long z_squared = x_squared + y * y;
            int digit = (int) z_squared % 10;

            if (digit == 3 || digit == 7)
            {
                continue;  // minimalist non-perfect square elimination
            }

            long z = (long) Math.Sqrt(z_squared);

            if (z * z == z_squared)
            {
                c++;
            }
        }
    }

    Console.WriteLine(c);
}

我遵循了这个趋势，并且像问题的提出者的代码中暗示的那样省略了"退化解决方案"，但并没有明确说明。最初的回答。

- cdlane