素数定理

高木貞治先生の近世数学史談に, Abelの手紙の日付が³√6064321219　とあり, それが1823年8月4日と解されている.

計算してみると, この数の立方根は
(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は合成数を示す.