https://mathtod.online/@mashiro/452539
Pythonがものすごく遅い原因は、あの極めて豊富なライブラリ群(数値計算用のライブラリはCやC++で書かれており十分高速)ではないので、Pythonの数値計算用のライブラリを使う分には速度面で失うものは少ないと思います。
Pythonはforループが滅茶苦茶遅い。
forループを使わないコードを書こうとすると、大きめの配列を確保するようなコードになってしまうことが多く、配列の確保も結構重いので、高速なライブラリ群を使っていても、Pythonだと遅くなってしまうわけです。(Juliaでも単に 1:n と書くべきところを [1:n;]とかcollect(1:n)と書くと遅くなる。)
Pythonでも高速なライブラリ群を活かして気楽に(←これ決定的に重要)高速なコードを書く方法があるなら教えて欲しいところ。
すでに公開済みの
https://gist.github.com/genkuroki/cde66e33bd6a2b9a7d185f836c30a0fa
では、Juliaで指数積分函数を使うための手段として
(1) GNU MPFR ライブラリをccall経由で使う。
(2) Python SciPyの特殊函数ライブラリを使う。
の2つの手段を試してみて、速度比較もしています。
GNU MPFRの側は256bit浮動小数点なのでさすがに遅かったです。
素数定理の確認などではそこまでの精度は必要ないので、より短時間で計算できる Python SciPyライブラリを使った方が便利です。
学生の人なんかが、Julia用に高速で便利な指数積分函数のコードをを書けば、ちょっとした自慢になると思います。
現時点では100万以下の素数について素数定理が数値的に成立していることを確認するときの律速段階は対数積分函数を計算するところです。
気楽に色々計算して遊びたい場合には、気楽にコードを書けることが必須。
だからもしもJuliaからPythonライブラリを使うことが気楽にできないなら、Pythonの豊富なライブラリを使用可能でも、かなりがっかりだということになります。
現実にはほとんどPythonを使っていることを意識せずにPythonライブラリを使えます。atマークを{at}と書くことにすると、
using PyCall
{at}pyimport scipy.special as s
と書くだけで、SciPyの特殊函数ライブラリをJuliaでそのまま使えるようになります。
しかもMPFRで高精度計算する選択肢を使うよりも速い。
こういうのはとてもありがたいです。
もしも、SciPyの特殊函数ライブラリよりも速い指数積分函数をJuliaで書けるなら、誰か書くとよいと思います。
#julialang #python #MPFR #SciPy #jupyter #gist
Juliaで対数指数函数を使うための二つの選択肢の比較のノートブックを更新しました。
https://gist.github.com/genkuroki/cde66e33bd6a2b9a7d185f836c30a0fa
GNU MPFRのBigFloatによる計算はsetprecision()で精度を変更できるので、SciPy側と精度を合わせて速度比較をしてみました。
それでも私のコードだと、MPFR etan()を使うよりも、SciPy側のscipy.special.expi()を使った方が倍くらい速いです。
私の現時点での実力では、PythonのSciPyライブラリを使うことは速度面においても十分に魅力的な選択肢です。
@genkuroki MPFR側は(精度を制限しても)多倍長計算で,一方SciPy側のEi(x)はFortranのdoubleなので,速度の比較をするのは酷な気がします.
https://github.com/scipy/scipy/blob/master/scipy/special/specfun/specfun.f#L5587
#julialang #MPFR #python #scipy
GNU Scientific Libraryのマニュアルの指数積分函数の部分へのリンク
https://www.gnu.org/software/gsl/doc/html/specfunc.html?#exponential-integrals
積分の書き方が $dt$ を左側に書くスタイルになっている:
\[
\mathrm{Ei}(x)
= -PV\left(\int_{-x}^\infty dt\,\exp(-t)/t\right)
\]
SciPy側の対応する函数のマニュアル
https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.expi.html#scipy.special.expi