; 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

#t
5
7
9
3
10
10
(7 6)