累積百ます計算,パスカルの三角形,関・ベルヌーイ数を計算する


百ます計算

百ます計算を解くプログラムを書いてみます.


masu100.hs

smpl :: Num a => ([a], [a]) -> [[a]]

smpl (cs, rs) = [[r+c | c <- cs] | r <- rs]

smpl は行と列のそれぞれの数列をリストで受け取り,答を二次元配列の代わりにリストのリストで返します.

ふたつのリストの直積の和を返すだけのシンプルな実装です.

以降で使う問題を用意します.


masu100.hs

cols = [8,2,4,6,1,8,9,3,1,7]

rows = [4,2,5,6,7,1,3,9,3,2]

smpl に渡して計算してみます.

λ> smpl (cols, rows)

[[12,6,8,10,5,12,13,7,5,11]
,[10,4,6,8,3,10,11,5,3,9]
,[13,7,9,11,6,13,14,8,6,12]
,[14,8,10,12,7,14,15,9,7,13]
,[15,9,11,13,8,15,16,10,8,14]
,[9,3,5,7,2,9,10,4,2,8]
,[11,5,7,9,4,11,12,6,4,10]
,[17,11,13,15,10,17,18,12,10,16]
,[11,5,7,9,4,11,12,6,4,10]
,[10,4,6,8,3,10,11,5,3,9]]
λ>

表示用のdraw'関数を実装してあれば以下のように分かりやすい表示もできます.

今回の記事ではdraw'関数についての詳細な説明は省略します.

ひとことだけ添えておくと最初の引数は表示における列幅指定になっています.

もちろん全部計算が終わってから最大値の十進表記における桁数を計算して自動で列幅を求めることも可能ですが

この記事ではあえて指定する方針でいきます.

ともかく実行するとこんな感じになりました.

λ> draw' 3 smpl cols rows

8 2 4 6 1 8 9 3 1 7
========================================
4| 12 6 8 10 5 12 13 7 5 11
2| 10 4 6 8 3 10 11 5 3 9
5| 13 7 9 11 6 13 14 8 6 12
6| 14 8 10 12 7 14 15 9 7 13
7| 15 9 11 13 8 15 16 10 8 14
1| 9 3 5 7 2 9 10 4 2 8
3| 11 5 7 9 4 11 12 6 4 10
9| 17 11 13 15 10 17 18 12 10 16
3| 11 5 7 9 4 11 12 6 4 10
2| 10 4 6 8 3 10 11 5 3 9


累積百ます計算

百ます計算はます目の値を求めるのに,所与の2つの列からi番目とj番目をそれぞれ取り出して和を求めました.

これを直左と直上のふたつの数の和でます目を埋めるようにルールを変更します.

たとえば5行8列目のます目の数は4行8列目のます目の数と5行7列目のます目の数の和で求めるのです.

つまりr行c列目の値は以下のような漸化式で定義されます.

val(r, c) = val(r-1, c) + val(r, c-1)

ここで0行目は所与の問題の行を表現する数列で,0列目は所与の問題の列を表現する数列になります.


ナイーブな実装

まずは漸化式をそのままナイーブに実装してみます.


naive.hs

naive :: ([Int], [Int]) -> [[Int]]

