Project:
 ```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 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 ``` ```; modular arithmetic (define rand (let* ((a 3141592653) (c 2718281829) (m (expt 2 35)) (x 5772156649) (next (lambda () (let ((x-prime (modulo (+ (* a x) c) m))) (set! x x-prime) x-prime))) (k 103) (v (list->vector (reverse (let loop ((i k) (vs (list x))) (if (= i 1) vs (loop (- i 1) (cons (next) vs))))))) (y (next)) (init (lambda (s) (set! x s) (vector-set! v 0 x) (do ((i 1 (+ i 1))) ((= i k)) (vector-set! v i (next)))))) (lambda seed (cond ((null? seed) (let* ((j (quotient (* k y) m)) (q (vector-ref v j))) (set! y q) (vector-set! v j (next)) (/ y m))) ((eq? (car seed) 'get) (list a c m x y k v)) ((eq? (car seed) 'set) (let ((state (cadr seed))) (set! a (list-ref state 0)) (set! c (list-ref state 1)) (set! m (list-ref state 2)) (set! x (list-ref state 3)) (set! y (list-ref state 4)) (set! k (list-ref state 5)) (set! v (list-ref state 6)))) (else (init (modulo (numerator (inexact->exact (car seed))) m)) (rand)))))) (define (randint . args) (cond ((null? (cdr args)) (inexact->exact (floor (* (rand) (car args))))) ((< (car args) (cadr args)) (+ (inexact->exact (floor (* (rand) (- (cadr args) (car args))))) (car args))) (else (+ (inexact->exact (ceiling (* (rand) (- (cadr args) (car args))))) (car args))))) (define (check? a n) (let loop ((r 0) (s (- n 1))) (if (even? s) (loop (+ r 1) (/ s 2)) (if (= (expm a s n) 1) #t (let loop ((j 0) (s s)) (cond ((= j r) #f) ((= (expm a s n) (- n 1)) #t) (else (loop (+ j 1) (* s 2))))))))) (define (prime? n) (cond ((< n 2) #f) ((= n 2) #t) ((even? n) #f) (else (let loop ((k 50)) (cond ((zero? k) #t) ((not (check? (randint 1 n) n)) #f) (else (loop (- k 1)))))))) (define (grand-daddy m n) (cond ((< m n) (grand-daddy m (- n m))) ((< n m) (grand-daddy (- m n) n)) (else m))) (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 (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 (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-syntax (with-modulus stx) (syntax-case stx () ((with-modulus e expr ...) (with-syntax ((modulus (datum->syntax-object (syntax with-modulus) 'modulus)) (== (datum->syntax-object (syntax with-modulus) '== )) (+ (datum->syntax-object (syntax with-modulus) '+ )) (- (datum->syntax-object (syntax with-modulus) '- )) (* (datum->syntax-object (syntax with-modulus) '* )) (/ (datum->syntax-object (syntax with-modulus) '/ )) (^ (datum->syntax-object (syntax with-modulus) '^ )) (sqrt (datum->syntax-object (syntax with-modulus) 'sqrt ))) (syntax (letrec ((fold (lambda (op base xs) (if (null? xs) base (fold op (op base (car xs)) (cdr xs)))))) (let* ((modulus e) (mod (lambda (x) (if (not (integer? x)) (error 'with-modulus "all arguments must be integers") (modulo x modulus)))) (== (lambda (x y) (= (mod x) (mod y)))) (+ (lambda xs (fold (lambda (x y) (mod (+ x (mod y)))) 0 xs))) (- (lambda (x . xs) (if (null? xs) (mod (- 0 x)) (fold (lambda (x y) (mod (- x (mod y)))) x xs)))) (* (lambda xs (fold (lambda (x y) (mod (* x (mod y)))) 1 xs))) (/ (lambda (x . xs) (if (null? xs) (inverse x e) (fold (lambda (x y) (* x (inverse y e))) x xs)))) (^ (lambda (base exp) (expm base exp e))) (sqrt (lambda (x) (mod-sqrt x e)))) expr ...))))))) (with-modulus 12 (display (== 17 5)) (newline) ; #t (display (+ 8 9)) (newline) ; 5 (display (- 4 9)) (newline) ; 7 (display (* 3 7)) (newline) ; 9 (display (/ 9 7)) (newline)) ; 3 (with-modulus 13 (display (^ 6 2)) (newline) ; 10 (display (^ 7 2)) (newline) ; 10 (display (sqrt 10)) (newline)) ; ±6 ```
 ```1 2 3 4 5 6 7 8 ``` ```#t 5 7 9 3 10 10 (7 6) ```