有限体をかじる中出てきたアルゴリズムをHaskell で実装してみました。
アルゴリズムおよびpseudo code は以下のwikipedia から
Extended Euclidean algorithm
ここでは基本的な集合論と代数構造の知識を仮定してます。
環$Z_n$
整数全体$\mathbb{Z}$ は環になりますが、更にある自然数$n$ を取ってきて
\mathbb{Z}/n\mathbb{Z}
なる商集合を考えてもこれまた環になります。
平たく言うと整数全体を$n$ で割った余りで同値類を作ってその同値関係で割ったものということですね。
\{0, \cdots, (n-1)\}
なる有限集合と考えてもよし、演算は全部
\mod n
で考えれば良いということになります。
体
ここで$p$ を素数とすると
\mathbb{Z}/p\mathbb{Z}
が体になることが示せます。
体というのは、環のうち、零でない元が掛け算に対する逆数を持つようなもののことです。
ここで$p$ が素数なことを思い出すと、
\{0, \cdots, (p-1)\}
の内$p$ と互いに素でないのは0だけで、それ以外は最大公約数が1 となることから分かります。
拡張されたユークリッドの互除法
ここで力を発揮するのがタイトルにもなったユークリッドの互除法の拡張版です、普通のユークリッドの互除法は最大公約数を求めるだけですが、このアルゴリズムだと積に対する逆元の候補も計算できるので嬉しいというわけです。
pseudo code
以下wikiからコピペです。
function extended_gcd(a, b) s := 0; old_s := 1 t := 1; old_t := 0 r := b; old_r := a while r ≠ 0 quotient := old_r div r (old_r, r) := (r, old_r - quotient * r) (old_s, s) := (s, old_s - quotient * s) (old_t, t) := (t, old_t - quotient * t) output "Bézout coefficients:", (old_s, old_t) output "greatest common divisor:", old_r output "quotients by the gcd:", (t, s)
整数a,b に対してこのアルゴリズムはそれらの最大公約数と以下の等号が成り立つ二整数x,y を返します：
\gcd(a,b) = a*x + b*y
従って例えば
\gcd(a,p) = a*x + b*p
を考えてみると、もしa とp が互いに素だと
1 = a*x + b*p
となります。
この右辺は$\mod p$ では$a*x (\mod p)$ となりこのx、正確には$x (\mod p)$ がa の逆数となるわけです。
Haskell での実装
なるほど、遅延リストを使えそうとまずあたりが付きます、すなわち例えばリストr は最初の要素２つが入力のa とb であとは相互に再帰的につけていったら良いという訳です。
アイディアとしては一個ずらしたリストをtail で作ってzipWith で必要な演算を挟んでくっつけていったという感じです。
終了条件をどうしようかなと迷ったので、素朴な高階関数を一つ定義して使っています。
こんないい加減で書き下しただけのコードでも動くのはうれしいですね。
> module Ffield where
> import Data.Ratio ((%))
> coprime :: Integral a => a -> a -> Bool
> coprime a b = (gcd a b) == 1
Consider a finite ring
Z_n := [0..(n-1)]
> haveInverse :: Integral a => a -> [Bool]
> haveInverse n = map (coprime n) [0..(n-1)]
If any non-zero element has its multiplication inverse, then the ring is a field:
> isField :: Integral a => a -> Bool
> isField n = and $ tail $ haveInverse n
If p is prime, then Z_p is a field:
zip [2..] $ map isField [2..13]
[(2,True),(3,True),(4,False),(5,True),(6,False),(7,True),(8,False),(9,False),(10,False),(11,True),(12,False),(13,True)]
Here we would like to implement the extended Euclidean algorithm.
See the algorithm, examples, and pseudo code at:
https://en.wikipedia.org/wiki/Extended_Euclidean_algorithm
> exGCD' :: Integral c => c -> c -> ([c], [c], [c], [c])
> exGCD' a b = (q, r, s, t)
> where
> q = zipWith quot r (tail r)
> r = takeBefore (==0) r'
> r' = a:b: (zipWith (-) r' (zipWith (*) q (tail r')))
> s = 1:0: (zipWith (-) s (zipWith (*) q (tail s)))
> t = 0:1: (zipWith (-) t (zipWith (*) q (tail t)))
>
> takeBefore :: (a -> Bool) -> [a] -> [a]
> takeBefore _ [] = []
> takeBefore p (l:ls)
> | p l = []
> | otherwise = l : (takeBefore p ls)
This example is from wikipedia:
*Ffield> exGCD' 240 46
([5,4,1,1,2],[240,46,10,6,4,2],[1,0,1,-4,5,-9,23],[0,1,-5,21,-26,47,-120])
*Ffield> gcd 240 46
2
*Ffield> 240*(-9) + 46*(47)
2
> -- a*x + b*y = gcd a b
> exGcd a b = (g, x, y)
> where
> (_,r,s,t) = exGCD' a b
> g = last r
> x = last . init $ s
> y = last . init $ t
*Ffield> exGcd 46 240
(2,47,-9)
*Ffield> 46*47 + 240*(-9)
2
*Ffield> gcd 46 240
2
解説
coprimeは互いに素かどうかの判定、組み込みのgcd を使ってます。
haveInverse は商環$\mathbb{Z}/n\mathbb{Z}$をリスト[0..(n-1)]で与えて、それらの元がn と互いに素かどうかで逆元が存在するかどうかを判定します。
isField は$\mathbb{Z}/n\mathbb{Z}$ のn を取って、その商環が体になるかどうかの判定です、素朴に0 以外がn と互いに素かどうかを見ています。
メインの前に補助関数としてexGCD'を作りました、これは拡張版ユークリッドの互除法のアルゴリズムに沿ってリストたちを返す関数です。
返り値のリストに$\gcd$ と欲しい二整数も入ってる寸法です。
where の中ですが、フィボナッチの時と同様に一つずらしてzipWith で演算とまとめるをバカの一つ覚えのように多用します。
あとは終了判定のための高階関数`takeBefore と遅延リストr', s,tです。
問題点
とりあえず動くは動きますがやや不満が残ります。
強Haskeller の皆様（凶でも狂でもこの際構いません）もしよろしければお力を貸してください、よろしくお願いします。
似たパターンで作られる遅延リスト
最初に思いついたのは高階関数を用意して使いまわせないかな？ということでした。
> exGCD2 :: Integral c => c -> c -> ([c], [c], [c], [c])
> exGCD2 a b = (q, r, s, t)
> where
> q = zipWith quot r (tail r)
> r = takeBefore (==0) r'
> r' = steps r' q a b
> s = steps s q 1 0
> t = steps t q 0 1
>
> -- The follwoing higher order function does NOT work.
> steps :: Integral a => [a] -> [a] -> a -> a -> [a]
> steps rr@(_:rs) qs a b = a:b: (zipWith (-) rr $ zipWith (*) qs rs)
一見良さそうですし、ghci に読ませても何も文句を言われないので型も合ってますが、走らせるとループになりますと怒られます。
*Ffield> exGCD2 13 7
(*** Exception: <<loop>>
*Ffield> exGCD' 13 7
([1,1,6],[13,7,6,1],[1,0,1,-1,7],[0,1,-1,2,-13])
恐らくこの高階関数の引数に結果として返したいリスト自身があてがわれてるせいで、必ずバインドされるわけですが、呼ばれた段階ではまだ中身が決まってないので、、、ということでしょうか？
fold で書き直せそうな高階関数takeBefore
最初takeUntil という名前でしたがどっちが良かったのでしょう？
さておき、なんとなくfilter に似てる気がするので畳み込めそうな気がします。
高階関数を使った追記
困ったときのfold(r) のチュートリアルを読んでそのまま書き直してみました、基本的にif then else を書きたくはないので同じ構造をそのまま持ってきましたが右から畳み込んでることを明示できるのでこちらの方が気分が楽な気がします。
> takeUntil :: (a -> Bool) -> [a] -> [a]
> takeUntil p = foldr func []
> where
> func x xs
> | p x = []
> | otherwise = x : xs
安全とは言えなんかいやなlast とか
パターンマッチでなんとか取り出せないかな？とも考えましたが、分からなかったので放置してます。
ここで安全というのは初期値が必ず二つ入ってるので、この操作は必ず何かの要素にヒットすることが保証されてる、ということです。
例$Z_{11}$
素数11 を使ってちょっと実験してみました。
ここで使ったのは新バージョン
> exGCD' :: Integral n => n -> n -> ([n], [n], [n], [n])
> exGCD' a b = (qs, rs, ss, ts)
> where
> qs = zipWith quot rs (tail rs)
> -- rs = takeBefore (==0) r'
> rs = takeUntil (==0) r'
> r' = steps a b
> ss = steps 1 0
> ts = steps 0 1
> steps a b = rr
> where rr@(_:rs) = a:b: zipWith (-) rr (zipWith (*) qs rs)
>
> takeUntil :: (a -> Bool) -> [a] -> [a]
> takeUntil p = foldr func []
> where
> func x xs
> | p x = []
> | otherwise = x : xs
さて以下が実験です、まずは$Z_{11}$ が体になるかどうか、すなわち素数かどうかです。
そのあと少々インタープリターの上でいじってみて、非零元の逆元を返すような関数を実装しています。
Example Z_{11}
*Ffield> isField 11
True
*Ffield> map (exGcd 11) [0..10]
[(11,1,0),(1,0,1),(1,1,-5),(1,-1,4),(1,-1,3),(1,1,-2),(1,-1,2),(1,2,-3),(1,3,-4),(1,-4,5),(1,1,-1)]
*Ffield> map ((`mod` 11) . (\(_,_,x)->x) . exGcd 11) [1..10]
[1,6,4,3,9,2,8,7,5,10]
*Ffield> zip [1..10] it
[(1,1),(2,6),(3,4),(4,3),(5,9),(6,2),(7,8),(8,7),(9,5),(10,10)]
> inverses :: Int -> Maybe [(Int, Int)]
> inverses n
> | isField n = Just lst
> | otherwise = Nothing
> where
> lst' = map ((`mod` n) . (\(_,_,c)->c) . exGcd n) [1..(n-1)]
> lst = zip [1..] lst'
*Ffield> inverses 11
Just [(1,1),(2,6),(3,4),(4,3),(5,9),(6,2),(7,8),(8,7),(9,5),(10,10)]