pi関数

TAOCPにpi関数というのが登場する. 真理値表が 1100 1001 0000 1111 となる4変数の関数である. これは, 種明かしがしてあるが, 円周率πを二進で表わしたものだ. ランダムな真理値表が欲しいとき, πや素数表を使う.

十進のπの値は良く知られている.

TAOCPの第1巻の終りの方に八進のπの値が書いてあり,

3.11037 55242 10264 30215 14230 63050 56006 70163 2112

これを二進にすれば

11.001 001 000 011 111 101 101 010 100 010 ...

だから, 始めの方はたしかに上の真理値表だ.

ところが512ビットの真理値表というのも現れる.

1100100100011...00101

だというのである. 本当かな?と思って計算してみた. ついでだから, 1024ビットまで求めて, 十六進で書いたのが下のものだ.


2行目の最後が512桁で, 45だから00101だった.

この十六進の文字の分布も計算した. 下のリストの左から0, 1, 2,... , fの出現数で, 合計は256になっている.

(apply + '(15 16 17 16 21 17 21 10 14 13 14 19 18 14 17 14))
=> 256

ところでこの計算法だが, Lispのbignumを活用する. 十進法のπの値はspigotアルゴリズムで計算した. もっともここの段階で十六進に出来るはずだが, 面倒だったので, 十進の値を出した. 1024ビットは十進では310桁くらいなので, 320桁を用意した.


小数点は下から319桁にあるから, 10の319乗を用意する.

(define pow10 (expt 10 319))

とする. また最初を1100にしたいから, piを4倍する.

(set! pi (* pi 4))

十六進の印字プログラムも必要

(define (hexprint n)
(define hex "0123456789abcdef")
(display (string-ref hex n)))

そして計算する.

(hexprint (quotient pi pow10))
(set! pi (modulo pi pow10))

ここで最初のcが出力される.

(do ((i 0 (+ i 1))) ((= i 256))
(hexprint (quotient (* pi 16) pow10))
(set! pi (modulo (* pi 16) pow10)))

ここから90fdaa...が得られた.

折角プログラムがあるのだから, 自然対数の底のeでも同じ計算をしよう.


TAOCPの八進の値は

2.55760 52130 50535 51246 52773 42542 00471 72363 6166

で, 十六進1024ビットは以下のとおりである.