; 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

1000
1000
inverse: divide by zero

 === context ===
Line 18:4: f
Line 16:0: crt