(use gauche.time)

(define (mersenne-prime? 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 (modulo (- (* s s) 2) m)))
        ((= i (+ p 1)) (= s 0)))))

(define (prime? number)
  "Determines whether p is prime."
  (call/cc
   (lambda (return)
     (do ((i 2 (+ i 1)))
         ((> i (sqrt number)) #t)
       (when (= (modulo number i) 0)
         (return #f))))))

(define (find-mersenne-primes max)
  "Lists all Mersenne primes less than 2^max - 1."
  (dotimes (i (- max 3))
    (let ((n (+ i 3)))
      (when (and (prime? n) (mersenne-prime? n))
        (format #t "2^~d - 1 is a prime.\n" n)))))

(define (find-mersenne-primes-verbose max)
  (dotimes (i (- max 3))
    (let ((n (+ i 3)))
      (when (prime? n)
        (format #t "~d " n) (flush-all-ports)
        (when (mersenne-prime? n)
          (newline)
          (format #t "2^~d - 1 is prime: ~d\n" n (- (expt 2 n) 1))))))
  (newline))

(define (find-primes max)
  (dotimes (i max)
    (when (prime? i)
      (format #t "~s " i)))
  (newline))

;; NOTE: 2 is also a mersenne prime, see mupad: numlib::mersenne()

(time (find-mersenne-primes 1000))

(define mersenne-primes-2000 '(3 5 7 13 17 19 31 61 89 107 127 521 607 1279))
;;(map mersenne-prime? '(3 5 7 13 17 19 31 61 89 107 127 521 607 1279))