Python - 完美数搜索的优化

Posted 2023-03-05

【中文标题】Python - 完美数搜索的优化

【英文标题】：Python - Optimisation of Perfect Number search

【发布时间】：2017-03-30 12:47:41

【问题描述】：

p = []
for x in range(1, 50000000):
    count = 0
    for y in range(1, x // 2 + 1):
        if (x % y == 0):
            count += y
    if (count == x):
        p.append(x)

这是我的代码，用于尝试查找源自 1 到 50000000 之间的所有完美数字。它对前 3 个数字工作正常，它们在 1 到 10000 之间。但是随着它的进展，它变得非常缓慢。就像可能每 10 秒通过 1000 个数字。然后最终每 5 秒通过 10 个数字。

现在有什么办法可以让这更快吗？我尝试包含一些较小的东西，例如将 x 减 2 以确保我们不会超过一半（不会是 x 的倍数）

【问题讨论】：

我认为你不应该用蛮力的方式来解决这个数学问题，我建议你阅读更多关于梅森素数的内容，然后你就会明白你的方法在全部:)

【参考方案1】：

你可以做得更好。考虑以下几点：

1) 考虑 36 的因式分解，例如：1x36、2x18、3x12、4x9、6x6。就是这样。走得更远不会增加任何新东西。下一个分解将是 9x4，但我们已经知道 (4x9)。这个优势越来越大（比较你必须检查的最后一个数字的根）

2) 没有奇完美数。这实际上是一个猜想（尚未证明），但他们尝试了低于 10^300 的所有内容，但没有找到。所以肯定没有

from math import ceil
p = []
for x in range(2, 50000000, 2):
    divisors = 1
    for y in range(2, ceil(x**0.5) + 1):
        if x % y == 0:
            divisors.update(y, (x//y))
    if sum(divisors) == x:
        print('-', x)
        p.append(x)
#- 6
#- 28
#- 496
#- 8128

这应该会快很多，但肯定可以做更多的事情。

【讨论】：

如果要跳过奇数，不妨利用 Euclid-Euler 定理，只计算 2^(k-1)·(2^k - 1)。

我会吃，但我正在吃午饭！不过还是谢谢你提醒我！！

确实，正如@spectras 所暗示的那样，距离这个数学问题的正确解决方案还很遥远

应该是 "x // 2 + 1" 而不是 "c​​eil(x**0.5)):"

@Reece 你不能是认真的。

【参考方案2】：

这是一个使用 Mersenne Primes 的一些预先计算值的解决方案：

mersenne_prime_powers = [2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127]

def perfect_numbers(n=10):
    return [(2**p - 1) * (2**(p - 1)) for p in mersenne_prime_powers[:n]]

print(perfect_numbers(n=5))

输出：

[6, 28, 496, 8128, 33550336]

【讨论】：

这里有一个更简单但也很有意义的方法：return [6, 28, 496, 8128, 33550336]。

@StefanPochmann 从那条评论中可以很清楚地看到“你根本不明白我的答案......”，一个线索，分析“梅森素数生成”的复杂性，然后回顾 OP 的原始算法:) 。另外，我建议你看看“今天的号码”here，差不多了；）

我想是你不明白。我的观点是，如果您要查找并使用您自己无法想到的东西，那么您不妨查找您所追求的东西。我认为你的方式没有优势。

好吧，公平地说，这看起来像是一项家庭作业，意义不是来自结果，而是来自思考过程。我什至会争辩说，在学习处理这类问题时，实施这种幼稚的解决方案并观察它的速度有多快，这本身就是一种体验。

@PM2Ring 确实如此 :) 。这适用于所有数学水平，事实上，通常当我开始解决“中等/困难”欧拉项目时，我倾向于从天真的显而易见的蛮力方法开始并从那里优化，直到我意识到我不会再进一步聪明的优化，当我达到那个意识点时，我就会开始更多地学习必要的数学知识，以“以正确的方式”解决问题。无论如何，我的回答只是试图激励 OP 更多地了解“梅森素数”，我根本不打算优化它的循环，我认为这会误导学习路径

【参考方案3】：

如前所述，没有找到奇完美数。根据Wikipedia article on perfect numbers，如果存在任何奇完美数，它们必须大于 10¹⁵⁰⁰。显然，寻找奇完美数需要复杂的方法和/或很多时间。 :)

