2012年11月5日月曜日

Piの1000桁

IEEE Annals of the History of Computing July-September 2012に, 1949年のLabor Dayの休暇にENIACで円周率を2035計算したという記事があった.

円周率(Ludolph's number)といえばWilliam Shanksが1873年に707桁計算し, 自分の墓に刻んたという話が有名であったが, その値が528桁目から下は違っていて, それを見付けたのはENIACの計算によるのかと思っていたが, ウェブで調べると見つけたのは同時代だが, 電動計算機だったらしい.

今から半世紀前, 日本でもあちこちの大学や研究所で計算機が作られはじめ, デモ用に自然対数の底e(Napier's number)を1000桁計算するプログラムが作られた. 私たちもパラメトロン計算機PC-1でプログラムを作ったのは 2012年5月9日のこのブログに書いた通りだ.

ENIACの記事によれば, 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·2395)-1/(7·2397)+...

見て分るようにこの式はeの式

e=1+1/1!+1/2!+1/3!+...

より面倒なので, πに手を出す人はあまりいなかった. しかしtan-1は交代級数なので, 計算は楽な筈である. 1000桁の計算ならtan-1(1/5)は1/(2n+1·52n+1)が1/101000より小さくなる辺りまで計算すればよい.

手始めに100桁でやってみる.
(define c (expt 10 100))
(do ((n 1 (+ n 2))) ((< (/ c n (expt 5 n)) 1) n)) => 141
(do ((n 1 (+ n 2))) ((< (/ c n (expt 239 n)) 1) n)) => 43
これで必要な項数は判ったので, 105桁のbignumで計算する. aは tan-1(1/5), bはtan-1(1/239)で, 最後にaの16倍からbの4倍を引き, 最後の5桁を捨てる. opはループで毎回0と1に切り替わりそれに従って次の項を足したり引いたりする.
(define c (expt 10 105))
(define op 0) (define a 0)
 (do ((n 1 (+ n 2))) ((= n 143))
  (set! a ((if (= op 0) + -) a
   (quotient (quotient c n) (expt 5 n))))
  (set! op (- 1 op)))
(define op 0) (define b 0)
(do ((n 1 (+ n 2))) ((= n 45))
 (set! b ((if (= op 0) + -) b
  (quotient (quotient c n) (expt 239 n))))
 (set! op (- 1 op)))
(quotient (- (* a 16) (* b 4)) 100000)
で結果は
314159265358979323846264338327950288419716939937510
 58209749445923078164062862089986280348253421170679
πの100桁, 1000桁の値はすぐに見付かる. 私は城, 牧之内「計算機械」にあったのを知っているので書棚から取り出し同じであることを確認した.

1000桁も同様に計算したが, 書くまでもあるまい.