CRT

考えるの面倒臭くなってmapとapplyで全て済ましたというメモ。

(define (ExGCD x y)
  ;;;Input: x, y in Z
  ;;;Output: a and b s.t. ax + by = gcd(a,b)
  (let loop ((r0 x) (r1 y) (a0 1) (a1 0) (b0 0) (b1 1))
    (if (= r1 0)
        (values r0 a0 b0)
      (receive (q r) (quotient&remainder r0 r1)
               (loop
                 r1 r
                 a1 (- a0 (* q a1))
                 b1 (- b0 (* q b1)))))))

(define (inverse-p x p)
  ;;;Input: x in Z
  ;;;Output: x^{-1} mod p
  ;;;find a s.t. ax + bp = 1 and return a mod p.
  (receive (r a b) (ExGCD x p)
           (modulo a p)))

(define (CRT lis)
  ;;;INPUT: ((b1 . d1) (b2 . d2) ...)
  ;;;Output: N s.t. N \equiv b_i mod d_i
  (let1 N (apply * (map cdr lis))
        (modulo 
         (apply +
                (map (lambda (c)
                       (apply *
                              (list (car c)
                                    (/ N (cdr c))
                                    (inverse-p (/ N (cdr c)) (cdr c)))))
                     lis))
         N)))