naive (cs, rs) = [[val (i, j) | j <- [0..c']] | i <- [0..r']]
where
(c', r') = (length cs - 1, length rs - 1)
val (0, 0) = rs !! 0 + cs !! 0
val (0, j) = val (0, j-1) + cs !! j
val (i, 0) = rs !! i + val (i-1, 0)
val (i, j) = val (i, j-1) + val (i-1, j)

このnaiveを使って同じ問題を解くと以下のようになります.

λ> draw' 6 naive cols rows

8 2 4 6 1 8 9 3 1 7
======================================================================
4| 12 14 18 24 25 33 42 45 46 53
2| 14 28 46 70 95 128 170 215 261 314
5| 19 47 93 163 258 386 556 771 1032 1346
6| 25 72 165 328 586 972 1528 2299 3331 4677
7| 32 104 269 597 1183 2155 3683 5982 9313 13990
1| 33 137 406 1003 2186 4341 8024 14006 23319 37309
3| 36 173 579 1582 3768 8109 16133 30139 53458 90767
9| 45 218 797 2379 6147 14256 30389 60528 113986 204753
3| 48 266 1063 3442 9589 23845 54234 114762 228748 433501
2| 50 316 1379 4821 14410 38255 92489 207251 435999 869500
λ>

出来ました.

正しく計算されていますが,このnaiveは行数や列数を増やすと急速に遅くなります.

例えば10行10列のどちらかを倍にしてdraw' 6 naive (cols++rows) rowsなどとすれば十分遅さを体感できるかと思います.

これは漸化式にしたがって順々により小さなます目の数を求めていて,しかもその小さなます目の数を求めるのにさらに小さなます目の数を求めていて,しかもそうして途中で求めた数はどこかに保持しているわけでもないので必要になるたびに同じます目に相当する数を何度も再計算しているからです.

とくに同じ行でも右の方にいくに従いどんどん遅くなることも体感してみてください.


累積百ます計算を高速化する

この累積百ます計算を高速化したい思います.

そもそもこの漸化式をなぞっただけの実装ではなく,おおきく構造を変更します.

つまり二分木を使います.

val(r, c) = val(r-1, c) + val(r, c-1)

このような漸化式だったので同じ構造の部分式を2つ子要素に持つ再帰構造をしているので,そう不自然な話ではないと思います.

つまり10x10のます目であれば,一番右下のます目である(10,10)の位置にあるます目をルートとした二分木を構成したいと思います.

ということで二分木を基本構造にしますが,その上で今回の高速化のポイントはふたつです.


  1. 各ます目について計算結果を保持できるようにする.

  2. 各ます目は直右のます目の子ノードとして,そして直下のます目の子ノードとして共有されるように構成する.

ポイントのひとつ目は一般にtabulationとして知られています.

いわゆるDynamic Programming(以下DPと略す)の具体的な実現方法のひとつです.

Bottom upなDPと見做すことができます.

一方,Top downなDPはもちろんmemo(r)izationです.

実際に手で百ます計算をする場合には左上のます目から数を埋めていき,すでに計算済みのます目の数をそのまま利用して隣接するます目を埋めていくと思いますが,まさにそのように解き進めようというアイディアです.

ポイントのふたつ目は一般にNexus(ネクサス)として知られています.

ただ単に途中の計算結果をアノテーションするだけで計算の重複を除去できるわけではありません.

上流にある複数の異なるノードから同じノードへ合流するノードがあれば,重複する計算を共有することができるので高速化の可能性がでてきます.

ちょっと具体的に説明しましょう.

先程の結果を再掲します.

λ> draw' 6 naive cols rows

8 2 4 6 1 8 9 3 1 7
======================================================================
4| 12 14 18 24 25 33 42 45 46 53
2| 14 28 46 70 95 128 170 215 261 314
5| 19 47 93 163 258 386 556 771 1032 1346
6| 25 72 165 328 586 972 1528 2299 3331 4677
7| 32 104 269 597 1183 2155 3683 5982 9313 13990
1| 33 137 406 1003 2186 4341 8024 14006 23319 37309
3| 36 173 579 1582 3768 8109 16133 30139 53458 90767
9| 45 218 797 2379 6147 14256 30389 60528 113986 204753
3| 48 266 1063 3442 9589 23845 54234 114762 228748 433501
2| 50 316 1379 4821 14410 38255 92489 207251 435999 869500
λ>

4行4列目の328という数が入っているます目に注目します.

この328は直上の163と直左の165の和で,それらを子要素にもつ二分木の親ノードになっています.

この163の子ノードのうちのひとつ93という数が3行3列目にありますが,この93は165の方の子要素でもありますね.

         |        |

....... 93 <--- 163
^ ^
| |
....... 165 <--- 328

完全二分木であればこの163から来る93と165から来る93が別々のノードになってしまいます.

そうすると結局2度別々に計算されてしまうので,やっぱり嬉しくありません.

この図のように共有した構造を構成したいわけです.

そうしておいて下から,この表で言えば1行めや1列目といった端の方から計算を進めてます目を埋めていけば

人間が手でます目を埋めていく時のように効率的に計算できるわけです.

こうすると最終的には10行10列目だけが求まって終りますが,それ以外のます目の値も欲しいので,それらのノードをつまんでおくようなリストのリストも同時に組み上げるようにします.

というわけで概略説明はおわり.

とっとと実装してみます.

まずはじめに二分木のもととなる台関手を定義します.


box100.hs

data TreeF a x = Tip a | Bin x x deriving (Show)

type Tree a = Fix (TreeF a)

tip :: a -> Tree a
tip = In . Tip

bin :: Tree a -> Tree a -> Tree a
bin l r = In (Bin l r)


二分木型はこのTreeFの不動点として定義されます.

ノードに値を保持できるようにします.


box100.hs

tip' :: a -> Cofree (TreeF a) a

tip' n = Cf (In (Hisx (n, Tip n)))

bin' :: Num a => (Cofree (TreeF t) a, Cofree (TreeF t) a) -> Cofree (TreeF t) a
bin' (l, r) = Cf (In (Hisx (extract l + extract r, Bin (unCf l) (unCf r))))


これで二分木のノードにそのノードの値を保持できるようになりました.

ですが上でも説明した通り,これだけだと左の枝の右の枝と,右の枝の左の枝とが同じマス目を意味するのに実際には別のノードになってしまいます.

結局別々に計算されてしまって意味がないのでうまくネクサスを構成するように組み上げます.


box100.hs

winder :: ((a, b) -> c) -> (b, [a]) -> Maybe (c, (c, [a]))

winder f (y, xxs) = case xxs of
[] -> Nothing
(x:xs) -> Just (y', (y', xs)) where y' = f (x, y)

