Mathtodon#gcc

Juliaとgccの比較について、円周率モンテカルロのような単純計算ではなく、実際に役に立つ数値計算で比較したらどうなるかに興味がある人は多いと思います。

インターネット上で見付かる指数積分函数 Ei(x) のCによる数値計算のコードをほとんどそのままJuliaに忠実に翻訳したものと、もとのCのコードをmingwのgccでコンパイルしたものとの速度比較を実行してみました。

結論:計算の速さはほぼ同じ。

gist.github.com/genkuroki/dabc

JuliaをCの代わりに用いても数値計算の速さ的には何も問題がないというのが結論。

github.com/stevengj/18S096-iap
のような最適化は行っていません。

mathtod.online/media/JfVZMKAwD

gcc連合軍とJulia帝国軍の戦いの記録を1枚のファイルにしました。

gist.github.com/genkuroki/c9ae

例の円周率モンテカルロによる比較です。

mathtod.online/media/76DAtpHMd

ローカルタイムラインにHaskellという文字列を見たので、てもとのWindows機にまだ入れていなかったことに気付き、Stack経由でインストール中。

stack setup が進行中。ダウンロードがなぜか滅茶苦茶遅い。23:20に作業をスタートしたのにまだ終わらない。😭

もしかして間違った道に入り込んでいる?

確か、GHC (Haskellコンパイラー)も激速だったはず。例の単純な円周率モンテカルロを試してみたい。誰かもうやった人っている?40億回のループに何秒かかりましたか?

私のパソコンだとgccで2秒台、Juliaで4秒台です。(Juliaでも擬似乱数生成を別のに変えれば3秒台にまで縮まる。)

我々のmathtodonにおける経験で判明したことは、 gcc でJuliaに速度面で勝つためには相当に非自明なコードを書かなければいけないということ。

そして、Cythonの速さは実質的に C の速さです。

非自明なコードを C で直接書かずに、Cython で単純なコードをコンパイルして使っても、Juliaに速度面で勝てないのは予想通りの結果でした。

ところで、

Julia側からPythonのライブラリを使うのは易しいです。ほとんどPythonの存在を意識せずに使えるくらい簡単。

Python側からJuliaを気軽に呼び出して使える機能はすでに実装されているのでしょうか?

ちょっと検索しても見つからない。キーワードがわからない。

(1) Juliaでは、gccを呼び出してCで書かれた函数をコンパイルしてJuliaの函数として利用できる。

(2) JuliaからPythonのライブラリを

using PyCall
@pyimport scipy.special as s

のように書いて利用できる。(SciPyの特殊函数ライブラリ(十分に高速)をこれで利用できる。)

matplotlib.pyplot をJuliaユーザーは最も気楽に使えるplotting backend として利用している。

SymPyもJuliaで利用できる。プログラムを書くときに、多項式やTaylor展開などの面倒な計算をしなければいけないことが結構あるので、これはとてもありがたい。

(3) JuliaではRも呼び出せる。Rのコードも

R"mod <- lm(y ~ X-1)"

のように書いて実行できる。

luiarthur.github.io/usingrcall

Rでやった方が楽なことはRですませるというようなことも簡単にできます。

github.com/stevengj/18S096-iap
の講義録は滅茶苦茶参考なりそうな内容でした。講義の主題はどうも「高速計算」のようです。

そうか、JuliaとgccとPythonの速度比較をやるには、すべてをJuliaから呼び出してやればよかったのか。

github.com/stevengj/18S096-iap
の最初の部分に、Juliaからgccを呼び出して函数をコンパイルして、それをそのままJuliaで利用する方法が書いてあります。

わざわざ Makefile を書いて試したのはちょっとアホだったかも。

Jupyter notebookの形式になっていた方が取り扱いが便利だった。

複数の言語を一つの言語から自由に利用できるのはやはりとてもよいですね。

問題:次の条件収束する交代級数の値を小数点以下第10桁まで求めよ:$$
\alpha = \sum_{n=1}^\infty \frac{(-1)^{n-1}}{\sqrt{n}}.
$$
すでに我々は で並列処理を行えば、軽く1分以内に $10^{10}$ 個の和を直接的に足し上げられそうなことを知っています。😊

実際にそれをやってみれば、$10^{2n}$ 項の和で小数点以下 $n$ 桁まで求まっていそうなことがわかります。😅

だから小数点以下10桁まで求めるためには $10^{20}$ 個の項の和を取らなければいけなさそうです。😥

さすがにそれは苦し過ぎ。それではどうすればよいかというのが問題です。

たぶん、答えを知っている人はたくさんいると思います。

ちなみにマクロ定義を除いたJulia側のコードは次の一行だけです。

findpi_macrosum(n::Int) = (4.0/n) * (@sum i 1:n ifelse(rand()^2 + rand()^2 <= 1.0, 1, 0))

この1行だけで円周率をモンテカルロで計算できる。

gist.github.com/genkuroki/45e4

OMP_NUM_THREADS の値によってどれだけ速さが変わるかに興味がある人は添付画像を見て下さい。

CPU: Intel(R) Core(TM) i7-4770HQ CPU @ 2.20GHz

mathtod.online/media/FNHimMxh0

mathtod.online/@tetu/413121

