; 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 (goldbach2 79)) (newline)
(display (goldbach2 28)) (newline)
(display (goldbach2 986332)) (newline)