Python – Optimización de búsqueda de número perfecto

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) 

Este es mi código para intentar encontrar todos los números perfectos que se originan entre 1 y 50000000. Funciona bien para los primeros 3 números, están entre 1 y 10000. Pero a medida que avanza, se vuelve extremadamente lento. Como tal vez pasar por 1000 números cada 10 segundos. Luego, finalmente, pasando por 10 números cada 5 segundos.

Ahora hay de todos modos podría hacer esto más rápido? Intenté incluir algunas cosas más pequeñas, como bucear x por 2 para asegurarme de que no superamos la mitad (no vamos a ser un múltiplo de x)

Como ya se ha mencionado, no se han encontrado números perfectos impares. Y de acuerdo con el artículo de Wikipedia sobre números perfectos , si existen números impares perfectos, deben ser mayores que 10 1500 . Claramente, buscar números impares perfectos requiere métodos sofisticados y / o mucho tiempo. 🙂

Como se indica en Wikipedia , Euclid demostró que incluso los números perfectos pueden producirse a partir de números primos de Mersenne, y Euler demostró que, por el contrario, todos los números incluso perfectos tienen esa forma.

Entonces podemos producir una lista de números incluso perfectos generando números primos de Mersenne. Y podemos (relativamente) probar rápidamente si un número es primo de Mersenne mediante la prueba de Lucas-Lehmer .

Aquí hay un progtwig corto que hace precisamente eso. La función de primes utilizada aquí fue escrita por Robert William Hanks, el rest del código fue escrito por mí hace unos minutos. 🙂

 ''' 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') 

salida

 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 

Esa salida se produjo en 4,5 segundos en una máquina de 2GHz de 32 bits. Puede producir fácilmente números primos de Mersenne más grandes y números perfectos, pero no espere que sean rápidos.

Puedes hacerlo mucho mejor. Considera lo siguiente:

1) Piense en la factorización de 36, por ejemplo: 1×36, 2×18, 3×12, 4×9, 6×6. Y eso es. Ir más allá no añade nada nuevo. La siguiente factorización sería 9×4 pero ya lo sabemos (4×9). Esta ventaja se hace cada vez más grande (compare la raíz del último número que debe comparar con su mitad)

2) No hay números perfectos impares. Esto es en realidad una conjetura (aún no probada) pero probaron todo por debajo de 10 ^ 300 y no encontraron ninguna. Entonces definitivamente no hay ninguno <50000000. Eso significa que puede omitir la mitad de los términos.

 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 

Esto debería ser mucho más rápido, pero definitivamente hay más cosas que uno puede hacer.

Aquí hay una solución que utiliza algunos valores precomputados de 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)) 

Salida:

 [6, 28, 496, 8128, 33550336]