Python - Optimisation of Perfect Number search

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)

This is my code to try and find all the perfect numbers that originate between 1 and 50000000. It works fine for the first 3 numbers, they are between 1 and 10000. But as it progresses it becomes extremely slow. Like maybe going through 1000 numbers every 10 seconds. Then eventually going through 10 numbers every 5 seconds.

Now is there anyway I could make this faster? I tried including some smaller things, like diving x by 2 to make sure we don't go over half (not going to be a multiple of x)

回答1:

As has already been mentioned, no odd perfect numbers have been found. And according to the Wikipedia article on perfect numbers, if any odd perfect numbers exist they must be greater than 10¹⁵⁰⁰. Clearly, looking for odd perfect numbers takes sophisticated methods &/or a lot of time. :)

As stated on Wikipedia, Euclid proved that even perfect numbers can be produced from Mersenne primes, and Euler proved that, conversely, all even perfect numbers have that form.

So we can produce a list of even perfect numbers by generating Mersenne primes. And we can (relatively) quickly test if a number is a Mersenne prime via the Lucas-Lehmer test.

Here's a short program that does just that. The primes function used here was written by Robert William Hanks, the rest of the code was just written by me a few minutes ago. :)

''' Mersenne primes and perfect numbers '''

def primes(n):
    """ Return a list of primes < n """
    # From http://stackoverflow.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')

output

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

That output was produced in 4.5 seconds on a 2GHz 32 bit machine. You can easily produce larger Mersenne primes and perfect numbers, but don't expect it to be fast.

回答2:

You can do a lot better. Consider the following:

1) Think of the factorization of 36 for example: 1x36, 2x18, 3x12, 4x9, 6x6. And that's it. Going further doesn't add anything new. The next factorization would be 9x4 but we already know that (4x9). This advantage gets progressively larger (compare the root of the last number you have to check against its half)

2) There are no odd perfect numbers. This is actually a conjecture (not proven yet) but they tried everything below 10^300 and didn't find any. So there are definately exactly none < 50000000. That means you can skip half the terms.

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

This should be a lot faster but there are definately more things one can do.

回答3:

Here's a solution using some precomputed values of 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))

Output: