(require 'stdlib)

(defun mersenne-prime-p (p)
  "Determines whether 2^p - 1 is a prime, assuming that p is prime."
  (let ((m (- (expt 2 p) 1)))
    (do ((i 3 (+ i 1))
	 (s 4 (mod (- (* s s) 2)
		   m)))
	((= i (+ p 1)) (= s 0)))))

(defun prime-p (number)
  "Determines whether p is prime."
  (catch 'return
    (do ((i 2 (+ i 1)))
        ((> i (sqrt number)) t)
      (when (= (mod number i) 0)
        (throw 'return nil)))))

(defun find-mersenne-primes (max)
  "Lists all Mersenne primes less than 2^max - 1."
  (dotimes i (- max 3)
    (let ((n (+ i 3)))
      (when (and (prime-p n) (mersenne-prime-p n))
	(printf "2^%d - 1 is a prime.\n" n)))))

;;(find-mersenne-primes 2000)