Gauß の公式


最初に断っておくが, この program は土居範久, 慶應義塾大学助教授 (当時)「PASCAL 入門」(電子計算機のプログラミング 9) 培風館, 1985, 第一部 第十二節「長桁計算」の例題 34 をそのまま porting したものである。従ってあんまり解説しすぎると著作権に触れる恐れがあるので, 適当に save して解説する。勿論, 色々余計な変数を使っているのを止めたりなどの改変はしてある。

私の円周率を解説した page に逆正接函数系と称して C.F. Gauß の公式が出ている。もう一度書くと
π = 48 Arctan(1/18) + 32 Arctan(1/57) - 20 Arctan(1/239)
である。これは逆正接函数系の公式の中では比較的効率の良い方なのでこれを採用する。問題は Arctan(x) をどう計算するかであるが, 良く知られているように次の無限級数に展開することができるので, これを (誤差が出なくなるような) 適当な桁数まで計算することによって (時間と memory の許す限り) 何処までも計算することが出来る。
Arctan x = Σn = 1[(-1)n-1x2n-1 / (2n -1)], 但し |x| ≦ 1.

例によって配列一つあたりは十進で 9 桁の数を格納することにする。大域的な配列を三つ用意し, pi[] に, 円周率の値を入れていく。助変数 work1[], work2[] には各々 x2n-1 と x2n-1 / (2n - 1) を計算しておく。つまり先ず work1 の方で x2n-1 を計算しておき, それを 2n - 1 で割った値を work2 に入れ, 順に pi[] に加えていけば良いのである。明らかに work1[], work2[] の最初の方は次第に 0 ばかりになってくるので, 計算を速くするために, 何処から 0 でなくなるのかというのを示す変数を受け渡ししている。(この変数を所謂「参照渡し」で渡してやらなければならないのだが, JAVA でそれを実現する為に要素 1 つの配列で渡してしまうという少々 tricky な方法を取らねばならなかった。ということは本当はきちんと class の設計をしなければならないということでもある。)

Sunday, 7th May, 2000.


BACK