計算してみると, この数の立方根は
(expt 6064321219 1/3)=>1823.5908275197114
であり, この小数部分が本当に8月4日であろうかという気分になる. そこでユリウス日を計算する関数jdを使い,
(define (abel m d)
(exact->inexact
(/ (- (jd 1823 m d) (jd 1823 1 1))
(- (jd 1824 1 1) (jd 1823 1 1)))))
(abel 8 4)=> .589041095890411
(abel 8 5)=> .5917808219178082
となり, 8月4日を確認した.
あるいは
(* .5908275197114 365) => 215.65204469466102
1月から3月は90日, 4月から6月は91日, 7月は31日だから, そこまでで212日. つまり212.なにがしなら8月1日であり, 215.なにがしは8月4日である.
実は日付はどうでもよく, 問題はその手紙にAbelがLegendreの素数定理の式を書いていることである.
π(x)=x/(log x - 1.08366) (a)
昔読んだ素数定理は, Gaussが考えたもので,
π(x)=x/log x (b)
という気持のいいものであるが, それはあまりよく合わないといわれている.
早速絵を描いてみる.
x=2から1000000まで, 2つの式で左辺/右辺をプロットした. 青が(b), 赤が(a)である.999999点を描いているのだが, 最近の計算機ではあっという間である.
Legendreのは流石によくあっている. 高尚な理論あるようだが, 最小自乗法で合わせたかと思うほどである.
π(x)は, Eratosthenesの篩を使い, 1000000までの素数の奇数のビットマップを作った. MITのSchemeにはbitstringという便利なものがある.
(define x 500000) (define ps '(2)) ;xは上限 psは素数の表
(define sieve (make-bit-string x #t)) ;篩を作る
(bit-string-clear! sieve 0) ;素数でない1の場所をクリアする
(do ((i 0 (+ i 1)) (c 1)) ((>= i x))
(if (bit-string-ref sieve i) ;iの場所が1なら2i+1が素数
(let ((p (+ i i 1))) ;pの値
(set! c (+ c 1)) (set! ps (cons p ps)) ;c,psを更新
(do ((q (/ (- (* p p) 1) 2) (+ q p))) ((>= q x))
(bit-string-clear! sieve q)))));(p*p-1)/2からp置きにクリア
(define port (open-output-file "bittable"))
(do ((i 0 (+ i 128))) ((>= i x));128bitごとに符号なし整数に
(display (substring (number->string (+
(bit-string->unsigned-integer (bit-substring sieve i
(min x (+ i 128)))) (expt 2 128)) 16) 1 33) port)
(newline port)) ;十六進で出力
(close-output-port port)
これで3907行の出力が得られる. 最初の8行は
2196820D864A4C32816D129A64B4CB6E
4A2882D129861144A48961205A0434C9
148A48844225064B0834992132424030
65048928125108A00B40B4086C304205
C02104C94112422180124496804C3098
220825B0826896810804490000982D32
69009244304340069004265940A28948
086122D22400C06012410DA408088210
先頭の語の最後の 6E は 01101110で, 右から3,5,7,11,13が素数. 1,9,15は合成数を示す.
0 件のコメント:
コメントを投稿