「再帰のこころ」「反復のこころ」そして「メモ化」少々

「n個の対応する括弧のパターン」(http://d.hatena.ne.jp/ymotongpoo/20101209/1291879286)についてあれこれを徒然に考えてみました.

再帰で考える

まずはパターン数を数えることにしましょう.この手の問題はまずは再帰で考えるのが基本です.

m個のカッコ(開き括弧)とn個コッカ(閉じ括弧)を並べるのであるが,m \le nであり,並べる途中においてもコッカの数がカッコの数を越えることがないようにするものとします.このときの並べ方の数を kakoka (m,n) としましょう.この kakoka (m,n) を構成できれば,n組の対応する括弧を含む列の数は parenseq n = kakoka (n,n) というわけです.

さて,この kakoka を考える最もシンプルな方法は再帰です.再帰というと難しいと思うプログラマもいるかもしれませんが,「再帰のこころ」さえ掴んでしまえば,もっとも簡単にものを考える手法になります.「再帰のこころ」は実に単純で,

一歩手前まで出来たとしたら,完成させるにはどうする?

です.まず一歩手前を考えましょう.

  • m == n すなわち,カッコとコッカが同じだけなら,最後に付け加えたのはコッカ.したがって,一歩手前,すなわち,そのコッカを付け加える前の状態の数 kakoka (m,n-1) と同じ.
  • m > n の場合は,最後に付け加えたのはカッコ,または,コッカ.カッコを付け加える前の状態の数は,kakoka (m-1,n),また,コッカを付け加える前の状態の数は,kakoka (m,n-1).したがって,kakoka (m,n) は kakoka (m-1,n) と kakoka (m,n-1)の和に等しい.

簡単ですよね.もちろん,一歩手前を考えるまでもなく自明の場合というのを押えておかなければなりません.

kakoka :: (Int,Int) -> Integer
kakoka (_,0)   = 1
kakoka (m,n) 
  | m >  n     = kakoka (m-1,n) + kakoka (m,n-1)
  | m == n     =                  kakoka (m,n-1)
  | otherwise  = error $ "(m,n) = " ++ show (m,n) ++ ": At least, 'm' must equal to 'n'"

この定義はわかりやすい定義です.問題そのものの構造が再帰構造で,それをそのままプログラムとして表現しているからです.この形状の再帰は枝分れしていくので樹状再帰(tree recursion)と呼ばれたりします.

すこし実行してみましょう.

ghci> let parenseq n = kakoka (n,n)
ghci> parenseq 10
16796
(0.06 secs, 4893428 bytes)
ghci> parenseq 11
58786
(0.17 secs, 16238636 bytes)
ghci> parenseq 12
208012
(0.57 secs, 56756540 bytes)
ghci> parenseq 13
742900
(1.99 secs, 202941288 bytes)
ghci> parenseq 14
2674440
(7.11 secs, 727493260 bytes)
ghci> parenseq 15
9694845
(25.81 secs, 2630816028 bytes)

おそろしく効率が悪そうです.これは枝分かれした再帰の先で重複した計算が何度もおこなわれているので効率が悪いのです.この手の問題に対処するにはいくつかの方法(まとめて「動的計画法」「Dynamic Programming」「DP」などと呼ばれる方法)があります.

メモ化

重複する計算を,メモ表検索に置き換える方法をメモ化(memoization メモワイゼーション)といいます.メモ化では関数適用の際に引数をキーにしてメモ表を引き,メモ表になければ,通常の関数適用計算を行い,引数をキーに結果をメモ表に登録します.メモ表に結果が登録されていれば,通常の関数適用計算は行わずにメモ表に登録されている計算結果を返します.

通常の関数をメモ化する関数 memoize があれば便利です.しかし,メモ化というのは関数の機能を変更するのではなく,関数適用の機能を変更することなので,再帰的に定義された関数では memoize が効きません.

「こういうときは不動点が kokako になるような汎関数 kakokaF を考えよ.」というのがHaskell一門の家訓です.:p (汎関数という用語はtwitterで教えていただきました)

kakokaF :: (Num a, Ord a, Num b) => ((a,a) -> b) -> (a,a) -> b
kakokaF _ (_,0) = 1
kakokaF f (m,n)
  | m >  n    = f (m-1, n) + f (m,n-1)
  | m == n    =              f (m,n-1)
  | otherwise  = error $ "(m,n) = " ++ show (m,n) ++ ": At least, 'm' must equal to 'n'"

kakokaMemo :: (Int,Int) -> Integer
kakokaMemo = memoize kakokaF
  where
    memoize   =  memo tableinfo
    tableinfo =  tableinfo :: (Memo M.Map) (Int,Int) Integer

汎関数の型が多相型になっているのは,通常再帰版とメモ化再帰版でこの汎関数を再利用できるようにしているためです.通常版は,たとえば,kakoka = fix kakokaF と定義できます.

tableinfoはメモ表の実装を示すためだけにあるり,必要なのは型情報だけです(値は\botで構いません).高階関数 memo は sampou.org 謹製:pのControl.Memoで定義されています.

 % darcs get http://darcs.sampou.org/memo/

で取得でき,cabal install でインストールできます.インストール後 Control.Memo をインポートしてください.このモジュールの詳細はソースコードをご覧ください.

メモ化の威力はどうでしょうか.

ghci> parenseq 15
9694845
(0.00 secs, 531092 bytes)
ghci> parenseq 150
620925183926009621146978506218967449531342090729015621989883130549504437230725772687824
(0.10 secs, 22565612 bytes)

なかなかのものですね.このようなメモ化はメモ化されていない版と汎関数を共有しますので,関数の意味が判りやすいのがメリットです.

当然のことながら,メモ表を引くのはタダではありませんし,検索がヒットしなければ,検索コスト以外に新しいエントリをメモ表に追加するコストがかかります.したがって,重複が少い場合にはメモ化すると逆効果ということもありえます.

この「括弧列問題」では問題になりませんが,メモ化した関数はメモ表のキーとなる引数に関して,常に eager evaluation になることにも注意してください.たとえば,竹内のたらい回し関数の場合,メモ化すると lazy evaluation ではなくなるので遅くなります.

ボトムアップ

メモ法は枝分かれする再帰計算で起こる冗長な重複計算をメモ表をつかって回避するというものでした.こんどは別の方向で考えてみましょう.関数 kakoka は「再帰の心」で素直に定義されたもので,ゴールからトップダウンに考えたものです.トップダウンで再帰的に分岐された計算はそれぞれ独立ですが,計算内容が重複しているので効率が悪いというわけです.

そこで,スタートからゴールに向けてボトムアップに計算を組み立ててみます.
ボトムアップで鍵になるのは,「反復の心」です.

これまで出来てたとしたら,一歩先はどうなる?

これが「反復のこころ」です.まぁ要するに「末尾再帰」というやつです.
蓄積変数を使って成果を次に伝えるというのが基本的な実装方法です.
もちろん,どこで終了するかを判定するための仕組みが必要です.

この図の右上向き矢印は'('の追加を表し,右下向き矢印は')'の追加を表しています.黒の経路を左から辿ると,(()(()(()))) となります.

1ステップでカッコまたはコッカを追加していくことを考えます.

全体は 2n ステップですが,半分の n ステップで進んだところで残りは対称になっているので,n ステップで到達できる位置での経路数を自乗して和をとれば全体の経路数になります.

というわけで,nステップまでの範囲で考えます.カッコを p 個,コッカを q 個 を使う場合,矢印に沿って到達する位置を (p,q) と表現することにします.

  • 現在 2i-1 ステップまで来たとすると,現在到達した位置は,(2i-1,0),(2i-2,1),..,(i+1,i) です.次の2iステップ目で到達できるのは,現在の括弧列の後ろへ"("を連結して,(2i,0),(2i-1,1),..,(i,i) の位置に来るか,現在の括弧列の後ろへ")"を連結して,(2i-1,1),(2i-2,2),..,(i,i)の位置に来る.
  • 現在 2i ステップまで来たとすると,現在到達した位置は,(2i,0),(2i-1,1),..,(i,i) です.次の2i+1ステップ目で到達できるのは,現在の括弧列の後ろへ"("を連結して,(2i+1,0),(2i,1),..,(i+1,i) の位置に来るか,現在の括弧列の後ろへ")"を連結して,(2i,1),(2i-1,2),..,(i+1,i)の位置に来る.

そういうわけで,奇数ステップと偶数ステップでは少しだけ異ります.これを同じ関数で表現することもできますが,交互に使えばよいだけなので,ここでは別の関数として定義します.

kkkso,kkkse :: Int -> [Integer] -> [Integer]
kkkso 0 xs  = xs
kkkso n xs  = kkkse (n-1) (zipWith (+) (xs ++ [ ]) ([0] ++ xs))

kkkse 0 xs  = xs
kkkse n xs  = kkkso (n-1) (zipWith (+) (xs ++ [0]) ([0] ++ xs))

全体では

parenseq n = sum $ map (^2) $ kkkso n [1]

となります.

ghci> let parenseq n = sum $ map (^2) $ kkkso n [1]
ghci> parenseq 15
9694845
(0.00 secs, 0 bytes)
ghci> parenseq 150
620925183926009621146978506218967449531342090729015621989883130549504437230725772687824
(0.01 secs, 1061312 bytes)

途中のkステップで辿り付ける「複数の」点への経路をすべて求めるのは,一見冗長に見えますが,結局これはどれも必要になるので冗長ではないわけです.

最速なのは

実はこのパターンの数は n に関する閉じた式として導けます.

C_n = \frac{(2n!)}{(n+1)!n!}

これは良く知られたCatalan数というものです.(実際にどう導くかはググって下さい.)
というわけで,これを使えば一瞬で計算が終ります.

(!) :: Int -> Integer
(!) 0 = 1
(!) n = toEnum n * ((n-1)!)

catalan :: Int -> Integer
catalan n = ((2*n)!) `div` ((n+1)!) `div` (n!)

括弧列を文字列として列挙する

文字列の列挙であっても,経路数を数え上げるのとプログラムの構造は同じです.

再帰版
kkakkokka :: (Int,Int) -> [String]
kkakkokka (m,0) = [replicate m '(']
kkakkokka (m,n)
  | m == n      =                                    map (++")") (kkakkokka (m,n-1))
  | m >  n      = map (++"(") (kkakkokka (m-1,n)) ++ map (++")") (kkakkokka (m,n-1))
  | otherwise   = error $ "(m,n) = " ++ show (m,n) ++ ": At least, 'm' must equal to 'n'"

計算してみましょう.

ghci> let parenseq n = kkakkokka (n,n)
ghci> parenseq 5
["()()()()()","(())()()()","()(())()()","(()())()()","((()))()()","()()(())()"
,"(())(())()","()(()())()","(()()())()","((())())()","()((()))()","(()(()))()"
,"((()()))()","(((())))()","()()()(())","(())()(())","()(())(())","(()())(())"
,"((()))(())","()()(()())","(())(()())","()(()()())","(()()()())","((())()())"
,"()((())())","(()(())())","((()())())","(((()))())","()()((()))","(())((()))"
,"()(()(()))","(()()(()))","((())(()))","()((()()))","(()(()()))","((()()()))"
,"(((())()))","()(((())))","(()((())))","((()(())))","(((()())))","((((()))))"]
(0.02 secs, 1072644 bytes)

括弧の数が多くなると列挙するのは大変なので,括弧列を数えてみましょう.

ghci> length $ parenseq 10
16796
(0.10 secs, 24669152 bytes)
ghci> length $ parenseq 11
58786
(0.40 secs, 91803232 bytes)
ghci> length $ parenseq 12
208012
(1.44 secs, 347797948 bytes)
ghci> length $ parenseq 13
742900
(5.19 secs, 1325709196 bytes)
ghci> length $ parenseq 14
2674440
(18.89 secs, 5080014448 bytes)
ghci> length $ parenseq 15
9694845
(69.17 secs, 19534267712 bytes)
メモ版

これも汎関数を使います.

kkakkokkaF :: (Monad f, Functor f, Monoid (f [String]))
           => ((Int, Int) -> f [String]) -> (Int, Int) -> f [String]
kkakkokkaF f (m,0) = return [replicate m '(']
kkakkokkaF f (m,n)
  | m == n         =                                          fmap (map (++")")) (f (m,n-1))
  | m >  n         = fmap (map (++"(")) (f (m-1,n)) `mappend` fmap (map (++")")) (f (m,n-1))
  | otherwise      = error $ "(m,n) = " ++ show (m,n) ++ ": At least, 'm' must equal to 'n'"

kkakkokkaMemo :: (Int,Int) -> [String]
kkakkokkaMemo = memoise kkakkokkaF
  where
    memoise   =  memo tableinfo
    tableinfo =  tableinfo :: (Memo M.Map) (Int,Int) [String]

実行してみましょう.

ghci> parenseq 5
["()()()()()","(())()()()","()(())()()","(()())()()","((()))()()","()()(())()"
,"(())(())()","()(()())()","(()()())()","((())())()","()((()))()","(()(()))()"
,"((()()))()","(((())))()","()()()(())","(())()(())","()(())(())","(()())(())"
,"((()))(())","()()(()())","(())(()())","()(()()())","(()()()())","((())()())"
,"()((())())","(()(())())","((()())())","(((()))())","()()((()))","(())((()))"
,"()(()(()))","(()()(()))","((())(()))","()((()()))","(()(()()))","((()()()))"
,"(((())()))","()(((())))","(()((())))","((()(())))","(((()())))","((((()))))"]
(0.02 secs, 1688036 bytes)
ghci> length $ parenseq 10
16796
(0.02 secs, 4725776 bytes)
ghci> length $ parenseq 11
58786
(0.04 secs, 15731016 bytes)
ghci> length $ parenseq 12
208012
(0.21 secs, 54009780 bytes)
ghci> length $ parenseq 13
742900
(0.71 secs, 193462952 bytes)
ghci> length $ parenseq 14
2674440
(2.52 secs, 697282508 bytes)
ghci> length $ parenseq 15
9694845
(9.19 secs, 2531698676 bytes)
ボトムアップ

こちらも,場合の数を数えるプログラムと同じ構造なので簡単に書けます.

kkksso,kkksse :: Int -> [[(String,String)]] -> [[(String,String)]]
kkksso 0 xs  = xs
kkksso n xs  = kkksse (n-1) (zipWith (++) (map oc xs ++ [  ]) ([[]] ++ map co xs))

kkksse 0 xs  = xs
kkksse n xs  = kkksso (n-1) (zipWith (++) (map oc xs ++ [[]]) ([[]] ++ map co xs))

oc,co :: [(String,String)] -> [(String,String)]
oc = map (('(':) *** (')':))
co = map ((')':) *** ('(':))

全体としては

parenseqIter n = concatMap (prod . unzip) $ kkksso n [[("","")]]
prod (xss,yss) = [foldl (flip (:)) ys xs | xs <- xss, ys <- yss]

これも実行してみましょう.

ghci> parenseqIter 5
["((((()))))","()((()))()","()((())())","()((()()))","()(((())))","(()(()))()"
,"(()(())())","(()(()()))","(()((())))","((()()))()","((()())())","((()()()))"
,"((()(())))","(((())))()","(((()))())","(((())()))","(((()())))","()()()()()"
,"()()()(())","()()(())()","()()(()())","()()((()))","(())()()()","(())()(())"
,"(())(())()","(())(()())","(())((()))","()(())()()","()(())(())","()(()())()"
,"()(()()())","()(()(()))","(()())()()","(()())(())","(()()())()","(()()()())"
,"(()()(()))","((()))()()","((()))(())","((())())()","((())()())","((())(()))"]
(0.01 secs, 1075948 bytes)

さきほどの parenseq 5 とはすこし違うように見えるかもしれませんが,同じ括弧列が列挙されているはずです.両方ともソートして比較してみれば同じであることがわかります.

ghci> sort (parenseq 5) == sort (parenseqIter 5)
True
(0.01 secs, 1052336 bytes)

効率を見ましょう.

ghci> length $ parenseqIter 10
16796
(0.02 secs, 2119860 bytes)
ghci> length $ parenseqIter 11
58786
(0.04 secs, 5835536 bytes)
ghci> length $ parenseqIter 12
208012
(0.10 secs, 19139136 bytes)
ghci> length $ parenseqIter 13
742900
(0.36 secs, 66548864 bytes)
ghci> length $ parenseqIter 14
2674440
(1.07 secs, 239042944 bytes)
ghci> length $ parenseqIter 15
9694845
(3.28 secs, 864143584 bytes)

列挙の効率

文字列列挙の場合にはメモ化の効果やボトムアップの効果が意外に少いと思った方,正直に手を上げなさい.:)

こういう解の列挙問題の計算量について,たまに思い違いをするプログラマの方がいらっしゃいますが,全数列挙するなら当然その計算量はどんなに頑張っても,列挙する対象の数のに比例する計算量になります.

引っ掛け問題に,「2分探索木に含まれる要素を列挙するのにかかる計算時間を見積もれ」というのがあるそうです.わかりますよね.

追記(18:00:35)

母関数を使う解法.解説が素晴しい.id:nineties さんの母関数を使って括弧のパターンを求めてみたも是非.