2021年6月7日月曜日

Piの1000桁

前回のブログのPiの1000桁をもう少し詳しく説明したい. Piの値の計算に Machinの式

π/4 = 4 tan-1(1/5)- tan-1(1/239)
tan-1(1/5)=1/5-1/(3·53)+1/(5·55) -1/(7·57)+...
tan-1(1/239)=1/239-1/(3·2393)+1/(5·239 5)-1/(7·2397)+...

を使った. これは絵を描いてみるとこういうことだ.
円弧の中に底辺:高さが5:1の直角三角形を描く. その中心角をαとする. つまり α=tan-1(1/5). αは11.3°くらいだ.

αの斜辺の上にまた今の直角三角形を描く. それを3回繰り返すと4αの 角度が得られ, 上の値を使えば4α=45.2°だ. 一方π/4は45°だから, 今度は底辺:高さ=239:1の直角三角形(赤線)を角度を引く方法に描くとその斜辺が 丁度45°になるというのがMachinの式であった. 計算すると次のようになる.
一番上はtanの加法公式である. これを使いαから2α, 4α を 計算すると4α=120/119が得られる. この120と119の差が1なのが 重要だ. 4αとπとの差の角度βのtanはまた加法公式を使うと1/239 になるわけだ. 一番下のようにMachinの式になる.

Machinが分ったから計算に進む. 最初のPiの1000桁のブログでやったように, 40桁計算する様子を見るには, 小数点を40桁右に移動して長い整数(bignum)で計算する. つまり1/5とせず, ss=1040/5 から始める.

それを1で割ってtsに足し, ss=ss/(5×5)にしてから ss/3をtsから引き, ss=ss/(5×5)にしてから ss/5をtsに足し, ... を計算する. しかし最初だけ 5で割り, 後は5×5で割るなら, 最初の分子を5倍しておき, いつも5×5 で割るのがプログラムは簡単だ.

さらにMachinの式はπ/4を計算するが, 我々はπが欲しいから, 最初のss は 4×4×5×1040とする.

従って
(set! ts 0)
(set! ss (* 4 4 5 (expt 10 40)))
続けて
(set! ss (quotient ss 25))
(set! ts (+ ts (quotient ss 1)))
(set! ss (quotient ss 25))
(set! ts (- ts (quotient ss 3)))
...
これをssが0になるまで実行する. 初期化の後 doループを使い
(set! ts 0)
(set! ss (* 4 4 5 (expt 10 40)))
(do ((n 1 (+ n 2)) (ff #t (not ff))) ((= ss 0) n)
  (set! ss (quotient ss 25))
  (set! ts ((if ff + -) ts (quotient ss n))))
doループは, nを1, ffを#tで始め, ループを回る度にnはn+2にし, ffは ¬ffにしてssが0になるまで, 本体を実行すると読む. つまりnは1,3,5,7,... ffは#t,#f,#t,#f,...となる. 終了条件s=0の後のnはdoループの返す値である. nをどこまで進めるとsが0になったかを知るためだ.

((if ff + -) ts (quotient ss n))は, ffが#tなら (+ ts (quotient...)), #fなら(- ts (quotient ...))のことである. つまり1回置きに加算と減算を繰り返す.

初期化したssとループでのssとtsを出力すると
ss 800000000000000000000000000000000000000000(
(n 1)
ss 32000000000000000000000000000000000000000
ts 32000000000000000000000000000000000000000
(n 3)
ss  1280000000000000000000000000000000000000
ts 31573333333333333333333333333333333333334
(n 5)
ss    51200000000000000000000000000000000000
ts 31583573333333333333333333333333333333334
(n 7)
ss     2048000000000000000000000000000000000
ts 31583280761904761904761904761904761904763

...

(n 57)
ss                                        23
ts 31583289575980921339207962431166446951613
(n 59)
ss                                         0
ts 31583289575980921339207962431166446951613
61
最後の61はdoループが返したnの値である. 下の方の23や0は元々は右端に 揃っていたのだが, ブログはそういう風には出してくれない.

同じようにtan-11/239も計算すると
(n 17)
ss                                         1
ts 31415926535897932384626433832795028841973
(n 19)
ss                                         0
ts 31415926535897932384626433832795028841973
21
と円周率が見えてくる.

余談だが, 私の恩師 故高橋秀俊先生は中学生の時円周率を筆算で50桁計算されたそうだ. 40桁でもn=57まで計算するのだから, 50桁はまだまだ先までで大変だっただろう.

次はこの計算を有限語長ので分割しながら実行するわけだが, それは次のブログにしよう.