なぜHaskellでタライ回しが速いのか、あるいは遅延評価とSTG


この記事は遅延評価を探る試みの一つです。遅延評価にはメリット, デメリット両方存在しており、この記事を根拠に遅延評価最高！と喧伝は出来ません。

以下、GHCは7.10.1です。


遅延評価と言えばタライ回しです。

http://ja.wikipedia.org/wiki/竹内関数

遅延評価を持つ言語でタライ回しはやたら速いと言われています。Haskellで書いてみましょう。

```haskell
import System.Environment (getArgs)

tarai :: Int -> Int -> Int -> Int
tarai x y z
    | x <= y    = y
    | otherwise = tarai
					(tarai (x-1) y z)
					(tarai (y-1) z x)
					(tarai (z-1) x y)
{-# NOINLINE tarai #-}

main :: IO ()
main = do
    args <- getArgs
    let (x:y:z:_) = fmap read args
	putStrLn $ "Result: " ++ (show $ tarai x y z)
```

`tarai`がINLINE化されてしまうことを防ぐために`NOINLINE`指定しています。後でSTGを見るためです。
特に最適化指定せずにコンパイルします。

```bash
$ ghc tarai.hs
```

実行します。

```bash
$ time ./tarai 12 6 0
Result: 12
./tarai 12 6 0  0.00s user 0.00s system 76% cpu 0.009 total
```

一瞬で終わります。タライ回しの言語比較は検索すれば幾つかヒットするのでここではやりません。
この異様なまでの速さがHaskellでのタライ回しです。


# タライは何回回されるのか？

Eagerな評価戦略を持つ言語だと、tarai関数は引数が(12, 6, 0)のときに12,604,860回呼ばれるようで、計算にもそれなりに時間が掛かります。
さてHaskellでは`tarai`は何回呼ばれているのでしょう。

カウントのために、`tarai`関数が呼ばれたときに標準エラー出力にログを吐くようにします。

```haskell
tarai :: Int -> Int -> Int -> Int
tarai x y z
    | trace ("(" ++ show x ++ "," ++ show y ++ ")") False = error "ここには来ない"
    | x <= y    = y
    | otherwise = tarai
					(tarai (x-1) y z)
					(tarai (y-1) z x)
					(tarai (z-1) x y)
```

繰り返しますが、この`trace`は`tarai`が呼ばれたときに、必ず実行されます(そうなるように仕込んでいます)。そのため、標準エラー出力に表示された行を数えれば`tarai`の呼び出し回数が分かるというわけです。
幸い`x`, `y`はすぐにthunkが潰されることが見えるため（`tarai`呼び出し直後に `x <= y` 比較しているためです。thunkを潰さないと値の比較は出来ません）、`x`, `y`を表示するようにします。
ここで`z`は表示出来ません。というのも、`z`を表示すると`z`のthunkが潰されてしまうからです。ログのために挙動が変わってしまうことは望ましくないので、ここでは表示しません。また、thunkの中を潰さずに覗く方法も恐らく現状はありません。
まあ`z`が表示出来ないと言っても呼び出しカウントには問題ありません。実行します。

```bash
$ ./tarai 12 6 0 |& wc -l
     110
```

どうやら110回程度しかタライは回されないようです(実際は1行カウントミスで、109回)。Eagerなときの12,604,860回に比べたらとても小さい数です。
数を大きくしてみます。

```bash
$ ./tarai 100 50 0 |& wc -l
    7502
```

7501回。それなりに大きくしたつもりですがこんなものか、といったところですね。
もちろん7500回程度の単純な計算など、現在のマシンならすぐに終わります。

```bash
$ time ./tarai 100 50 0
Result: 100
./tarai 100 50 0  0.00s user 0.00s system 70% cpu 0.011 total
```

ところでタライ回しを引数を変えて色々実行して見ると、呼び出し回数に規則が見えてきました。

