[ create a new paste ] login | about

Project: programmingpraxis
Link: http://programmingpraxis.codepad.org/ZnZycmvr    [ raw code | output | fork ]

programmingpraxis - Scheme, pasted on Apr 1:
; 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) #f
          (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 loop ((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))
                        (loop (/ n d) d (cons d fs)))
                      (else (loop n (+ d 2) fs)))))))))

(define prime?
  (let ((ps (primes 100)))
    (lambda (n)
      (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
            (if (member n ps) #t
              (do ((ps ps (cdr ps)))
                  ((or (null? ps) (not (spsp? n (car ps))))
                    (null? ps)))))))))

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


Output:
1
2
3
4
5
6
7
(2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)
78498
#f
(71 839 1471 6857)
#f
#t
(71 839 1471 6857)


Create a new paste based on this one


Comments: