; prime -- a modest library
 
; primes n -- list of primes not greater than n in ascending order
(define (primes n) ; assumes n is an integer greater than one
  (let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
    (let loop ((i 0) (p 3) (ps (list 2))) ; sieve of eratosthenes
      (cond ((< n (* p p))
              (do ((i i (+ i 1)) (p p (+ p 2))
                   (ps ps (if (vector-ref bits i) (cons p ps) ps)))
                  ((= i len) (reverse ps))))
            ((vector-ref bits i)
              (do ((j (+ (* 2 i i) (* 6 i) 3) (+ j p)))
                  ((<= len j) (loop (+ i 1) (+ p 2) (cons p ps)))
                (vector-set! bits j #f)))
            (else (loop (+ i 1) (+ p 2) ps))))))
 
; prime? n -- #f if provably composite, else #t if probably prime
(define prime? ; strong pseudoprime to prime bases less than 100
  (let ((ps (primes 100))) ; assumes n an integer greater than one
    (lambda (n)
      (define (expm b e m)
        (define (times p q) (modulo (* p q) m))
          (let loop ((b b) (e e) (x 1))
            (if (zero? e) x
              (loop (times b b) (quotient e 2)
                    (if (odd? e) (times b x) x)))))
      (define (spsp? n a) ; #t if n is a strong pseudoprime base a
        (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
            ((odd? d)
              (if (= (expm a d n) 1) #t
                (do ((r 0 (+ r 1)))
                    ((or (= (expm a (* d (expt 2 r)) n) (- n 1))
                         (= r s))
                      (< r s)))))))
      (if (member n ps) #t
        (do ((ps ps (cdr ps)))
            ((or (null? ps) (not (spsp? n (car ps))))
              (null? ps)))))))
 
; factors n -- list of prime factors of n in ascending order
(define (factors n) ; assumes n is an integer, may be negative
  (if (<= -1 n 1) (list n) ; pollard's rho algorithm
    (if (< n 0) (cons -1 (factors (- n)))
      (let fact ((n n) (c 1) (fs (list)))
        (define (f x) (modulo (+ (* x x) c) n))
        (if (even? n) (fact (/ n 2) c (cons 2 fs))
          (if (= n 1) fs
            (if (prime? n) (sort (cons n fs) <)
              (let loop ((t 2) (h 2) (d 1))
                (cond ((= d 1) (let ((t (f t)) (h (f (f h))))
                                 (loop t h (gcd (- t h) n))))
                      ((= d n) (fact n (+ c 1) fs)) ; cyclic
                      ((prime? d)
                        (fact (/ n d) (+ c 1) (cons d fs)))
                      (else (fact n (+ c 1) fs)))))))))))
 
(display (primes 50)) (newline)
 
(do ((n 1 (+ n 1)))
    ((< 70 n))
  (display n)
  (display ": ")
  (for-each
    (lambda (x)
      (display " ")
      (display x))
    (factors (- (expt 2 n) 1)))
  (newline))

OyBwcmltZSAtLSBhIG1vZGVzdCBsaWJyYXJ5Cgo7IHByaW1lcyBuIC0tIGxpc3Qgb2YgcHJpbWVzIG5vdCBncmVhdGVyIHRoYW4gbiBpbiBhc2NlbmRpbmcgb3JkZXIKKGRlZmluZSAocHJpbWVzIG4pIDsgYXNzdW1lcyBuIGlzIGFuIGludGVnZXIgZ3JlYXRlciB0aGFuIG9uZQogIChsZXQqICgobGVuIChxdW90aWVudCAoLSBuIDEpIDIpKSAoYml0cyAobWFrZS12ZWN0b3IgbGVuICN0KSkpCiAgICAobGV0IGxvb3AgKChpIDApIChwIDMpIChwcyAobGlzdCAyKSkpIDsgc2lldmUgb2YgZXJhdG9zdGhlbmVzCiAgICAgIChjb25kICgoPCBuICgqIHAgcCkpCiAgICAgICAgICAgICAgKGRvICgoaSBpICgrIGkgMSkpIChwIHAgKCsgcCAyKSkKICAgICAgICAgICAgICAgICAgIChwcyBwcyAoaWYgKHZlY3Rvci1yZWYgYml0cyBpKSAoY29ucyBwIHBzKSBwcykpKQogICAgICAgICAgICAgICAgICAoKD0gaSBsZW4pIChyZXZlcnNlIHBzKSkpKQogICAgICAgICAgICAoKHZlY3Rvci1yZWYgYml0cyBpKQogICAgICAgICAgICAgIChkbyAoKGogKCsgKCogMiBpIGkpICgqIDYgaSkgMykgKCsgaiBwKSkpCiAgICAgICAgICAgICAgICAgICgoPD0gbGVuIGopIChsb29wICgrIGkgMSkgKCsgcCAyKSAoY29ucyBwIHBzKSkpCiAgICAgICAgICAgICAgICAodmVjdG9yLXNldCEgYml0cyBqICNmKSkpCiAgICAgICAgICAgIChlbHNlIChsb29wICgrIGkgMSkgKCsgcCAyKSBwcykpKSkpKQoKOyBwcmltZT8gbiAtLSAjZiBpZiBwcm92YWJseSBjb21wb3NpdGUsIGVsc2UgI3QgaWYgcHJvYmFibHkgcHJpbWUKKGRlZmluZSBwcmltZT8gOyBzdHJvbmcgcHNldWRvcHJpbWUgdG8gcHJpbWUgYmFzZXMgbGVzcyB0aGFuIDEwMAogIChsZXQgKChwcyAocHJpbWVzIDEwMCkpKSA7IGFzc3VtZXMgbiBhbiBpbnRlZ2VyIGdyZWF0ZXIgdGhhbiBvbmUKICAgIChsYW1iZGEgKG4pCiAgICAgIChkZWZpbmUgKGV4cG0gYiBlIG0pCiAgICAgICAgKGRlZmluZSAodGltZXMgcCBxKSAobW9kdWxvICgqIHAgcSkgbSkpCiAgICAgICAgICAobGV0IGxvb3AgKChiIGIpIChlIGUpICh4IDEpKQogICAgICAgICAgICAoaWYgKHplcm8/IGUpIHgKICAgICAgICAgICAgICAobG9vcCAodGltZXMgYiBiKSAocXVvdGllbnQgZSAyKQogICAgICAgICAgICAgICAgICAgIChpZiAob2RkPyBlKSAodGltZXMgYiB4KSB4KSkpKSkKICAgICAgKGRlZmluZSAoc3BzcD8gbiBhKSA7ICN0IGlmIG4gaXMgYSBzdHJvbmcgcHNldWRvcHJpbWUgYmFzZSBhCiAgICAgICAgKGRvICgoZCAoLSBuIDEpICgvIGQgMikpIChzIDAgKCsgcyAxKSkpCiAgICAgICAgICAgICgob2RkPyBkKQogICAgICAgICAgICAgIChpZiAoPSAoZXhwbSBhIGQgbikgMSkgI3QKICAgICAgICAgICAgICAgIChkbyAoKHIgMCAoKyByIDEpKSkKICAgICAgICAgICAgICAgICAgICAoKG9yICg9IChleHBtIGEgKCogZCAoZXhwdCAyIHIpKSBuKSAoLSBuIDEpKQogICAgICAgICAgICAgICAgICAgICAgICAgKD0gciBzKSkKICAgICAgICAgICAgICAgICAgICAgICg8IHIgcykpKSkpKSkKICAgICAgKGlmIChtZW1iZXIgbiBwcykgI3QKICAgICAgICAoZG8gKChwcyBwcyAoY2RyIHBzKSkpCiAgICAgICAgICAgICgob3IgKG51bGw/IHBzKSAobm90IChzcHNwPyBuIChjYXIgcHMpKSkpCiAgICAgICAgICAgICAgKG51bGw/IHBzKSkpKSkpKQoKOyBmYWN0b3JzIG4gLS0gbGlzdCBvZiBwcmltZSBmYWN0b3JzIG9mIG4gaW4gYXNjZW5kaW5nIG9yZGVyCihkZWZpbmUgKGZhY3RvcnMgbikgOyBhc3N1bWVzIG4gaXMgYW4gaW50ZWdlciwgbWF5IGJlIG5lZ2F0aXZlCiAgKGlmICg8PSAtMSBuIDEpIChsaXN0IG4pIDsgcG9sbGFyZCdzIHJobyBhbGdvcml0aG0KICAgIChpZiAoPCBuIDApIChjb25zIC0xIChmYWN0b3JzICgtIG4pKSkKICAgICAgKGxldCBmYWN0ICgobiBuKSAoYyAxKSAoZnMgKGxpc3QpKSkKICAgICAgICAoZGVmaW5lIChmIHgpIChtb2R1bG8gKCsgKCogeCB4KSBjKSBuKSkKICAgICAgICAoaWYgKGV2ZW4/IG4pIChmYWN0ICgvIG4gMikgYyAoY29ucyAyIGZzKSkKICAgICAgICAgIChpZiAoPSBuIDEpIGZzCiAgICAgICAgICAgIChpZiAocHJpbWU/IG4pIChzb3J0IChjb25zIG4gZnMpIDwpCiAgICAgICAgICAgICAgKGxldCBsb29wICgodCAyKSAoaCAyKSAoZCAxKSkKICAgICAgICAgICAgICAgIChjb25kICgoPSBkIDEpIChsZXQgKCh0IChmIHQpKSAoaCAoZiAoZiBoKSkpKQogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAobG9vcCB0IGggKGdjZCAoLSB0IGgpIG4pKSkpCiAgICAgICAgICAgICAgICAgICAgICAoKD0gZCBuKSAoZmFjdCBuICgrIGMgMSkgZnMpKSA7IGN5Y2xpYwogICAgICAgICAgICAgICAgICAgICAgKChwcmltZT8gZCkKICAgICAgICAgICAgICAgICAgICAgICAgKGZhY3QgKC8gbiBkKSAoKyBjIDEpIChjb25zIGQgZnMpKSkKICAgICAgICAgICAgICAgICAgICAgIChlbHNlIChmYWN0IG4gKCsgYyAxKSBmcykpKSkpKSkpKSkpCiAgICAgICAgICAgICAgICAgICAgICAKKGRpc3BsYXkgKHByaW1lcyA1MCkpIChuZXdsaW5lKQoKKGRvICgobiAxICgrIG4gMSkpKQogICAgKCg8IDcwIG4pKQogIChkaXNwbGF5IG4pCiAgKGRpc3BsYXkgIjogIikKICAoZm9yLWVhY2gKICAgIChsYW1iZGEgKHgpCiAgICAgIChkaXNwbGF5ICIgIikKICAgICAgKGRpc3BsYXkgeCkpCiAgICAoZmFjdG9ycyAoLSAoZXhwdCAyIG4pIDEpKSkKICAobmV3bGluZSkp