素数定理

2009年5月7日の素数定理のブログに, 青赤の線で素数分布の絵を描いた. それを描いた自分自身でも, その図で意外と思うのは, xが5万くらいのところで, どちらの線でも振幅が急に大きくなることだ. 素数の分布が変わるか? と一瞬思ったが, (まさか!) 青と赤の線で振幅の変わるxの位置が違うというのも怪しげなことである.

その原因を調べようと, 違う範囲の図を描いてみたら, またまた妙なことに, この振幅の変わる位置が, そのたびに違っていた. これはどうもPostscript(のシミュレータgs)に原因があるのではないか, というので理由を考えている. (まだ分かっていない.)

とりあえず, 描画の線幅を1から, 0.2に変えてみたら, お! 実にちゃんとした図が描けたのである. 0.5だと振幅の変化が相変わらず見られたが, うまくいった0.4で描いた図を下に示す. 縦の倍率は前回の5倍にしてある.

振幅の変化もなく, 赤い線が美しく1に接近しているのは感激だ.

さて, 傾向がよく見えるのはGaussの青線の方で, xが素数だとπ(x)が増えるので, 線は右上がりになり, 合成数だと増えないから, 右下がりになるのがよく分かる. つまり2から始まって, 3へは上昇, 4へは下降, 5で上昇, 6で下降, 7で上昇の後, 8,9,10で下降するのが見える. その後多少顕著なのは, x=100の先で, 青も赤も急上昇する. これは前回のQ₀の左の方の816d, 二進では1000 0001 0110 1101.
このビットは右から97, 99, 101, 103, 105, 107, 109, 111, 113, 115, 117, 119, 121, 123, 125, 127に対応し, 赤字が素数で, 97から113で上昇が見られるのである.

素数分布で話題になるのは, 大きい方での双子の素数と, ギャップである. 偶数が邪魔するので, 2つの続く整数が素数のことはないが, 1つ置きになら存在でき, それがどの辺まで分かっているかとか, どの辺にずーっと素数のない区間があるかとか, 関心が持たれている. TAOCPのMMIXでの素数表の解答のところには, このようなプログラムを使い, 418032645936712127とその次の素数の間に, 1370のギャップを見つけたと書いてある.

素数定理

The Art of Computer Programming, Volume 4, Pre-Fascicle 1Aに奇数に対する素数表がある. (5ページ)



早くいえば, 前回の話題の最後の方の十六進のビットマップを二進で表示したものである.

Q₀のビット列の右からの位置, i=0,1,...の各ビットは整数p=2i+1に対応し, pが素数なら1, 合成数なら0になっている. Q₀の左端はi=63, p=127, これは2⁷-1なので, 有名なMersenne素数であり, ビットも1になっている.

各Qは64ビット, 8語全体で512ビットなので, 最後(Q₇の左端)はi=511, p=1023に相当する. 1023は, 2¹⁰-1だが, 冪の10が合成数なので, 1023ももちろん合成数である. 1023=3*11*31

TAOCPは最近はMMIXというRISCマシンで記述されている. MMIXは1語最大64ビットである. Eratosthenesの篩を使った, 上の表の作り方を, 今日はC言語で書いて説明したい. C言語は, 私にとっては第一外国語みたいなもので, 大体はOKだが, 多少怪しいという言語だ. これはTAOCPの演習問題7.1.3-24であり, 同書にMMIXALによるプログラムも載っている.

