codepad
[
create a new paste
]
login
|
about
Language:
C
C++
D
Haskell
Lua
OCaml
PHP
Perl
Plain Text
Python
Ruby
Scheme
Tcl
; combined n +/- 1 primality prover ; 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) (let loop ((b b) (e e) (x 1)) (if (zero? e) x (loop (modulo (* b b) m) (quotient e 2) (if (odd? e) (modulo (* b x) m) 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 (< n 2) #f (if (member n ps) #t (do ((ps ps (cdr ps))) ((or (null? ps) (not (spsp? n (car ps)))) (null? ps)))))))) ; factors n -- list of prime factors of n in ascending order (define (factors n) ; assumes n is an integer, may be negative (if (<= -1 n 1) (list n) ; pollard's rho algorithm (if (< n 0) (cons -1 (factors (- n))) (let fact ((n n) (c 1) (fs (list))) (define (f x) (modulo (+ (* x x) c) n)) (if (even? n) (fact (/ n 2) c (cons 2 fs)) (if (= n 1) fs (if (prime? n) (sort < (cons n fs)) (let loop ((t 2) (h 2) (d 1)) (cond ((= d 1) (let ((t (f t)) (h (f (f h)))) (loop t h (gcd (- t h) n)))) ((= d n) (fact n (+ c 1) fs)) ; cyclic ((prime? d) (fact (/ n d) (+ c 1) (cons d fs))) (else (fact n (+ c 1) fs))))))))))) (define next-prime ; uses 2,3,5,7-wheel (let ((gap (vector 1 10 9 8 7 6 5 4 3 2 1 2 1 4 3 2 1 2 1 4 3 2 1 6 5 4 3 2 1 2 1 6 5 4 3 2 1 4 3 2 1 2 1 4 3 2 1 6 5 4 3 2 1 6 5 4 3 2 1 2 1 6 5 4 3 2 1 4 3 2 1 2 1 6 5 4 3 2 1 4 3 2 1 6 5 4 3 2 1 8 7 6 5 4 3 2 1 4 3 2 1 2 1 4 3 2 1 2 1 4 3 2 1 8 7 6 5 4 3 2 1 6 5 4 3 2 1 4 3 2 1 6 5 4 3 2 1 2 1 4 3 2 1 6 5 4 3 2 1 2 1 6 5 4 3 2 1 6 5 4 3 2 1 4 3 2 1 2 1 4 3 2 1 6 5 4 3 2 1 2 1 6 5 4 3 2 1 4 3 2 1 2 1 4 3 2 1 2 1 10 9 8 7 6 5 4 3 2 1 2))) (lambda (n) (if (< n 2) 2 (if (< n 3) 3 (if (< n 5) 5 (let loop ((n (+ n (if (even? n) 1 2)))) (if (prime? n) n (loop (+ n (vector-ref gap (modulo n 210)))))))))))) (define (expm b e m) ; modular exponentiation (let loop ((b b) (e e) (x 1)) (if (zero? e) x (loop (modulo (* b b) m) (quotient e 2) (if (odd? e) (modulo (* b x) m) x))))) (define (fermat-test? n ps) (let loop ((ps ps) (a 1)) (if (null? ps) #t (if (and (= (expm a (- n 1) n) 1) (= (gcd (expt a (/ (- n 1) (car ps))) n) 1)) (loop (cdr ps) 1) (loop ps (+ a 1)))))) (define (jacobi a m) ; jacobi symbol (let loop1 ((a (modulo a m)) (m m) (t 1)) (if (zero? a) (if (= m 1) t 0) (let ((z (if (member (modulo m 8) (list 3 5)) -1 1))) (let loop2 ((a a) (t t)) (if (even? a) (loop2 (/ a 2) (* t z)) (loop1 (modulo m a) a (if (and (= (modulo a 4) 3) (= (modulo m 4) 3)) (- t) t)))))))) (define (selfridge n) ; initialize lucas sequence (let loop ((d-abs 5) (sign 1)) (let ((d (* d-abs sign))) (cond ((< 1 (gcd d n)) (values d 0 0)) ((= (jacobi d n) -1) (values d 1 (/ (- 1 d) 4))) (else (loop (+ d-abs 2) (- sign))))))) (define (lucas p q m n) ; lucas sequences u[n] and v[n] and q^n (mod m) (define (even e o) (if (even? n) e o)) (define (mod n) (if (zero? m) n (modulo n m))) (let ((d (- (* p p) (* 4 q)))) (let loop ((un 1) (vn p) (qn q) (n (quotient n 2)) (u (even 0 1)) (v (even 2 p)) (k (even 1 q))) (if (zero? n) (values u v k) (let ((u2 (mod (* un vn))) (v2 (mod (- (* vn vn) (* 2 qn)))) (q2 (mod (* qn qn))) (n2 (quotient n 2))) (if (even? n) (loop u2 v2 q2 n2 u v k) (let* ((uu (+ (* u v2) (* u2 v))) (vv (+ (* v v2) (* d u u2))) (uu (if (and (positive? m) (odd? uu)) (+ uu m) uu)) (vv (if (and (positive? m) (odd? vv)) (+ vv m) vv)) (uu (mod (/ uu 2))) (vv (mod (/ vv 2)))) (loop u2 v2 q2 n2 uu vv (* k q2))))))))) (define (lucas-test? n ps) (call-with-values (lambda () (selfridge n)) (lambda (d pp qq) (let loop ((ps ps) (p pp) (q qq)) (if (null? ps) #t (call-with-values (lambda () (lucas p q n (+ n 1))) (lambda (un+1 vn+1 qn+1) (call-with-values (lambda () (lucas p q n (/ (+ n 1) (car ps)))) (lambda (un+1/p vn+1/p qn+1/p) (if (and (zero? un+1) (= (gcd un+1/p n) 1)) (loop (cdr ps) pp qq) (loop ps (+ p 2) (+ p q 1)))))))))))) (define-syntax while (syntax-rules () ((while pred? body ...) (do () ((not pred?)) body ...)))) (define (bls-prime? n) (define (last-pair xs) (if (null? (cdr xs)) xs (last-pair (cdr xs)))) (define (cycle xs) (set-cdr! (last-pair xs) xs) xs) (define wheel (list 2 4 2 4 6 2 6 4 2 4 6 6 2 6 4 2 6 4 6 8 4 2 4 2 4 8 6 4 6 2 4 6 2 6 6 4 2 4 6 2 6 4 2 4 2 10 2 10)) (define (rho n c b) (define (f y) (modulo (+ (* y y) c) n)) (define (g p x y) (modulo (* p (abs (- x y))) n)) (let stage1 ((x 2) (y (+ 4 c)) (z (+ 4 c)) (j 1) (q 2) (p 1)) (if (= j b) #f ; timeout (if (= x y) (rho n (+ c 1) (- b j)) ; cycle (if (= j q) (let ((t (f y))) (stage1 y (f y) z (+ j 1) (* q 2) (g p y t))) (if (positive? (modulo j 100)) (stage1 x (f y) z (+ j 1) q (g p x y)) (let ((d (gcd p n))) (if (= d 1) (stage1 x (f y) y (+ j 1) q (g p x y)) (if (and (< 1 d n) (bls-prime? d)) d ; factor (let stage2 ((k 1) (z (f z))) (if (= k 100) (rho n (+ c 1) (- b j)) (let ((d (gcd (- z x) n))) (if (= d 1) (stage2 (+ k 1) (f z)) (if (= d n) (rho n (+ c 1) (- b j)) (if (bls-prime? d) d ; factor (rho n (+ c 1) (- b j))))))))))))))))) (define (remove-twos n) (let loop ((f 1) (r n)) (if (odd? r) (values f (list 2) r) (loop (* f 2) (/ r 2))))) (define (join x xs) (if (member x xs) xs (cons x xs))) (define (enough? n p f1 f2) (< n (* (max (+ (* p f1) 1) (- (* p f2) 1)) (+ (* p p f1 f2 1/2) 1)))) (if (not (prime? n)) #f ; sanity check (let-values (((f1 f1s r1) (remove-twos (- n 1))) ((f2 f2s r2) (remove-twos (+ n 1)))) (let loop ((p 3) (ws (cons 2 (cons 2 (cons 4 (cycle wheel)))))) (cond ((and (< p 10000000) (not (enough? n p f1 f2))) (case (modulo (+ n 1) p) ((0) (while (zero? (modulo r2 p)) (set! f2 (* f2 p)) (set! f2s (join p f2s)) (set! r2 (/ r2 p)))) ((2) (while (zero? (modulo r1 p)) (set! f1 (* f1 p)) (set! f1s (join p f1s)) (set! r1 (/ r1 p))))) (when (< r1 (* p p)) (set! f1 (- n 1)) (set! f1s (join r1 f1s)) (set! r1 1)) (when (< r2 (* p p)) (set! f2 (+ n 1)) (set! f2s (join r2 f2s)) (set! r2 1)) (loop (+ p (car ws)) (cdr ws))) (else (set! p (min p 10000000)) (when verbose? (display (list 'wheel n p f1 f1s r1 f2 f2s r2)) (newline)) (let loop ((more1? #t) (more2? #t) (more3? #t)) (cond ((and (not (enough? n p f1 f2)) (or more1? more2?)) (if more1? (let ((f (rho r1 1 10000000))) (if (not f) (loop #f #t #t) (begin (set! f1 (* f1 f)) (set! f1s (join f f1s)) (set! r1 (/ r1 f)) (loop #t #t#t)))) (let ((f (rho r2 1 10000000))) (if (not f) (loop #f #f #t) (begin (set! f2 (* f2 f)) (set! f2s (join f f2s)) (set! r2 (/ r2 f)) (loop #f #t #t)))))) (more3? (when verbose? (display (list 'rho n p f1 f1s r1 f2 f2s r2)) (newline)) (loop #f #f #f)) ((not (enough? n p f1 f2)) #f) ; failure to prove ((< r1 f1) (when verbose? (display "fermat only") (newline)) (fermat-test? n f1s)) ; condition 1 only ((< r2 f2) (when verbose? (display "lucas only") (newline)) (lucas-test? n f2s)) ; condition 2 only (else (when verbose? (display "fermat and lucas") (newline)) (and (fermat-test? n (cons r1 f1s)) ; 1 and 3 (lucas-test? n (cons r2 f2s)))))))))))) ; 2&4 (define rand #f) (define randint #f) (let ((two31 #x80000000) (a (make-vector 56 -1)) (fptr #f)) (define (mod-diff x y) (modulo (- x y) two31)) ; generic version ; (define (mod-diff x y) (logand (- x y) #x7FFFFFFF)) ; fast version (define (flip-cycle) (do ((ii 1 (+ ii 1)) (jj 32 (+ jj 1))) ((< 55 jj)) (vector-set! a ii (mod-diff (vector-ref a ii) (vector-ref a jj)))) (do ((ii 25 (+ ii 1)) (jj 1 (+ jj 1))) ((< 55 ii)) (vector-set! a ii (mod-diff (vector-ref a ii) (vector-ref a jj)))) (set! fptr 54) (vector-ref a 55)) (define (init-rand seed) (let* ((seed (mod-diff seed 0)) (prev seed) (next 1)) (vector-set! a 55 prev) (do ((i 21 (modulo (+ i 21) 55))) ((zero? i)) (vector-set! a i next) (set! next (mod-diff prev next)) (set! seed (+ (quotient seed 2) (if (odd? seed) #x40000000 0))) (set! next (mod-diff next seed)) (set! prev (vector-ref a i))) (flip-cycle) (flip-cycle) (flip-cycle) (flip-cycle) (flip-cycle))) (define (next-rand) (if (negative? (vector-ref a fptr)) (flip-cycle) (let ((next (vector-ref a fptr))) (set! fptr (- fptr 1)) next))) (define (unif-rand m) (let ((t (- two31 (modulo two31 m)))) (let loop ((r (next-rand))) (if (<= t r) (loop (next-rand)) (modulo r m))))) (init-rand 19380110) ; happy birthday donald e knuth (set! rand (lambda seed (cond ((null? seed) (/ (next-rand) two31)) ((eq? (car seed) 'get) (cons fptr (vector->list a))) ((eq? (car seed) 'set) (set! fptr (caadr seed)) (set! a (list->vector (cdadr seed)))) (else (/ (init-rand (modulo (numerator (inexact->exact (car seed))) two31)) two31))))) (set! randint (lambda args (cond ((null? (cdr args)) (if (< (car args) two31) (unif-rand (car args)) (floor (* (next-rand) (car args))))) ((< (car args) (cadr args)) (let ((span (- (cadr args) (car args)))) (+ (car args) (if (< span two31) (unif-rand span) (floor (* (next-rand) span)))))) (else (let ((span (- (car args) (cadr args)))) (- (car args) (if (< span two31) (unif-rand span) (floor (* (next-rand) span)))))))))) (define (rand-prime n . args) (let ((b (if (pair? args) (car args) 10))) (let loop ((r (randint 1 b)) (n n)) (if (= n 1) (next-prime r) (loop (+ (* r b) (randint b)) (- n 1)))))) (define verbose? #t) (let ((n (rand-prime 42))) (display "proving primality of ") (display n) (newline) (display (bls-prime? n)) (newline))
Private
[
?
]
Run code
Submit