codepad
[
create a new paste
]
login
|
about
Language:
C
C++
D
Haskell
Lua
OCaml
PHP
Perl
Plain Text
Python
Ruby
Scheme
Tcl
; aks primality prover (define (split n xs) (let loop ((n n) (xs xs) (zs '())) (if (or (zero? n) (null? xs)) (values (reverse zs) xs) (loop (- n 1) (cdr xs) (cons (car xs) zs))))) (define (make-list n x) (let loop ((n n) (xs '())) (if (zero? n) xs (loop (- n 1) (cons x xs))))) (define (square x) (* x x)) (define (log2 n) (/ (log n) (log 2))) (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 (ilog b n) (let loop1 ((lo 0) (b^lo 1) (hi 1) (b^hi b)) (if (< b^hi n) (loop1 hi b^hi (* hi 2) (* b^hi b^hi)) (let loop2 ((lo lo) (b^lo b^lo) (hi hi) (b^hi b^hi)) (if (<= (- hi lo) 1) (if (= b^hi n) hi lo) (let* ((mid (quotient (+ lo hi) 2)) (b^mid (* b^lo (expt b (- mid lo))))) (cond ((< n b^mid) (loop2 lo b^lo mid b^mid)) ((< b^mid n) (loop2 mid b^mid hi b^hi)) (else mid)))))))) (define (td-prime? n) (if (even? n) (= n 2) (let loop ((d 3)) (cond ((< n (* d d)) #t) ((zero? (modulo n d)) #f) (else (loop (+ d 2))))))) (define (uniq-factors n) (define (uniq-cons x xs) (if (null? xs) (list x) (if (= x (car xs)) xs (cons x xs)))) (let twos ((n n) (fs '())) (if (even? n) (twos (/ n 2) (uniq-cons 2 fs)) (if (= n 1) fs (let odds ((n n) (d 3) (fs fs)) (cond ((< n (* d d)) (reverse (uniq-cons n fs))) ((zero? (modulo n d)) (odds (/ n d) d (uniq-cons d fs))) (else (odds n (+ d 2) fs)))))))) (define (prime-power? n) (if (even? n) (if (= (expt 2 (ilog 2 n)) n) 2 #f) (let loop ((a 2) (n n)) (let* ((b (expm a n n)) (p (gcd (- b a) n))) (cond ((or (= p 1) (= a b)) #f) ((prime? p) (if (= (expt p (ilog p n)) n) p #f)) (else (loop (+ a 1) p))))))) (define (ord r n) (do ((k 2 (+ k 1))) ((= (expm n k r) 1) k))) (define (compute-r n) (let ((target (* 4 (square (log2 n))))) (let loop ((r 3)) (if (not (= (gcd r n) 1)) (loop (+ r 1)) (if (< target (ord r n)) r (loop (+ r 1))))))) (define (phi n) (let loop ((fs (uniq-factors n)) (t n)) (if (null? fs) t (loop (cdr fs) (* t (- 1 (/ (car fs)))))))) (define (poly-mult-mod xs ys r n) (define (times x) (lambda (y) (* x y))) (define (plus xs ys) (let loop ((xs xs) (ys ys) (zs (list))) (cond ((null? xs) (reverse (append (reverse ys) zs))) ((null? ys) (reverse (append (reverse xs) zs))) (else (loop (cdr xs) (cdr ys) (cons (+ (car xs) (car ys)) zs)))))) (define (mod-poly xs) (let-values (((hs ts) (split r (reverse xs)))) (reverse (plus hs ts)))) (define (mod-n x) (modulo x n)) (let ((xs (reverse xs))) (let loop ((xs (cdr xs)) (zs (map (times (car xs)) ys))) (if (null? xs) (map mod-n (mod-poly zs)) (loop (cdr xs) (plus (cons 0 zs) (map (times (car xs)) ys))))))) (define (poly-power-mod bs e r n) (let loop ((bs bs) (e e) (rs (list 1))) (if (zero? e) rs (loop (poly-mult-mod bs bs r n) (quotient e 2) (if (even? e) rs (poly-mult-mod rs bs r n)))))) (define (binomial-test? a r n) (not (equal? (poly-power-mod (list 1 a) n r n) (append (list 1) (make-list (- r 1) 0) (list a))))) (define (aks-prime? n) (if (prime-power? n) #f (let* ((r (compute-r n)) (phi-r (phi r)) (sqrt-phi-r (sqrt phi-r)) (log2-n (log2 n)) (sqrt-phi-r-log2-n (* sqrt-phi-r log2-n))) (let loop ((a 1)) (if (<= a r) (if (< 1 (gcd a n) n) #f (loop (+ a 1))) (if (<= n r) #t (let loop ((a 1)) (if (<= a sqrt-phi-r-log2-n) (if (binomial-test? a r n) (loop (+ a 1)) #f) #t)))))))) (display (compute-r 89)) (newline) (display (aks-prime? 89)) (newline)
Private
[
?
]
Run code
Submit