让我们N是一个数字(10<=N<=10^5) .

我必须把它分成 3 个数字 (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());

我的问题:

我的解决方案是花时间获得更大的数字(我认为这是 for 循环)，我该如何改进它？

是否有更好的方法？

最佳答案

这里有一个枚举三元组的方法，而不是详尽地测试它们，使用这里描述的数论: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。

阅读 Sum of Squares Function 的好处发现我们可以排除任何候选 z^2 - 1，它包含形式为 4k + 3 且具有奇次幂的素数。仅使用该检查，我就能够使用下面的 Rosetta 质因数分解代码将 10^5 上的 Prune 循环从 214 秒减少到 19 秒(在 repl.it 上)。

这里的实现只是一个演示。它没有限制 x 和 y 的处理或优化。相反，它只是在进行时进行枚举。一起玩吧 here .

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://stackoverflow.com/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

关于c# - 如何枚举 x^2 + y^2 = z^2 - 1(有附加约束)，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/54110614/