Re:Haskellで書いてみたらC++の10倍遅かった 5倍程度になりました
10倍は遅すぎませんか?
本の虫: Haskellで書いてみたらC++の10倍遅かった
を読みました.
10倍は差が出過ぎなのではないかと思いました.
C++ソースコード
まずC++のソースコードが完全なものではないので補完しました.
#include <algorithm>
#include <array>
#include <iostream>
#include <random>
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";
}
}
2018-02-06T12:05:17 ncaq@strawberry/pts/4(1) ~/Documents/archive/2018-02
% g++ -O2 -std=c++17 coin.cpp
2018-02-06T12:05:30 ncaq@strawberry/pts/4(0) ~/Documents/archive/2018-02
% time ./a.out
100 : 14%
1000 : 8.3%
10000 : 8.8%
100000 : 8.557%
./a.out 0.21s user 0.00s system 99% cpu 0.217 total
Haskellソースコードにhlintをかけたもの
Haskellソースコードを写して, 元のソースコードが読みにくいのでhlintに従って整形しました.
import Data.List
import System.Random
coinSeq :: (RandomGen g) => g -> [Bool]
coinSeq = randoms
splitN :: Int -> [a] -> [[a]]
splitN n s = take n s : splitN n (drop n s)
hasContiguousElems :: (Eq a) => Int -> [a] -> Bool
hasContiguousElems len s = any (>= len) (map length (group s))
coinToss :: Int -> [Bool] -> Bool
coinToss = hasContiguousElems
monteCarlo :: Int -> [Bool] -> Double
monteCarlo try_count s = (fromIntegral n / fromIntegral try_count) * 100.0
where
seq_n = take try_count (splitN 100 s)
n = length . filter id $ map (coinToss 10) seq_n
main :: IO ()
main = do
gen <- getStdGen
let s = coinSeq gen
mapM_ (\n -> putStrLn (show n ++ "\t : " ++ show (monteCarlo n s) ++ "%") )
[100, 1000, 10000, 100000]
2018-02-06T13:43:57 ncaq@strawberry/pts/4(130) ~/Documents/archive/2018-02
% stack ghc -- -O2 oldCoin.hs
2018-02-06T13:43:59 ncaq@strawberry/pts/4(0) ~/Documents/archive/2018-02
% time ./oldCoin
100 : 10.0%
1000 : 8.1%
10000 : 8.780000000000001%
100000 : 8.559999999999999%
./oldCoin 3.23s user 0.03s system 99% cpu 3.268 total
なるほど確かに10倍程度遅いようですね.
プロファイリング
ところでghcにはプロファイリングオプションがあります.
計測せよって偉い人も言ってました.
2018-02-06T12:11:55 ncaq@strawberry/pts/6(0) ~/Documents/archive/2018-02
% stack ghc -- -prof -fprof-auto -rtsopts coin.hs
[1 of 1] Compiling Main ( coin.hs, coin.o )
Linking coin ...
2018-02-06T12:13:11 ncaq@strawberry/pts/6(0) ~/Documents/archive/2018-02
% time ./coin +RTS -p
100 : 91.0%
1000 : 91.2%
10000 : 91.43%
100000 : 91.259%
./coin +RTS -p 10.76s user 0.11s system 99% cpu 10.882 total
Tue Feb 6 12:13 2018 Time and Allocation Profiling Report (Final)
coin +RTS -p -RTS
total time = 10.15 secs (10151 ticks @ 1000 us, 1 processor)
total alloc = 11,164,362,552 bytes (excludes profiling overheads)
COST CENTRE MODULE SRC %time %alloc
hasContiguousElems Main coin.hs:11:1-68 21.2 25.6
randomIvalInteger System.Random System/Random.hs:(468,1)-(489,76) 21.0 26.5
stdNext System.Random System/Random.hs:(518,1)-(528,64) 13.8 15.8
randomIvalInteger.f System.Random System/Random.hs:(486,8)-(489,76) 9.3 2.1
randomIvalInteger.f.v' System.Random System/Random.hs:489:25-76 5.3 4.3
randomIvalInteger.b System.Random System/Random.hs:473:8-54 4.0 5.7
randomIvalInteger.(...) System.Random System/Random.hs:472:8-36 3.8 0.0
randomIvalInteger.magtgt System.Random System/Random.hs:483:8-21 3.7 1.4
randomR System.Random System/Random.hs:(387,3)-(397,30) 2.9 0.0
next System.Random System/Random.hs:218:3-17 2.8 2.1
splitN Main coin.hs:8:1-43 2.7 5.4
randoms.\ System.Random System/Random.hs:316:42-67 2.4 5.0
randomIvalInteger.f.(...) System.Random System/Random.hs:488:25-39 2.3 0.0
randomIvalInteger.k System.Random System/Random.hs:482:8-20 1.8 1.4
random System.Random System/Random.hs:399:3-49 1.1 0.0
stdNext.z System.Random System/Random.hs:520:17-34 0.6 1.4
stdNext.s1'' System.Random System/Random.hs:524:17-64 0.3 1.4
stdNext.s2'' System.Random System/Random.hs:528:17-64 0.2 1.4
やっぱり乱数生成が遅い
あーやっぱりRandomが遅い. そんな気はしていました.
Haskellの乱数事情 - Qiita で触れられてるとおり, System.Randomは遅いのです.
乱数を速度を重視して真面目に考えたことがないので何を選択すれば良いのか全くわからない. 先程の記事で推奨されていたmwc-randomを使ってみます.
一番速いものを選ぼうかと思ったのですが, 純Haskellで実装されたものでないと「Haskellの速度」を測るレギュレーションに違反しているのではないかと思ったので, 純Haskellで実装されたものを選びました.
System.Random.MWCでの再実装
import Data.List
import qualified Data.Vector as V
import System.Random.MWC
coinSeq :: Int -> IO (V.Vector Bool)
coinSeq l = withSystemRandom . asGenST $ \gen -> uniformVector gen l
splitN :: Int -> [a] -> [[a]]
splitN n s = take n s : splitN n (drop n s)
hasContiguousElems :: (Eq a) => Int -> [a] -> Bool
hasContiguousElems len s = any (>= len) (map length (group s))
coinToss :: Int -> [Bool] -> Bool
coinToss = hasContiguousElems
monteCarlo :: Int -> [Bool] -> Double
monteCarlo try_count s = (fromIntegral n / fromIntegral try_count) * 100.0
where
seq_n = take try_count (splitN 100 s)
n = length . filter id $ map (coinToss 10) seq_n
main :: IO ()
main = do
s <- V.toList <$> coinSeq 100000
mapM_ (\n -> putStrLn (show n ++ "\t : " ++ show (monteCarlo n s) ++ "%") )
[100, 1000, 10000, 100000]
2018-02-06T13:08:09 ncaq@strawberry/pts/6(0) ~/Documents/archive/2018-02
% stack ghc -- -O2 coin.hs
2018-02-06T13:08:13 ncaq@strawberry/pts/6(0) ~/Documents/archive/2018-02
% time ./coin
100 : 9.0%
1000 : 9.2%
10000 : 0.9199999999999999%
100000 : 9.2e-2%
./coin 0.03s user 0.00s system 98% cpu 0.038 total
これでC++の2倍程度の遅さに収まりましたね. 良かったですね.
と言いたいところですが, 乱数が偏りすぎていてまともな答えが帰ってきていない. これではコードがちゃんと動いたとは言い難いです.
まともな偏りを持った純Haskell製乱数生成ライブラリは少しだけ調べた結果見つかりませんでした. 実用的には他言語へのバインディングライブラリ使えば良いんでしょうけど, 今回のケースを解決するものは見つかりませんでした. 残念です.
答えとしては片手落ちですね. 誰か知ってる人が居たら教えてください.
そもそも無限リストなんかじゃなくてRandomIO
とか使って必要なところに乱数を使えばそこそこ早くなるんじゃないのと思ったのですがとても眠いので大幅な書き換えをやる気が起きません.
指摘してもらいました
MWC乱数を使ったコードにバグがあるようです。準備する乱数が、100000個しかないので、100回のコイントスは1000回までしか試行
— Nobuo Yamashita (@nobsun) 2018年2月6日
できません。1000回まではいいのですが、10000回だと 1/10 の値、100000回だと 1/100 の値が出てしまいます。
なるほど.
元のアルゴリズムを全く読まずに適当な書き換えをしたツケですね.
乱数に偏りがあったという私の理解は誤りでした.
修正版
import Data.List
import qualified Data.Vector as V
import System.Random.MWC
coinSeq :: Int -> IO (V.Vector Bool)
coinSeq l = withSystemRandom . asGenST $ \gen -> uniformVector gen l
splitN :: Int -> [a] -> [[a]]
splitN n s = take n s : splitN n (drop n s)
hasContiguousElems :: (Eq a) => Int -> [a] -> Bool
hasContiguousElems len s = any (>= len) (map length (group s))
coinToss :: Int -> [Bool] -> Bool
coinToss = hasContiguousElems
monteCarlo :: Int -> [Bool] -> Double
monteCarlo try_count s = (fromIntegral n / fromIntegral try_count) * 100.0
where
seq_n = take try_count (splitN 100 s)
n = length . filter id $ map (coinToss 10) seq_n
main :: IO ()
main = do
s <- V.toList <$> coinSeq (100000 * 100)
mapM_ (\n -> putStrLn (show n ++ "\t : " ++ show (monteCarlo n s) ++ "%") )
[100, 1000, 10000, 100000]
2018-02-06T16:54:40 ncaq@strawberry/pts/9(0) ~/Documents/archive/2018-02
% time ./coin
100 : 8.0%
1000 : 7.7%
10000 : 8.959999999999999%
100000 : 8.698%
./coin 1.09s user 0.04s system 99% cpu 1.133 total
これで1.133 / 0.217 = 5.221198156682028というわけで, 5倍程度の速さ(遅さ)になりました.
リストじゃなくてVectorのまま取り扱ったらもっと早くなるだろうって? それは読者の宿題としようと思います(めんどくさい).