#julialang #python #scipy #numpy #exponentialintegral #gist
指数積分函数のうちの一つである Ei 函数をJuliaで書いて試してみました。Numerical Recipes in C のコードのほぼ忠実なパクりで、Julia特有のスピードアップのためのチューニングは一切していません。
Python上のscipy.special.expi()とほぼ同じ速さで計算してくれました。何も工夫していないのに！
scipyの特殊函数自体は速いのですが、Juliaからの呼び出しに律速段階があるらしく、scipyの特殊函数をJuliaから速度重視の問題で使うのは止めた方が良さそうです。この点についておかしなことを言っていたので訂正しておきます。ごめんなさい。
GNU MPFRと256bit精度で比較してみましたが、これはMPFRのeint() の方が明確に速かったです。
詳しくは次のリンク先を参照
https://gist.github.com/genkuroki/99854018c45f3adde4b3cf2945b97800
$$
\mathrm{Ei}(x)
=PV\int_{-\infty}^x \frac{e^t}{t}\,dt
$$指数積分函数を数値計算するプログラムを書くことは、解析学の演習としても極めて優れものだと思いました。
こういう函数が Ei(log(x))=li(x) を通して、素数定理なんかの数値的確認のために役に立つという話もあって、色々教育的な感じ。
すでに触れましたが、Dirichlet積分$$
\int_0^\infty\frac{\sin x}{x}\,dx
$$の計算の話も関係してきます。
手元にある本では、
一松信著
『特殊関数入門』
森北出版1999
が色々なことがコンパクトにまとまっていて良い本でした。この本は買っても損がないと思いました。(電話帳のような本はすすめにくい。)
https://www.amazon.co.jp/dp/4627038291
https://www.amazon.co.jp/dp/4627038216
#julialang #numeticalanalysis #jupyter #gist
https://gist.github.com/genkuroki/99854018c45f3adde4b3cf2945b97800
Ei函数の別のアルゴリズムを採用したら1割ほど速くなりました。
Python SciPy との比較：
SciPy (FORTRAN) ～ 290μsec
Julia ～ 220μsec
大した違いではないとも言える。
数値計算用の函数のコードも Julia で書いて十分に実用になるんですね。
こういうのはとてもありがたいと思います。
採用したのは
http://www.mymathlib.com/c_source/functions/exponential_integrals/exponential_integral_Ei.c
と同じ計算法。
先の Numerical Recipes in C との違いは $x\ne 0$ が $0$ に近いところで数値計算を $x$ に近い整数点でのTaylor展開で行っていることです。
高速計算の専門家によるexp integralの計算では、
"Pretty good! We are 6 times faster than the SciPy routine for complex arguments, and 5 times faster for real arguments, even though SciPy internally calls an optimized Fortran routine."
なのだそうです。
Juliaで書かれたコードは、SciPyで採用されているFortranのルーチンよりも、実函数で5倍、複素函数で6倍速く計算できると言っています。
Fortranの5～6倍！
マジかよ！！！
https://github.com/stevengj/18S096-iap17/blob/master/pset3/pset3-solutions.ipynb https://mathtod.online/media/zMz7NZ4FcbowX1mhkfg
https://github.com/stevengj/18S096-iap17/blob/master/pset3/pset3-solutions.ipynb
からの引用
"It's actually not unusual for special functions written in Julia to beat older optimized C and Fortran implementations. ～The reason is that the inlining performed by macros like ＠evalpoly and ＠e₁_cf64 is hard to replicate in lower-level languages — hardly anyone bothers to write the code-generation programs that would be required — but it is fairly easy in Julia."
おそらく、Julia側でやっていることのすべてをCやFortranでもやればCやFortranの側でも相当に速くなるのだと思います。
しかし、プログラムは人間が書くものなので、面倒な作業が追加されるとやる気が失せるし、誤りも増えるし、他の仕事もできなくなる。
しかし、現代的な高級言語であれば、プログラムを楽に書くための仕組みが最初から言語に組み込まれているので、高速化のための技法を楽に実行できると。
https://github.com/stevengj/18S096-iap17/blob/master/pset3/pset3-solutions.ipynb
をダウンロードし、Juliaのバージョンアップのせいで必要になった一部手直しをした後に、手元の環境で実行してみました。
確かにJuliaの方がSciPy(Fortran)より5～6倍速いです。
添付画像は私が走らせた結果です。
なるほど、アルゴリズムの工夫による数値計算の高速化は偉大な知恵だな。
原理的に「こうすれば速くなる」とわかっていてもそういうコードを楽に書ける仕組みがないと、人間はそういうコードをほとんど書かなくなるものだと思う。
正直、マクロを使いまくりのコードにまだ着いて行けていないのですが、かなり参考になりそうです。
この水準の数値解析の講義であれば、受講している人たちは相当に楽しいと思う。
https://github.com/stevengj/18S096-iap17/blob/master/pset3/pset3-solutions.ipynb
#julialang #numericalanalysis
Juliaでは積分指数函数を使う標準的な方法がまだないのですが、GNU MPFR は高精度過ぎるし、Juliaから scipy.special.expi() を呼びだすのも遅いです。
だから、Juliaで書いた指数積分函数を私はこれから使うことになりそうです。
現時点で私が知っている選択肢の中ではこれが一番速い。
li(x)=Ei(log(x)) の形で使う。
大量の数を $x$ に入れて計算したりする可能性があるので、計算速度はかなり大事。