yasuabe blog

2012年6月13日水曜日

ピタゴラス数の小ネタ

今までプログラマやってきて特に気にすることもなかったけど、project euler を始めてみるとピタゴラス数に関連する問題に出会うことがある。で、いろいろ調べてみると、ちょっと面白いネタというか Tips が見つかった。

■ 原始ピタゴラス数の無限リスト

このサイトで原始ピタゴラス数（primitive Pythagorean triple (PPT)）を次々に生成する行列が紹介されている。こんな行列。

（ -1　 -2　 2）
（ -2 　-1　 2）
（ -2　 -2　 3）

この行列に符号を3通りに変えたピタゴラス数を掛けると、子のピタゴラス数が3つ生成される。例えば(3,4,5)から始めると以下のようになる。

(-1 -2 2) (-3)   ( 5)
(-2 -1 2) ( 4) = (12) 
(-2 -2 3) ( 5)   (13)

(-1 -2 2) ( 3)   (15)
(-2 -1 2) (-4) = ( 8) 
(-2 -2 3) ( 5)   (17)

(-1 -2 2) (-3)   (21)
(-2 -1 2) (-4) = (20) 
(-2 -2 3) ( 5)   (29)

算出された子のピタゴラス数から、さらに３つずつの孫を作ることができ、これを再帰的に繰り返すとピタゴラス数の無限リストを作ることができる。Haskell で書くとこんな感じになる。

pyt1 (a,b,c) = ((-1)*a+(-2)*b+2*c,(-2)*a+(-1)*b+2*c,(-2)*a+(-2)*b+3*c)
pyt2 (a,b,c) = [pyt1(-a,b,c), pyt1(-a,-b,c), pyt1(a,-b,c)]
pyt3 ts = ts ++ pyt3 (concatMap pyt2 ts)

以下のような結果になる。

ghci< pyt3 [(3,4,5)]
[(3,4,5),(5,12,13),(21,20,29),(15,8,17),(7,24,25),(55,48,73),(45,28,53),(39,80,89),(119,120,169),(77,36,85),(33,56,65),(65,72,97),(35,12,37),(9,40,41),(105,88,137),(91,60,109),(105,208,233),(297,304,425),(187,84,205),(95,168,193),(207,224,305),(117,44,125),(57,176,185),(377,336,505),(299,180,349),(217,456,505),(697,696,985),(459,220,509),(175,288,337),(319,360,481),(165,52,173),(51,140,149),(275,252,373),(209,120,241),(115,252,277),(403,396,565),(273,136,305),(85,132,157),(133,156,205),(63,16,65),(11,60,61),(171,140,221),(153,104,185),(203,396,445),(555,572,797),・・・

■ 高さ＝底辺の二等辺三角形に漸近する漸化式

このサイトでは、ピタゴラス数を直角二等辺三角形に近づけていく漸化式が紹介されていた。直角を構成する辺 a と b の差を1にこんな式。

a_n+2 = 6*a_n+1 - a_n + 2  a₁ = 3, a₂ = 20
c_n+2 = 6*b_n+1 - b_n - 2  b₁ = 4, b₂ = 21
c_n+2 = 6*b_n+1 - c_n      c₁ = 5, c₂ = 29

で、ここから類推して高さ＝底辺の二等辺三角形（の片側の三角形）に漸近する漸化式もきっとあるだろうなと思って探してみた（きっかけは euler 138）。

まず直角を構成する2辺のうち長い辺と短い方を二倍したものとの差が1であるものを何個か探す。これは上述の無限リストに対して、条件に合うものをフィルタリングした。以下のようなものが見つかった（これ以上は数が大きくなりすぎて難しい）。

(15,8,17)
(273,136,305)
(4895,2448,5473)
(87841,43920,98209)
(1576239,788120,1762289)
(28284465,14142232,31622993)
(507544127,253772064,567451585)
(9107509825,4553754912,10182505537)

で見つかったピタゴラス数から連立方程式を立ててこれを解いた。手で計算するのは面倒なので、ネット上の連立方程式を解くプログラムを使った。例えば、斜辺ならこんな式を解けば良い。

305*x + 17y + z = 5473
5473*x + 305*y + z= 98209
98209*x + 5473*y + z =1762289

解）x=18, y=-1, z=0

で、c₁ = 17, c₂ = 305, c_n+2 = 18*c_n+1 - c_n という漸化式が得られる。Haskell で書くとこうなる。

l = 17:305:zipWith (\a b->18*b-a) l (tail l)

同様に、直角を構成する2辺についても漸化式が得られる。 a1 = 15, a2 = 273, a_n+3 = 17*a_n+2 + 17*a_n+1 - a_n (bも同様)

試しに20番目の要素を見てみると、a= 10151021471800938910964641, b = 5075510735900469455482320, c = 11349187026003431978487841)で、a² + b² = c², 1 = |2b-a| が確かに成立する。

