比您的速度快数百倍，而且没有浮点问题。
比 kaya3 的 O(n²) 解决方案快数千倍。
我运行它直到 MAX_NUM = 4,000,000 并没有发现任何结果。花了大约 12 分钟。
利用特殊数字。
这不只是一个普通的 3SUM。这些数字很特殊，我们可以利用它。它们的形式为 ab/c²，其中 (a,b,c) 是原始毕达哥拉斯三元组。
假设我们有一个数字 x=ab/c²，我们想找到另外两个这样的数字加起来为 x：
取消后，分母 c² 和 (fi)² 变为 c²/k 和 (fi)²/m（对于某些整数 k 和 m），我们有 c²/k = (fi) ²/米。令 p 为 c²/k 的最大素因子。然后 p 也除以 (fi)²/m，从而得到 f 或 i。所以至少有一个数字 de/f² 和 gh/i² 有一个可以被 p 整除的分母。我们称它为 y，而另一个称为 z。
那么对于某个 x，我们如何找到拟合的 y 和 z？我们不必为 y 和 z 尝试所有个数字。对于 y，我们只尝试那些分母可以被 p 整除的。对于z？我们将其计算为 x-y 并检查我们是否有该数字（在哈希集中）。
它有多大帮助？如果您天真地尝试所有（小于 x）数字，我的解决方案会计算有多少 y 候选人，以及我的方式有多少 y 候选人以及少多少：
MAX_NUM naive mine % less
--------------------------------------------------
10,000 1,268,028 17,686 98.61
100,000 126,699,321 725,147 99.43
500,000 3,166,607,571 9,926,863 99.69
1,000,000 12,662,531,091 30,842,188 99.76
2,000,000 50,663,652,040 96,536,552 99.81
4,000,000 202,640,284,036 303,159,038 99.85
伪代码
以上代码形式的说明：
h = hashset(numbers)
for x in the numbers:
p = the largest prime factor in the denominator of x
for y in the numbers whose denominator is divisible by p:
z = x - y
if z is in h:
output (x, y, z)
基准测试
各种 MAX_NUM 的时间（以秒为单位）及其结果 n：
MAX_NUM: 10,000 100,000 500,000 1,000,000 2,000,000 4,000,000
=> n: 1,593 15,919 79,582 159,139 318,320 636,617
--------------------------------------------------------------------------------
Original solution 1.6 222.3 - - - -
My solution 0.05 1.6 22.1 71.0 228.0 735.5
kaya3's solution 29.1 2927.1 - - - -
复杂性
这是 O(n²)，实际上可能更好。我不太了解数字的性质，无法对它们进行推理，但上述基准确实使它看起来比 O(n²) 好得多。对于二次运行时，从 n=318,320 到 n=636,617，您预计运行时间增加因子 (636,617/318,320)² ≈ 4.00，但实际增加仅为 735.5/228.0 ≈ 3.23。
我没有针对所有尺寸运行您的，但由于您的增长至少是二次方，因此在 MAX_NUM=4,000,000 时，您的解决方案至少需要 222.3 * (636,617/15,919)² = 355,520 秒，比我的慢 483 倍.同样，kaya3 的速度会比我的慢 6365 倍。
用这个奇怪的把戏浪费时间
Python 的 Fraction 类很简洁，但也很慢。特别是它的散列。转换为元组并散列该元组的速度大约快 34 倍：
>set SETUP="import fractions; f = fractions.Fraction(31459, 271828)"
>python -m timeit -s %SETUP% -n 100000 "hash(f)"
100000 loops, best of 5: 19.8 usec per loop
>python -m timeit -s %SETUP% -n 100000 "hash((f.numerator, f.denominator))"
100000 loops, best of 5: 581 nsec per loop
Its code 说：
[...] 这种方法很昂贵 [...] 为了确保 Fraction 的哈希值与数字相等的整数、浮点数或小数实例的哈希值一致，我们遵循数字哈希的规则文档中概述。
其他操作也有点慢，所以除了输出我不使用Fraction。我改用 (numerator, denominator) 元组。
解决方案代码
from math import gcd
def solve_stefan(triples):
# Prime factorization stuff
largest_prime_factor = [0] * (MAX_NUM + 1)
for i in range(2, MAX_NUM+1):
if not largest_prime_factor[i]:
for m in range(i, MAX_NUM+1, i):
largest_prime_factor[m] = i
def prime_factors(k):
while k > 1:
p = largest_prime_factor[k]
yield p
while k % p == 0:
k //= p
# Lightweight fractions, represented as tuple (numerator, denominator)
def frac(num, den):
g = gcd(num, den)
return num // g, den // g
def sub(frac1, frac2):
a, b = frac1
c, d = frac2
return frac(a*d - b*c, b*d)
class Key:
def __init__(self, triple):
a, b, c = map(int, triple)
self.frac = frac(a*b, c*c)
def __lt__(self, other):
a, b = self.frac
c, d = other.frac
return a*d < b*c
# The search. See notes under the code.
seen = set()
supers = [[] for _ in range(MAX_NUM + 1)]
for triple in sorted(triples, key=Key):
a, b, c = map(int, triple)
x = frac(a*b, c*c)
denominator_primes = [p for p in prime_factors(c) if x[1] % p == 0]
for y in supers[denominator_primes[0]]:
z = sub(x, y)
if z in seen:
yield tuple(sorted(Fraction(*frac) for frac in (x, y, z)))
seen.add(x)
for p in denominator_primes:
supers[p].append(x)
注意事项：
- 我通过增加分数值的三元组，即增加 x 值。
- 我的
denominator_primes 是x 分母的质因数列表。请记住，这是 c²/k，因此它的质因数也必须是 c 的质因数。但是 k 可能已经取消了一些，所以我检查了 c 的质因数并检查它们是否除分母。为什么如此“复杂”而不是仅仅查找 c²/k 的素数？因为这可能非常大。
-
denominator_primes 是递减的，所以 p 就是 denominator_primes[0]。顺便说一句，为什么要使用最大的？因为更大意味着更稀有意味着更少的 y 候选意味着更快。
-
supers[p] 列出了其分母可以被 p 整除的数字。它用于获取 y 候选。
- 用完 x 后，我使用
denominator_primes 将 x 放入 supers 列表中，因此它可以作为未来 x 值的 y。
- 我在循环期间（而不是之前）构建
seen 和supers 以保持它们的小。毕竟，对于 x=y+z 的正数，y 和 z 必须小于 x，因此寻找更大的将是浪费。
验证
如果没有结果，您如何验证结果？据我所知，我们的解决方案都没有找到。所以没有什么可以比较的，除了虚无，这并不完全令人信服。好吧，我的解决方案不依赖于毕达哥拉斯，所以我创建了一组原始三元组并检查了我的解决方案的结果。它计算了与参考实现相同的 25,336 个结果：
def solve_reference(triples):
fractions = {Fraction(int(a) * int(b), int(c)**2)
for a, b, c in triples}
for x, y in combinations_with_replacement(sorted(fractions), 2):
z = x + y
if z in fractions:
yield x, y, z
MIN_NUM = 2
MAX_NUM = 25
def triples():
return list((a, b, c)
for a, b, c in combinations(range(MIN_NUM, MAX_NUM+1), 3)
if gcd(a, gcd(b, c)) == 1)
print(len(triples()), 'input triples')
expect = set(solve_reference(triples()))
print(len(expect), 'results')
output = set(solve_stefan(triples()))
print('output is', ('wrong', 'correct')[output == expect])
输出：
1741 input triples
25336 results
output is correct