[ create a new paste ] login | about

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

programmingpraxis - Scheme, pasted on Apr 14:
; proving primality

(define (expm b e m)
  (define (m* x y) (modulo (* x y) m))
  (cond ((zero? e) 1)
        ((even? e) (expm (m* b b) (/ e 2) m))
        (else (m* b (expm (m* b b) (/ (- e 1) 2) m)))))

(define (td-factors n)
  (let loop ((n n) (x 2) (fs '()))
    (cond ((< n (* x x)) (reverse (cons n fs)))
          ((zero? (modulo n x)) (loop (/ n x) x (cons x fs)))
          (else (loop n (+ x 1) fs)))))

(define (prime? n)
  (if (even? n) (= n 2)
    (let* ((n-1 (- n 1)) (fs (td-factors n-1)))
      (let loop1 ((b 2))
        (cond ((= b n) #f)
              ((= (expm b n-1 n) 1)
                (let loop2 ((qs fs))
                  (cond ((null? qs) #t)
                        ((= (expm b (/ n-1 (car qs)) n) 1)
                          (loop1 (+ b 1)))
                        (else (loop2 (cdr qs))))))
              (else (loop1 (+ b 1))))))))

(display (prime? (- (expt 2 89) 1))) (newline)

(define (primes n)
  (let* ((max-index (quotient (- n 3) 2))
         (v (make-vector (+ 1 max-index) #t)))
    (let loop ((i 0) (ps '(2)))
      (let ((p (+ i i 3)) (startj (+ (* 2 i i) (* 6 i) 3)))
        (cond ((>= (* p p) n)
               (let loop ((j i) (ps ps))
                  (cond ((> j max-index) (reverse ps))
                        ((vector-ref v j)
                          (loop (+ j 1) (cons (+ j j 3) ps)))
                        (else (loop (+ j 1) ps)))))
              ((vector-ref v i)
                (let loop ((j startj))
                  (if (<= j max-index)
                      (begin (vector-set! v j #f)
                             (loop (+ j p)))))
                      (loop (+ 1 i) (cons p ps)))
              (else (loop (+ 1 i) ps)))))))

(define (test-prime? n)
  (let ((ps (primes n)))
    (do ((i 2 (+ i 1))) ((= i n))
      (let ((p? (prime? i)))
        (when (and p? (not (member i ps)))
          (display i) (display " found to be prime") (newline))
        (when (and (not p?) (member i ps))
          (display i) (display " found to be composite") (newline))))))

(test-prime? 200) ; no news is good news


Output:
1
#t


Create a new paste based on this one


Comments: