Haskell で動的計画法を書くための3つの方針

Haskell

アルゴリズムの代表っぽい存在とも言えるDPですが，Haskellは参照透明なので書きにくいと思われがちです．
しかし，実際は，C言語やSTLなしのC++より遥かに簡単に動的計画法が書けます．

リストを用いる

最初に知るであろう方法．
フィボナッチ数列の第100項だと，

let f = 0 : 1 : zipWith (+) a (tail a) in f!!100

です．
「ずらして足したものを後ろにつなげる」と言えばいいのでしょうか．
これについては，他でよく解説されているので詳しくは説明しません．

利点：

  • リストの知識のみでよい．
  • 無限リストの恩恵が受けられる．
  • importが不要．

欠点：

  • 使えるケースが限られる．
  • 書くときに，混乱することも（上の例だと，F_{k+2}=F_k+F_{k+1}と考える必要がある．
  • ランダムアクセスができないので，場合によってはO(n)倍の時間がかかる*1．

添字を自由に用いることができないので，ただの「Haskellでフィボナッチ数が楽にかける」という宣伝ぐらいにしか使えないかもしれません．

Arrayを用いる

添字を用いて参照したいときに役立つ方法．

import Array
let f = array (0,100) ((0,0):(1,1):[(k, f!(k-2)+f!(k-1))| k<-[2..100]]) in f!100

と書けるのですが，要は．F_0=0,F_1=1,F_k=F_{k-2}+F_{k-1}(2\leq k\leq 100)と定義を並べるだけでよいということです．
変数を101個作り，依存関係はHaskellに任せられると考えるといいのかも知れません．


遅延評価のおかげで，

import Array
let f = array (0,200) ((0,0):(2,1):[(k, f!(k-4)+f!(k-2))| k<-[4,6..200]]) in f!200

と使わない添字があってもちゃんと動作します*2し，

import Array
let f = array (0,100) ((1,1):(0,0):reverse [(k, f!(k-2)+f!(k-1))| k<-[2..100]]) in f!100

のように，各添字の値の定義を書く順番を変えても動作します*3．

利点：

  • 普通のDPならほとんど書ける．
  • 定義を並べるだけで書ける．

欠点：

  • Arrayがいる．
  • 計算に必要な添字の範囲を知っている必要がある．

Stateモナドを用いる

Arrayを用いた場合は添字の範囲をあらかじめ定めなければなりません*4．
これは手続き型言語で書いた時に配列の大きさを決めなければならないという問題に対応し，こういう時は連想配列を使うことが定石です．
さて，Haskellでどのように書けばよいかというと，計算中に答えを蓄えた連想配列を持ちまわるという「状態つき計算」を行うことになるので，Stateモナドを用いること楽に書けます．
参考URL：

Stateモナドがghc6から標準では付属していないようなので，日記の最後に適当な私の実装を載せておきます．

まず，DPをするためのライブラリはこんな感じです（ghc向け）．

-- DP.hs
module DP (
  DP
  , dp
  , evalDP, evalDPAll
  ) where

import qualified Data.Map
import State
-- ここは環境に応じて適切に書き換えてください．

type Memo a b = Data.Map.Map a b
type DP a b = a -> State (Memo a b) b

emptyMemo :: Memo a b
emptyMemo = Data.Map.empty
lookupMemo :: Ord a => a -> Memo a b -> Maybe b
lookupMemo = Data.Map.lookup
insertMemo :: Ord a => a -> b -> Memo a b -> Memo a b
insertMemo = Data.Map.insert

dp :: Ord a => DP a b -> DP a b
dp f x = do
  memo <- gets (lookupMemo x)
  case memo of
       Just y  -> return y
       Nothing -> do
         y <- f x
         modify (insertMemo x y)
         return y

evalDP :: DP a b -> a -> b
evalDP f x = evalState (f x) emptyMemo

evalDPAll :: DP a b -> [a] -> [b]
evalDPAll f xs = evalState (sequence (map f xs)) emptyMemo

細部を省いて説明すると，DP a bというのが，型aから型bへのメモ化再帰をする関数の型です．
メモ用の連想配列にはData.Mapというバランス木を用い，dpは，関数を受け取り「引数をlookupしてメモにあればそれを返し，なければ，計算して結果をinsertしつつ返す」という新しい関数を返します．
evalDPは空のメモから実行して結果を返す関数，
evalDPAllはその複数の引数でメモを使い回して実行するバージョンです．

使い方は，

fib :: Int -> Integer
fib = evalDP fib'
fib' :: DP Int Integer
fib' = dp $ \ n ->
  if n <= 1
     then return (toInteger n)
     else do
       a <- fib' (n-2)
       b <- fib' (n-1)
       return (a+b)

のように，関数定義にdpを挟むだけ．
というと，少し言い過ぎで，モナドなので，returnとかdoとか*5書く必要があります．
実際，再帰呼び出しの順序が異なると異なる計算になる*6ので，この制約は回避不可能です．


利点：

  • とりあえず万能．

欠点：

  • 毎回ライブラリを手で書くのは面倒*7．
例題

Problem 14 - Project Euler
「100万未満の数から始めて，コラッツ予想の数列を作るとき，1になるまでのステップが最も長いのは？」
の解答：

import DP

main = print . solve $ 10^6 - 1

solve nMax = snd . maximum . flip zip [1..] . evalDPAll collatzLength $ [1..nMax]

collatzLength :: DP Integer Int
collatzLength = dp $ \ n ->
  if n == 1
     then return 0
     else do
       l <- collatzLength (collatzIter n)
       return $ l + 1

collatzIter :: Integer -> Integer
collatzIter n =
  if even n
     then n`div`2
     else 3*n+1
Stateモナド実装例
-- State.hs
module State (
  State
  , runState
  , get, put
  , modify, gets
  , evalState, execState
  ) where

newtype State s a = State { runState :: (s -> (a,s)) }

instance Monad (State s) where
  return a = State $ \ s -> (a,s)
  (State x) >>= f = State $ \ s ->
    let (v,s') = x s in runState (f v) s'

get :: (State s) s
put :: s -> (State s) ()
get = State $ \ s -> (s,s)
put s = State $ \ _ -> ((),s)

modify :: (s -> s) -> (State s) ()
modify f = get >>= put . f

gets :: (s -> a) -> (State s) a
gets f = get >>= return . f

evalState :: State s a -> s -> a
evalState = (fst .) . runState
execState :: State s a -> s -> s
execState = (snd .) . runState

*1:時間がかかっても良いから書けるかと言われると少し自信がありません

*2:ただし，メモリは消費します

*3:ちなみに，同じ添字を2回以上書くとghcでは最後のを用いるが，本当はundefinedらしい

*4:リストだとそもそも自分より大きい添字を参照できない

*5:liftM2を用いた方が楽に書けるという人はそれでもOK

*6:メモされる順番が異なり，木の細部が異なってくる

*7:つまり，「持ち込み可の試験」に限られる