(x, y, z) = (2N, N, 0)のときに、ちょうど(3 * N^2 + 1)回のようです。
(2N, N, 0)を指定している理由は、そのときにそこそこ回数が大きくなっているように思われるからです。
まあ計算量としては恐らくO(x^2)なわけです。これに意味はあるのでしょうか。


# タライはどのように回されるのか？

もう少しログ情報を加えてみます。具体的には `x <= y` で"DONE"と表示します。タライ回しが終わって`y`が返されることを表す印です。

```haskell
tarai :: Int -> Int -> Int -> Int
tarai x y z
	| x <= y && trace ("(" ++ show x ++ "," ++ show y ++ ") DONE") False = undefined
	| x > y && trace ("(" ++ show x ++ "," ++ show y ++ ")") False = undefined
    | x <= y    = y
    | otherwise = tarai
					(tarai (x-1) y z)
					(tarai (y-1) z x)
					(tarai (z-1) x y)
```

小さな数で実行します。

```bash
$ ./tarai 4 2 0
(4,2)
(3,2)
(2,2) DONE
(1,0)
(0,0) DONE
(-1,3) DONE
(0,3) DONE
(2,3) DONE
(1,0)
(0,0) DONE
(-1,4) DONE
(0,4) DONE
(3,4) DONE
Result: 4
```

とりあえず表示された数は`z-1`と`x`の間に収まっているようですね。
重複がないとかだと分かりやすいのですけどそういうわけでもないようです。
少し長くなりますが、もう少し大きい例を出してみます。

```bash
$ ./tarai 8 4 0
(8,4)
(7,4)
(6,4)
(5,4)
(4,4) DONE
(3,0)
(2,0)
(1,0)
(0,0) DONE
(-1,5) DONE
(0,5) DONE
(-1,5) DONE
(5,5) DONE
(-1,5) DONE
(5,5) DONE
(4,5) DONE
(3,0)
(2,0)
(1,0)
(0,0) DONE
(-1,6) DONE
(0,6) DONE
(-1,6) DONE
(6,6) DONE
(-1,6) DONE
(6,6) DONE
(5,6) DONE
(3,0)
(2,0)
(1,0)
(0,0) DONE
(-1,7) DONE
(0,7) DONE
(-1,7) DONE
(7,7) DONE
(-1,7) DONE
(7,7) DONE
(6,7) DONE
(3,0)
(2,0)
(1,0)
(0,0) DONE
(-1,8) DONE
(0,8) DONE
(-1,8) DONE
(8,8) DONE
(-1,8) DONE
(8,8) DONE
(7,8) DONE
Result: 8
```

DONEが定期的に連続して現れています。その連続回数は8で、`x`に一致するようです。
DONEのところだけ追っていくと、5...6...7...8と`y`が増えていって、つまり`y`が`x/2`から始まり`x`まで到達して終わる、という形になっているようです。先ほどの計算量O(x^2)ともなんとなく符号してますね。
あ、計算量の証明はしませんよ。


# サンク(Thunk)

Haskellは未評価の値を暗黙にサンク(thunk; thinkの過去分詞から)として扱っています。
これはEagerな評価戦略に慣れているプログラマにとっては大きな重荷になっているようです。

遅延評価はこのサンクによって実現しているわけですが、その実態はどうなっているのでしょうか。


# サンクが明示されるSTG

サンクがどのような表現でどのように使われているのか、それを確認するためにSTG(Spineless Tagless G-machine)を見てみます。GHCにおいて、コンパイル時にCoreの後に生成される形式です。
STGは実際の挙動との対応がだいぶ明確になっています。

* caseはevaluationと一対一対応
* letはheap allocationと一対一対応
* 全てのthunkが明示
* 型情報はおおよそ捨てられて、コード生成時に必要な分だけ残っている
* ClousureとThunkは共にlambda-formで表現される(引数のないlambda-formがthunk)
* 引数が埋まっている(saturated)関数は部分適用関数とは別に扱われる(Spineless)

参考:
https://ghc.haskell.org/trac/ghc/wiki/Commentary/Compiler/StgSynType
http://dev.stephendiehl.com/hask/#stg
http://www.scs.stanford.edu/11au-cs240h/notes/ghc.html