==== ==== ==== ====

さっき problem 138 を解いて、euler project を150問解いた事になる。簡単な問題を選んでやってきたが、そろそろ難しくなってきた。かなり頑張ったけど、まだ半分にも至ってないのか…

2012年5月2日水曜日

Pell方程式を解くための実装メモ

Project Euler をやってると、Pell方程式 X²-DY²=1 に関連する問題をたまに見かける。これを解くとき、連分数を使うやり方がよく使われるようだけど、半分位の手数で解が得られるやり方がここで紹介されていたので、Haskell で実装してみた。

■ やり方

だいたいこんな感じ

g_-1=0	g₀=0	g_n+1 = -g_n + k_nh_n
h_-1=0	h₀=1	h_n+1 = (D - g_n+1²)/h_n
k_-1=0	k₀=[√D]	k_n+1 = [(k₀ + g_n)/h_n]
x_-1=0	x₀=1	x_n+1 = g_ny_n + h_ny_n-1
y_-1=0	y₀=0	y_n+1 = y_n-1 + k_ny_n

で、g_n とg_n+1 で同じ値が連続したら、

X = (x_n² + Dy_n²)/h_n

Y = 2x_ny_n/h_n

h_n とh_n+1 で同じ値が連続したら、

X' = (x_nx_n+1 + Dy_ny_n+1)/h_n

Y' = (x_ny_n+1 + x_n+1y_n)/h_n

X = X'² + DY'²

Y = 2X'Y'

■ 例

D = 23の場合

n	g	h	k	x	y
-1	0	0	0	0	1
0	0	1	4	1	0
1	4	7	1	4	1
2	3	2	3	5	1
3	3	7	1	19	4

24 = (5² + 23*1²)/2

5 = 2*5*1/2

D = 29の場合

n	g	h	k	x	y
-1	0	0	0	0	1
0	0	1	5	1	0
1	5	4	2	5	1
2	3	5	1	11	2
3	2	5	1	16	13

70 = (11*16 + 29*2*3)/5

13 = (11*3 + 16*2)/5

9801 = 70² + 29*13²

1820 = 2*70*13

■ コード

pell d=pell' d (0,0,0,0,1) (0,1,floor(sqrt(fromIntegral d)),1,0) 
  where 
   pell' d (g,h,k,x,y) n@(gn,hn,kn,xn,yn)
    | g'== gn    = let xr = div (xn^2 + d*yn^2) hn
                       yr = div (2*xn*yn) hn
                   in (xr, yr)
    | h'== hn    = let xr = div (x'*xn + d*y'*yn) hn
                       yr = div (x'*yn + xn*y') hn
                   in (xr^2 + d*yr^2, 2*xr*yr)
    | otherwise  = pell' d n (g',h',k',x',y')
    where k0 = floor$sqrt$fromIntegral d
          g' = (-gn) + kn*hn
          h' = div (d - g'^2) hn
          k' = div (k0 + g') h'
          y' = y + kn*yn
          x' = g'*y' + h'*yn

■ 結果

ghci> let isNotSquare n= (round . sqrt $ fromIntegral n) ^ 2 /= n
ghci> mapM_ (\n->print$(n,pell n)) $filter isNotSquare [2..15]
(2,(3,2))
(3,(2,1))
(5,(9,4))
(6,(5,2))
(7,(8,3))
(8,(3,1))
(10,(19,6))
(11,(10,3))
(12,(7,2))
(13,(649,180))
(14,(15,4))
(15,(4,1))
ghci> maximumBy (\a b->(on compare (snd.snd))a b )$ map (\n->(n,pell n)) $filter isNotSquare [2..2000]
(1621,(6298101812493732343034974500091457815529942308667051412857352310169665125001,156429324369979112128445583345098338627552043874824108399177922442751050500))

Problem 66 で、D≦1000でXが最大になるものを求める問題があるけど、D≦10000でもすぐに答えが出てくる。

2012年3月27日火曜日

ブルートフォース癖

ガウスがちっちゃい頃の逸話で、先生が算数の時間に、生徒達に１から100までの合計を計算させて、その間に雑用か何かを片付けようとしてたら、ガウスが101 * 50 = 5050と瞬時に答えをだしたので先生が驚いたってのを聞いたことがある。

プログラマたるもの大いに見習わなきゃならんなと常々思っていたけど、最近、Project Euler を始めてみて、かなり自分にブルートフォース癖があるのに気づいて反省。

例えば Problem 169 の以下のような問題がある（それほど難しくないと言われている）。

Define f(0)=1 and f(n) to be the number of different ways n can be expressed as a sum of integer powers of 2 using each power no more than twice.

For example, f(10)=5 since there are five different ways to express 10:

    1 + 1 + 8
    1 + 1 + 4 + 4
    1 + 1 + 2 + 2 + 4
    2 + 4 + 4
    2 + 8

