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.

About these ads

Pages: 1 2

13 Responses to “Mersenne Primes”

  1. 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]
    
  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)
    
  13. Claude said

    3single dot com is a great site. It has live chat so you can talk to other people online. Good luck.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Follow

Get every new post delivered to your Inbox.

Join 611 other followers

%d bloggers like this: