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

Follow

Get every new post delivered to your Inbox.

Join 598 other followers

%d bloggers like this: