正木です。


|[ruby-math:00753] [Q]Gamma & Zeta functions
|From: Masaaki Sakano <mas / star.le.ac.uk>

|今、Ruby で、Γ関数と Riemann の ζ関数とを使いたいと思っています。
|今の目的では、精度も速度も不要です :-)
|# 10秒で3桁が出ればOKです。

|正木さんの過去投稿の中で、Real クラスとして、両関数が含まれたものを
|公開されていたのを見つけました。、が、それら過去投稿の修正を
|いろいろ加えた Real クラスを試してみたところ、残念ながら
|同クラスの定義ファイルを require している途中でこけます。

|正木さんの手元には最終(最新)版 Real class があると思うのですが、
|それをどこかで公開してらっしゃたりはしないでしょうか?
|もしくは ruby-math かあるいは直接メイルで頂けたら嬉しいのですが…。

Real Class の修正が重なったのと、自分で使っていて、ある程度実用的になった
と思うので、ここで公開して、皆さんの御批判及び advice を受けたいと思ってい
るんですが、興味を持つ人がいなさそうなのと、あまり長い script をここに流す
のも迷惑かも知れないと思って躊躇しています。
坂野さんには直接 mail でお送りします。
もし他にも興味のある人がいれば mail を下さい。

Real Class の基本的な idea は次の通りです。
無限列を表す class Sequence を作る。
Rational(Integer を含む),Complex(Rational,Rational) の Cauchy 列 s と、
その近似値を表出するための補助的な 0 に収束する非負数列(誤差項) err を
全ての n に対して |(limit s) - s[n]| <= err[n]
を満たすように選んで Real.new(s,err) と定義する。

四則演算と関数を上の条件が自動的に満たされるように定義する。

今の所使える関数は
Sqrt,Exp,Log,Cos,Sin,Tan,Atan,Gamma,Zeta,Power
ですが algorithm さえ分かっていれば関数を定義するのは簡単です。

Numeric Class としては == の判定ができないと言う原理的な問題があります。
例えば

Exp[Pi*Complex::I]==-1
のような判定はいつまでも終らない。

1/(Exp[Pi*Complex::I]+1)
は error にならずに、memory や stack を使い果たすまで計算し続ける。

計算速度に関してはまだ充分改善の余地があると思っています。参考のために
test program の出力結果を以下に付記します。(Pentium  1.6GHz)

EulerConstant=0.5772156649015328606065120900824024310421
user	0.69 sec
ZetaA[1]=0.5772156649015328606065120900824024310421
(ZetaA[s]=Zeta[s]-1/s)
user	0.56 sec
Pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679
user	0.48 sec
1/Pi=0.3183098861837906715377675267450287240689192914809128974953346881177935952684530701802276055325061719
user	0.42 sec
Zeta[Complex(1/2,1/2)]=Complex(-0.45930,-0.96125)
user	2.11 sec
Zeta[2]=1.6449340668482264364724151666460251892189
user	0.05 sec
Gamma[1/2]=1.77245385
user	0.71 sec
Gamma[Complex(1/2,1/2)]=Complex(0.81816,-0.76331)
user	3.15 sec
Log[Complex(0,1)]=Complex(0,1.57079632679489661923)
user	0.57 sec
Sqrt[Complex(0,1)]=Complex(0.70710678118654752440,0.70710678118654752440)
user	0.02 sec
Exp[Pi*Complex::I]=Complex(-1.00000000,0.00000000)
user	3.72 sec
Zeta[Complex(0.5,14)]=Complex(0.0222,-0.1032)
user	0.02 sec