以下のように-ddump-stgオプションをつけて生成します。

```bash
$ ghc -ddump-stg tarai.hs > tarai.stg
```

`tarai`のSTGは以下です。

```haskell
tarai_rlW
  :: GHC.Types.Int -> GHC.Types.Int -> GHC.Types.Int -> GHC.Types.Int
[GblId, Arity=3, Str=DmdType, Unf=OtherCon []] =
    sat-only \r srt:SRT:[rlW :-> tarai_rlW,
                         ro4 :-> GHC.Classes.$fOrdInt, rAU :-> GHC.Num.$fNumInt] [x_s1Cx
                                                                                  y_s1Cy
                                                                                  z_s1Cz]
        case
            GHC.Classes.<= GHC.Classes.$fOrdInt x_s1Cx y_s1Cy
        of
        _ [Occ=Dead]
        { GHC.Types.False ->
              let {
                sat_s1CJ [Occ=Once] :: GHC.Types.Int
                [LclId, Str=DmdType] =
                    \u srt:SRT:[rlW :-> tarai_rlW, rAU :-> GHC.Num.$fNumInt] []
                        let {
                          sat_s1CI [Occ=Once] :: GHC.Types.Int
                          [LclId, Str=DmdType] =
                              \u srt:SRT:[rAU :-> GHC.Num.$fNumInt] []
                                  let {
                                    sat_s1CH [Occ=Once] :: GHC.Types.Int
                                    [LclId, Str=DmdType] =
                                        NO_CCS GHC.Types.I#! [1];
                                  } in  GHC.Num.- GHC.Num.$fNumInt z_s1Cz sat_s1CH;
                        } in  tarai_rlW sat_s1CI x_s1Cx y_s1Cy; } in
              let {
                sat_s1CG [Occ=Once] :: GHC.Types.Int
                [LclId, Str=DmdType] =
                    \u srt:SRT:[rlW :-> tarai_rlW, rAU :-> GHC.Num.$fNumInt] []
                        let {
                          sat_s1CF [Occ=Once] :: GHC.Types.Int
                          [LclId, Str=DmdType] =
                              \u srt:SRT:[rAU :-> GHC.Num.$fNumInt] []
                                  let {
                                    sat_s1CE [Occ=Once] :: GHC.Types.Int
                                    [LclId, Str=DmdType] =
                                        NO_CCS GHC.Types.I#! [1];
                                  } in  GHC.Num.- GHC.Num.$fNumInt y_s1Cy sat_s1CE;
                        } in  tarai_rlW sat_s1CF z_s1Cz x_s1Cx; } in
              let {
                sat_s1CD [Occ=Once] :: GHC.Types.Int
                [LclId, Str=DmdType] =
                    \u srt:SRT:[rlW :-> tarai_rlW, rAU :-> GHC.Num.$fNumInt] []
                        let {
                          sat_s1CC [Occ=Once] :: GHC.Types.Int
                          [LclId, Str=DmdType] =
                              \u srt:SRT:[rAU :-> GHC.Num.$fNumInt] []
                                  let {
                                    sat_s1CB [Occ=Once] :: GHC.Types.Int
                                    [LclId, Str=DmdType] =
                                        NO_CCS GHC.Types.I#! [1];
                                  } in  GHC.Num.- GHC.Num.$fNumInt x_s1Cx sat_s1CB;
                        } in  tarai_rlW sat_s1CC y_s1Cy z_s1Cz;
              } in  tarai_rlW sat_s1CD sat_s1CG sat_s1CJ;
          GHC.Types.True -> y_s1Cy;
        };
```

lambda-formを確認しておきます。lambda-formは \ (backslash)で表されています。
おおよそ以下です。

```haskell
\  -- Represent Lambda-form
	r -- UpdateFlag (u | r | s)
	srt:SRT:[...] -- (Static Reference Table)
	[x_s1Cx y_s1Cy z_s1Cz] -- (Arguments)
	let = ...  -- body
```

