Trabajando con primos grandes en Python

¿Cuál es una forma eficiente de trabajar con números primos grandes con Python? Busca aquí o en google, y encuentra muchos métodos diferentes para hacerlo … tamices, algoritmos de prueba de primalidad … ¿En qué formas funcionan los primos más grandes?

Para determinar si un número es primo, hay pruebas de tamices y primalidad.

 # for large numbers, xrange will throw an error. # OverflowError: Python int too large to convert to C long # to get over this: def mrange(start, stop, step): while start < stop: yield start start += step # benchmarked on an old single-core system with 2GB RAM. from math import sqrt def is_prime(num): if num == 2: return True if (num < 2) or (num % 2 == 0): return False return all(num % i for i in mrange(3, int(sqrt(num)) + 1, 2)) # benchmark is_prime(100**10-1) using mrange # 10000 calls, 53191 per second. # 60006 function calls in 0.190 seconds. 

Esto parece ser el más rápido. Hay otra versión usando not any que veas,

 def is_prime(num) # ... return not any(num % i == 0 for i in mrange(3, int(sqrt(num)) + 1, 2)) 

Sin embargo, en los puntos de referencia obtuve 70006 function calls in 0.272 seconds. sobre el uso de all 60006 function calls in 0.190 seconds. durante la prueba si 100**10-1 era primo.

Si necesita encontrar el próximo primo más alto, este método no funcionará para usted. Necesitas ir con una prueba de primalidad, he encontrado que el algoritmo de Miller-Rabin es una buena opción. Es un poco más lento que el método Fermat , pero más preciso contra las pseudoprimes. El uso del método mencionado lleva más de 5 minutos en este sistema.

Algoritmo de Miller-Rabin :

 from random import randrange def is_prime(n, k=10): if n == 2: return True if not n & 1: return False def check(a, s, d, n): x = pow(a, d, n) if x == 1: return True for i in xrange(s - 1): if x == n - 1: return True x = pow(x, 2, n) return x == n - 1 s = 0 d = n - 1 while d % 2 == 0: d >>= 1 s += 1 for i in xrange(k): a = randrange(2, n - 1) if not check(a, s, d, n): return False return True 

Algoritmo Fermat :

 def is_prime(num): if num == 2: return True if not num & 1: return False return pow(2, num-1, num) == 1 

Para obtener la siguiente prima más alta:

 def next_prime(num): if (not num & 1) and (num != 2): num += 1 if is_prime(num): num += 2 while True: if is_prime(num): break num += 2 return num print next_prime(100**10-1) # returns `100000000000000000039` # benchmark next_prime(100**10-1) using Miller-Rabin algorithm. 1000 calls, 337 per second. 258669 function calls in 2.971 seconds 

Al utilizar la prueba de Fermat , obtuvimos una referencia de 45006 function calls in 0.885 seconds. , pero tienes una mayor probabilidad de pseudoprimes.

Por lo tanto, si solo necesita comprobar si un número es primo o no, el primer método para is_prime funciona bien. Es el más rápido, si usas el método mrange con él.

Idealmente, querría almacenar los números primos generados por next_prime y simplemente leerlos.

Por ejemplo, usando next_prime con el algoritmo de Miller-Rabin :

 print next_prime(10^301) # prints in 2.9s on the old single-core system, opposed to fermat's 2.8s 1000000000000000000000000000000000000000000000000000000000000000000000000000 0000000000000000000000000000000000000000000000000000000000000000000000000000 0000000000000000000000000000000000000000000000000000000000000000000000000000 00000000000000000000000000000000000000000000000000000000000000000000000531 

No podría hacer esto con return all(num % i for i in mrange(3, int(sqrt(num)) + 1, 2)) de manera oportuna. Ni siquiera puedo hacerlo en este viejo sistema.

Y para estar seguros de que next_prime(10^301) y Miller-Rabin arrojan un valor correcto, esto también se probó utilizando los algoritmos de Fermat y Solovay-Strassen .

Consulte: fermat.py , miller_rabin.py y solovay_strassen.py en gist.github.com

Editar: Se corrigió un error en next_prime

He escrito dos artículos sobre esto, junto con puntos de referencia, para ver qué métodos son más rápidos.

Los números primos con Python fueron escritos antes del conocimiento de la prueba de primalidad Baillie-PSW .

Los números primos con Python v2 se escribieron después, comparando las Lucas pseudoprimes y la prueba de primalidad de Baillie-PSW .

En respuesta a la posible imprecisión de math.sqrt he math.sqrt dos métodos diferentes para realizar una isqrt(n) . isqrt_2(n) viene de este artículo y este código C

El método más común visto:

 def isqrt_1(n): x = n while True: y = (n // x + x) // 2 if x <= y: return x x = y cProfile.run('isqrt_2(10**308)') 

Resultados de referencia:

 isqrt_1 at 10000 iterations: 12.25 Can perform 816 calls per second. 10006 function calls in 12.904 seconds Ordered by: standard name ncalls tottime percall cumtime percall filename:lineno(function) 1 0.000 0.000 12.904 12.904 :1() 1 0.690 0.690 12.904 12.904 math.py:10(func) 10000 12.213 0.001 12.213 0.001 math.py:24(isqrt_1) 1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects} 1 0.000 0.000 0.000 0.000 {range} 2 0.000 0.000 0.000 0.000 {time.time} 

Este método es increíblemente lento. Así que probamos el siguiente método:

 def isqrt_2(n): if n < 0: raise ValueError('Square root is not defined for negative numbers.') x = int(n) if x == 0: return 0 a, b = divmod(x.bit_length(), 2) n = 2 ** (a + b) while True: y = (n + x // n) >> 1 if y >= n: return n n = y cProfile.run('isqrt_2(10**308)') 

Resultados de referencia:

 isqrt_2 at 10000 iterations: 0.391000032425 Can perform 25575 calls per second. 30006 function calls in 1.059 seconds Ordered by: standard name ncalls tottime percall cumtime percall filename:lineno(function) 1 0.000 0.000 1.059 1.059 :1() 1 0.687 0.687 1.059 1.059 math.py:10(func) 10000 0.348 0.000 0.372 0.000 math.py:34(isqrt_2) 10000 0.013 0.000 0.013 0.000 {divmod} 10000 0.011 0.000 0.011 0.000 {method 'bit_length' of 'long' objects} 1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects} 1 0.000 0.000 0.000 0.000 {range} 2 0.000 0.000 0.000 0.000 {time.time} 

Como puede ver, la diferencia entre isqrt_1(n) e isqrt_2(n) es un sorprendente 11.858999967575 seconds más rápido.

Puedes ver esto en acción en Ideone.com u obtener el código

nota: Ideone.com resultó en un tiempo de espera de ejecución para isqrt_1(n) por lo que el índice de referencia se redujo a 10**200