;;; perfect.scm --- compute mersenne primes and perfect numbers

;; Author: Noah Friedman <friedman@splode.com>
;; Created: 1997-07-05
;; Public domain

;; $Id: perfect.scm,v 1.4 1997/07/08 11:10:38 friedman Exp $

;;; Commentary:

;; This implementation is probably specific to MIT Scheme.
;; In particular, the use of bitstrings may need to be substituted for
;; vectors in other implementations; but allocating large bitstrings is
;; orders of magnitude faster than vectors of equivalent size in MIT Scheme.

;; To run, compile this code, load it, then evaluate one of
;; (print-mersenne-expts) or (print-perfect-numbers).

;; TODO: dynamically increase primes table when end of table is reached.

;;; Code:

;; Allows open coding of some primitive procedures in MIT Scheme.
(declare (usual-integrations))

;; If not false, print all mersenne numbers as they are tested for primality.
(define perfect-verbose #t)

(define primes-table #f)
(define primes-table-length 0)

(define (log2 n)
  (/ (log n) (log 2)))

(define (prime-initial-table-size)
  ;; 2^21 == 2097152, which is larger than the largest known mersenne expt.
  ;; (log2 n) estimates the number of binary bits required to
  ;; represent n.  Taking the ceiling roughly doubles the size of the table.
  ;; But always allocated at least 2^19==524288.
  (max (ceiling->exact (log2 (highest-computed-mersenne-expt)))
       (expt 2 19)))

(define (initialize)
  (set! primes-table (make-primes-table (prime-initial-table-size)))
  (set! primes-table-length (bit-string-length primes-table)))

(define (make-primes-table to)
  (define tbl (make-bit-string to #t))
  (define (cross-out-multiple! i index)
    (if (< i to)
        (begin
          (bit-string-clear! tbl i)
          (cross-out-multiple! (+ i index) index))))
  (define (cross-out! index)
    (if (< index to)
        (begin
          (if (bit-string-ref tbl index)
              (cross-out-multiple! (* 2 index) index))
          (cross-out! (+ 1 index)))))
  (bit-string-clear! tbl 0)
  (bit-string-clear! tbl 1)
  (cross-out! 2)
  tbl)

(define (next-prime-after n)
  (define (next i)
    (if (bit-string-ref primes-table i)
        i
        (next (+ 1 i))))
  (or primes-table
      (initialize))
  (next (+ 1 n)))

(define (primes-table->list tbl)
  (do ((i 0 (+ 1 i))
       (bsl (bit-string-length tbl))
       (l '()))
      ((= i bsl)
       (reverse! l))
    (and (bit-string-ref tbl i)
         (set! l (cons i l)))))

(define (prime? n)
  (bit-string-ref primes-table n))


;; Note: as of 1997-07-08, no one has verified that there are no mersenne
;; primes between M(859433) and M(1257787).  The highest mersenne number
;; that has been checked to be composite is M(1042100).
;; (define mersenne-expts
;;   '(2 3 5 7 13 17 19 31 61 89 107 127 521 607 1279 2203 2281 3217
;;     4253 4423 9689 9941 11213 19937 21701 23209 44497 86243 110503
;;     132049 216091 756839 859433 1257787 1398269))
(define mersenne-expts '(2))

(define mersenne-expts-tail
  (list-tail mersenne-expts (- (length mersenne-expts) 1)))

(define (highest-computed-mersenne-expt)
  (car mersenne-expts-tail))

(define (first-computed-mersenne-expt)
  (car mersenne-expts))

;; This is used for aborted searches: if you quit the search, you can
;; inspect this variable, then restart the computation later by calling
;; next-mersenne-expt-after with this value.
(define last-mersenne-expt-tried (highest-computed-mersenne-expt))

;; The mersenne number corresponding to last-mersenne-expt-tried.
;; This is used by extend-mersenne-expts! to avoid recomputing 2^n.
;; Instead, we use the fact that 2^(i+j)==(2^i)(2^j) to do incremental
;; multiplication.
(define last-mersenne-tried+1 (expt 2 last-mersenne-expt-tried))

;; Returns m such that (2^m)-1 is the next mersenne prime after (2^n)-1.
;; (2^n)-1 need not itself be a mersenne prime.
(define (next-mersenne-expt-after n)
  (cond ((>= n (highest-computed-mersenne-expt))
         (do ()
             ((< n (highest-computed-mersenne-expt))
              (highest-computed-mersenne-expt))
           (extend-mersenne-expts!)))
        (else
         (do ((m mersenne-expts (cdr m)))
             ((> (car m) n)
              (car m))))))

(define (extend-mersenne-expts!)
  (let* ((next (next-prime-after last-mersenne-expt-tried))
         (c (* last-mersenne-tried+1
               (expt 2 (- next last-mersenne-expt-tried)))))
    (output-verbose ";Trying:")
    (do ((n next next)
         (m (- c 1) (- c 1)))
        ((mersenne-prime? m n)
         (output-verbose "\n")
         (set! last-mersenne-expt-tried n)
         (set! last-mersenne-tried+1 c)
         (set-cdr! mersenne-expts-tail (cons n '()))
         (set! mersenne-expts-tail (cdr mersenne-expts-tail)))
      (set! last-mersenne-expt-tried n)
      (set! last-mersenne-tried+1 c)
      (set! next (next-prime-after n))
      (set! c (* c (expt 2 (- next n)))))))

;; Lucas-Lehmer test for mersenne prime candidates.
;; p is the exponent base, n is (2^p)-1
;; (pass p as an arg to avoid needing to decompose (2^p)-1)
(define (mersenne-prime? n p)
  (define (lucas-lehmer? s i)
    (if (zero? i)
        (zero? s)
        (lucas-lehmer? (modulo (- (expt s 2) 2) n) (- i 1))))
  (output-verbose " " p)
  (lucas-lehmer? 4 (- p 2)))


(define (perfect m)
  (let ((c (expt 2 (- m 1))))
    (* c (- (* 2 c) 1))))

(define (output . seqs)
  (define (output-iter s)
    (cond ((null? s))
          (else
           (display (car s))
           (output-iter (cdr s)))))
  (output-iter seqs))

(define (output-verbose . seqs)
  (and perfect-verbose
       (apply output seqs)))

(define (print-perfect-numbers)
  (define (iter n i)
    (output "n[" i "] = " n ", p(n) = " (perfect n) "\n")
    (iter (next-mersenne-expt-after n) (+ 1 i)))
  (output "\n")
  (iter (first-computed-mersenne-expt) 0))

(define (print-mersenne-expts)
  (define (iter n i)
    (output "n[" i "] = " n "\n")
    (iter (next-mersenne-expt-after n) (+ 1 i)))
  (output "\n")
  (iter (first-computed-mersenne-expt) 0))

;;; perfect.scm ends here.