Mastodon#gcc

#numericalanalysis

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

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

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

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

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

#julialang #gcc #dsfmt #openmp #fma

ちなみにマクロ定義を除いた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

#julialang #gcc #dsfmt #openmp #fma

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

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

mathtod.online/media/FNHimMxh0

mathtod.online/@tetu/413121

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

#julialang #gcc #dSFMT #OpenMP #FMA

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

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

Julia帝国軍 → 4.2秒

gcc dSFMT OpenMP FMA 連合軍 → 1.9秒

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

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

mathtod.online/media/DM-mn4Tze

#julialang #mingw64 #gcc #dSFMT #OpenMP #SSE2

二大怪獣最終決戦?😁

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

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

結果

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

julia @paralle → 4.3秒

Julia @parallel 強いです。

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

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

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

#julialang #python #numpy #gcc #gpp #openmp #dsfmt #mersennetwister

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

gist.github.com/genkuroki/740f

にあり、Python と Julia の Jupyter notebooks が

gist.github.com/genkuroki/45e4

にあります。

勉強になっています。

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

#julialang #python #numpy #gcc #gpp #openmp

【更新】私の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を凌駕できました!

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

#julialang #gcc #dsfmt #openmp

#julialang #python #numpy #gcc #gpp #openmp

【更新】私の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

#julialang #gcc #dsfmt

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

#julialang #python #numpy #gcc #gpp #openmp

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

Python + numpy → 2.3秒台

mingw64 gcc + Mersenne Twister → 0.8秒台

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

Julia (非並列化) → 0.4秒台

Julia (並列化) → 0.12秒

#julialang #gcc #mingw64 #pi #montecarlo #mersennetwister #gist

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

gist.github.com/genkuroki/d256

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

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

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

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

#julialang #gcc #mingw64 #pi #montecarlo #mersennetwister

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-