How to calculate perfect numbers quickly?

Question

I am trying to perform an exercise to show the perfect numbers present within a certain range, but I can only accomplish such a feat until the perfect fourth number. If I raise the range, it takes too long and doesn’t respond:

def perfeito(n_max):
    lista = []
    for n in range(1, n_max + 1):
        soma = 0
        for div in range(1, n):
            if n % div == 0:
                soma += div
        if n == soma:
            lista.append(n)
    return lista

Here we define the range, in this case 500, and it will return the first 3:

print(perfeito(500)) 

Set perfect numbers: {6, 28, 496, 8128, 33550336, 8589869056, ...}

After the 4th number. I can no longer display. Any optimization tips for my code?

Answer 1

There is a direct relationship between the perfect numbers and the prime numbers of Mersenne. A prime number of Mersenne is nothing more than a prime number that can be written in the form M_n = 2ⁿ - 1, for data n integer, and its relation to perfect numbers is a power of 2. It is worth emphasizing in the answer that the concepts applied here apply to, only, even perfect numbers, however the solution remains valid given the fact that it is not yet known by mathematics no odd perfect number - the day they discover I edit the solution, I promise.

Thus, given a prime number of Mersenne M_n = 2ⁿ - 1, we can get the perfect number by doing P_n = (2^n-1) (2ⁿ - 1).

• P_n=1 = (2^1-1) (2¹ - 1) = 1
• P_n=2 = (2^2-1) (2² - 1) = 6
• P_n=3 = (2^3-1) (2³ - 1) = 28
• P_n=4 = (2^4-1) (2⁴ - 1) = 120
• P_n=5 = (2^5-1) (2⁵ - 1) = 496
• ...

Note that when 2ⁿ – 1 is not prime, in cases of n=1 and n=4, the result will not be a perfect number. So the first challenge of this approach is to make sure we have a 2ⁿ - 1 a prime number.

And this challenge is very easy to solve. Mathematically, we can define that 2ⁿ - 1 is a prime number and, ready, we solve our problem. However, solving one, another appears: but 2ⁿ - 1 is a cousin of Mersenne?

To ensure that the prime number P_n = 2ⁿ - 1 is a Mersenne prime there must be an integer value n that validates the expression. That is, if there is a value n such that n = log₂(P_n+1) be whole, then P_n will be a cousin of Mersenne.

The approach will then be to scroll through the list of prime numbers, check whether it is a Mersenne prime and, when it is, calculate its perfect number. With luck, we already know how to calculate prime numbers in an efficient way

• How to generate 200,000 primes as fast as possible in Python?

Then just go through these values and check which are Mersenne primes by calculating the perfect number:

# Percorre a lista de números primos menores que N:
for prime in sieve_of_eratosthene(N):

    # Calcula o valor de n que define o Pn:
    n = math.log2(prime+1)

    # Verifica se n é inteiro, sendo um primo de Mersenne:
    is_mersenne = n.is_integer()

    # Se for um primo de Mersenne, calcula o número perfeito:
    if is_mersenne:
        print(2**(n-1) * prime)

^{See working on Repl.it | Ideone | Github GIST}

So, testing here, it took about 34 seconds to calculate to the perfect number 137438691328. Above that I started to have problems with memory, that I will try to solve soon.

Additional references

Two videos I recommend watching on the subject are:

• Perfect Numbers and Mersenne Primes - Numberphile
• Perfect Number Proof - Numberphile

---

Seeking to optimize the code, I was able to do without using Sieve. Basically what the code does is to define a generator that returns all Mersenne prime numbers by checking if the number is prime and using the Lucas-Lehmer primality test.

def is_prime(number):
    # Condição para number=2 ignorada propositalmente, visto que a menor entrada será 3
    i = 3
    while i**2 <= number:
        if number % i == 0:
            return False
        i += 2
    return True

def lucas_lehmer(p):
    s = 4
    M = 2**p - 1

    for _ in range(p - 2):
        s = ((s * s) - 2) % M
    return s == 0

def mersenne_primes():
    p = 3
    while True:
        if is_prime(p) and lucas_lehmer(p):
            yield (p, 2**p - 1)
        p += 2

To use it, for example, searching for the first 10 perfect numbers, just do:

numbers = mersenne_primes()

for _ in range(10):
    p, mersenne = next(numbers)
    perfect = 2**(p-1) * mersenne
    print(perfect)

^{See working on Repl.it | Ideone | Github GIST}

Some results I got:

• To the top ten perfect numbers: 0.0003993511199951172s
• To the 15 first perfect numbers: 1.9178969860076904s
• To the Top 16 perfect numbers: 6.957615613937378s
• To the Top 17 perfect numbers: 28.416791677474976s
• To the 18 first perfect numbers: 78.89492082595825s
• To the 19 first perfect numbers: 93.7487268447876s

For 20 I was too lazy to wait.

Answer 2

Your program is just taking a really long time to calculate.

Notice that you have one for inside each other. The outside loop will rotate more than 33 million times to find the 5th number. The inside loop will traverse in each iteration the same value as the external loop number minus 1.

If the outer loop rotates n_max times, this means that the loop inside will end up running a number of times equal n_max * (n_max - 1). If n_max for 33550336, then the internal loop will have to rotate 1.125.625.012.162.560 times until the fifth number is found. To reach this number of more than a quadrillion iterations, it is possible that the program will take a few days.

Already to find the 8589869056, it would probably take years.

It is possible to make some simplifications:

• Cut the internal loop when soma overtake n.

• When you find that n % div == 0, you find two dividers, not just one. One of them is div and the other is n / div. With this, you can decrease the upper limit of the iteration and only need to go to the square root of n.

• Already do soma start with 1 and start testing from 2. It may seem silly, but without this the two previous optimizations does not work because when finding 1, he would also find the n and cut the ribbon immediately.

With this, it is possible to optimize the internal loop more or less like this:

soma = 1
for div in range(2, math.ceil(math.sqrt(n))):
    if n % div == 0:
        soma += div
        div2 = int(n / div)
        if (div != div2):
            soma += div2
        if soma > n:
            break

You can also add a parameter n_min and instead of running from 1 to n_max, traverse of n_min until n_max. With these changes, by searching the perfect numbers in the range 33550337 and 33550338, he finds 33550336 in the blink of an eye.

The whole code goes like this:

import math

def perfeito(n_min, n_max):
    lista = []
    for n in range(n_min, n_max + 1):
        soma = 1
        for div in range(2, math.ceil(math.sqrt(n))):
            if n % div == 0:
                soma += div
                div2 = int(n / div)
                if (div != div2):
                    soma += div2
                if soma > n:
                    break
        if n == soma:
            lista.append(n)
    return lista

print(perfeito(33550335, 33550337))

See here working on ideone.

I tried using print(perfeito(2, 33550337)). The numbers 6, 28, 496 and 8128 appeared in milliseconds (I put a print before the append to see this). However, it still takes several hours for the 33550336 to appear.

Answer 3

Adding to the already published answer, one point to note is that every perfect number is a hexagonal number. So instead of testing one by one of the numbers, you can scroll through only the hexagonal numbers. Ex:

def gen_hexagonal():
    n = 1
    while True:
        yield (2 * n * (2 * n - 1)) // 2
        n += 1

for i, num in gen_hexagonal():
    print(num)

# 1, 6, 15, 28, 45, 66, 91, 120, 153, 190, 231, 276, 325, 378, 435, ...

Working example