; programming with prime numbers (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 (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)