How to calculate perfect numbers quickly?

10

I'm trying to perform an exercise to show the perfect numbers present within a given range, but I can only do this to the perfect fourth number. If I increase the range, it takes a long time and does not 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?

    
asked by anonymous 24.09.2018 / 19:56

3 answers

7

There is a direct relationship between perfect numbers and the Mersenne prime numbers . A prime number of Mersenne is nothing more than a prime number that can be written in the form M n = 2 n - 1, for given n integer, and the relation of the perfect numbers is a power of 2. It is worth emphasizing in the answer that the concepts applied here are valid for only perfect even numbers, but the solution holds true given the fact that no perfect perfect number is known yet by mathematics - on the day they find out I edit the solution, I promise.

Thus, given a prime number of Mersenne M n = 2 n -1, we can obtain the perfect number by making P n = (2 n-1 ) (2 n -1).

  • P n = 1 = (2 1-1 ) (2 1 - 1) = 1
  • (2 2 - 1) = 6 (2 3 - 1) = 28 (2 4 - 1) = 120 (2 5 - 1) = 496
  • ...

Notice that when 2 n - 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 ensure that we have a prime number as 2 n - 1.

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

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

The approach will then be to go through the list of prime numbers, check if it is a Mersenne cousin, and when it does, calculate its perfect number. With luck, we already know how to calculate primes efficiently

So just go through those values and see which are Mersenne's cousins, 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 at Repl.it | Ideone | GitHub GIST

So by testing here, it took about 34 seconds to figure out the perfect number 137438691328. Above this I started having memory problems, which I will try to solve soon.

Additional references

Two videos I recommend watching on the subject are:

Looking to optimize the code, I was able to do without using Sieve. Basically what the code does is define a generator that returns all the prime numbers of Mersenne 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 at Repl.it | Ideone | GitHub GIST

Some results I got:

  • For the top 10 perfect numbers: 0.0003993511199951172s
  • For 15 first perfect numbers: 1.9178969860076904s
  • For the 16 first perfect numbers: 6.957615613937378s
  • For 17 first perfect numbers: 28.416791677474976s
  • For 18 first perfect numbers: 78.89492082595825s
  • For 19 first perfect numbers: 93.7487268447876s

For 20 I was too lazy to wait.

    
24.09.2018 / 20:42
6

Your program is only taking a very long time to compute.

Note that you have for inside the other. The outer loop will rotate more than 33 million times to find the 5th number. The inner loop will run at each iteration the same value as the outer loop number minus 1.

If the outer loop runs n_max times, this means that the loop inside will end up executing a number of times equal n_max * (n_max - 1) . If n_max is 33550336, then the inner 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 one quadrillion iterations, it may take a few days.

Already to find 8589869056, it would probably take years.

Some simplifications are possible:

  • Cut the inner loop when soma exceeds n .

  • When you find that n % div == 0 , you find two divisors, 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 just go to the square root of n .

  • Now having soma start with 1 and start testing from 2. It may seem silly, but without that the two previous optimizations will not work because when you find 1, it would also find n e would cut the loop immediately.

With this, you can 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 going from 1 to n_max , go from n_min to n_max . With these changes, by searching for the perfect numbers in the range 33550337 and 33550338, he finds 33550336 in the blink of an eye.

The complete code looks 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)) . Numbers 6, 28, 496, and 8128 have appeared in milliseconds (I put print before append to see this). However, it still takes several hours for the 33550336 to appear.

    
24.09.2018 / 20:07
3

Adding to the answers already posted, one point to note is that any perfect number is a hexagonal number hexagonal number < a>. 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, ...

Example working

    
24.09.2018 / 20:57