; cornacchia's algorithm
; primes n -- list of primes not greater than n in ascending order
(define (primes n) ; assumes n is an integer greater than one
(let* ((len (quotient (- n 1) 2)) (bits (make-vector len #t)))
(let loop ((i 0) (p 3) (ps (list 2))) ; sieve of eratosthenes
(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))))))
; prime? n -- #f if provably composite, else #t if probably prime
(define prime? ; strong pseudoprime to prime bases less than 100
(let ((ps (primes 100))) ; assumes n an integer greater than one
(lambda (n)
(define (expm b e m)
(define (times p q) (modulo (* p q) m))
(let loop ((b b) (e e) (x 1))
(if (zero? e) x
(loop (times b b) (quotient e 2)
(if (odd? e) (times b x) x)))))
(define (spsp? n a) ; #t if n is a strong pseudoprime base a
(do ((d (- n 1) (/ d 2)) (s 0 (+ s 1)))
((odd? d)
(if (= (expm a d n) 1) #t
(do ((r 0 (+ r 1)))
((or (= (expm a (* d (expt 2 r)) n) (- n 1))
(= r s))
(< r s)))))))
(if (member n ps) #t
(do ((ps ps (cdr ps)))
((or (null? ps) (not (spsp? n (car ps))))
(null? ps)))))))
; expm b e m -- modular exponentiation: b^e (mod m)
(define (expm b e m)
(define (times p q) (modulo (* p q) m))
(let loop ((b b) (e e) (x 1))
(if (zero? e) x
(loop (times b b) (quotient e 2)
(if (odd? e) (times b x) x)))))
; jacobi a m -- jacobi symbol
(define (jacobi a m)
(let loop ((a (modulo a m)) (m m) (t 1))
(if (= a 0) (if (= m 1) t 0)
(let twos ((a a) (t t))
(if (even? a)
(twos (quotient a 2)
(if (member (modulo m 8) '(3 5)) (- t) t))
(loop (modulo m a) a
(if (= 3 (modulo a 4) (modulo m 4)) (- t) t)))))))
; legendre a n -- legendre symbol
(define (legendre a n)
(if (and (< 2 n) (prime? n)) (jacobi a n)
(error 'legendre "modulus must be odd prime")))
; mod-sqrt a p -- modular square root: x such that x^2 == a (mod p)
(define (mod-sqrt a p) ; tonelli-shanks algorithm
(define (both n) (list n (- p n)))
(when (not (and (odd? p) (prime? p)))
(error 'mod-sqrt "modulus must be an odd prime"))
(when (not (= (legendre a p) 1))
(error 'mod-sqrt "must be a quadratic residual"))
(let ((a (modulo a p)))
(case (modulo p 8)
((3 7) (both (expm a (/ (+ p 1) 4) p)))
((5) (let* ((x (expm a (/ (+ p 3) 8) p))
(c (expm x 2 p)))
(if (= a c) (both x)
(both (modulo (* x (expm 2 (/ (- p 1) 4) p)) p)))))
(else (let* ((d (let loop ((d 2))
(if (= (jacobi d p) -1) d
(loop (+ d 1)))))
(s (let loop ((p (- p 1)) (s 0))
(if (odd? p) s
(loop (quotient p 2) (+ s 1)))))
(t (quotient (- p 1) (expt 2 s)))
(big-a (expm a t p))
(big-d (expm d t p))
(m (let loop ((i 0) (m 0))
(cond ((= i s) m)
((= (- p 1)
(expm (* big-a (expm big-d m p))
(expt 2 (- s 1 i)) p))
(loop (+ i 1) (+ m (expt 2 i))))
(else (loop (+ i 1) m))))))
(both (modulo (* (expm a (/ (+ t 1) 2) p)
(expm big-d (/ m 2) p)) p)))))))
; isqrt n -- largest integer x such that x*x <= n
(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)))))
; divides? d n -- #t if n/d is integral, else #f
(define (divides? d n) (zero? (modulo n d)))
; square? n -- sqrt(n) if n is a perfect square, else #f
(define square?
(let ((q11 (make-vector 11 #f))
(q63 (make-vector 63 #f))
(q64 (make-vector 64 #f))
(q65 (make-vector 65 #f)))
(do ((k 0 (+ k 1))) ((< 5 k))
(vector-set! q11 (modulo (* k k) 11) #t))
(do ((k 0 (+ k 1))) ((< 31 k))
(vector-set! q63 (modulo (* k k) 63) #t))
(do ((k 0 (+ k 1))) ((< 31 k))
(vector-set! q64 (modulo (* k k) 64) #t))
(do ((k 0 (+ k 1))) ((< 32 k))
(vector-set! q65 (modulo (* k k) 65) #t))
(lambda (n)
(if (not (integer? n)) (error 'square? "must be integer")
(if (< n 1) #f
(if (not (vector-ref q64 (modulo n 64))) #f
(let ((r (modulo n 45045)))
(if (not (vector-ref q63 (modulo r 63))) #f
(if (not (vector-ref q65 (modulo r 65))) #f
(if (not (vector-ref q11 (modulo r 11))) #f
(let ((q (isqrt n)))
(if (= (* q q) n) q #f))))))))))))
(define-syntax fold-of
(syntax-rules (range in is)
((_ "z" f b e) (set! b (f b e)))
((_ "z" f b e (v range fst pst stp) c ...)
(let* ((x fst) (p pst) (s stp)
(le? (if (positive? s) <= >=)))
(do ((v x (+ v s))) ((le? p v) b)
(fold-of "z" f b e c ...))))
((_ "z" f b e (v range fst pst) c ...)
(let* ((x fst) (p pst) (s (if (< x p) 1 -1)))
(fold-of "z" f b e (v range x p s) c ...)))
((_ "z" f b e (v range pst) c ...)
(fold-of "z" f b e (v range 0 pst) c ...))
((_ "z" f b e (x in xs) c ...)
(do ((t xs (cdr t))) ((null? t) b)
(let ((x (car t)))
(fold-of "z" f b e c ...))))
((_ "z" f b e (x is y) c ...)
(let ((x y)) (fold-of "z" f b e c ...)))
((_ "z" f b e p? c ...)
(if p? (fold-of "z" f b e c ...)))
((_ f i e c ...)
(let ((b i)) (fold-of "z" f b e c ...)))))
(define-syntax list-of (syntax-rules ()
((_ arg ...) (reverse (fold-of
(lambda (d a) (cons a d)) '() arg ...)))))
(define (cornacchia d p)
(if (negative? (jacobi (- d) p)) #f
(let loop ((a p) (b (apply max (mod-sqrt (- d) p))))
(if (< p (* b b)) (loop b (modulo a b))
(let* ((n (- p (* b b))) (c (/ n d)))
(cond ((not (divides? d n)) #f)
((not (square? c)) #f)
(else (list b (isqrt c)))))))))
(display (cornacchia 4 1733)) (newline)
(display (cornacchia 58 3031201)) (newline)
(display
(list-of (cons p (cornacchia 1 p))
(p in (primes 1000))
(= (modulo p 4) 1)))
(newline)
programmingpraxis
-
Scheme,
pasted
on Apr 2: