加快欧拉计划问题75的解决方案

from functools import reduce
import math

primes=set([2,3,5,7,11,13,17,19,23])
def isPrime(n):
    n=abs(n)
    if n in primes:
       return True
    if n<2:
       return False
    if n==2:
       return True
    if n%2==0:
       return False
    for x in range(3, int(n**0.5)+1, 2):
        if n % x == 0:
            return False
    primes.add(n)
    return True

def aFactors(n):
    if isPrime(n):
        return [1,n]
    return set(reduce(list.__add__,([i,n//i] for i in range(1,int(math.sqrt(n))+1) if n%i==0)))


count=0
number=12
while number<=1500000:
    p=number/2
    f=aFactors(p)
    triangles=[]
    pairs=[(i, int(p/i)) for i in f]
    add=triangles.append
    for i in pairs:
        mList=aFactors(i[0])
        pairsOfmc=[(k,int(i[0]/k)) for k in mList]
        for j in pairsOfmc:
            add((2*i[1]*i[0]-i[1]*i[1]*j[0],2*i[0]*i[1]-2*i[0]*j[1],2*i[0]*j[1]+i[1]*i[1]*j[0]-2*i[0]*i[1]))
    r=0
    while r<len(triangles):
        if any(triangles[r][i]<=0 for i in range(len(triangles[r]))):
            del triangles[r]
        else:
            l=list(triangles[r])
            l.sort()
            triangles[r]=tuple(l)
            r+=1
    trianglesFinal=list(set(triangles))
    for i in trianglesFinal:
        print(number, i)
    if len(trianglesFinal)==1:
        count+=1
    number+=2
print(count)

#1 楼

作为记录，这是欧拉计划问题75。

1。告诫

您写道：


请注意，我并不是在寻找另一种计算方法（我相信有一种计算方法，但是对我来说，Project欧拉（Euler）打算寻找自己的方法，对我来说，使用您的方法会作弊。）


放弃这种态度！编程的关键部分（与编码同样重要）正在开发各种算法和技术，可将其应用于所面临的问题。尽管为了您自己的知识满意度，最好是自己发现算法，但是您应该始终将自己的解决方案与其他人发现的最佳解决方案进行比较，以便您可以改善下一个问题的库。 />特别是，当您想加速程序时，说“我不想实现不同的算法，我只是想加快编写的代码”会适得其反。最大的加速来自找到更好的算法。

您说您的程序“运行良好”，但这似乎并不正确。欧拉计画说：


每个问题都是按照“一分钟规则”设计的，这意味着尽管设计成功的算法可能会花费几个小时，但遇到的困难却更多，一个有效的实现将允许在不到一分钟的时间内在功率适中的计算机上获得解决方案。 （并且我的笔记本电脑运行得很热），所以我将其杀死。

因此，在下面的第2部分中，我将使用较小的测试大小，即100,000，而不是1,500,000，以使运行时可管理。 />
2。分段优化



您的程序很难在交互式解释器中进行测试，因为它具有顶级代码。最好将主程序放在函数中，以便您可以从解释器中调用它，然后添加一个if __name__ == '__main__':部分，以便在需要时可以将其作为脚本运行。

函数名称problem75，它接受一个参数，它是问题中导线长度的最大值。现在：

>>> from timeit import timeit
>>> timeit(lambda:problem75(10**5),number=1)
12 (3, 4, 5)
24 (6, 8, 10)
30 (5, 12, 13)
[... much output deleted ...]
33.87288845295552

不需要在主循环末尾使用print语句。这样可以节省大约3秒钟，使时间减少到31.0秒。

您的函数aFactors对isPrime进行了初始调用，采用\ $ O（\ sqrt n）\ $，以避免循环这也需要\ $ O（\ sqrt n）\ $。测试成本与节省的成本一样多，因此不值得。删除对isPrime的调用，然后像这样重写函数：并缓存了result.add的值，因此无需在每次循环时都进行查找。这样可以节省大约4秒钟；现在的时间为27.7秒。


由于您要处理的一组因素中唯一要做的就是将其变成成对然后进行迭代，所以为什么不生成这样的对，如下所示： br />

def factors(n):
    """Return the set of factors of n."""
    result = set([1, n])
    add = result.add
    for i in range(2, int(math.sqrt(n)) + 1):
        if n % i == 0:
            add(i)
            add(n // i)
    return result

，然后像这样处理它们：

def product_pairs(n):
    """Generate pairs (i, j) such that i * j = n."""
    yield 1, n
    yield n, 1
    for i in range(2, int(math.sqrt(n)) + 1):
        if n % i == 0:
            j = n // i
            yield i, j
            yield j, i

这避免了在内存中构造中间列表。这样可以节省大约3秒钟；现在的时间为24.3秒。

而不是编写for i in product_pairs(...):然后查找i[0]和i[1]，而是在循环中将其分配为两个变量时将其分成两个变量，如下所示：

triangles = []
append = triangles.append
for i in product_pairs(number // 2):
    for j in product_pairs(i[0]):
        append((2*i[1]*i[0]-i[1]*i[1]*j[0],2*i[0]*i[1]-2*i[0]*j[1],2*i[0]*j[1]+i[1]*i[1]*j[0]-2*i[0]*i[1]))

这样可以避免序列查找。这样可以节省大约3秒钟；现在时间21.2秒。


您构造一个三角形列表，其中一些三角形的边长为负或零。然后，您遍历列表，并检查无效的三角形。但是，列表上的del可能是一项昂贵的操作：它必须复制列表的其余部分以保持其连续性。有关Python内置数据结构上基本操作的成本，请参见Python Wiki上的TimeComplexity页面，在该页面上您可以看到列表上的“删除项目”操作花费\ $ O（n）\ $。

所以不要首先将这些无效的三角形添加到列表中。实际上，根本不用理会列表，只需直接构造三角形集即可：

现在时间15.1秒。


通过命名del，i，j和k的各个乘积可以避免某些重复的乘法： >
这节省了大约1秒钟；时间现在为13.9秒。


我恐怕还不及对代码进行分段加速。总共提高了约60％。修改后的程序现在可以在我的笔记本电脑上大约八分钟内获得解决欧拉计画问题的完整方法：

要大幅提高速度，您确实需要...

3。更好的算法

（如果需要的话，请睁开眼睛。）

基本见解是，我们可以遍历周界并找到具有该周界的勾股三元组，而不必以更方便的顺序遍历所有勾股三叉戟，并计算在每个周长处发现多少个三角形。

我们可以通过仅生成原始毕达哥拉斯三元组，然后将原始三元组乘以连续的数字1、2、3，...来生成所需范围内的其余毕达哥拉斯三元组，从而进一步加快处理速度。实际上，我们只需要乘以它们的周长即可。

Euclid的公式可用于生成所有原始勾股三元组。给定共质数正整数\ $ m \ $和\ $ n \ $，其中\ $ m> n \ $，并且恰好是\ $ m \ $和\ $ n \ $之一，$$ \ eqalign {a＆= m ^ 2 − n ^ 2 \\ b＆= 2mn \\ c＆= m ^ 2 + n ^ 2} $$是一个原始的勾股三元组。不会按周长的长度顺序生成三元组，因此需要某种方法来确定何时停止。我在以下代码中的策略是迭代\ $ m \ $的升序值。三元组的周长为\ $ a + b + c = 2m ^ 2 + 2mn \ $，至少为\ $ 2m（m +1）\ $，因为\ $ n≥1 \ $。因此\ $ 2m（m +1）\ $是迄今为止所有三元组都已生成的下界。当超出限制时，我们可以停止搜索：在所需范围内找不到更多的三元组。

for i, j in product_pairs(number // 2):
    for k, l in product_pairs(i):
        append((2*j*i - j*j*k, 2*i*j - 2*i*l, 2*i*l + j*j*k - 2*i*j))

更合理的时间：

triangles = set()
add = triangles.add
for i, j in product_pairs(number // 2):
    for k, l in product_pairs(i):
        triangle = 2*j*i - j*j*k, 2*i*j - 2*i*l, 2*i*l + j*j*k - 2*i*j
        if all(side > 0 for side in triangle):
            add(tuple(sorted(triangle)))

#2 楼

不必构造所有三角形来确定是否恰好有一个三角形：找到两个三角形后就可以停下来。

有些小事情可以清理代码，从而可以稍微提高速度也可以：

而不是索引(triangles[r][i]<=0 for i in range(len(triangles[r])))，而是直接迭代：可以避免以下问题：

(x <= 0 for x in triangles[r])

可以代替以下reduce(list.__add__, ...外部循环使用：

(j for i in range(1,int(math.sqrt(n))+1) if n%i==0 for j in (i,n//i))

（这是针对Python 3的，应在Python 2中使用while）