[ create a new paste ] login | about

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

programmingpraxis - Scheme, pasted on Feb 26:
; goldbach's conjecture

(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 (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 (digits n . args)
  (let ((b (if (null? args) 10 (car args))))
    (let loop ((n n) (d '()))
      (if (zero? n) d
          (loop (quotient n b)
                (cons (modulo n b) d))))))

(define (isqrt n)
  (let loop ((x n) (y (quotient (+ n 1) 2)))
    (if (<= 0 (- y x) 1) x
      (loop y (quotient (+ y (quotient n y)) 2)))))

(define (square? n)
  (let ((n2 (isqrt n)))
    (= n (* n2 n2))))

(define (jacobi a n)
  (if (not (and (integer? a) (integer? n) (positive? n) (odd? n)))
      (error 'jacobi "modulus must be positive odd integer")
      (let jacobi ((a a) (n n))
        (cond ((= a 0) 0)
              ((= a 1) 1)
              ((= a 2) (case (modulo n 8) ((1 7) 1) ((3 5) -1)))
              ((even? a) (* (jacobi 2 n) (jacobi (quotient a 2) n)))
              ((< n a) (jacobi (modulo a n) n))
              ((and (= (modulo a 4) 3) (= (modulo n 4) 3)) (- (jacobi n a)))
              (else (jacobi n a))))))

(define legendre jacobi)

(define (inverse x n)
  (let loop ((x (modulo x n)) (a 1))
    (cond ((zero? x) (error 'inverse "division by zero"))
          ((= x 1) a)
          (else (let ((q (- (quotient n x))))
                  (loop (+ n (* q x)) (modulo (* q a) n)))))))

(define (miller? n a)
  (let loop ((r 0) (s (- n 1)))
    (if (even? s) (loop (+ r 1) (/ s 2))
      (if (= (expm a s n) 1) #t
        (let loop ((r r) (s s))
          (cond ((zero? r) #f)
                ((= (expm a s n) (- n 1)) #t)
                (else (loop (- r 1) (* s 2)))))))))

(define (chain m f g x0 x1)
  (let loop ((ms (digits m 2)) (u x0) (v x1))
    (cond ((null? ms) (values u v))
          ((zero? (car ms)) (loop (cdr ms) (f u) (g u v)))
          (else (loop (cdr ms) (g u v) (f v))))))

(define (lucas? n)
  (let loop ((a 11) (b 7))
    (let ((d (- (* a a) (* 4 b))))
      (cond ((square? d) (loop (+ a 2) (+ b 1)))
            ((not (= (gcd n (* 2 a b d)) 1))
              (loop (+ a 2) (+ b 2)))
            (else (let* ((x1 (modulo (- (* a a (inverse b n)) 2) n))
                         (m (quotient (- n (legendre d n)) 2))
                         (f (lambda (u) (modulo (- (* u u) 2) n)))
                         (g (lambda (u v) (modulo (- (* u v) x1) n))))
                    (let-values (((xm xm1) (chain m f g 2 x1)))
                      (zero? (modulo (- (* x1 xm) (* 2 xm1)) n)))))))))

(define (prime? n)
  (cond ((or (not (integer? n)) (< n 2))
          (error 'prime? "must be integer greater than one"))
        ((even? n) (= n 2)) ((zero? (modulo n 3)) (= n 3))
        (else (and (miller? n 2) (miller? n 3) (lucas? n)))))

(define (goldbach1 n)
  (let loop ((ps (primes n)))
    (if (member (- n (car ps)) ps)
        (list (car ps) (- n (car ps)))
        (loop (cdr ps)))))

(define (goldbach2 n)
  (let ((ps '(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 101 103 107 109 113 127 131 137
         139 149 151 157 163 167 173 179 181 191 193 197 199)))
    (let loop ((ps ps))
      (cond ((null? ps)
              (let loop ((p 211))
                (cond ((< n p) (error 'goldbach "conjecture fails"))
                      ((and (prime? p) (prime? (- n p))) (list p (- n p)))
                      (else (loop (+ p 2))))))
            ((prime? (- n (car ps)))
              (list (car ps) (- n (car ps))))
            (else (loop (cdr ps)))))))

(display (goldbach1 28)) (newline)
(display (goldbach2 28)) (newline)
(display (goldbach2 986332)) (newline)


Output:
1
2
3
(5 23)
(5 23)
(353 985979)


Create a new paste based on this one


Comments: