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 #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