\直後にuとかrとかがあって、SRTがあって、引数が並んでて、その後にbodyが続くと。

さて上記taraiのSTGは良く見ると、ほぼ`tarai`そのままであることが分かります。
分かりやすく簡略化して、変数名も変えてみましょう。

```haskell
tarai_rlW :: Int -> Int -> Int -> Int [...] =
    sat-only \r srt:SRT:[...] [arg_x arg_y arg_z]
        case
            <= $fOrdInt arg_x arg_y -- x <= y (STGでは演算子は中置ではなく前置)
        of
        _ [Occ=Dead]
        { False ->
              let {
                thunk1 [Occ=Once] :: Int [...] =
                    \u srt:SRT:[...] [] -- lambda-formで引数がないので、これはthunk
                        let { thunk1' = ... } in  tarai_rlW thunk1' arg_x arg_y; } in
              let {
                thunk2 [Occ=Once] :: Int [...] =
                    \u srt:SRT:[...] []
                        let { thunk2' = ... } in  tarai_rlW thunk2' arg_z arg_x; } in
              let {
                thunk3 [Occ=Once] :: Int [...] =
                    \u srt:SRT:[...] []
                        let { thunk3' = ... } in  tarai_rlW thunk3' arg_y arg_z;
              } in  tarai_rlW thunk3 thunk2 thunk1; -- thunkを渡しつつ再帰

          True -> arg_y; -- yを返すだけ
        };
```

はじめに引数の`x`と`y`を比較して、その結果がTrueなら`y`をそのまま返す。
Falseなら`thunk1`, `thunk2`, `thunk3`を作ってからそれらthunkを`tarai`に渡す。`tarai`そのままですね。

注目する点は、関数適用時にthunkを作って、それを関数の返り値の代わりに扱っていることでしょう。
そしてそのthunkが潰されるのはcaseの時ですが、それに関してはSTGではわからないようです。まあthunkか否かの違いはSTG上では分からないようなので仕方がないですかね。グラフ簡約もしますし、そちらはランタイムの仕事ということなのでしょう。


# サンクの正体と扱い方

STGを確認して、ThunkがClosureの特別な場合であることを確認しました。これはGHCの場合の実装であって、一般的な実装ではないかもしれませんが。
ここまででthunkに関して幾つかの挙動が確認出来たのでまとめておきます。

* 関数適用時はすぐに評価せずにthunkを作る。関数適用した値を関数に渡すときは、thunkを作ってそのthunkを渡す。
* 関数にthunkを渡すということは、関数が受け取る引数はthunkが積んである可能性があるということ。
* thunkはcase時に評価される。評価はWHNFまで、つまりはConstructorかLamdaが見えるまで。分岐でどちらへ行けばいいか判断ができるまで。

細かい挙動についてはまだまだ抑えるべき点はあるかもしれませんが、とりあえずはこれだけわかっておけばいいのではないでしょうか。

実際にサンクの正体を知っておくと、遅延評価を扱う際に安心できるのではないでしょうか。もし詳しく知りたければ今回のようにSTGを見ればいいでしょう。何回か見れば慣れるものです。

# 結局ナンデ速かったのか

さて`u`, `r`といったUpdateFlagですが、今回は関係ありません。`u`(updatable)のthunkは他の同じ計算が終わると、その計算結果でthunkが置き換わります。なのでupdatableのthunkを大量に作りつつ、一つが終わったら大量のthunkがぶわーっと更新されるのかなーとか少し期待しないでもなかったのですが、そんなことはありませんでした。
実際にいくらか手動でタライを回してみたところ、`tarai`の速さは引数`z`上に積まれているthunkを一切計算しなくても答えが出ることが原因のようです。まあ`x <= y`調べて`y`返しているだけですからね。

結局、`z`をeagerに計算するか否かでこれだけの呼び出し回数の差が出ていたということでした。
言ってみれば、「`tarai`が遅延評価で速いのは、遅延評価で速くなるような問題だから」ということです。「タライ回しが速いから遅延評価最高！」とはならないというのはそういうことです。