; carmichael numbers
(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 (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)))))10
(loop (+ 1 i) (cons p ps)))
(else (loop (+ 1 i) ps)))))))
(define (fermat-pseudo-prime? a n)
(and (= (gcd a n) 1)
(= (expm a (- n 1) n) 1)))
(display (fermat-pseudo-prime? 2 561)) (newline)
(display (fermat-pseudo-prime? 3 561)) (newline)
(define (strong-pseudo-prime? a n)
(let loop ((r 0) (s (- n 1)))
(if (even? s) (loop (+ r 1) (/ s 2))
(if (= (expm a s n) 1) #t
(let loop ((j 0) (s s))
(cond ((= j r) #f)
((= (expm a s n) (- n 1)) #t)
(else (loop (+ j 1) (* s 2)))))))))
(display (strong-pseudo-prime? 2 561)) (newline)
(display (strong-pseudo-prime? 3 561)) (newline)
(define pseudo-prime?
(let ((as (primes 200)))
(lambda (method? n)
(if (member n as) #t
(let loop ((as as))
(cond ((null? as) #t)
((not (method? (car as) n)) #f)
(else (loop (cdr as)))))))))
(define (fermat-prime? n)
(pseudo-prime? fermat-pseudo-prime? n))
(define (strong-prime? n)
(pseudo-prime? strong-pseudo-prime? n))
(display (fermat-prime? 561)) (newline)
(display (strong-prime? 561)) (newline)
(define (carmichael? n)
(let ((ps (primes n)))
(if (= (car (reverse ps)) n) #f
(let loop ((ps (primes n)))
(cond ((null? ps) #t)
((not (= (gcd (car ps) n) 1)) (loop (cdr ps)))
((not (fermat-pseudo-prime? (car ps) n)) #f)
(else (loop (cdr ps))))))))
(display
(let loop ((c 4000) (cs '()))
(cond ((= c 2) cs)
((carmichael? c)
(loop (- c 1) (cons c cs)))
(else (loop (- c 1) cs)))))
(newline)
(define (factors n)
(if (even? n) (cons 2 (factors (/ n 2)))
(let loop ((n n) (f 3) (fs '()))
(cond ((< n (* f f)) (reverse (cons n fs)))
((zero? (modulo n f)) (loop (/ n f) f (cons f fs)))
(else (loop n (+ f 2) fs))))))
(define (carmichael? n)
(let ((fs (factors n)))
(and (odd? n) (pair? (cdr fs)) (apply < fs)
(let loop ((fs fs))
(cond ((null? fs) #t)
((not (zero? (modulo (- n 1) (- (car fs) 1)))) #f)
(else (loop (cdr fs))))))))
(display
(let loop ((c 10000) (cs '()))
(cond ((= c 2) cs)
((carmichael? c)
(loop (- c 1) (cons c cs)))
(else (loop (- c 1) cs)))))