2012年10月14日日曜日

多面体描画道楽

このタイトルのブログを最後にアップロードしたのは今年の9月9日で, 大二十面体の6枚の面に色づけした.

その頃から試みたかったのは, ある方向からの光線で陰影をつけることであった.

その前に見える面のすべてに色づけした図を描いてみるとこのようになった. 色の選び方には自信がないが.



これが出来たので, いよいよ陰影の図に取りかかる. 光の来る方向と各面の法線とのスカラー積を計算し, +1は光線に向いているから白に, -1は黒にすればよいのではと考えて次のような図が出来た.



しかし灰色が似ていてやはり面の識別は難しいというのが感想である. 枠となる正二十面体の頂点0,2,3のところだけ取り出したのがこの図である. 緑の正五角形の中が凹み, そこに5/2角錐が立っているのが見てとれればよいが, これもなかなか難しい.



何回か前のブログに書いたように, やはり難解な星形多面体である.

2012年10月8日月曜日

EDSACのプログラム技法

私が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."とコメントしている.