-- EEA
eea :: (Integral m) => m -> m -> (m, m, m)
eea f g = loop f g 1 0 0 1
  where
    loop r0 r1 s0 s1 t0 t1
      | r2 == 0   = if r1 < 0 then ((-r1), (-s1), (-t1)) else (r1, s1, t1)
      | otherwise = loop r1 r2 s1 s2 t1 t2
      where
        (q, r2) = r0 `divMod` r1
        s2      = s0 - s1 * q
        t2      = t0 - t1 * q

-- 途中経過も出力
eeaList :: (Integral m) => m -> m -> [(m, m, m)]
eeaList f g = loop f g 1 0 0 1
  where
    loop r0 r1 s0 s1 t0 t1
      | r2 == 0   = [(r0, s0, t0)] ++ case r1 < 0 of True -> [((-r1), (-s1), (-t1))]
                                                     _    -> [(r1, s1, t1)]
      | otherwise = [(r0, s0, t0)] ++ loop r1 r2 s1 s2 t1 t2
      where
        (q, r2) = r0 `divMod` r1
        s2      = s0 - s1 * q
        t2      = t0 - t1 * q

-- 余りの絶対値を小さくとる割り算
absDivMod :: (Integral m) => m -> m -> (m, m)
absDivMod a b
  | abs b >= abs (2 * r) = (q, r)
  | otherwise            = (q + 1, r - b)
  where
    (q, r) = a `divMod` b

-- EEA with absDivMod
absEea :: (Integral m) => m -> m -> (m, m, m)
absEea f g = loop f g 1 0 0 1
  where
    loop r0 r1 s0 s1 t0 t1
      | r2 == 0   = if r1 < 0 then ((-r1), (-s1), (-t1)) else (r1, s1, t1)
      | otherwise = loop r1 r2 s1 s2 t1 t2
      where
        (q, r2) = r0 `absDivMod` r1
        s2      = s0 - s1 * q
        t2      = t0 - t1 * q

absEeaList :: (Integral m) => m -> m -> [(m, m, m)]
absEeaList f g = loop f g 1 0 0 1
  where
    loop r0 r1 s0 s1 t0 t1
      | r2 == 0   = [(r0, s0, t0)] ++ case r1 < 0 of True -> [((-r1), (-s1), (-t1))]
                                                     _    -> [(r1, s1, t1)]
      | otherwise = [(r0, s0, t0)] ++ loop r1 r2 s1 s2 t1 t2
      where
        (q, r2) = r0 `absDivMod` r1
        s2      = s0 - s1 * q
        t2      = t0 - t1 * q

-- EEA の検算
checkEea :: (Integral m, Random m) => (m -> m -> (m, m, m)) -> Bool
checkEea f
  | False `elem` checkedeList = False
  | otherwise                 = True
  where
    checkedeList = check f <$> as <*> bs
    as = map (fst . random . mkStdGen) [0..99]
    bs = map (fst . random . mkStdGen) [0..(-99)]
    check f a b
      | s * a + t * b == gcd a b = True
      | otherwise                = False
      where
        (_, s, t) = f a b