tetuさん、お世話になっています。
ついに、ついに、ついに、………

Julia と gcc の死闘、ついに終止符が打たれるのか?!?!?!

40億回のモンテカルロループにかかった時間

Julia帝国軍 → 4.2秒

gcc dSFMT OpenMP FMA 連合軍 → 1.9秒

ついに、gcc連合軍の圧倒的勝利!

しかし、たくさんの血が流された……………

mathtod.online/media/DM-mn4Tze

二大怪獣最終決戦?😁

10億回のループでは短い感じなので、40億回にしてみました。

適当にあいだをあけながら二三十回試して一番短い方の時間をピックアップしました。

結果

mingw64 gcc with dSFMT, SSE2, and OpenMP → 5.3秒

julia @paralle → 4.3秒

Julia @parallel 強いです。

なかなか倒せない。数字がばらつくのですが、安定してJuliaの方が速い感じです。

統計処理をしてみようかと思ったのですが、色々面倒なのでやめました。

CPUの温度が計算にかかる時間のゆらぎに関係しているっぽい雰囲気。 mathtod.online/media/ga6fhfgzu

使用したCとC++のコードは

gist.github.com/genkuroki/740f

にあり、Python と Julia の Jupyter notebooks が

gist.github.com/genkuroki/45e4

にあります。

勉強になっています。

みなさん、本当にどうもありがとうございます。

【更新】私のWindows環境での現時点での最高速のまとめ【更新】

ループを1億から10億に増やしました。

モンテカルロ法による円周率の計算でループを10億回まわした場合

g++ → 26秒台

gcc default rand() → 約25秒

Python & numpy → 約22.5秒

gcc Mersenne Twister → 約8.5秒

g++ OpenMPで並列化 → 約4秒

Julia → 約4秒

gcc dSFMT → 約3.2秒
ついにJuliaを凌駕した!!!

gcc dSFMT OpenMPで並列化 → 1.3秒台
Juliaの並列化のケースとほぼ互角だと言ってよいでしょう!

Juliaで並列化 → 1.2秒台

今回の件で、モンテカルロシミュレーションで乱数生成が高速であることの重要性を初めて理解できました。

@tetu

すんばらしいです!!!とても速いです!

それでも微妙に並列処理版のJuliaよりちょっと遅い。

そこで、非OpenMP版も作って試してみました。(普段使うときには並列処理に書き直す手間をかけない。)

Juliaの非並列処理版よりも速かったです!!!

ついに普段使っている状況でJuliaを凌駕できました!

もう少ししたら、返答を外して、詳しい情報を公開します。

【更新】私のWindows環境での現時点での最高速のまとめ【更新】

円周率を求めるモンテカルロの1億回ループにかかる時間

g++ → 2.6秒台

Python & numpy → 2.3秒台

g++ & OpenMPで並列化 → 0.4秒台

gcc & dSFMT → 0.4秒台

Julia → 0.4秒台

Juliaで並列化 → 0.12秒

mathtod.online/@waidotto/40984

y.さんのおかげで、gccでJuliaと同じ速さ、円周率モンテカルロをできるようになりました!

gcc で並列化を使わずに0.4秒台が出ました!

私が追試するときに使ったコードとコンパイルの仕方(Makefile)は次の場所にあります。

gist.github.com/genkuroki/67a5

findpi_dSFMT.c のコンパイルのためには
math.sci.hiroshima-u.ac.jp/~m-
から dSFMT-src-2.2.3 をダウンロードして展開しておく必要があります。

mathtod.online/media/-8HAEg24d

私のWindows環境での最高速をまとめておきます。円周率のモンテカルロの1億回ループにかかる時間です。

Python + numpy → 2.3秒台

mingw64 gcc + Mersenne Twister → 0.8秒台

mingw64 g++ + OpenMP (並列化) → 0.4秒台

Julia (非並列化) → 0.4秒台

Julia (並列化) → 0.12秒

mingw64 gcc -O2 用のMakefile, *.c などを次の場所に置いておきました。

gist.github.com/genkuroki/d256

C コンパイラーを使った円周率モンテカルロで Julia よりも速く計算できた人がいればそのやり方を教えて下さい。(できれば誰でも無料で使える方法を希望。Juliaは無料なので。)

計算の速さが欲しければ、gccなどを使ったりせずに、Julia言語を使う方がずっと楽そうですね。

必要なら並列処理のコードも簡単に書けるし。

コードを書くのが楽なのは実際に使ってみると滅茶苦茶ありがたいです。

Juliaでは軽く1秒を切る(私の環境では0.4秒台)のに、Cではそれよりずっと遅い。

私も手元のWindows機の mingw64 gcc -O2 で試してみました。2.5秒以上かかる。

for ループ内では整数しか使わないようにしても2秒をやっと切る程度。

Juliaの0.4秒台に適わない。

そこでMersennne Twisterのコードを

math.sci.hiroshima-u.ac.jp/~m-

からコピペして試してみました。

すると、0.8秒台まで計算時間が縮まりました。Juliaの半分程度の速さ。

Cの側がJuliaとオーダーが違うほど遅くなっていた原因は、rand() 函数が滅茶苦茶遅いことだったようです。

finfpi_mt がメルセンヌツイスター版

mathtod.online/media/WZv0SGBb-