codepad
[
create a new paste
]
login
|
about
Language:
C
C++
D
Haskell
Lua
OCaml
PHP
Perl
Plain Text
Python
Ruby
Scheme
Tcl
; modular square root with non-prime modulus (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? n) ; spsp(2), spsp(3), and lucas pseudoprime (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 (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 (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))))) (define (square? n) (let ((n2 (isqrt n))) (= n (* n2 n2)))) (define (jacobi a n) (if (not (and (integer? a) (integer? n) (positive? n) (odd? n))) (error 'jacobi "modulus must be positive odd integer") (let jacobi ((a a) (n n)) (cond ((= a 0) 0) ((= a 1) 1) ((= a 2) (case (modulo n 8) ((1 7) 1) ((3 5) -1))) ((even? a) (* (jacobi 2 n) (jacobi (quotient a 2) n))) ((< n a) (jacobi (modulo a n) n)) ((and (= (modulo a 4) 3) (= (modulo n 4) 3)) (- (jacobi n a))) (else (jacobi n a)))))) (define legendre jacobi) (define (inverse x n) (let loop ((x (modulo x n)) (a 1)) (cond ((zero? x) (error 'inverse "division by zero")) ((= x 1) a) (else (let ((q (- (quotient n x)))) (loop (+ n (* q x)) (modulo (* q a) n))))))) (define (miller? n a) (let loop ((r 0) (s (- n 1))) (if (even? s) (loop (+ r 1) (/ s 2)) (if (= (expm a s n) 1) #t (let loop ((r r) (s s)) (cond ((zero? r) #f) ((= (expm a s n) (- n 1)) #t) (else (loop (- r 1) (* s 2))))))))) (define (chain m f g x0 x1) (let loop ((ms (digits m 2)) (u x0) (v x1)) (cond ((null? ms) (values u v)) ((zero? (car ms)) (loop (cdr ms) (f u) (g u v))) (else (loop (cdr ms) (g u v) (f v)))))) (define (lucas? n) (let loop ((a 11) (b 7)) (let ((d (- (* a a) (* 4 b)))) (cond ((square? d) (loop (+ a 2) (+ b 1))) ((not (= (gcd n (* 2 a b d)) 1)) (loop (+ a 2) (+ b 2))) (else (let* ((x1 (modulo (- (* a a (inverse b n)) 2) n)) (m (quotient (- n (legendre d n)) 2)) (f (lambda (u) (modulo (- (* u u) 2) n))) (g (lambda (u v) (modulo (- (* u v) x1) n)))) (let-values (((xm xm1) (chain m f g 2 x1))) (zero? (modulo (- (* x1 xm) (* 2 xm1)) n))))))))) (cond ((or (not (integer? n)) (< n 2)) (error 'prime? "must be integer greater than one")) ((even? n) (= n 2)) ((zero? (modulo n 3)) (= n 3)) (else (and (miller? n 2) (miller? n 3) (lucas? n))))) (define (euclid x y) (let loop ((a 1) (b 0) (g x) (u 0) (v 1) (w y)) (if (zero? w) (values a b g) (let ((q (quotient g w))) (loop u v w (- a (* q u)) (- b (* q v)) (- g (* q w))))))) (define (inverse x m) (if (not (= (gcd x m) 1)) (error 'inverse "divisor must be coprime to modulus") (call-with-values (lambda () (euclid x m)) (lambda (a b g) (modulo a m))))) (define (jacobi a n) (if (not (and (integer? a) (integer? n) (positive? n) (odd? n))) (error 'jacobi "modulus must be positive odd integer") (let jacobi ((a a) (n n)) (cond ((= a 0) 0) ((= a 1) 1) ((= a 2) (case (modulo n 8) ((1 7) 1) ((3 5) -1))) ((even? a) (* (jacobi 2 n) (jacobi (quotient a 2) n))) ((< n a) (jacobi (modulo a n) n)) ((and (= (modulo a 4) 3) (= (modulo n 4) 3)) (- (jacobi n a))) (else (jacobi n a)))))) (define (mod-sqrt a p) (define (both n) (list n (- p n))) (cond ((not (and (odd? p) (prime? p))) (error 'mod-sqrt "modulus must be an odd prime")) ((not (= (jacobi a p) 1)) (error 'mod-sqrt "must be a quadratic residual")) (else (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))))))))) > (define (lift-root n p) (let* ((r (apply min (mod-sqrt n p))) (2r (+ r r)) (t (modulo (- (inverse 2r p)) p))) (+ r (* t p)))) (define (lift-root n p) (let* ((r (apply min (mod-sqrt n p))) (t (modulo (- (inverse (+ r r) p)) p))) (+ r (* t p)))) (display (lift-root 2 7)) (newline) (display (lift-root 13290059 127)) (newline)
Private
[
?
]
Run code
Submit