## Mersenne Primes

### June 3, 2011

Here is our version of the Lucas-Lehmer test for Mersenne primes:

```(define (mersenne-prime? p)   (let ((m (- (expt 2 p) 1)))     (if (not (prime? m)) #f       (let loop ((i 3) (s 4))         (if (< p i) (zero? s)           (loop (+ i 1)             (modulo (- (* s s) 2) m)))))))```

We made two changes to the algorithm stated on the previous page: First, the S sequence is accumulated modulo 2p−1, which keeps the numbers smaller. Second, the divisibility test has become a test for zero-congruence modulo 2p−1, which is necessary because we no longer have the true value of S. The indices of the Mersenne primes less than 256 are shown below; we had to handle M2 specially, since the Lucas-Lehmer test specifies only odd primes:

```> (cons 2 (filter mersenne-prime? (range 256))) (2 3 5 7 13 17 19 31 61 89 107 127)```

We used range and filter from the Standard Prelude, and `prime?` from a previous exercise. You can run the program at http://programmingpraxis.codepad.org/VEUVtatQ.

Pages: 1 2

### 12 Responses to “Mersenne Primes”

```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]
```
2. Graham said

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.

3. @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.

4. Graham said

@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."

5. @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.

6. Graham said

@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)

7. programmingpraxis said

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

8. Graham said

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

9. Graham said

I reused a revised version of the “Python Cookbook’s” Sieve of Eratosthenes recipe for `primes()`. Also, my original code computed the function `s` on 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)]
```
10. David said

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 }

11. Sam Kennedy said

I managed to find 2^44497 – 1 was prime before giving up :)

12. David said

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)
```