Lucas-Lehmer, corroborando primos de Mersenne

En el post anterior sobre los números primos, hablabamos de los números primos de Mersenne, y una de las cosas que surgieron es que con números medianamente grandes (ok, estamos hablando del órden de 2²⁹, o sea, grandes dependiendo para que), esto podía tardar una eternidad. Sabía que tenían que haber métodos para hacer esto más rápido. Buscando, me topé con el test de Lucas-Lehmer, que en teoría es más rápido. Ahora, veamos de que se trata.

El test

Supongamos la siguiente sucesión (ok, la imágen la saqué de Wikipedia, se hace lo que se puede):

Vemos que los primeros términos de esta sucesión son (ver A003010): 4, 14, 194, 37634, 1416317954, 2005956546822746114, 4023861667741036022825635656102100994 y 16191462721115671781777559070120513664958590125499158514329308740975788034.

Ahora bien, con esto en mente, lo que dice el test de Lucas-Lehmer es que si tenemos un número de Mersenne Mp, entonces, si s_p-2 Mod Mp=0 el número es primo. Si no, es compuesto.

Por ejemplo, pensemos en el número de Mersenne 31 (2⁵-1). dado que p=5, entonces agarramos s₃ que es 37634. 37634/31=1214 y el módulo es 0, por lo tanto es primo.

Ahora intentemos hacer esto en Python.

El primer paso es hacer que dado un número p, nos devuelva s_p-2. Esta parte parece sencilla:

def lucasLehmer(numero):
    s=4
    for i in range(1,numero-1):
        s=(s**2)-2
    return s

n=input("Ingrese p: ")
print lucasLehmer(n)

El segundo paso es obtener el número de Mersenne a partir de p. Eso ya lo teníamos, de hecho no hace falta más que hacer la cuenta, y por último hacer el módulo. Luego si el módulo==0 es primo, si no, no lo es.

from datetime import *

def lucasLehmer(numero):
    s=4
    for i in range(1,numero-1):
        s=(s**2)-2
    return s

p=input("Ingrese p: ")
print datetime.today()
ll=lucasLehmer(p)
mersenne=(2**p)-1
if ll%mersenne==0:
    print str(mersenne)+" es primo"
else:
    print str(mersenne)+" no es primo"
print datetime.today()

En principio funciona, por ejemplo para p=19 (524287)

2013-07-01 12:09:22.821000
524287 es primo
2013-07-01 12:09:22.832000

son 11 centésimas contra las 14 del método anterior, así que algo avanzamos (también lo hice correr en 10 centésimas, vamos mejorando). Con p=31 al menos no colgaba la máquina, pero ni idea cuanto puede tardar, lo podría dejar andando toda la noche a ver que pasa y les cuento.

Si quieren la demostración aquí la tienen.

Pronto tendremos más tests de primalidad, hay mucho en este area para calcular.