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

Solution 1:

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 inrange(2, 50000000, 2):
    divisors = {1}
    for y inrange(2, ceil(x**0.5) + 1):
        if x % y == 0:
            divisors.update({y, (x//y)})
    ifsum(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.

Solution 2:

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 '''defprimes(n):
    """ Return a list of primes < n """# From http://stackoverflow.com/a/3035188/4014959
    sieve = [True] * (n//2)
    for i inrange(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 + 1for i inrange(1, n//2) if sieve[i]]

deflucas_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 = 4for i inrange(p - 2):
        s = (s * s - 2) % m
    return s == 0and m or0#Lucas-Lehmer doesn't work on 2 because it's even, so we just print it manually :)print('1 2\n3\n6\n')
count = 1for 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

12362372835314964712781285138191335503366171310718589869056719524287137438691328831214748364723058430081399521289612305843009213693951265845599156983174465469261595384217610896189700196426901374495621111915619426082361072947933780843036381309973215481692161110716225927682921336339157801028812713164036458569648337239753460458722910223472318386943117783728128121271701411834604692317316873037158841057271447401115466452442794637312608598848157367749147483588906635434913119915212813521686479766013060971498190079908139321726943530014330540939446345918554318339765605212255964066145455497729631139148085803712198799971664381257402829111505715123562723457267347065789548996709904988477547858392600710143027597506337283178622239730365539602600561360255566462503270175052892578043215543382498428777152427010394496918664028644534128033831439790236838624033171435922356643219703101720713163527487298747400647801939587165936401087419375649057918549492160555646976146075311379928167670986895882065524686273295931177270319231994441382004035598608522427391625022652292856688893294862465010153465793376527072394095199787665873519438312708353932190317281271410537837067120690632079580860631898814867435147156678388386759999548677426523801141041933290376902515619505687098293271640877243663700871167312681593136524874506524398058772962072974467232951666582288469268077866528701889208678794514783645693139220603706950647360735723786951764730552668262532848863837150729743244638353000531384294602965751433680655707595373281281512791040793219466439908192524032736408553861526224726670480531911235040360805967336029801223944173232418484242161395428100779138356624832346490813990660567732076292412950938922034577318334966158355047295942054768981121169367714754847886696250138443826029173234888531116082853841658502825560466622483189091880184706822220314052102669843548873295802887805086973618690071472071055570316872908754162526284365847412654465374391316140856490539031695784603920818387206994158534859198999921056719921919057390080263646159280013827605439746262788903057303445505827028395139475207769044924431494861729435113126280837904930462740681717960465867348720992572190569465545299629919823431031092624244463547789635441481391719816441605586788092147886677321398756661624714551726964302217554281784254817319611951659855553573937788923405146222324506715979193757372820860878214322052227584537552897476256179395176624426314480313446935085203657584798247536021172880403783048602873621259313789994900336673941503747224966984028240806042108690077670395259231894666273615212775603535764707952250173858305171028603021234896647851363949928904973292145107505979911456221519899345764984291328

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.

Solution 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]

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

print(perfect_numbers(n=5))

Output:

[6, 28, 496, 8128, 33550336]