2009年5月8日金曜日

素数定理

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



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

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

各Qは64ビット, 8語全体で512ビットなので, 最後(Q7の左端)はi=511, p=1023に相当する. 1023は, 210-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のプログラムはもう少し凝っているが, 趣旨は殆んど同じである.

0 件のコメント: