codepad
[
create a new paste
]
login
|
about
Language:
C
C++
D
Haskell
Lua
OCaml
PHP
Perl
Plain Text
Python
Ruby
Scheme
Tcl
(define (expm b e m) (define (times x y) (modulo (* x y) m)) (let loop ((b b) (e e) (r 1)) (if (zero? e) r (loop (times b b) (quotient e 2) (if (odd? e) (times b r) r))))) (define (iroot k n) (let loop ((u n)) (let* ((s u) (t (+ (* (- k 1) s) (quotient n (expt s (- k 1))))) (u (quotient t k))) (if (< u s) (loop u) s)))) (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 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 (iroot 2 n))) (if (= (* q q) n) q #f)))))))))))) (define (spsp? n a) (do ((d (- n 1) (/ d 2)) (s 0 (+ s 1))) ((odd? d) (let ((t (expm a d n))) (if (or (= t 1) (= t (- n 1))) #t (do ((s (- s 1) (- s 1)) (t (expm t 2 n) (expm t 2 n))) ((or (zero? s) (= t (- n 1))) (positive? s)))))))) (define (pseudoprime? n) (and (spsp? n 2) (spsp? n 3) (spsp? n 5))) (define (td n limit) (let twos ((n n) (fs (list))) (cond ((= n 1) fs) ((even? n) (twos (/ n 2) (cons 2 fs))) (else (let odds ((n n) (f 3) (fs fs)) (cond ((< limit f) fs) ((< n (* f f)) (reverse (cons n fs))) ((zero? (modulo n f)) (odds (/ n f) f (cons f fs))) (else (odds n (+ f 2) fs)))))))) (define (rho n limit) (let loop ((t 2) (h 2) (d 1) (c 1) (limit limit)) (define (f x) (modulo (+ (* x x) c) n)) (cond ((zero? limit) #f) ((= d 1) (let ((t (f t)) (h (f (f h)))) (loop t h (gcd (- t h) n) c (- limit 1)))) ((= d n) (loop 2 2 1 (+ c 1) (- limit 1))) ((pocklington? d) d) (else (rho d (- limit 1)))))) (define (pocklington? n) (if (not (pseudoprime? n)) #f (if (< n #e1e12) (= (car (td n #e1e6)) n) (let* ((n-1 (- n 1)) (root2 (iroot 2 n)) (root3 (iroot 3 n)) (fs (td n-1 1000)) (f (apply * fs))) (let f-loop ((f f) (fs fs)) (if (< f root3) (let ((p (rho (/ n-1 f) #e1e6))) (if p (f-loop (* f p) (cons p fs)) 'unknown)) (let a-loop ((a 2)) (if (not (= (expm a n-1 n) 1)) #f (let loop ((fs fs)) (let ((g (gcd (- (expm a (/ n-1 (car fs)) n) 1) n))) (cond ((< 1 g n) #f) ((= g n) (a-loop (+ a 1))) ((pair? (cdr fs)) (loop (cdr fs))) ((< root2 f) #t) (else (let* ((ds (digits n f)) (c2 (car ds)) (c1 (cadr ds)) (s (- (* c1 c1) (* 4 c2)))) (not (square? s))))))))))))))) (define (mersenne n) (- (expt 2 n) 1)) (display (pocklington? (mersenne 89))) (newline) (display (pocklington? (mersenne 521))) (newline) (display (pseudoprime? 3825123056546413051)) (newline) (display (pocklington? 3825123056546413051)) (newline)
Private
[
?
]
Run code
Submit