codepad
[
create a new paste
]
login
|
about
Language:
C
C++
D
Haskell
Lua
OCaml
PHP
Perl
Plain Text
Python
Ruby
Scheme
Tcl
; 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)
Private
[
?
]
Run code
Submit