; proving primality
(define (expm b e m)
(define (m* x y) (modulo (* x y) m))
(cond ((zero? e) 1)
((even? e) (expm (m* b b) (/ e 2) m))
(else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))
(define (td-factors n)
(let loop ((n n) (x 2) (fs '()))
(cond ((< n (* x x)) (reverse (cons n fs)))
((zero? (modulo n x)) (loop (/ n x) x (cons x fs)))
(else (loop n (+ x 1) fs)))))
(define (prime? n)
(if (even? n) (= n 2)
(let* ((n-1 (- n 1)) (fs (td-factors n-1)))
(let loop1 ((b 2))
(cond ((= b n) #f)
((= (expm b n-1 n) 1)
(let loop2 ((qs fs))
(cond ((null? qs) #t)
((= (expm b (/ n-1 (car qs)) n) 1)
(loop1 (+ b 1)))
(else (loop2 (cdr qs))))))
(else (loop1 (+ b 1))))))))
(display (prime? (- (expt 2 89) 1))) (newline)
(define (primes n)
(let* ((max-index (quotient (- n 3) 2))
(v (make-vector (+ 1 max-index) #t)))
(let loop ((i 0) (ps '(2)))
(let ((p (+ i i 3)) (startj (+ (* 2 i i) (* 6 i) 3)))
(cond ((>= (* p p) n)
(let loop ((j i) (ps ps))
(cond ((> j max-index) (reverse ps))
((vector-ref v j)
(loop (+ j 1) (cons (+ j j 3) ps)))
(else (loop (+ j 1) ps)))))
((vector-ref v i)
(let loop ((j startj))
(if (<= j max-index)
(begin (vector-set! v j #f)
(loop (+ j p)))))
(loop (+ 1 i) (cons p ps)))
(else (loop (+ 1 i) ps)))))))
(define (test-prime? n)
(let ((ps (primes n)))
(do ((i 2 (+ i 1))) ((= i n))
(let ((p? (prime? i)))
(when (and p? (not (member i ps)))
(display i) (display " found to be prime") (newline))
(when (and (not p?) (member i ps))
(display i) (display " found to be composite") (newline))))))
(test-prime? 200) ; no news is good news