#include <stdio.h>
main (){
  unsigned long long int qs [8];
  unsigned long long int
    x33=0x9249249249249249llu|0x4008010020040080llu;
  unsigned long long int
    x35=0x8421084210842108llu|0x0408102040810204llu;
  unsigned long long q;
  int j,pp=1,p=13,m;
  qs[0]=0x816d129a64b4cb6eull;
  for(j=1;j<8;j++){
    qs[j]=(x33|x35)^0xffffffffffffffffull;
    x33=x33<<2|x33>>31;x35=x35<<6|x35>>29;}
  for(j=0;j<8;j++)printf("%16llx\n",qs[j]);printf("\n");
  q=qs[0]>>6;
  do{
    for(m=p*p/2;m<8*64;m=m+p){qs[m/64]=qs[m/64]&
      (0x00000000000000001llu<<(m%64)
      ^0xffffffffffffffffull);}
    q=q>>1;p=p+2;
    while(q==0){q=qs[pp];pp=pp+1;p=(p|0x7f)+2;}
    while((q&1)==0){;q=q>>1;p=p+2;}
  }while(p<1024);
  for(j=0;j<8;j++)printf("%16llx\n",qs[j]);}
}

14行目と24行目で, qsを十六進で出力している.

816d129a64b4cb6e
2996c20d865a4c32
b49961225a2534c9
4a6982d12d86134c
 c36996132ca45b0
948a6894d325264b
4b48b6186d30d225
65a4c92916d128a6

816d129a64b4cb6e
2196820d864a4c32
a48961205a0434c9
4a2882d129861144
 834992132424030
148a48844225064b
 b40b4086c304205
65048928125108a0

この後半の十六進表記を二進にしたのが, 最初の表である.

qs[0]は127までの素数表で, これは手作業で作る. 100までの素数は25個しかなく, 127は13*13より小さいから, なんでもない. qs[1]からqs[7]までは, 篩だから, 通常はすべて1にするが, このプログラムでは, 3か5か7か11かで割れる数(に対応するビット)は, 最初から0にしておき, 13から篩い始める.

初期パターンを作るため, x33とx35を用意する. x33のビット毎ORの左は, 129から始めて3で割れる位置にビットをたててある. 129は3で割れるから, 右端は1, そのあと左へ向って001001と続ければ良い. 同じ語の右は11で割れるもののパターンで, 121が割りきれるから, 次に割りきれる奇数は143になる. これは129を0ビット目とすると, 7ビット目なので, そこから左へ1ビットめ毎に1をたてる. x35は同様なことを5と7についてやっている.

nで割れるパターンは長さnで繰り返す. 従って3と11で割れるパターンは33で繰り返し, 5と7で割れるパターンは35で繰り返す. MMIXは1語64ビットなので, 3と11, 5と7の繰り返しのパターンは1語に収まり, つまり合成できる. このことは私がKnuthに提案したので, TAOCPの解答には, We can combine 3 with 11 and 5 with 7, as suggested by E. Wada. と書いてある.

これらのパターンで1が立っているところは, 初期パターンでは0にするので, ビットの否定など, 多少注意が必要だ. またx33もx35も, 次の初期パターンに対しては, ローテートする必要がある.

とにかくこうしてqs[1]からqs[7]が設定出来た. 実行結果の上半分は, この設定直後のqsを示す.

次は篩だ. qにqs[0]を入れ, 11まで済んだので, 6ビット右シフとし, pも13に設定する. そしてqs[1]から篩始める. 最初の場所はp*p/2から. pの場所は素数を示すから, 消してはいけない. 次のpの倍数で奇数の場所はpだけ先である. しかし, 例えば5で篩う場合, 15はすでに3で篩われているから, 25から篩うので充分である. 25はi=25/2のビットに対応する. (小数部は棄てる.)

ある素数で篩ったら, qを右シフトし, pを2増やしながら, 次のビット位置を調べる. そうこうしている内に, qが0になったらどうするか. qを64ビットシフトしたか確認したくなるのが人情である. このプログラムではそういうことはせず, qsの次の語を取ってくる. 前のqsの語が終ったとき, ビット位置は64の整数倍引く1で, pは128の整数倍引く1の筈である. 従ってpは0x7fをビット毎ORして次のqsに対応する値にしている. この技もちょっとしたもである.

MMIXALのプログラムはもう少し凝っているが, 趣旨は殆んど同じである.

素数定理

高木貞治先生の近世数学史談に, 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は合成数を示す.