module ModArithmetic where

modmult :: Integer -> Integer -> Integer -> Integer
modmult m x y = (x * y) `mod` m

modadd :: Integer -> Integer -> Integer -> Integer
modadd m x y = (x + y) `mod` m  

modexp :: Integer -> Integer -> Integer -> Integer
modexp m a x = case compare x 1 of
				LT -> if x == 0 then 1 else modexp m a (x `mod` m)
				EQ -> a `mod` m
				GT ->  if even x
						then let p = (modexp m a (x `div` 2)) `mod` m
							 in  (p^2) `mod` m 
						else (a * modexp m a (x-1)) `mod` m

modinv :: Integer -> Integer -> Integer
modinv m a
		| g /= 1    = error "No inverse exists"
	    | otherwise = x `mod` m 
           		where (g,x,_) = gcd_e a m

gcd_e :: Integer -> Integer -> (Integer,Integer,Integer)
gcd_e a b = if d < 0 then (-d,-x,-y) else (d,x,y)
             	where (d,x,y) = euclidean_e (a,1,0) (b,0,1)

-- euclidean_e implements the extended euclidean algorithm
-- usually used to find the gcd
euclidean_e :: (Integral t) => (t, t, t) -> (t, t, t) -> (t, t, t) 
euclidean_e (r1, x1, y1) (r2, x2, y2)
              | (r2 == 0) = (r1, x1, y1) 
              | otherwise = euclidean_e (r2, x2, y2) (r, x1 - (q * x2), y1 - (q * y2))
           		where (q,r) = divMod r1 r2

-- chineseRemainder solves a set of simultaneous congruences if possible
chineseRemainder :: [Integer] -> [Integer] -> Integer
chineseRemainder xs ms = let parameters = zip3 xs ms (map (div m) ms)
						 in  foldl1 (modadd m) [m_i' * (modmult m_i x_i (modinv m_i m_i')) | (x_i, m_i, m_i') <- parameters ]
							where m = product ms;

-- Alternative formulation of the chinese remainder theorem
--
-- chineseRemainder' xs ms = fst (foldl1 chineseRemainder2 (zip xs ms))
-- chineseRemainder2 (x1,m1) (x2,m2) = ((x1 + (x2 - x1) * (m1 * m1')) `mod` m1m2, m1m2)
--		where { m1m2 = m1 * m2;
--				m1'  = modinv m2 m1; }