What is f(10^25)?

これを解くのに、可能な手をそれぞれ試行してからルールに合わないものをふるい落とすなんてやり方（試行回数を減らす工夫はしてるつもりだけど）だと、いつまで経っても終わらない。

途中の計算結果をキャッシュでもしてみようかとも思ったけど、もうそんなのは止めておこうと思い至り、初めから再考する事にした。

で、改めてよく考えてみると、「1」を使わない部分の合計を x としたとき、n が奇数ならば必ず「x + 1」の形になる。ここで x 部分のバリーション数は f (n div 2) で、これが即ち f(n)と等しいということが分かる。

次に n が偶数のときを考えると、表記法のバリエーションは「x」か「x + 1 + 1」のいずれかのパターンに分類され、前者の表記法の数は f (n div 2) で、後者は f ((n-2) div 2)となるから、f (n) = f (n div 2) + f ((n-2) div 2)となる。

ここまでを問題文にある f(0)=1 を考慮してまとめると以下の様になる。

n = 0     ⇒ f(n) = 1
n ∈ odd  ⇒ f(n) = f(n')
n ∈ even ⇒ f(n) = f(n') + f(n'-1)
※ n' = n div 2

うん、本当は、この辺りまで公式化してから、さて実装はどうしようかと考えるべきなんだよなあ。

改めて「分析」ってやっぱ大事だなと思う。実際の開発現場では、重量プロセスで必死になって大量にドキュメントを書いてるくせして、「分析」が全然なされてなされてなかったりするし、下手糞なアジャイルもどきの開発でも、「後でなんぼでもリファクタできるから」的なノリで実装を進めて、変な方向にはまり込んで引き返せなくなったりもする。ちゃんと分析しようやと。

話を戻すと、実は上の公式をそのまんま実装しても、二重再帰が含まれてるのでものすごく遅い。だから、ここからが実装の話になるけど、これを防ぐには、一回の f(n) の計算で f(n) とf(n-1)の両方を返すようにすれば良い。以下の表を観察すると規則性が見えてくる。

	1	2	3	4	5	6	7	8	9	10	11	12
f(n)	1	2	1	3	2	3	1	4	3	5	2	5
f(n-1)	1	1	2	1	3	2	3	1	4	3	5	2

n が奇数のとき f(n-1) は f(n')+f(n'-1) になっていて、偶数ならば f(n-1) は f(n'-1)になっている。例えば、f(9-1)=f(8)=f(4)+f(4-1)で、f(10-1)=f(9)=f(5)+f(5-1)。ここまで分かれば、f(n) と f(n-1)を同時に返して末尾再帰させるコードを簡単に書ける。自分書いた4行ばかりの Haskellコードは、実行すると瞬時に正答が算出された（二重再帰版は未だに帰ってこない）。

2012年3月12日月曜日

Monoid をちょっと調べてみる自習

Monoid というのが、なんか地味っぽいけどよく使われているみたいなので、調べてみる。

◆ Data.Monoid

ghci で :i Monoidとやると、Data.Monoid で定義されているインスタンスがたくさん表示されてくるので、ソースを見ながら一個ずつ調べてみる。（「⇨」の左辺に ghci への入力、右辺に ghci の出力を書いてみた）

■ Any、All

Bool値をとる Monoid。

mempty::Any ⇨ Any {getAny = False}
Any False `mappend` Any False `mappend` Any True ⇨ Any {getAny = True}
mconcat [Any False, Any True, Any False] ⇨ Any {getAny = True}

All はこれの反対で、類推は容易。

■ Ordering

これがちょっとおもしろい。LT、EQ、GTの三つの値があって、EQ の優先度が低く、LT と GTでは早い者勝ちという感じ。

mempty::Ordering ⇨ EQ
mconcat [LT,LT] ⇨ LT, mconcat [LT,EQ] ⇨ LT, mconcat [LT,GT] ⇨ LT
mconcat [EQ,GT] ⇨ GT, mconcat [EQ,EQ] ⇨ EQ, mconcat [EQ,GT] ⇨ GT
mconcat [GT,LT] ⇨ GT, mconcat [GT,EQ] ⇨ GT, mconcat [GT,GT] ⇨ GT
mconcat [EQ,EQ,GT,LT] ⇨ GT

■ ()

() は、mempty、mappend、mconcat ともに ()

■ [a]

これも大体想像通り。mconcat では空リストが無くなる。

mempty::[Int] ⇨ []
[3] `mappend` [1,4] ⇨ [3,1,4]
mconcat [[3,1],[],[4]] ⇨ [3,1,4]

■ First a、Last a

First は Maybe a 型の値を持つ Monoid で、Nothing 以外の値で早い者勝ち。Last はその逆。

mempty::First Int ⇨ First {getFirst = Nothing}
mconcat[mempty, First Nothing, First (Just 27), mempty, First (Just 18)]
  ⇨ First {getFirst = Just 27}
mempty::Last Int ⇨ Last {getLast = Nothing}
Last (Just 27) `mappend` Last (Just 18) ⇨ Last {getLast = Just 18}
Last (Just 27) `mappend` mempty `mappend` Last (Just 18) `mappend` Last Nothing
  ⇨ Last {getLast = Just 18}

■ Endo a

なんだか endomorphism というものに関係があるらしい。wikipedia を引くと自己準同型とか自己射なるものの事らしいけど、正直、今のところよくからない。まあソースを読めば関数の合成に関係がある事はなんとなくわかるけど・・・

instance Monoid (Endo a) where
    mempty = Endo id
    Endo f `mappend` Endo g = Endo (f . g)
appEndo (Endo (*2) `mappend` mempty `mappend` Endo (+1)) 3 ⇨ 8
appEndo (mconcat [Endo (+3), Endo (+1), Endo (+4)]) 3 ⇨ 11

Foldable のソースでも使われていたりする。

■ Sum a、Product a …（a は Num）

文字通り、Sum が足し算で、Product が掛け算。

mempty::Sum Int ⇨ Sum {getSum = 0}
Sum 3 `mappend` Sum 1 `mappend` Sum 4 ⇨ Sum {getSum = 8}
getSum $mconcat [Sum 3, Sum 1, Sum 4] ⇨ 8
mempty::Product Int ⇨ Product {getProduct = 1}
Product 3 `mappend` Product 1 `mappend` Product 4 ⇨ Product {getProduct = 12}
mconcat [Product 3, Product 1, Product 4] ⇨ Product {getProduct = 12}

■ Maybe a …（a は Monoid）

中身が Monoid の Maybe が Monoid になるつう事なんだろうか。

Just [2,7] `mappend` Just [1,8] ⇨ Just [2,7,1,8]
Just EQ `mappend` Just LT ⇨ Just LT
Nothing `mappend` (Just $Sum 1) ⇨ Just (Sum {getSum = 1})
Nothing `mappend` Just GT ⇨ Just GT

■ Dual a …（a は Monoid）

なんかひっくり返すらしい。「双対」というものに関係がありそうだけど、詳細はよくわからない。Foldable の foldl の実装で使われてた。

mempty::Dual [Int] ⇨ Dual {getDual = []}
Dual [3,1] `mappend` Dual [4,1] ⇨ Dual {getDual = [4,1,3,1]}
mconcat [Dual [3], Dual [1], Dual [4]] ⇨ Dual {getDual = [4,1,3]}

■ (a -> b) …（b が Monoid）

Monoid を返す関数を Monoid として扱うって事か。

:t mempty::Monoid b =>(Int->b)
   ⇨ mempty::Monoid b =>(Int->b) :: Monoid b => Int -> b
 
(mempty::(a->String)) 333 ⇨ ""
(mempty::(a->())) "hello" ⇨ ()
((:"!") `mappend` (:"?")) '@' = "@!@?"
((compare 1) `mappend` (compare 10)) 5 ⇨ LT
(((\a b->[a + b]) `mappend` (\a b->[a * b])) 3) 10 ⇨ [13,30]

■ (a, b...) …（a から最大 e まで Monoid）

Monoid を値に持つタプルを Monoid として見て、縦に並べて計算する感じ。
5個分まで定義。

mempty:: (String, Maybe (Sum Int)) ⇨ ("",Nothing)
("Hello, ", Sum 0.3) `mappend` ("World!", Sum 0.014)
   ⇨ ("Hello, World!",Sum {getSum = 0.314})
(GT, Sum 0.3) `mappend` mempty ⇨ (GT, Sum {getSum = 0.3})

◆ 試しに書いてみる

mconcat はデフォルト実装があるので、mempty と mappend だけ書けばいいらしい（参考URL）。

ただし Monoid law に気をつける必要がある。

mempty `mappend` x = x
x `mappend` mempty = x
(x `mappend` y) `mappend` z = x `mappend` (y `mappend` z)

とりあえず「範囲」を Monoidとして書いてみた。

mappend は両方の範囲を含む最小の範囲とした。左辺が右辺より大きくなる Range も作れてしまうが、両辺ともに mempty でない場合の mappend で正規化するようにした。（mempty をふくむ mappend でこれをやると Monoid law が成立しなくなる）

data Range a = NullRange | Range (a, a) deriving (Show, Read, Eq)

instance Ord a => Monoid (Range a) where
    mempty                        = NullRange
    NullRange `mappend` range     = range
    range     `mappend` NullRange = range
    range1    `mappend` range2    =
        let Range (l1, r1) = canonicalize range1
            Range (l2, r2) = canonicalize range2
        in  Range (min l1 l2, max r1 r2)
        where canonicalize (Range (l,r)) =
          Range (min l r, max l r)

GHCI に読ませてちょっと試してみる。

ghci> mempty::Range Int
NullRange
ghci> Range (3, 6) `mappend` Range (4, 10)
Range (3,10)
ghci> mconcat [Range (2,8), NullRange, Range (13,5), Range (5,5)]
Range (2,13)

でさらに、Writer に Monoid の制約があるので、これと組み合わせて使ってみる。

printRange :: (Ord a, Show a) => (Range a) ->WriterT (Range a) IO ()
printRange range = do
    tell range
    liftIO $ print range

main = runWriterT $ do
    printRange $ Range ('g', 'u')
    printRange $ Range ('l', 'r')
    printRange $ Range ('t', 'w')

printRange は、受け取った範囲を標準出力に書き出しながら、Writer にも出力する関数。これを連続的に実行すると、最終的に「ログ」の部分は全ての範囲を含む最小の範囲になる。

ghci> :main
Range ('g','u')
Range ('l','r')
Range ('t','w')
((),Range ('g','w'))

◆ 雑感

Endo Monoid のところで、どうやら自己射／自己準同型というものに関係があるらしいと書いたけど、そもそも Monoid からして（もちろん Monad も）圏論周辺の数学用語らしいので、普通に Haskell を使う分には数学は不要とは言うけど、やっぱりある程度理解しておく必要はあるのかなと思う。

2012年3月3日土曜日

Haskell で三目並べ (２)

前回のつづきで三目並べ。
マシン側の指し手の選択を、少しだけ知的にしてみる。

以前の記事で、ルール言語の Drools を使って、下記のようなルールで三目並べを書いた。

後一手で列が完成する盤面なら、その列を完成して終了。
ユーザが後一手で一列完成できる盤面なら、それをブロックする。
後一手で二方向の王手が作れる盤面なら、それを作る。
ユーザが後一手で二方向の王手が作れる盤面なら、ブロックする。
中央が空いていたら、そこに置く。
相手が角を取っていて対角が空いていたら、そこに置く。
角が空いていたらそこに置く。
空いているマスがあったら、そこに置く。

これを Haskell で書いてみたい。（参考wiki）

◆ 方針

前回の暫定実装した、マシンの指し手を選択する関数をこんな風に書き換えてみる。

machinesTurn :: StateT Board IO Int
machinesTurn = do
    board <- get
    let idx = fromJust $ select board
    updateBoard Machine idx
    return idx
  where
    select = 
        findWinPos       <||>
        findBlockPos     <||>
        findForkPos      <||>
        findForkBlockPos <||>
        tryCenter        <||>
        findAnyBlank
    f1 <||> f2 = mplus <$> f1 <*> f2

select 以下に並べられた findWinPos から findAnyBlank までの関数は、盤面 Board を取って、指し手を Maybe Intとして返すもの。

これらの関数そのものを「(->) r」というタイプの Applicative とみなして合成する演算子を、<||> として定義した。MonadPlus としての Maybe の振る舞いを利用している。

結果として合成された select 関数は、ある関数に盤面を与えて、指し手が得られたらそれを返し、ダメなら次の関数に盤面を与えるという計算を数珠つなぎにしたものになる。 ※((->) r)を Functor や Applicative として扱う手法はこのあたりの記事に詳しい。Applicative 一般の入門としても、超絶わかりやすく解説されている。

◆ コード全体

findWinPos 〜 findAnyBlank の実装は以下のようなものになる（前回書いた部分も含めて全部晒してしまうことにした）。

import Control.Monad.State
import Control.Applicative
import Data.List
import Data.Functor.Identity
import Control.Monad.List
import Data.Maybe
import qualified Data.List.Utils as LU
import System

data Mark = User|Machine|Blank deriving Eq
instance Show Mark where
  show User    = "U"
  show Machine = "M"
  show Blank   = " "

type Board = [Mark]

main :: IO ()
main = do
    args <- getArgs
    let machineFirst = 0 < length args && "m"==args!!0
    evalStateT (play machineFirst) initialBoard >>= putStrLn

initialBoard :: Board
initialBoard = take 9 $ repeat Blank

printBoard :: Board -> IO ()
printBoard board =
    putStrLn $ unlines $ LU.join ["-+-+-"] $ map showRow [0..2]
  where
    showRow row = [LU.join "|" $ map show $ cells row]
    cells   row = take 3 $ snd $ splitAt (row * 3) board

play :: Bool -> StateT Board IO String
play isMachinesTurn = do
    if isMachinesTurn then machinesTurn
                      else usersTurn
    board <- get
    liftIO $ printBoard board
    case judge board of
        Nothing     -> (play $ not isMachinesTurn)
        Just result -> return result

judge :: Board -> Maybe String
judge board
    | wonBy User board                 = return "won by user."
    | wonBy Machine board              = return "won by machine."
    | Nothing == elemIndex Blank board = return "drawn."
    | otherwise                        = Nothing

wonBy :: Mark -> Board -> Bool
wonBy mark board = any threeInRow rows
  where
    threeInRow indices = all (mark==) $ rowStates indices
    rowStates indices  = map (board!!) indices

rows :: [[Int]]
rows = [[0,3,6],[1,4,7],[2,5,8],
        [0,1,2],[3,4,5],[6,7,8],
        [0,4,8],[2,4,6]]

updateBoard :: Mark -> Int -> StateT Board IO ()
updateBoard mark pos = get >>= put.(update mark pos)

update :: Mark -> Int -> Board -> Board
update mark pos board =
    let (f, (x:xs)) = splitAt pos board in f ++ (mark:xs)

usersTurn :: StateT Board IO Int
usersTurn = do
    position <- liftIO getDigit
    board    <- get
    if Blank /= board!!position
        then (liftIO $putStrLn "wrong position") >> usersTurn
        else updateBoard User position >> return position
  where
    getDigit = do
        ch <- getChar
        putStrLn ""
        if elem ch ['0'..'8']
            then return (read [ch])
            else putStrLn "wrong input" >> getDigit

machinesTurn :: StateT Board IO Int
machinesTurn = do
    board <- get
    let idx = fromJust $ select board
    updateBoard Machine idx
    return idx
  where
    select =
      findWinPos       <||>
      findBlockPos     <||>
      findForkPos      <||>
      findForkBlockPos <||>
      tryCenter        <||>
      findAnyBlank
    f1 <||> f2 = mplus <$> f1 <*> f2

findWinPos :: Board -> Maybe Int
findWinPos   = findWinPos' Machine

findBlockPos :: Board -> Maybe Int
findBlockPos = findWinPos' User

findWinPos' :: Mark -> Board -> Maybe Int
findWinPos' mark board =
    ff ([mark, mark, Blank], 2) `mplus`
    ff ([mark, Blank, mark], 1) `mplus`
    ff ([Blank, mark, mark], 0)
  where
    ff (css, idx) = do findPattern css board >>= return.(!!idx)

findForkPos      :: Board -> Maybe Int
findForkPos      = findForkPos' Machine

findForkBlockPos :: Board -> Maybe Int
findForkBlockPos = findForkPos' User

findForkPos':: Mark -> Board -> Maybe Int
findForkPos' mark brd = do
    let posibles = elemIndices Blank brd
    idx <- findIndex (\pos->containsFork (update mark pos brd)) posibles
    return $ posibles!!idx
  where
    containsFork  b = (2<=) $ length $ findTwoInARow b
    findTwoInARow b = findIndices (matchTwoInARow b) rows
    matchTwoInARow b row =
        let stts = map (b!!) row
        in 2 == (countState mark stts) && 1 == (countState Blank stts)

findOppositeCorner :: Board -> Maybe Int
findOppositeCorner brd = do
    let pairs = [(0,8), (8,0), (2,6), (6,2)]
    r' <- findIndex (\(f,s)-> brd!!f==User && brd!!s==Blank) pairs
    return $ snd $ pairs!!r'

findBlankCorner :: Board -> Maybe Int
findBlankCorner = findBlankCell [0, 2, 8, 6]

findBlankCell   :: [Int] -> Board -> Maybe Int
findBlankCell candidate board =
    findIndex ((Blank==).(board!!)) candidate
    >>= return.(candidate!!)

findAnyBlank :: Board -> Maybe Int
findAnyBlank board = elemIndex Blank board

tryCenter     :: Board -> Maybe Int
tryCenter brd = if Blank == brd!!4 then Just 4 else Nothing

findPattern :: [Mark] -> Board -> Maybe [Int]
findPattern pattern board = do
    idx <- elemIndex pattern $ map (\row->map (\idx->board!!idx) row) rows
    return $ rows!!idx

countState :: Mark -> [Mark] -> Int
countState state cells = length $elemIndices state cells

◆ まとめ

実は、上のアルゴリズムではユーザ先手で 0, 8, 6 と指すと、ユーザが勝ってしまう。三目並べはちゃんとやると引き分けになるんだけど、ちょっと手抜きをしている（ユーザのフォークをマシンが防ぐ辺り）。三目並べ自体が目的ではないので、割愛した。
少し書いては直し、また書いては直しと常時リファクタしているが、コンパイルが通るだけで、感覚的には JUnit で 6割くらいのコード・カバレッジがあるくらいの安定感がある。この 150行弱のコードだと、デグレ的な事は一回も起きなかった。
上で少し触れたサイト、「Learn You a Haskell for Great Good!」の、Functor 〜 Applicative 〜 Monado の解説(この順序で説明されている）が、ものすごくわかりやすくて感動。

Haskell で三目並べ (1)

Haskell の Monad transformer が大体わかってきたので、State と IO を組み合わせる練習コードを書いてみる。題材として、三目並べ（Tic-Tac-Toe）を選んでみた。

こんな仕様でやってみる。

コンソールアプリとする
実行時引数で"m"を渡したらマシン、でなければユーザの先手
0から8の文字で位置を指定する（左上が0、Enter不要）。
不正入力時、その旨表示して再入力を促す。
マシン／ユーザ双方の指し手の後に盤面を表示する。
決着が着いたら、結果を表示する

====

◆ 盤面の表現

盤面は九つのマスから成り、それぞれ U(ユーザ)か M(マシン)か空のいずれかの状態

data Mark = User|Machine|Blank deriving Eq
instance Show Mark where
  show User    = "U"
  show Machine = "M"
  show Blank   = " "

type Board = [Mark]

initialBoard :: Board
initialBoard = take 9 $ repeat Blank

下記のような関数で、盤面を文字列化してコンソールに表示できる。（LU は Data.List.Utilsのエイリアス）

printBoard :: Board -> IO ()
printBoard board =
    putStrLn $ unlines $ LU.join ["-+-+-"] $ map showRow [0..2]
  where
    showRow row = [LU.join "|" $ map show $ cells row]
    cells   row = take 3 $ snd $ splitAt (row * 3) board

タテ、ヨコ、ナナメの列を表すためのインデクスの配列は以下のように書いた。計算で導出することもできるけど、却ってわかりにくかったの明示的な定数として書き直した。

rows :: [[Int]]
rows = [[0,3,6],[1,4,7],[2,5,8],
        [0,1,2],[3,4,5],[6,7,8],
        [0,4,8],[2,4,6]]

◆ ゲームの流れの表現

main は、先手を判別してゲームを開始し、結果を表示する。

main :: IO ()
main = do
    args <- getArgs
    let machineFirst = 0 < length args && "m"==args!!0
    evalStateT (play machineFirst) initialBoard >>= putStrLn

play が実際のゲームを進行する関数で、状態としての盤面を保持しながらユーザ入力を受け付けたりする。というわけで型としては State Board と IO を組み合わせた StateTになっている。

ユーザの指し手とマシンの指し手を交互に取得しながら、盤面の状態を更新・再表示した上で勝敗を判定。決着が着けばこれを結果とし、着かなければ play に再帰するという形になる。

play :: Bool -> StateT Board IO String
play isMachinesTurn = do
    if isMachinesTurn then machinesTurn
                      else usersTurn
    board <- get
    liftIO $ printBoard board
    case judge board of
        Nothing     -> (play $ not isMachinesTurn)
        Just result -> return result

judge :: Board -> Maybe String
judge board
    | wonBy User board                 = return "won by user."
    | wonBy Machine board              = return "won by machine."
    | Nothing == elemIndex Blank board = return "drawn."
    | otherwise                        = Nothing

wonBy :: Mark -> Board -> Bool
wonBy mark board = any threeInRow rows
  where
    threeInRow indices = all (mark==) $ rowStates indices
    rowStates indices  = map (board!!) indices

◆ ユーザの指し手の取得

標準入力から一文字取得して、指し手を得る。既にマークのあるマスに置こうとしたり、0〜8以外のキーを押下するとエラー。問題なければ盤面を更新して、得られた指し手を結果とする。

usersTurn :: StateT Board IO Int
usersTurn = do
    position <- liftIO getDigit
    board    <- get
    if Blank /= board!!position
        then (liftIO $putStrLn "wrong position") >> usersTurn
        else updateBoard User position >> return position
  where
    getDigit = do 
        ch <- getChar
        putStrLn ""
        if elem ch ['0'..'8']
            then return (read [ch])
            else putStrLn "wrong input" >> getDigit

盤面の更新はこんな関数を使う。update 関数は、後で非Monad な文脈でも使えるように切り出した。

updateBoard :: Mark -> Int -> StateT Board IO ()
updateBoard mark pos = get >>= put.(update mark pos)

update :: Mark -> Int -> Board -> Board
update mark pos board =
    let (f, (x:xs)) = splitAt pos board in f ++ (mark:xs)

◆ マシンの指し手の決定（暫定）

長くなってきたので、詳細は後で別に書くことにして、盤上左上から開いているマスを探して見つかり次第返すだけの暫定実装で、とりあえず動かしてみる。

machinesTurn :: StateT Board IO Int
machinesTurn = do
    board <- get
    let idx = fromJust $ elemIndex Blank board
    updateBoard Machine idx
    return idx

ここまでのコードを、適切に import 文を追加した上で一個のファイルにまとめて、ghci に読ませると、プロンプトから下のように実行できる。

ghci> :main
0
U| | 
-+-+-
 | | 
-+-+-
 | | 

U|M| 
-+-+-
 | | 
-+-+-
 | | 

4
U|M| 
-+-+-
 |U| 
-+-+-
 | | 

U|M|M
-+-+-
 |U| 
-+-+-
 | | 

8
U|M|M
-+-+-
 |U| 
-+-+-
 | |U

won by user.
ghci>

◆ ここまでのまとめ

まあ、マシン側の指し手が全く勝つ気も守る気もない暫定実装だけど、ここまででも一応、対話的に入力を受け付けながら状態を管理していくやり方が、だいたいわかった。もっと上達した後に見直したら、いろいろツッコミどころがあるんだろうけど、とりあえず先に進んでみる。

2012年2月18日土曜日

手計算で monad を理解してみる

去年の秋頃から、暇なときにちょこちょこ Haskell をいじっているけど実に面白い。

やっぱり一番おいしいところは、Monad をはじめとする計算のやり方なんだろうけど、一行に収まるくらいの短い式を書いて、手書きで計算してみると意外と理解に役立つ事に気がついた。

例えば、runCont (return "hello") (++"!") を計算すると "hello!" になるけど、たぶん慣れてる人には自明過ぎるようなこんな計算でも、初めての時は何だか分かったような分かってないようなあやふやな感じがする。

これをこんな風に手計算してみる。（関数等の定義は「All About Monads」のここを参考にした。）

runCont (return "hello") (++"!")
= runCont(Cont $ \k->k "hello")) (++"!") … ①
= (\k->k "hello")) (++"!") … runCont の定義より
= (++"!") "hello"
= "hello!" 

① return の定義、return a = Cont $ \k -> k a より

個人的には、これであやふやな感じが解消してスッキリする。

続けて更に Cont から別の例を出してみると、callCC なんかも下の定義だけだは、何がどうなるのかよく分からない。

callCC f = Cont $ \k -> runCont (f (\a -> Cont $ \_ -> k a)) k

この辺りのサンプルコードを読むと、使い方が何となく分かるが、どういう仕組みでそうなるのか腑に落ちずモヤモヤする。

例えば「runCont (callCC (\ex -> do {ex "bye"; return "hi"})) (++"!") 」みたいな式で、なんで "hi!" にならないのか不思議だったりする。

こんなのも手で計算してみると、細かいところがはっきりする。

  runCont (callCC (\ex -> do {ex "bye"; return "hi"})) (++"!")
= runCont (callCC (\ex -> ex "bye" >>= (\f -> return "hi"))) (++"!")
= runCont (Cont $ \k ->　(\_ -> (k "bye")) k) (++"!") … (a)
= (\k ->　(\_ -> (k "bye")) k) (++"!")
= (\_ -> ((++"!") "bye")) (++"!")
= (++"!") "bye"
= "bye!"

(a)
callCC f = Cont $ \k->　runCont (f(\a-> Cont $ \_->k a)) k より

  callCC (\ex->ex "bye">>=(\f->return "hi"))
= Cont $ \k->　runCont ((\ex -> ex "bye">>=(\_ -> return "hi")) (\a-> Cont $ \_ -> k a)) k
= Cont $ \k->　runCont (((\a -> Cont $ \_ -> k a) "bye") >>= (\_ -> return "hi")) k
= Cont $ \k->　runCont ((Cont $ \_ -> k "bye") >>= (\_ -> return "hi")) k
= Cont $ \k->　runCont (Cont $ \_ -> (k "bye")) k … (b)
= Cont $ \k->　(\_ -> (k "bye")) k

(b)
(Cont c) >>= f = Cont $ \k' -> c (\a -> runCont (f a) k') より

  (Cont $ \_->k "bye") >>= (\_ -> return "hi")
= Cont $ \k' -> (\_ -> k "bye") (\a -> runCont (f a) k')
= Cont $ \_ -> (k "bye")

実は最初、遅延評価に変にこだわってしまって手こずっていたんだけど、参照透明なんだから好きな順序で計算しても結果は変わらないわけで、手計算のときは余り気にする必要が無いと開き直って簡単になった。

2012年6月13日水曜日

ピタゴラス数の小ネタ

2012年5月2日水曜日

Pell方程式を解くための実装メモ

2012年3月27日火曜日

ブルートフォース癖

2012年3月12日月曜日

Monoid をちょっと調べてみる自習

2012年3月3日土曜日

Haskell で三目並べ (２)

Haskell で三目並べ (1)

2012年2月18日土曜日

手計算で monad を理解してみる

yasuabe blog

自己紹介

math

Link List A

このブログを検索

myscript

ラベル

ブログアーカイブ

Twitter

2012年6月13日水曜日

2012年5月2日水曜日

2012年3月27日火曜日

2012年3月12日月曜日

2012年3月3日土曜日

2012年2月18日土曜日

自己紹介

math

Link List A

このブログを検索

myscript

ラベル

ブログ アーカイブ

Twitter

ブログアーカイブ