; chinese remainder theorem
(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 (crt as ms)
(let ((p (apply * ms)))
(define (f a m)
(let* ((t (/ p m)) (b (inverse t m)))
(* a b t)))
(modulo (apply + (map f as ms)) p)))
(display (crt '(10 4 12) '(11 12 13))) (newline)
(define (inverse x n)
(let loop ((x (modulo x n)) (a 1))
(cond ((zero? x) (error 'inverse "divide by zero"))
((= x 1) a)
(else (let ((q (- (quotient n x))))
(loop (+ n (* q x)) (modulo (* q a) n)))))))
(display (crt '(10 4 12) '(11 12 13))) (newline)
(display (crt '(4 7 3) '(5 8 9))) (newline) ; should be 39, but fails