本の虫

著者:江添亮
ブログ: http://cpplover.blogspot.jp/
メール: boostcpp@gmail.com
Twitter: https://twitter.com/EzoeRyou
GitHub: https://github.com/EzoeRyou

アマゾンの江添のほしい物リストを著者に送るとブログ記事のネタになる

筆者にブログのネタになる品物を直接送りたい場合、住所をメールで質問してください。

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

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

コインを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より早いのではないかとも思ったが、まだ試していない。