正如Wikipedia 中所述，欧几里得证明了即使是完美数也可以从梅森素数中产生，而欧拉证明了相反，所有偶数都具有这种形式。

因此，我们可以通过生成梅森素数来生成一个偶数列表。我们可以（相对）通过Lucas-Lehmer test 快速测试一个数字是否是梅森素数。

这是一个可以做到这一点的简短程序。这里使用的primes 函数是由 Robert William Hanks 编写的，其余代码是我几分钟前刚刚编写的。 :)

''' Mersenne primes and perfect numbers '''

def primes(n):
    """ Return a list of primes < n """
    # From http://***.com/a/3035188/4014959
    sieve = [True] * (n//2)
    for i in range(3, int(n**0.5) + 1, 2):
        if sieve[i//2]:
            sieve[i*i//2::i] = [False] * ((n - i*i - 1) // (2*i) + 1)
    return [2] + [2*i + 1 for i in range(1, n//2) if sieve[i]]

def lucas_lehmer(p):
    ''' The Lucas-Lehmer primality test for Mersenne primes.
        See https://en.wikipedia.org/wiki/Mersenne_prime#Searching_for_Mersenne_primes
    '''
    m = (1 << p) - 1
    s = 4
    for i in range(p - 2):
        s = (s * s - 2) % m
    return s == 0 and m or 0

#Lucas-Lehmer doesn't work on 2 because it's even, so we just print it manually :)
print('1 2\n3\n6\n')
count = 1
for p in primes(1280):
    m = lucas_lehmer(p)
    if m:
        count += 1
        n = m << (p - 1)
        print(count, p)
        print(m)
        print(n, '\n')

输出

1 2
3
6

2 3
7
28 

3 5
31
496 

4 7
127
8128 

5 13
8191
33550336 

6 17
131071
8589869056 

7 19
524287
137438691328 

8 31
2147483647
2305843008139952128 

9 61
2305843009213693951
2658455991569831744654692615953842176 

10 89
618970019642690137449562111
191561942608236107294793378084303638130997321548169216 

11 107
162259276829213363391578010288127
13164036458569648337239753460458722910223472318386943117783728128 

12 127
170141183460469231731687303715884105727
14474011154664524427946373126085988481573677491474835889066354349131199152128 

13 521
6864797660130609714981900799081393217269435300143305409394463459185543183397656052122559640661454554977296311391480858037121987999716643812574028291115057151
23562723457267347065789548996709904988477547858392600710143027597506337283178622239730365539602600561360255566462503270175052892578043215543382498428777152427010394496918664028644534128033831439790236838624033171435922356643219703101720713163527487298747400647801939587165936401087419375649057918549492160555646976 

14 607
531137992816767098689588206552468627329593117727031923199444138200403559860852242739162502265229285668889329486246501015346579337652707239409519978766587351943831270835393219031728127
141053783706712069063207958086063189881486743514715667838838675999954867742652380114104193329037690251561950568709829327164087724366370087116731268159313652487450652439805877296207297446723295166658228846926807786652870188920867879451478364569313922060370695064736073572378695176473055266826253284886383715072974324463835300053138429460296575143368065570759537328128 

15 1279
10407932194664399081925240327364085538615262247266704805319112350403608059673360298012239441732324184842421613954281007791383566248323464908139906605677320762924129509389220345773183349661583550472959420547689811211693677147548478866962501384438260291732348885311160828538416585028255604666224831890918801847068222203140521026698435488732958028878050869736186900714720710555703168729087
54162526284365847412654465374391316140856490539031695784603920818387206994158534859198999921056719921919057390080263646159280013827605439746262788903057303445505827028395139475207769044924431494861729435113126280837904930462740681717960465867348720992572190569465545299629919823431031092624244463547789635441481391719816441605586788092147886677321398756661624714551726964302217554281784254817319611951659855553573937788923405146222324506715979193757372820860878214322052227584537552897476256179395176624426314480313446935085203657584798247536021172880403783048602873621259313789994900336673941503747224966984028240806042108690077670395259231894666273615212775603535764707952250173858305171028603021234896647851363949928904973292145107505979911456221519899345764984291328 

该输出是在 2GHz 32 位机器上在 4.5 秒内产生的。您可以轻松生成更大的梅森素数和完美数，但不要指望它会很快。