Haskellで書いてみたらC++の１０倍遅かった

エレガントな解法、エレファントな解法 〜モンテカルロ法を添えて〜

コインを100回投げたうちに、表もしくは裏が連続して10回でる確率をモンテカルロ法で計算する。

なにはともかくC++で実装してみた。

template < typename Random >
bool coin_toss( const unsigned int try_count, const unsigned int length, Random & r )
{
    unsigned int count{} ;
    int prev{} ;

    std::uniform_int_distribution<> d( 0, 1 ) ;

    for ( unsigned int i = 0 ; i != try_count ; ++i )
    {
        int result = d(r) ;
        if ( result == prev )
        {
            ++count ;
            if ( count == length )
                return true ;
        }
        else
        {
            prev = result ;
            count = 1 ;
        }
    }

    return false ;
}

template < typename Random >
double monte_carlo( unsigned int try_count, Random & r )
{
    unsigned int count{} ;
    for ( unsigned int i = 0 ; i != try_count ; ++i )
    {
        if ( coin_toss( 100, 10, r ) )
            ++count ;
    }

    return double(count) / double(try_count) ;
}

int main()
{
    std::array<unsigned int, sizeof(std::mt19937)/sizeof(unsigned int)> c ; ;
    std::generate( std::begin(c), std::end(c), std::random_device{} ) ;
    std::seed_seq s( std::begin(c), std::end(c) ) ;
    std::mt19937 engine(s) ;

    for ( unsigned int i = 100 ; i != 1000000 ; i *= 10 )
    {
        auto r = engine ;
        std::cout << i << "\t : " << monte_carlo( i, r ) * 100.0 << "%\n" ;
    } 
}

C++11から高度な乱数ライブラリが使えるようになったのだが、現在、実行時に適切に初期化する方法に面倒なボイラープレートコードが必要だ。実行時に乱数で初期化してくれる機能が現在提案されている。

ところで、最近Haskellを学んでいるので、このコードをHaskellでも書いてみようと思った。

まずは乱数だ。コインの表裏はBoolで表現するとして、Haskellらしく無限の乱数リストを生成しよう。

coin_seq :: (RandomGen g) => g -> [Bool]
coin_seq gen = randoms gen

あとはこれをtake 100して一回分の試行とする。1000回試行したければ、drop 100したのを再帰的にtake 100すればよい。しかし、もっといい方法がある。要素数100の[Bool]が入ったリスト、つまり[[Bool]]を無限に作ればいい。そのリストをtake 1000すれば1000回の試行分になる。

split_n :: Int -> [a] -> [[a]]
split_n n seq = take n seq : split_n n (drop n seq)

さて、seq_n = take 100 $ split_n 100 (coin_seq gen)のようなリストseq_n :: [[Bool]]は作れた。あとは[Bool]に連続した要素がn個あるかどうかを調べるcoin_toss nと、それをリストのすべての要素に対して行うmonte_carloを書けばいいだけだ。

coin_tossが結果をBoolで返すとすると、とりあえずseq_nにcoin_tossをmapしてTrueだけfilterして、その結果のリストの要素数を数えれば、成功した試行数がわかる。あとは試行数で割って確率を求めればいいだけだ。

monte_carlo :: Int -> [Bool] -> Double
monte_carlo try_count seq = ((fromIntegral n) / (fromIntegral try_count)) * 100.0
    where
        seq_n = take try_count (split_n 100 seq)
        n = length . filter id $ map (coin_toss 10) seq_n

coin_toss len seqはInt -> [Bool] -> Boolな関数で、seqの中にlen個の連続したTrueもしくはFalseがあるかどうかを数えればよい。

いろいろ考えたが、まずgroup seqして連続した同じ要素を分割して[[Bool]]を作り、それぞれの[Bool]の要素数をいれたリスト[Int]をmapし、len以上の要素だけをfilterして、結果のリストに要素があればTrueを返す計算をすればいいのではないか。つまり、

has_contiguous_elems :: (Eq a) => Int -> [a] -> Bool 
has_contiguous_elems len seq = not $ null $ filter (>= len) $ map length (group seq)

coin_toss :: Int -> [Bool] -> Bool
coint_toss _ [] = 0
coin_toss len seq = has_contiguous_elems len seq

となる。全体的にはこうだ。

import System.Random
import Data.List

coin_seq :: (RandomGen g) => g -> [Bool]
coin_seq gen = randoms gen

split_n :: Int -> [a] -> [[a]]
split_n n seq = take n seq : split_n n (drop n seq)

has_contiguous_elems :: (Eq a) => Int -> [a] -> Bool 
has_contiguous_elems len seq = not $ null $ filter (>= len) $ map length (group seq)

coin_toss :: Int -> [Bool] -> Bool
coint_toss _ [] = 0
coin_toss len seq = has_contiguous_elems len seq

monte_carlo :: Int -> [Bool] -> Double
monte_carlo try_count seq = ((fromIntegral n) / (fromIntegral try_count)) * 100.0
    where
        seq_n = take try_count (split_n 100 seq)
        n = length . filter id $ map (coin_toss 10) seq_n

main = do
    gen <- getStdGen
    let s = coin_seq gen
        in mapM (\n -> putStrLn ((show n) ++ "\t : " ++ (show (monte_carlo n s)) ++ "%") )
                [100, 1000, 10000, 100000] 

このコードは動いた。ただし、最適化オプションを使っても、C++の10倍遅かった。

やはりlen個の連続した要素が存在するかどうかを調べるhas_contiguous_elemsが遅いのではないかと思い、これを再帰で実装することにした。

has_contiguous_elems' :: (Eq a) => Int -> [a] -> Bool 
has_contiguous_elems' _ [] = False
has_contiguous_elems' len (x:xs) =
    if count >= len
        then True
        else has_contiguous_elems' len $ dropWhile (== x) xs
    where count =  1 + (length $ takeWhile (== x) xs)

この実装に差し替えたところ、実行速度はC++の9倍遅かった。

すると、他の部分もリストからリストを生成するという記述をやめて、再帰で実装すればもう少し早くなるのだろうか。

そして思うのは、同様のリストからリストを清々するような処理は、所詮リストのメモリーサイズがキャッシュに収まる程度なので、C++で書くとHaskellより早いのではないかとも思ったが、まだ試していない。