windCol :: Num a => (Cofree (TreeF t) a, [Cofree (TreeF t) a]) -> [Cofree (TreeF t) a]
windCol = unfoldr (winder bin')

nexus :: Num a => ([a], [a]) -> [[Cofree (TreeF a) a]]
nexus = unfoldr psi . tupply (map tip')
where
psi (cs, []) = Nothing
psi (cs, r:rs) = Just (ps, (ps, rs)) where ps = windCol (r, cs)

calc :: Num a => ([a], [a]) -> [[a]]
calc = map (map extract) . nexus


これで実装できました.

Cofree (TreeF t) aTreeF tという代数にa型の値をアノテーションとして付与した型にあたります.

実行して効果のほどを見てみましょう.

perfCheck act = do

s <- getCurrentTime
act
e <- getCurrentTime
return (diffUTCTime e s)

こんな関数を用意して時間を計測してみます.

λ> perfCheck (draw' 6 naive cols rows)

8 2 4 6 1 8 9 3 1 7
======================================================================
4| 12 14 18 24 25 33 42 45 46 53
2| 14 28 46 70 95 128 170 215 261 314
5| 19 47 93 163 258 386 556 771 1032 1346
6| 25 72 165 328 586 972 1528 2299 3331 4677
7| 32 104 269 597 1183 2155 3683 5982 9313 13990
1| 33 137 406 1003 2186 4341 8024 14006 23319 37309
3| 36 173 579 1582 3768 8109 16133 30139 53458 90767
9| 45 218 797 2379 6147 14256 30389 60528 113986 204753
3| 48 266 1063 3442 9589 23845 54234 114762 228748 433501
2| 50 316 1379 4821 14410 38255 92489 207251 435999 869500
0.312482933s

naiveな方はおよそ0.3秒でした.

λ> perfCheck (draw' 6 calc cols rows)

8 2 4 6 1 8 9 3 1 7
======================================================================
4| 12 14 18 24 25 33 42 45 46 53
2| 14 28 46 70 95 128 170 215 261 314
5| 19 47 93 163 258 386 556 771 1032 1346
6| 25 72 165 328 586 972 1528 2299 3331 4677
7| 32 104 269 597 1183 2155 3683 5982 9313 13990
1| 33 137 406 1003 2186 4341 8024 14006 23319 37309
3| 36 173 579 1582 3768 8109 16133 30139 53458 90767
9| 45 218 797 2379 6147 14256 30389 60528 113986 204753
3| 48 266 1063 3442 9589 23845 54234 114762 228748 433501
2| 50 316 1379 4821 14410 38255 92489 207251 435999 869500
0.003775219s
λ>

高速化した方は0.003秒でこのサイズだと100倍のパフォーマンス改善になりますが,

サイズを20x20とか増やしていけばそれ以上の効果となります.


関・ベルヌーイ数

関・ベルヌーイ数$B_{k}$は以下の等式を満たすような数列です.

\sum_{k=0}^{n} {}_{n+1}C_{k} B_{k} = n + 1

関・ベルヌーイ数を計算するのに,histomorphismを使うことにします.

アノテーションしておくのはより小さなすべての関・ベルヌーイ数です.


パスカルの三角形

累積百ます計算で実装したcalcを使ってパスカルの三角形のn段目の列を計算できます.


combination.hs

combs :: Integer -> [Integer]

combs n = map f (genIndex (fromIntegral n)) where f (x, y) = fromInteger (pascalTriangle !! x !! y)

pascalTriangle :: [[Integer]]
pascalTriangle = ones:(map (1:) (calc (ones, ones))) where ones = 1:ones



関・ベルヌーイ数の計算

上で実装したパスカルの三角形を計算する関数を使えば高速に関・ベルヌーイ数を計算できます.


bernoulli.hs

bernoulli :: Nat -> [Ratio Integer]

bernoulli = snd . histo phi
where
phi Z = (0, [1%1])
phi (S r) = (n+1, bs++[bn])
where
(n, bs) = extract r
ts = map fromInteger $ init $ combs (n+2)
bn = (fromInteger(n+2)-sum(zipWith (*) ts bs))*recip(ts!!fromInteger(n+1))

histomorphismの実装などはここではあらためて紹介しません.

コードの全体は以下です.

百ます計算

関・ベルヌーイ

表の表示関数

不動点となんとかもるふぃずむ群

これらのコードがimportしているものはaopに置いていますがちょくちょく手が入っていくのでなんかのタイミングではこの記事のコードと違ってる可能性もありますのでご容赦を.