Haksellで加速互除法・拡張ユークリッドの互除法 - これは圏です（はてな使ったら負けだとおもっていた）

……を、書いてみました。
ユークリッドの互除法といえば、a, b 二数が与えられたときに、その最大公約数を求めるアルゴリズムです。
世界で一番古いアルゴリズムとか云われていたような。次のような手順で求めます。

a > b > 0 としたとき、
(1) b == 0 なら a を返して終了。
(2) a を b で割った余り(剰余)をcとする。
(3) a = b
(4) b = c 
(5) また(1)から繰り返し

素因数分解するより遙かに手っ取り早いのだそうです。そりゃそうだわな。

で、加速互除法とはなんぞやと。素数入門―計算しながら理解できる (ブルーバックス)に載っていた、ユークリッドの互除法の改良版です。
どんなのかってえますと、普通、剰余は最も０に近い正の数ですが、正の数に限らず負の数でも良いじゃないかと。つまり、

(4) b = max( c, c-b )

とするわけです。こうすると、a, bの減り幅が増えて、少ない回数で求められるようになります。

で、拡張ユークリッドの互除法ってのは、最大公約数 r を求めるついでに、ax + by = r って成るようなa, bも求めちゃおうって云う、ちょっと欲張りなアルゴリズムです。
……説明はめんどくさいので、まあコード参照。

では、ソースです。

-- 剰余の絶対値が少ない方を取るmodとdiv
div1 a b  =  ((mod1 a b) - a) `div` b
mod1 a b  = if (abs m1) > (abs m2) then m2 else m1
            where m1 = a `mod` b
                  m2 = m1 - b

-- 加速互除法
euclid a 0 = (abs a)
euclid a b = euclid b (mod1 a b)

-- 拡張ユークリッドの互除法
euclidEx a b  = f a b 1 0 0 1
  where f r 0  a  _  b  _ = ((a, b), r)
        f x y a0 a1 b0 b1 = f y r a1 (a0-q*a1) b1 (b0-q*b1)
          where (q, r) = divMod x y

--回数計測用 (最大公約数, 呼び出し回数)を返す
gcd' x y = f 0 x y
  where f n x 0 = ((abs x), n+1)
        f n x y = f (n+1) y (x `mod` y)

euclid' x y = f 0 x y
  where f n x 0 = ((abs x), n+1)
        f n x y = f (n+1) y (x `mod1` y)

実行結果

Main> gcd 123456 123456543
3
Main> euclid 123456 123456543
3

Main> gcd' 123456 123456543
(3,11)
Main> euclid' 123456 123456543
(3,8)

Main> euclidEx 123456 123456543
((8867039,-8867),3)
Main> 8867039 * 123456 + (-8867) * 123456543
3

考えていること

拡張ユークリッドの互除法と加速互除法を合体できないかなあ……
説明が長すぎるなあ……