Lucas-Lehmer test using Python not working for large numbers

Question

I am trying to write a function that should take in a Mersenne number and return whether it is a prime number or not using the Lucas-Lehmer primality test. I am trying to return the last number generated in the Lucas-Lehmer sequence which should be 0 if it is a prime number

I have written the following function to do the above

def lucas_lehmer_modified(p):
   M=2**p-1
   for i in range (0,p-1):
     if i == 0:
       x = 4
     else : 
        x = (prev**2-2)%(M)
     prev=x
   return x

My problem is that this code works for small numbers such as 127 but not for large numbers like 2305843009213693951 and even for 524287. My Python Jupyter notebook hangs up. Any advice on how I can get a function that takes in a Mersenne prime as an input and returns whether it is a prime number or not using the Lucas Lehmer test. I need to make it work for at least 2^65-1

Answer 1

I wrote a Lucas Lehmer test from the psuedocode from the wikipedia entry, here is what i came up with ( p being the mersenne power number ):

def LucasLehmer(p):
   if p == 2:
     return True
   s = 4
   M = pow(2, p) - 1
   for x in range (1, (p-2)+1):
      s = ((s * s) - 2) % M
   if s == 0: return True
   else: return False

And some answers:

In [388]: LucasLehmer(9689)                                                                                                                                                                     
Out[388]: True

In [389]: LucasLehmer(19937)                                                                                                                                                                    
Out[389]: True

In [390]: LucasLehmer(500)                                                                                                                                                                      
Out[390]: False