very simple Mersenne prime search

lucas_lehmer.py

# very simple Mersenne prime search.
# AW2016-01-23

import time
import math

# http://stackoverflow.com/questions/16004407/a-fast-prime-number-sieve-in-python
def prime_sieve(n):
    size = n//2
    sieve = [1]*size
    limit = int(n**0.5)
    for i in range(1,limit):
        if sieve[i]:
            val = 2*i+1
            tmp = ((size-1) - i)//val 
            sieve[i+val::val] = [0]*tmp
    return [2] + [i*2+1 for i, v in enumerate(sieve) if v and i>0]

# https://en.wikipedia.org/wiki/Trial_division
def trial_division(n):
    """Return a list of the prime factors for a natural number."""
    if n < 2:
        return []
    prime_factors = []
    for p in prime_sieve(int(n**0.5) + 1):
        if p*p > n: break
        while n % p == 0:
            prime_factors.append(p)
            n //= p
    if n > 1:
        prime_factors.append(n)
    return prime_factors

def is_prime_trial_division(n):
    return len(trial_division(n)) == 1

# https://en.wikipedia.org/wiki/Lucas%E2%80%93Lehmer_primality_test
def is_prime_lucas_lehmer(p):
    s = 4
    M = pow(2,p) -1
    for n in range(p-2):
        s = (pow(s,2)-2) % M
    return s==0

# http://www.mersenne.org/primes/
# year of discovery
history = [-500, -500, -275, -275, 1456, 1588, 1588, 1772, 1883, 1911, 1914, 1876, 1952, 1952, 1952, 1952, 1952, 1957, 1961, 1961, 1963,1963,1963]
n=2 
t0 = time.time()
print "#\tp\ttime\tyear"
log_max_p=15 
for p in range(1,pow(2,log_max_p)):
    if is_prime_trial_division(p):
        if is_prime_lucas_lehmer(p):
            elapsed = time.time() - t0
            print "%4d\t%4d\t%.1f s\t%d" % (n, p, elapsed,history[n-1])
            n=n+1
"""
Typical output: (i7 laptop)
#	p	time	year
   2	   3	0.0 s	-500
   3	   5	0.0 s	-275
   4	   7	0.0 s	-275
   5	  13	0.0 s	1456
   6	  17	0.0 s	1588
   7	  19	0.0 s	1588
   8	  31	0.0 s	1772
   9	  61	0.0 s	1883
  10	  89	0.0 s	1911
  11	 107	0.0 s	1914
  12	 127	0.0 s	1876
  13	 521	0.1 s	1952
  14	 607	0.1 s	1952
  15	1279	0.8 s	1952
  16	2203	4.5 s	1952
  17	2281	5.1 s	1952
  18	3217	17.3 s	1957
  19	4253	49.2 s	1961
  20	4423	56.0 s	1961  
  21	9689	959.8 s	1963
  22	9941	1072.3 s	1963
  
"""