```1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 ``` ```; 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) ```
 ```1 2 ``` ```10 4711 ```