`; Programming With Prime Numbers`

; http://programmingpraxis.files.wordpress.com/

; 2012/09/programmingwithprimenumbers.pdf

```
```(define (sieve n)

(let ((bits (make-vector (+ n 1) #t)))

(let loop ((p 2) (ps '()))

(cond ((< n p) (reverse ps))

((vector-ref bits p)

(do ((i (+ p p) (+ i p))) ((< n i))

(vector-set! bits i #f))

(loop (+ p 1) (cons p ps)))

(else (loop (+ p 1) ps))))))

(define (primes n)

(if (or (not (integer? n)) (< n 2))

(error 'primes "must be integer greater than one")

(let* ((len (quotient (- n 1) 2))

(bits (make-vector len #t)))

(let loop ((i 0) (p 3) (ps (list 2)))

(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)))))))

(define (td-prime? n . args)

(if (or (not (integer? n)) (< n 2))

(error 'td-prime? "must be integer greater than one")

(let ((limit (if (pair? args) (car args) 1000000)))

(if (even? n) (= n 2)

(let loop ((d 3))

(cond ((< limit d)

(error 'td-prime? "limit exceeded"))

((< n (* d d)) #t)

((zero? (modulo n d)) #f)

(else (loop (+ d 2)))))))))

(define (td-factors n . args)

(if (or (not (integer? n)) (< n 2))

(error 'td-factors "must be integer greater than one")

(let ((limit (if (pair? args) (car args) 1000000)))

(let twos ((n n) (fs '()))

(if (even? n)

(twos (/ n 2) (cons 2 fs))

(let odds ((n n) (d 3) (fs fs))

(cond ((< limit d)

(error 'td-factors "limit exceeded"))

((< n (* d d))

(reverse (cons n fs)))

((zero? (modulo n d))

(odds (/ n d) d (cons d fs)))

(else (odds n (+ d 2) fs)))))))))

(define prime?

(let ((seed 3141592654))

(lambda (n)

(define (rand)

(set! seed (modulo (+ (* 69069 seed) 1234567) 4294967296))

(+ (quotient (* seed (- n 2)) 4294967296) 2))

(define (expm b e m)

(define (times x y) (modulo (* x y) m))

(let loop ((b b) (e e) (r 1))

(if (zero? e) r

(loop (times b b) (quotient e 2)

(if (odd? e) (times b r) r)))))

(define (spsp? n a)

(do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))

((odd? d)

(let ((t (expm a d n)))

(if (or (= t 1) (= t (- n 1))) #t

(do ((s (- s 1) (- s 1))

(t (expm t 2 n) (expm t 2 n)))

((or (zero? s) (= t (- n 1)))

(positive? s))))))))

(if (not (integer? n))

(error 'prime? "must be integer")

(if (< n 2) #f

(do ((a (rand) (rand)) (k 25 (- k 1)))

((or (zero? k) (not (spsp? n a)))

(zero? k))))))))

(define (rho-factors n . args)

(define (cons< x xs)

(cond ((null? xs) (list x))

((< x (car xs)) (cons x xs))

(else (cons (car xs) (cons< x (cdr xs))))))

(define (rho n limit)

(let loop ((t 2) (h 2) (d 1) (c 1) (limit limit))

(define (f x) (modulo (+ (* x x) c) n))

(cond ((zero? limit) (error 'rho-factors "limit exceeded"))

((= d 1) (let ((t (f t)) (h (f (f h))))

(loop t h (gcd (- t h) n) c (- limit 1))))

((= d n) (loop 2 2 1 (+ c 1) (- limit 1)))

((prime? d) d)

(else (rho d (- limit 1))))))

(if (not (integer? n))

(error 'rho-factors "must be integer")

(let ((limit (if (pair? args) (car args) 1000)))

(cond ((<= -1 n 1) (list n))

((negative? n) (cons -1 (rho-factors (- n) limit)))

((even? n)

(if (= n 2) (list 2)

(cons 2 (rho-factors (/ n 2) limit))))

(else (let loop ((n n) (fs '()))

(if (prime? n)

(cons< n fs)

(let ((f (rho n limit)))

(loop (/ n f) (cons< f fs))))))))))

`(display (sieve 100)) (newline)`

(display (primes 100)) (newline)

(display (length (primes 1000000))) (newline)

(display (td-prime? 600851475143)) (newline)

(display (td-factors 600851475143)) (newline)

(display (prime? 600851475143)) (newline)

(display (prime? 2305843009213693951)) (newline)

(display (rho-factors 600851475143)) (newline)