改进解决欧拉问题357的Python脚本

Question

我知道这个问题听起来有点奇怪和烦人，但我希望圣诞的气氛能让工作变得更轻松。简单来说，我尝试写一个Python脚本来解决欧拉问题357，内容是：“找出所有不超过100000000的正整数n，使得对于每个n的因子d，d+n/d是质数。”所以我写了这个Python脚本，并多次尝试改进它，甚至在网上搜索过。最终，我让它快了十倍以上，但即使运行了24小时，它还是无法处理1亿个数字（没错，我试过）。所以这是我的脚本：我想知道有没有什么问题，或者是什么让它运行得这么慢。请忽略注释，因为其中一些是我自己语言写的。

#Find the sum of all positive integers n not exceeding 100000000
#such that for every divisor d of n, d+n/d is prime.

import time
maximum=100000000

def isprime(n):
    # range starts with 2 and only needs to go up the squareroot of n
    for x in range(2,int(n**0.5)+1):
        if n%x==0:
            return False
    return True

#main()
start=time.time()
final=list()
for n in range(1,maximum):
    counter,d=1,1
    divisors=list()
    while d<=(n**0.5): #Mettendo d<=n, e poi sotto range(len(divisori)/2) ci si mette il DECUPLO del tempo
        if n%d==0:
            divisors.append(d)
        d+=1
    for divisors_index in range(len(divisors)):
        prime=divisors[divisors_index]+n/divisors[divisors_index]
        if isprime(prime)==True:
            counter+=1
        elif isprime(prime)==False:
            counter=0
            break
    if counter==(len(divisors))+1:
        #print "%d:%s--->%d"%(n,divisors,len(divisors))
        final.append(n)   
end=time.time()-start

print "Results: %d. Time: %f seconds"%(len(final),end)

补充一下：我在网上找到一个用Haskell写的解决方案，然后转成了我稍微懂一点的C语言。事实上，这个解决方案看起来和我的非常相似；不过我没有仔细看，因为我不想提前知道答案（其实现在我想知道了）。谢谢你，祝你圣诞快乐。

Answer 1

#!/usr/bin/env python                                                       

def primes(ubound):
    ubound = ubound + 1
    a = [False] * ubound
    p = 2
    primes = []
    while p < ubound:
        primes.append(p)
        for n in range(p, ubound, p):
            a[n] = True
        p = p + 1
        while p < ubound and a[p]:
            p = p + 1
    return primes

print primes(100000000)


$ time erastothenes.py > results.txt
real    1m32.441s
user    1m14.866s
sys     0m4.588s

这个程序还是有点慢，可能还有进一步改进的空间；我对Python不是很精通。也许可以用原生数组来提高速度。但我不知道怎么在不先创建列表的情况下直接创建数组，这样做感觉很浪费时间。

(是的，我知道我在浪费前两个数组的位置。:p )

编辑 根据@DanielFisher提到的算法优化（并且不浪费空间）：

#!/usr/bin/env python                                                       

def primes(ubound):
    size = (ubound - 3) / 2
    a = [False] * size
    s = 0
    primes = []
    while s < size:
        t = 2 * s
        p = t + 3
        primes.append(p)
        print p
        for n in range(t * (s + 3) + 3, size, p):
            a[n] = True
        s = s + 1
        while s < size and a[s]:
            s = s + 1
    return primes

primes(100000000)

$ time erastothenes.py > results.txt
real    0m36.737s
user    0m35.725s
sys     0m0.864s

Answer 2

我觉得你在你的 isprime 函数里花了很多时间。也许可以用一些额外的内存，提前用一个质数筛选的方法，把每个数字是否是质数的信息存储在一个列表里，这样会稍微有点帮助。

你甚至可以进一步优化，使用一种和筛选方法很相似的算法，同时生成质数和它们的因子列表。

Answer 3

不要先列出所有的因数再去测试它们。你可以在找到因数的时候就进行测试，一旦发现有不符合条件的，就立刻跳过这个数字。

换个说法就是：“找出所有不超过100000000的整数n，使得没有因数d满足n/d + d是合成数。”我们可以用Python更简单地表达这个问题，首先创建一个函数来测试一个因数是否会导致合成数的出现，然后利用Python内置的函数式编程工具。这样做通常比直接遍历列表要高效一些，而且还会自动实现“提前退出”的逻辑（也就是说，一旦我们找到一个会导致合成数的因数，就立刻停止对这个数字的测试）。

import time
maximum = 100000000

def isprime(n):
    # We can use those same tools in the prime tester...
    return not any(n % x == 0 for x in range(2, int(n**0.5)+1))

def divisor_with_composite_sum(n, d):
    return n % d == 0 and not isprime(d + n / d)

start = time.time()
result = sum(
    n for n in range(1, maximum)
    if not any(
        divisor_with_composite_sum(n, d)
        for d in range(1, int(n**0.5) + 1)
    )
)

现在只需要优化一下我们找到质数和因数的方法就可以了 :)

---

补充一下：经过进一步思考。这是一个欧拉问题；你应该有比蛮力更聪明的解决办法。

每个符合条件的数字在它的质因数分解中，每个质数的指数都是1。用反证法来证明：假设某个质数p，n能被p^2（或者更高的p的幂）整除。那么p显然能被p整除，而n/p也能被p整除（根据假设），所以它们的和也能被p整除，因此是合成数。此外，所有符合条件的数字都是质数减1，因为1总是一个因数。

所以：列出所有不超过100000000的质数（可以参考Amadan的回答，看看你能否应用上述技巧 ;)）。对于每个质数，减去1，然后检查结果是否有重复的质因数。如果没有，那么你就可以进行更详细的测试了。