如何枚举x ^ 2 + y ^ 2 = z ^ 2

本文介绍了如何枚举x ^ 2 + y ^ 2 = z ^ 2-1-1（具有其他约束）的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！ 问题描述

让 N 为一个数字（10 <== N <== 10 ^ 5）。

我必须将其分解为3个数字（x，y，z），以便验证以下条件。

I have to break it into 3 numbers (x,y,z) such that it validates the following conditions.

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

我必须找出多少组合我可以通过一种方法从给定的数字中得到。

I have to find how many combinations I can get from the given numbers in a method.

我尝试了以下方法，但是要花费更多的时间来获取更大的数字会导致超时。 p>

I have tried as follows but it's taking so much time for a higher number and resulting in a timeout..

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

我的问题：

我的解决方案需要更多时间（我认为这是的循环次数），我该如何改善呢？

My solution is taking time for a bigger number (I think it's the for loops), how can I improve it?

有没有

推荐答案

这里有一种枚举三元组的方法，而不是穷举测试它们，使用数论作为在此处进行描述：> https：// mathoverflow / questions / 29644 /枚举将整数分解为两个平方和的方法

Here's a method that enumerates the triples, rather than exhaustively testing for them, using number theory as described here: mathoverflow/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 。

Since the math took me a while to comprehend and a while to implement (gathering some code that's credited above it), and since I don't feel much of an authority on the subject, I'll leave it for the reader to research. This is based on expressing numbers as Gaussian integer conjugates. (a + bi)*(a - bi) = a^2 + b^2. We first factor the number, z^2 - 1, into primes, decompose the primes into Gaussian conjugates and find different expressions that we expand and simplify to get a + bi, which can be then raised, a^2 + b^2.

对平方和功能的阅读使人振奋，发现我们可以排除包含素数形式为 4k + 3 且质数为奇数的素数的任何候选项 z ^ 2-1-。仅使用该检查，就可以使用下面的Rosetta素数分解代码将10 ^ 5上的Prune循环从214秒减少到19秒（在repl.it上）。

A perk of reading about the Sum of Squares Function is discovering that we can rule out any candidate z^2 - 1 that contains a prime of form 4k + 3 with an odd power. Using that check alone, I was able to reduce Prune's loop on 10^5 from 214 seconds to 19 seconds (on repl.it) using the Rosetta prime factoring code below.

这里的实现只是一个演示。它没有限制 x 和 y 的处理或优化。相反，它只是枚举。在此处播放。

The implementation here is just a demonstration. It does not have handling or optimisation for limiting x and y. Rather, it just enumerates as it goes. Play with it here.

Python代码：

# math.stackexchange/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)) # rosettacode/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 stackoverflow/questions/54110614/ 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