Mersenne Primes
June 3, 2011
Numbers of the form Mn = 2n−1 are known as Mersenne numbers, named for the French monk Marin Mersenne who studied them early in the seventeenth century. Of the infinite set of Mersenne numbers, 47 are currently known to be prime; see Sloane’s A000043 for a list of their indices. Mersenne primes can be identified by the Lucas-Lehmer test, devised by Édouard Lucas in the 1870s and cast into its modern form by Derrick Lehmer in 1930:
For p an odd prime, the Mersenne number 2p−1 is prime if and only if 2p−1 divides Sp−1 where Sn+1 = Sn2−2 and S1 = 4.
The special form of Mersenne primes makes them easy to identify, and for many years the largest known prime has been a Mersenne prime. A cooperative project on the internet, the Great Internet Mersenne Prime Search (GIMPS), has found all of the recent Mersenne primes, because the numbers have grown so large that a single computer can’t handle the workload.
Your task is to use the Lucas-Lehmer test to find the Mersenne primes through M256. When you are finished, you are welcome to read or run a suggested solution, or to post your own solution or discuss the exercise in the comments below.
Programming Praxis 
My Haskell solution (see http://bonsaicode.wordpress.com/2011/06/03/programming-praxis-mersenne-primes/ for a version with comments):
import Data.Numbers.Primes mersennes :: [Int] mersennes = [p | p <- primes, iterate (\n -> mod (n^2 - 2) (2^p - 1)) 4 !! p-2 == 0] main :: IO () main = print $ takeWhile (<= 256) mersennes == [2,3,5,7,13,17,19,31,61,89,107,127]Perhaps it’s too early for me to try these without coffee, but I think there may be some mistakes. M_256 = 2^256 – 1, while not prime, is a 78 digit number. Moreover, 2 is not one less than a power of 2, so it’s not a Mersenne Prime.
@Graham: For your first point, both version give the exponents. The actual primes are trivially obtained by calculating 2^p – 1. For your second point, since we’re calculating the exponents, the actual value is 2^2 – 1 = 3, which is prime.
@Remco: Thanks, but it seems I’m still not getting it. For the second point, yes, 2^2-1 = 3 is prime, but that’s M_2. However, there is no n such that M_n = 2, since 2^n – 1 is an odd number. Hence, 2 is not a Mersenne Number, and therefore not a Mersenne Prime. I’m not quite sure what you mean for the first point. Basically, it seems to me that the question as written is not “Find all primes p less than 256 such that 2^p – 1 is prime,” but “Find all numbers of the form M_n = 2^n – 1 where M_n is prime for 2 < n < 256."
@Graham: From the wikipedia article on Mersenne primes: “Some definitions of Mersenne numbers require that the exponent p be prime, since the associated number must be composite if p is.” Testing only primes is an optimization since if p is not prime, 2^p – 1 won’t be. The definitions “Find all primes p less than 256 such that 2^p – 1 is prime” and “Find all numbers of the form M_n = 2^n – 1 where M_n is prime for 2 < n < 256." are therefore equivalent, save for the fact that the former returns p and the second 2^p – 1, which is the trivial transformation I mentioned.
@Remco: Thanks for the point about only testing primes from 2 to 256 instead of all integers. I hadn’t seen that. Appreciate you setting me straight! (and sorry for being so obtuse)
Graham: If you were confused about definitions, that’s my fault, since I set the terms of the exercise.
Mersenne numbers are numbers of the form 2^n-1, with n>0: 1, 3, 7, 15, 31, …. (Sometimes the sequence is extended to n=0, so the first term is 0, but we’ll stick with the standard definition.) A Mersenne number is either prime or composite; thus, the Mersenne number 7 is prime, and the Mersenne number 15=3*5 is composite.
Since Mersenne numbers grow so quickly, they are generally referred to by their index. For instance, 162259276829213363391578010288127 is a Mersenne number, but it’s easier to say M107. And it’s also prime, so it’s a Mersenne prime. And, as a special added bonus, we also know that since M107 is prime, so is 107, since it is a necessary (but not sufficient) condition that the index of a Mersenne prime is also a prime.
Sometimes when mathematicians talk about Mersenne numbers they really mean Mersenne numbers such that the index is prime, whether or not the number itself is prime, which is a ready source of confusion; that sequence is 3, 7, 31, 127, 2047, 8191, …, and it’s usually written 2^p-1 instead of 2^n-1. You may notice that Wikipedia uses the first definition, while Sloane uses the second. There is no “official” definition, you just have to sort out for yourself which meaning is intended.
The indices of the Mersenne primes are 2, 3, 5, 7, 13, 17, 19, 31, 61, 89, 107, 127, … and the Mersenne primes themselves are 3, 7, 31, 127, 8191, 131071, 524287, 2147483647, 2305843009213693951, 618970019642690137449562111, 162259276829213363391578010288127, 170141183460469231731687303715884105727, …; notice that M11=2047 is not a Mersenne prime, even though it is a Mersenne number under either definition.
If you want to know more about Mersenne numbers, of all kinds, a good place to start is OEIS (which used to be called “Sloane’s” after its founder N. J. A. Sloane, but he is now retired); just type “mersenne” in the search bar and click “Search.” After you’re finished being amazed, you can read some of the entries there, and follow the references given in each entry, and soon enough your entire evening will be gone. Sloane’s is one of my favorite web sites to sit and browse.
The program you were asked to write was to find the indices of the Mersenne primes to M256. The solution first tests the index, which must be prime, and then performs a Lucas-Lehmer test to determine if the Mersenne number is also prime. The list of indices and the corresponding numbers are shown in the prior paragraph.
Hope this helps,
Phil
Thanks Phil, I must have just been particularly slow yesterday. I appreciate the background, though, and look forward to learning more. Interestingly, the math software I use (Sage) has an interface to Sloane’s built in.
Thanks again to you and Remco for putting up with me and setting me straight, and thank you for this great blog!
Graham
I reused a revised version of the “Python Cookbook’s” Sieve of Eratosthenes recipe for
primes(). Also, my original code computed the functionson its own, but the code was incredibly slow (lots of recursive calls with very large numbers?). I opted for a translation of the Scheme solution.def lucas_lehmer(p): """Determines whether prime p is a Mersenne Prime via Lucas-Lehmer.""" if p == 2: return True m, i, s = pow(2, p) - 1, 3, 4 while i <= p: i += 1 s = (pow(s, 2) - 2) % m return s == 0 if __name__ == "__main__": print [p for p in primes(256) if lucas_lehmer(p)]Factor code. the sequence S1, S2, etc. is a lazily computed list. primes word is sieve of eratosthenes from previous exercise.
USING: kernel math lists lists.lazy sieve fry sequences ; IN: mersenne : lucas-lehmer-seq ( m -- list ) '[ sq 2 - _ mod ] 4 swap lfrom-by ; : lucas-lehmer? ( n -- ? ) dup 2^ 1 - lucas-lehmer-seq swap 2 - [ cdr ] times car 0 = ; : mersennes ( n -- list ) primes [ lucas-lehmer? ] filter 2 prefix ;( scratchpad ) 256 mersennes .
V{ 2 3 5 7 13 17 19 31 61 89 107 127 }
I managed to find 2^44497 – 1 was prime before giving up :)
Revisit in Clojure…
(load "next-prime") (defn square [x] (* x x)) (defn expt [x n] (cond (zero? n) 1 (even? n) (square (expt x (/ n 2))) :else (* x (expt x (dec n))))) (defn lucas-lehmer [p] (iterate (fn [x] (mod (- (square x) 2) p)) 4)) (defn mersenne-prime? [n] (let [p (dec (expt 2N n))] (= (first (drop (- n 2) (lucas-lehmer p))) 0))) (def mersennes (cons 2 (filter mersenne-prime? primes))) user=> (take 20 mersennes) (2 3 5 7 13 17 19 31 61 89 107 127 521 607 1279 2203 2281 3217 4253 4423)