ChineseRemainderTheorem.hs

{-# LANGUAGE BangPatterns #-}

import Control.Exception

-- crt [(3,2), (5,3), (7,2)] == 23
crt :: [(Integer, Integer)] -> Integer
crt xs = assert (all (\(n, a) -> ret `mod` n == a) xs) $ ret
  where
    ret = foldl add 0 [m' `mul` a `mul` t | (n, a) <- xs, let m' = m `div` n, let (d, t, _) = exgcd m' (- n), assert (d == 1) True]
    m = product [n | (n, a) <- xs]
    add x y = (x + y) `mod` m
    mul x y = (x * y) `mod` m

crt2 :: (Integer, Integer) -> (Integer, Integer) -> Integer
crt2 (m, a) (n, b) =
  case exgcd m n of
    (d, u, v) ->
      let ret = a `mul` n `mul` v + b `mul` m `mul` u
       in assert (d == 1 && ret `mod` m == a && ret `mod` n == b) ret
  where
    mn = m * n
    add x y = (x + y) `mod` mn
    mul x y = (x * y) `mod` mn

exgcd :: Integral a => a -> a -> (a, a, a)
exgcd f1 f2 = f $ go f1 f2 1 0 0 1
  where
    go !r0 !r1 !s0 !s1 !t0 !t1
      | r1 == 0   = (r0, s0, t0)
      | otherwise = go r1 r2 s1 s2 t1 t2
      where
        (q, r2) = r0 `divMod` r1
        s2 = s0 - q*s1
        t2 = t0 - q*t1
    f (g,u,v)
      | g < 0 = (-g, -u, -v)
      | otherwise = (g,u,v)