; 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)