python - 提高此搜索的效率，以检查此列表中是否有两个数字相加？

Question

我正在尝试找到最有效的方法，使用Python检查此列表中的任何两个数字是否等于该列表中的另一个数字。我已决定添加一些上下文，以使其更清晰并可能更易于优化。这是我的代码：

import numpy as np
from collections import Counter
from collections import deque


def gen_prim_pyth_trips(limit=None):
    u = np.mat(' 1  2  2; -2 -1 -2; 2 2 3')
    a = np.mat(' 1  2  2;  2  1  2; 2 2 3')
    d = np.mat('-1 -2 -2;  2  1  2; 2 2 3')
    uad = np.array([u, a, d])
    m = np.array([3, 4, 5])
    while m.size:
        m = m.reshape(-1, 3)
        if limit:
            m = m[m[:, 2] <= limit]
        yield from m
        m = np.dot(m, uad)

def find_target(values, target):

    dq = deque(sorted([(val, idx) for idx, val in enumerate(values)]))

    while True:
        if len(dq) < 2:
            return -1

        s =  dq[0][0] + dq[-1][0]

        if s > target:
            dq.pop()
        elif s < target:
            dq.popleft()
        else:
            break
    return dq[0], dq[-1]


ratioList = []

MAX_NUM = 500000

for i in list(gen_prim_pyth_trips(MAX_NUM)):
    ratioList.append((i[0]*i[1])/i[2]**2)
    if find_target(ratioList, (i[0]*i[1])/i[2]**2) != -1:
        print(find_target(ratioList, (i[0]*i[1])/i[2]**2))

gen_prim_pyth_trips()函数来自 here。 “慢”部分出现在三元组生成之后。 find_target来自 here。

当前它运行良好，但是我正在尝试找到一种方法来加快此速度，或者找到一种更快的全新方法。

人们在评论中说这是3SUM问题的一种变体，根据Wikipedia页面，它可以在O（n ^ 2）中完成，其中n是数字数（即我的比率数）。我还没有找到一种在一般和Python中实现此方法的方法。

任何加速都将有所帮助；它并不一定只是一个更好的算法（库等）。我认为目前这比O（n ^ 3）好吗？

另外，对于MAX_NUM = 100,000，还算不错（大约4分钟），但是对于500,000来说，这很不好（尚未停止运行）。

最终，我想做MAX_NUM = 1,000,000或更多。

编辑

我希望看到一个更快的算法，例如O（n ^ 2），或者速度大大提高。

Answer 1

比您的速度快数百倍，而且没有浮点问题。
比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）²/ m。令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。我改用（分子，分母）元组。

解决方案代码

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，因此寻找较大的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