私が2000年ころ書いたEDSAC用Eratosthenesの篩の説明 少し長いが興味と元気があればフォローして欲しい.
まずプログラムの先頭はこうなっている.
これはすべて定数か作業場所である.
以下のが本体.
本体のプログラムは先頭の T 96 K G K で以下のプログラムを96番地から格納する. 最後のE 96 K P F で96番地から実行を開始する. 以下相対番地を前回のように ' で示す.
0'の T F でアキュムレータをクリア. S 850 Dで850,851の P D, P F つまり長語の1を引き35ビットオール1を作る. T 992 Dで 992,993番地へ入れる. 2'番地をアキュムレータに置き, A 849 Fでアドレス部に2を足して2'へ戻す. 848の T 1024 D を引き, 負なら0'へ. このループで長語の992から1022まで16長語をオール1にして篩を用意する.
8'からのループは A 1022 Dでオール1にし, 850 D の1を引き, 111...10(負のストロブ)を作って長語852に入れる. アキュムレータはクリアされていてそれに 850 D の1を足し, 000...01 (正のストロブ)を作って長語922に入れる. これを左1ビットシフトし, 000...10 にして 850 D へ戻す.
16'から22'までは11'と13'の格納命令を2ずつ増やし, 結局長語852-920には負のストロブ, 922-990には正のストロブを用意する.
23', 24', 25'でFIGS, CR, LFを出力. 26', 27', 28'で2を作り, 29'で0番地へ格納して30', 31'のWheelerリンケージで出力ルーチン P6 へサブルーチンジャンプして最初の素数2を印字する. 33'からがいよいよ篩だ.
以下の表で左端の852の列は長語の番地 次の0の列はストロブや篩の長語の番号 その右の二進表示が長語の内容である.
852 0 11...110 負ストロブ
854 1 11...101
856 2 11...011
...
920 34 01...111
922 0 00...001 正ストロブ
924 1 00...010
926 2 00...100
...
990 34 10...000
9753←対応する奇数の素数
992 0 11...111 篩
994 1 11...111
996 2 11...111
...
1022 15 11...111
0番の篩の長語の各ビットはその上に書いてあるように右から3,5,7,9,...に対応する.従って最後の15番の篩の左端のビットは1121に対応する. EDSACのメモリー容量からすれば篩はもっと大きく出来るが, 篩の様子を水銀タンクで眺められるように, ちょうど1つのタンクの大きさにした.
篩は正のストロブを順に使い, 篩の右から順のビットとandをとって, 素数か合成数かをみる. 合成数なら次のストロブで次の篩のビットを調べに進む. 素数ならその素数をp とすると, 左の図のようにpビットごとにpの倍数があるから, 負のストロブを使っていま調べたビットからpビットおきに篩のビットを0に変えていく. 手順としてはこれだけだが, ストロブが下まで来たときに, 篩を次の長語にする; ストロブを先頭に戻す作業が必要である.
33'で正のストロブを乗数レジスタに置き34'で篩とandをとる. 844Dの長語の1を引き負なら合成数だったので75'へ. そうでないなら37'から素数出力; 841に素数の4倍があるので42'で右2ビットシフトし, 44',45'でP6へいく.
46'からは素数の倍数の篩を0にする, つまり篩う仕事を開始.
以下抜き書きしたプログラムは左から命令のある相対番地, 命令, 番地の示す場所の内容, つまりオペランド, 演算結果のアキュムレータ, 格納命令の場合は格納番地と内容を示す.
命令 オペランド アキュムレータ
47 A 34 @ C 992 D C 992 D
48 U 71 @ 71 C 992 D
49 S 843 F L F T 992 D
50 T 72 @ 72 T 992 D
51 A 33 @ H 922 D H 922 D
52 S 842 F P 70 F H 852 D
53 A 841 F P 6 F H 858 D
54 U 70 @ 70 H 998 D
55 S 840 F H 992 D P 6 F
56 G 69 @
47'で今素数を調べた篩を取り出し71'へ入れ, ファンクション部をTにした命令を72'へ入れる. 51'は素数を調べた正ストロブを乗数レジスタに入れる命令で, 篩を0にするには対応する負ストロブが必要である. それは7番地若い方にあるからまず70を引き, 素数のビット数だけ先から篩うので841番地にある素数の4倍の値を足し, 消すストロブの
番地を作る. 55'は負ストロブの範囲を超えたかのチェック.
範囲を超えていなければ69'へ来る. 以下のプログラムの70',71',72'は先ほどのプログラムが設定したものだ.
69 T F
70 H 998 D
71 C 992 D
72 T 992 D
73 A 70 @ H 998 D H 998 D
74 G 53 @
53 A 841 F P 6 F H1004 D
54 U 70 @ 70 H1004 D
73'からは篩をクリアする命令を更新する. H 998 Dをもって53'へ. そこで再び素数の数だけストロボのアドレスを増やす. こうしている内に負ストロブの範囲を飛び出し57'へ進む. 58'で篩の番地を取り出し, アドレスを2増やし71'へ戻す. 72'のT命令も増やす. しかし篩の方が範囲を超える心配があるので, 63'でT1024 Dを引き, 篩の範囲が終わっていなければストロブの番地をもとへ戻す. それが66'からでストロブを取り出し, 番地から70を引いて1周戻し, 54'へいって次のサイクルに入る
篩の最後まで来たら, 合成数の時と同じく75'へ行く.
75'からは次の奇数の素数テストになる. 76'からは素数の4倍をもっていたカウンタ841を4増やす. 79'からは正ストロプを乗数レジスタに置く33'の命令を2増やし, 82'で正ストロブの範囲をチェックし, 範囲内なら32'へ行って左隣のビットをテストする. 正ストロブが範囲を超えているなら, 84'で正ストロブを先頭に戻し, 86'から篩の長語を次にすすめ, 89'で篩も範囲を超えたなら91'でプログラムを終止する.
初期の機械語命令の時代はこういうプログラムは普通であった.
このプログラムが出来た時, EDSACのシミュレータを公開している英国Warick大学のMartin Campbell-Kelly君に送ったらシミュレータのホームページにWadaSieveとして公開してくれた.
彼は"The program is beautifully written and very fast. A little master work."とコメントしている.