## Programming With Prime Numbers: Source Code In Scheme

```; 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)```