黒木玄 Gen Kuroki is a user on mathtod.online. You can follow them or interact with them if you have an account anywhere in the fediverse. If you don't, you can sign up here.

mathtod.online/@mashiro/452539

Pythonがものすごく遅い原因は、あの極めて豊富なライブラリ群(数値計算用のライブラリはCやC++で書かれており十分高速)ではないので、Pythonの数値計算用のライブラリを使う分には速度面で失うものは少ないと思います。

Pythonはforループが滅茶苦茶遅い。

forループを使わないコードを書こうとすると、大きめの配列を確保するようなコードになってしまうことが多く、配列の確保も結構重いので、高速なライブラリ群を使っていても、Pythonだと遅くなってしまうわけです。(Juliaでも単に 1:n と書くべきところを [1:n;]とかcollect(1:n)と書くと遅くなる。)

Pythonでも高速なライブラリ群を活かして気楽に(←これ決定的に重要)高速なコードを書く方法があるなら教えて欲しいところ。

すでに公開済みの

gist.github.com/genkuroki/cde6

では、Juliaで指数積分函数を使うための手段として

(1) GNU MPFR ライブラリをccall経由で使う。

(2) Python SciPyの特殊函数ライブラリを使う。

の2つの手段を試してみて、速度比較もしています。

GNU MPFRの側は256bit浮動小数点なのでさすがに遅かったです。

素数定理の確認などではそこまでの精度は必要ないので、より短時間で計算できる Python SciPyライブラリを使った方が便利です。

学生の人なんかが、Julia用に高速で便利な指数積分函数のコードをを書けば、ちょっとした自慢になると思います。

現時点では100万以下の素数について素数定理が数値的に成立していることを確認するときの律速段階は対数積分函数を計算するところです。

mathtod.online/media/Lsw7CIGfe

気楽に色々計算して遊びたい場合には、気楽にコードを書けることが必須。

だからもしもJuliaからPythonライブラリを使うことが気楽にできないなら、Pythonの豊富なライブラリを使用可能でも、かなりがっかりだということになります。

現実にはほとんどPythonを使っていることを意識せずにPythonライブラリを使えます。atマークを{at}と書くことにすると、

using PyCall
{at}pyimport scipy.special as s

と書くだけで、SciPyの特殊函数ライブラリをJuliaでそのまま使えるようになります。

しかもMPFRで高精度計算する選択肢を使うよりも速い。

こういうのはとてもありがたいです。

もしも、SciPyの特殊函数ライブラリよりも速い指数積分函数をJuliaで書けるなら、誰か書くとよいと思います。

計算の速さは実際にコードを書いて実行してみないとわからないことが多いと思う。

がたがた言う前に実際にコードを書いて実行してみるべきだと思う。そして使ったコードは公開する。

Juliaで対数指数函数を使うための二つの選択肢の比較のノートブックを更新しました。

gist.github.com/genkuroki/cde6

GNU MPFRのBigFloatによる計算はsetprecision()で精度を変更できるので、SciPy側と精度を合わせて速度比較をしてみました。

それでも私のコードだと、MPFR etan()を使うよりも、SciPy側のscipy.special.expi()を使った方が倍くらい速いです。

私の現時点での実力では、PythonのSciPyライブラリを使うことは速度面においても十分に魅力的な選択肢です。

mathtod.online/media/tVGxfEs5x

黒木玄 Gen Kuroki @genkuroki

GNU Scientific Libraryのマニュアルの指数積分函数の部分へのリンク

gnu.org/software/gsl/doc/html/

積分の書き方が $dt$ を左側に書くスタイルになっている:
\[
\mathrm{Ei}(x)
= -PV\left(\int_{-x}^\infty dt\,\exp(-t)/t\right)
\]

SciPy側の対応する函数のマニュアル

docs.scipy.org/doc/scipy/refer

· Web · 0 · 1