codepad
[
create a new paste
]
login
|
about
Language:
C
C++
D
Haskell
Lua
OCaml
PHP
Perl
Plain Text
Python
Ruby
Scheme
Tcl
; tonelli-shanks algorithm (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 prime? (let ((seed 3141592654)) (lambda (n) (define (rand) (set! seed (modulo (+ (* 69069 seed) 1234567) 4294967296)) (+ (quotient (* seed (- n 2)) 4294967296) 2)) (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 (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)))))))) (if (not (integer? n)) (error 'prime? "must be integer") (if (< n 2) #f (do ((a (rand) (rand)) (k 25 (- k 1))) ((or (zero? k) (not (spsp? n a))) (zero? k)))))))) (define (tonelli a p) (define (expp b e) (expm b e p)) (define (modp n) (modulo n p)) (let* ((p-1 (- p 1)) (p-1/2 (/ p-1 2))) (if (even? p) (error 'tonelli "must be odd")) (if (not (prime? p)) (error 'tonelli "must be prime")) (if (< p 3) (error 'tonelli "must be greater than two")) (if (< 1 (gcd a p)) (error 'tonelli "must be coprime")) (if (= (expp a p-1/2) -1) (error 'tonelli "must be quadratic residue")) (let loop ((s p-1) (e 0)) (if (even? s) (loop (/ s 2) (+ e 1)) (let loop ((n 2)) (if (not (= (expp n p-1/2) p-1)) (loop (+ n 1)) (let loop1 ((x (expp a (/ (+ s 1) 2))) (b (expp a s)) (g (expp n s)) (r e)) (let loop2 ((m 0) (mm 1)) (if (not (= (expp b mm) 1)) (loop2 (+ m 1) (* mm 2)) (if (zero? m) (list x (- p x)) ; found result (loop1 (modp (* x (expp g (expt 2 (- r m 1))))) (modp (* b (expp g (expt 2 (- r m))))) (expp g (expt 2 (- r m))) m))))))))))) (display (tonelli 2 113)) (newline) (display (tonelli 5 40961)) (newline) (display (tonelli 2 360027784083079948259017962255826129)) (newline)
Private
[
?
]
Run code
Submit