2012年12月14日金曜日

生後n日目

最近 還暦, 古希, 喜寿, 傘寿などの祝事で兄弟があつまる事が少くない. そういう折に私が「生れて何万何千日たった」と発言したので, まわりが似たような計算をしたがった.

私自身はJulian Date(DJ)が計算できる例の個人用電卓を携帯しているので, 「今日の日付を入力, JDに変換, 生れた日付を入力, JDに変換, 引き算」で簡単に得られる. もっともこれでは生れた日が第0日になり, 生れた日を第1日にする常識的な勘定と違うが, 満n日といえば良いかもしれない.

しかし一方, 「ちょうどn日になるのはいつかなぁ」にはJDだけでは片付かない. そこで誕生日とnを入力すると, 生後n日目になる日付を計算するウェブアプリを作ってみようと考えた.

アルゴリズムは私のブログ, 2011年5月23日のunix timeにあるのを使えばよい. 現在生きている人が対象なら, Julian Calendarは考えなくてよい. Gregorian Calendarが昔まで続いていたとするFixed Date(Rata Die)というのがあるから, 年月日からFixed Dateへの変換とその逆変換があれば, 出来たようなものである.

私はhtmlのウェブページは何回も書いたことがあるが, 計算するにはJavaScriptでも使うことになるであろう. O'reillyのJavaScriptの本は持っている. その例をみながらやって見ることにした.

<html>
<head>
<title>nth day after birth</title>
<script language="JavaScript">

この辺まではお作法のうちだ. JavaScriptでは整数の除算が見当たらなかったので, 整数除算の関数を作る.
function quotient(a,b){
var r=a%b;
var q=(a-r)/b;
return q;}

つぎは私のカレンダーのプログラムによく登場する前月までの日数の和のリストで, 平年用と閏年用を用意する.
var mon0=[0,31,59,90,120,151,181,212,243,273,304,334];
var mon1=[0,31,60,91,121,152,182,213,244,274,305,335];

次は年月日からFixed Dateを計算する関数rdだ.
function rd(y,m,d){
var a=(y-1)*365;
var b=quotient(y+3,4);
var c=quotient(y,400)-quotient(y,100);
var f=(((y%100)==0)?((y%400)==0):((y%4)==0))?mon1[m-1]:
  mon0[m-1];
return a+b+c+f+d-1;}

ご覧のように簡単だ.

反対にFixed Dateから年月日を求める関数grだ. ブログにあったアルゴリズムはfloor関数を使うが, JavaScriptにはないらしいので, 1で整数除算をすることにした.
function gr(rd){
var d0=rd+366-1;
var y400=quotient(d0,146097);
var d1=d0%146097;
var y100=quotient(d1/36524.25,1);
var d2=d1-146097+quotient(36524.25*(4-y100),1);
var y4=(y100==0)?quotient(d2,1461):
  (24-quotient((36523-d2)/1460.96,1));
var d3=(y100==0)?(d2%1461):(d2-quotient(y4*1460.96,1));
var y1=((y100>0)&&(y4==0))?quotient(d3,365):
  quotient(d3/365.25,1);
var d4=((y100>0)&&(y4==0))?d3%365:
  (d3-1461+quotient(365.25*(4-y1),1));
var y=y400*400+y100*100+y4*4+y1;
var leap=((y%100)==0)?((y%400)==0):((y%4)==0);
var e=leap?mon1:mon0;
var mon=0;
for(i=0;e[i]<=d4;i++){mon=mon+1;}
var d=d4-e[mon-1]+1;
document.write(y.toString(10)+" "+mon.toString(10)+" "
  +d.toString(10));
document.write("<br>");}

これでテストしてみるとうまく行く. 結果は関数の中でdocument.writeする.

ところでy,m,d,nの入力はどうするか. 読み込んだものは文字列らしい. 解説書で使えそうなものを探すと, データ間はカンマで区切り, カンマの位置をindexOfで探し, substringを取れば, それぞれのデータが得られることが分った. またそれだけだと文字列であるが, 0を引くと整数になるという奥の手も知った. こうして書いたのが下の<body>部分である.
</script>
</head>
<body>
<script language="JavaScript">
document.write("<br>");
var s,i,y,m,d,n;
s = prompt("Type in your birthyear,birthmonth,birthday and n
  as `2000,1,1,1000' without space", "");
i = s.indexOf(',');
y = s.substring(0,i);
s = s.substring(i+1,s.length);
i = s.indexOf(',');
m = s.substring(0,i);
s = s.substring(i+1,s.length);
i = s.indexOf(',');
d = s.substring(0,i);
n = s.substring(i+1,s.length);
document.write("y = " + y);document.write("<br>");
document.write("m = " + m);document.write("<br>");
document.write("d = " + d);document.write("<br>");
document.write("n = " + n);document.write("<br>");
y=y-0;
m=m-0;
d=d-0;
n=n-0;
var r=rd(y,m,d);
document.write("Your "+ n.toString(10) +"th day is ");
gr(r+n);document.write("<br>");
</script>
</body>
</html>

さて上のプログラムを呼出すと下のような画面が現れる.





例えばAlan Turingは1912年6月23日生まれ. その滿1万5千日目をみるには

1912,6,23,15000

と入力すると
y = 1912
m = 6
d = 23
n = 15000
Your 15000th day is 1953 7 18

と得られる. Turingは1954年6月7日に他界したから, 1万5千日とすこししか生きていなかったわけである.






2012年11月5日月曜日

Piの1000桁

IEEE Annals of the History of Computing July-September 2012に, 1949年のLabor Dayの休暇にENIACで円周率を2035計算したという記事があった.

円周率(Ludolph's number)といえばWilliam Shanksが1873年に707桁計算し, 自分の墓に刻んたという話が有名であったが, その値が528桁目から下は違っていて, それを見付けたのはENIACの計算によるのかと思っていたが, ウェブで調べると見つけたのは同時代だが, 電動計算機だったらしい.

今から半世紀前, 日本でもあちこちの大学や研究所で計算機が作られはじめ, デモ用に自然対数の底e(Napier's number)を1000桁計算するプログラムが作られた. 私たちもパラメトロン計算機PC-1でプログラムを作ったのは 2012年5月9日のこのブログに書いた通りだ.

ENIACの記事によれば, Machinの式で計算しそうだ. つまり

π/4 = 4 tan-1(1/5)- tan-1(1/239)

tan-1(1/5)=1/5-1/(3·53)+1/(5·55)-1/(7·57)+...

tan-1(1/239)=1/239-1/(3·2393)+1/(5·2395)-1/(7·2397)+...

見て分るようにこの式はeの式

e=1+1/1!+1/2!+1/3!+...

より面倒なので, πに手を出す人はあまりいなかった. しかしtan-1は交代級数なので, 計算は楽な筈である. 1000桁の計算ならtan-1(1/5)は1/(2n+1·52n+1)が1/101000より小さくなる辺りまで計算すればよい.

手始めに100桁でやってみる.
(define c (expt 10 100))
(do ((n 1 (+ n 2))) ((< (/ c n (expt 5 n)) 1) n)) => 141
(do ((n 1 (+ n 2))) ((< (/ c n (expt 239 n)) 1) n)) => 43
これで必要な項数は判ったので, 105桁のbignumで計算する. aは tan-1(1/5), bはtan-1(1/239)で, 最後にaの16倍からbの4倍を引き, 最後の5桁を捨てる. opはループで毎回0と1に切り替わりそれに従って次の項を足したり引いたりする.
(define c (expt 10 105))
(define op 0) (define a 0)
 (do ((n 1 (+ n 2))) ((= n 143))
  (set! a ((if (= op 0) + -) a
   (quotient (quotient c n) (expt 5 n))))
  (set! op (- 1 op)))
(define op 0) (define b 0)
(do ((n 1 (+ n 2))) ((= n 45))
 (set! b ((if (= op 0) + -) b
  (quotient (quotient c n) (expt 239 n))))
 (set! op (- 1 op)))
(quotient (- (* a 16) (* b 4)) 100000)
で結果は
314159265358979323846264338327950288419716939937510
 58209749445923078164062862089986280348253421170679
πの100桁, 1000桁の値はすぐに見付かる. 私は城, 牧之内「計算機械」にあったのを知っているので書棚から取り出し同じであることを確認した.

1000桁も同様に計算したが, 書くまでもあるまい.

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

2012年9月28日金曜日

EDSACのプログラム技法

EDSACのプログラムの読み書きに重要なのが文字コードである. コードの知識なしでは仕事にならない. それ以前の計算機での入力法はよくわからないが, EDSACではプログラムを紙テープにパンチし, それをイニシアルオーダーで読み込むことが例の本で公開されたので, イニシアルオーダーを解読するためにも文字コードに関心をもつことになった.

EDSACは大学で作った計算機なので, 入出力は市販の機器を用いることになる. 当時はもちろんテレタイプは存在していたから, 当然それらを使うわけだ.

その頃 欧米で使われていたテレタイプは, 英大文字と数字と若干の記号のもので, 5単位テープが標準であった. このコードはInternational Telegraphic Alphabet No.2 (ITA2)という. (規格はここのFreely available itemsのEnglishから得られる) この種のテレタイプは昔は国際電電にいくと見られたがとおの昔に姿を消した.

テレタイプとしてはこのような形であった.



なかなか可愛らしくて1台ほしい. 数年前にアメリカのある空港で耳の不自由な人のための似たようなキーボードを見たことがあった.

コード表は次のとおり.



Figure shift(FIGS)を打つとその後は右の数字と記号を送受信することになる. Letter shift(LTRS)で英大文字になる.

このコードの特徴は, 使用頻度の高い文字は1が, つまりテープの孔が少ないことだ. Wikipediaでletter frequencyを見て, 頻度の順に孔の数をならべると
e t a i n o s r l d h c u m f p y g w v b k x j q z
1 1 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 4 3 4 4 3 4 2
孔が少ないとテープが丈夫なのか, ごみを減らそうとしたのか知らないが, Morseコードと同様な精神で出来ている. ITA2の特殊記号も欧文Morseコードと同じである.

さてEDSACではこういう機器を使おうとしたが, 数字のコードが
P 0 01101
Q 1 11101
W 2 11001
E 3 10000
R 4 01010
T 5 00001
Y 6 10101
U 7 11100
I 8 01100
O 9 00011
と数値に関係がない. そこで最上段の文字のコードを00000から01001に変更した. ただそうするとP 0のコードはまったく孔が開かず, ブランクテープと同じになってしまうので, テープでは一番左の位置はコードが0の時に孔が開くようにした.

変更は最小にしたらしいが, T→F→Y→T, O→D→W→O, R→X→U→S→R, I→A→I, P→B→E→P, Z→Q→Zと交換した. ITA2と同じなのはH, N, M, L, G, C, V, J, Kである.

これがEDSACの文字コードだ.



左の方がPerforator, テープ鑽孔機で, そのすぐ右がTeleprinter, 出力用タイプライターである. 図形文字でない機能鍵のFigure shiftもプログラムでは文字として使いたいというので, 鑽孔機にはπの文字がついているという具合である.

EDSACの最初の本のコードでは, 鑽孔機の記号の位置はタイプライターのそれとずいぶん違っていた. その辺は上の表では省略した.

EDSACのプログラムでは命令を A 3 F のように書くが, 数字の前にFigure shiftを打つかというとそうではない. 命令の最初は文字と思い, そのあとPQWERTYUIOJが続くあいだは数字として扱う. Jは10として使え, そう使うプログラムを存在する.

一方数値だけの疑命令でも, 先頭にはPを書かなければならない. 命令はコード表のFからVまでの文字で終わる. π, S, Z, Kには別の機能があった.

EDSACはこのようにコードを変更したが, 我々のパラメトロン計算機では, コードは市販の機器のままで, あとはプログラムで挑戦するという方法をとった. 日本国内では6単位が標準だったので文字数も多く, プログラムが見やすいシステムが作れた. その話はまたいつか.

2012年9月22日土曜日

EDSACのプログラム技法

英国ケンブリッジ大学のEDSAC(Electronic Delay Storage Automatic Calculator)が稼働し始めたのは1949年5月6日であった.

EDSACのメンバーはプログラムの解説書(Maurice V. Wilkes, David J. Wheeler, Stanley Gill: The Preparation of Programs for an Electronic Digital Computer)を早々に出版した. 世界で最初のプログラムの本であり, 多くの人がそれでプログラミングの楽しさを知った. 私もその1人である.

遠い昔のテクニックはどんどん忘れ去られる. 最近EDSACのプログラムを眺める機会があったので, その頃のプログラムを解説したくなった.

当時のメモリーは水銀遅延線かブラウン管であった. EDSACはそのDelay Storageから分かるように, 水銀遅延線の記憶装置を持つ. つまり水銀タンクの一方からPiezo効果を使い音波を送り, 反対側で受けた音波を逆Piezo効果で電 気信号にもどす. その遅れを記憶装置として利用する. (水銀遅延線の写真はここに)

EDSACのアーキテクチャは簡単である. レジスタとしては71ビットのアキュムレータと35ビットの乗算レジスタ. 記憶装置は17ビットの短語が1024語. 2n番地と2n+1番地の短語を2語つなぐと35ビットの長語になる. 不思議な計算だがその秘密は, 水銀遅延線での短語は18ビット分あり, 語と語の切替えに1ビットの時間がかかるので, 長語は36-1ビット, 短語は18-1ビット, アキュムレータは72-1ビットということだ.



命令語は左端の5ビットが英文字1字で表すファンクション部, 次の11ビットがアドレス部, 最後の1ビットがL/S(long or short)部で, 0だとアドレス部の短語, 1だと長語を示す.

数値は2の補数表示で, -1≤x<1の範囲の値を持つ. 左端のビットが1なら負数である. 具体的には左端の方だけを示すと, 0.012は0.25, 0.12は0.5, 0.112は0.75, 1.002は-1, 1.12は-0.5, 1.012は-0.75である.

EDSACの命令は A n F のように書く. 先頭の英字がファンクションで, つづいて十進でアドレスを書き(0なら何も書かない), 最後のコードレターという英字を書く. コードレターはアドレスの終りを示す他, FはL/Sビットが0, Dは1を示す.

A n はn番地の内容をアキュムレータに足す; S n は引く. T n はアキュムレータをn番地に格納してアキュムレータをクリア. U n は格納するがクリアせず. E n Fはアキュムレータが正ならn 番地へジャンプ, G n F は負ならジャンプで, 無条件ジャンプはなかった.

乗算は乗数を H 命令で乗算レジスタに置き, V n はn番地と乗算レジスタを掛けて積をアキュムレータに足す; N n は積を引く.

EDSACに除算命令はない.

基本的な説明は以上で終る. サブルーチンジャンプの説明も要ろう. n番地からm番地にあるサブルーチンを呼ぶには, 引数を指定された場所に入れ(T 命令で入れるからアキュムレータはクリアされている)n番地に A n F の命令を置く. n+1番地に E m F を置く. 英文字Aは左端のビットが1なので, アキュムレータは負, 従って E 命令でm番地へジャンプする.

m番地のサブルーチンに来たときは, アキュムレータに A n F があるから, これに U 2 F を足す. EDSACの文字コードでは A+U=E だから和は E n+2 F の命令になり, これをサブルーチンの最後に格納する. サブルーチンからはこの命令を実行してn+2番地へ戻る. これをWheelerリンケージという.

次は二進法の小数を十進に変換して出力する方法.

0.00012は十進では0.0625である. これを10倍する. 二進の方は0.10102, 十進の方は0.625 この時の整数部が元の小数の1桁目だが, 今は0だから0を出力.

また10倍する. 二進は110.01002, 十進は6.25. 整数部は6だから6を出力. 整数部を捨ててまた10倍する. 二進は10.12, 十進は2.5. そこで2を出力. さらに10倍して5を出力. という次第で0625が出来る.

EDSACのプリントルーチンP6は次のように書いてある.



左の0から31は相対行番号である. 欄外の24→9のようなのは, 24番地から9番地へジャンプして來るの意味. 25行目の(E F)のかっこは, この命令は変更されること. その下の横線は, 無条件ジャンプを示す. 29番地から31番地の左の2本の縦線は, 偽命令, つまり命令の形はしているが実行されず, 定数であることを示す.

最初の G K はKで終っているから, 制御指令で, 次の命令が格納される番地を相対番地の原点に設定する. つまり命令がθ(短語)やπ(長語)で終っている時は, A 3 F 命令の場所を0 とする. 以下では相対番地を'で示す.

0',1'番地はWheelerリンケージである. 帰り命令は25'番地に格納される. 2'番地は29'番地の偽命令を乗算レジスタに置く. Jは01010, 995は二進では11111000112で, 最後がFつまり0だから, 命令の形としては
01010 01111100011 0
であり, 数値的には
(/ #b01010011111000110 65536.0) => .655364990234375
だから, 216/105を上に丸めたものである.

0番地には短語の範囲に5桁の整数があり, それを5桁で出力するのだが, 最大の99999と最小の1とを見ると,
1 ]=> (* (/ 99999 65536) .655364990234375)

;Value: .9999976144172251

1 ]=> (* (/ 1 65536) .655364990234375)

;Value: 1.00000761449337e-5
だから, 純小数にして5桁出力する方針である.

4'番地の T 4 D でこうして出来た長語を4,5番地に置く. 5',6'番地は3'番地の命令V F を0番地に入れる. Vは11111なので, コメントのように-1/16を置くことになる. 0番地の負は, 出力の先頭の0をスペースにする目印である. 7'番地は乗算レジスタに10/16を置き, 10倍の準備をする. 8'番地は6'番地の命令を0から引き, Tが00101, 5なので1番地のカウンタを-5/16にする. 10'番地は4,5の長語に10/16を掛ける, すなわち10倍して小数点を4ビット右へ. つまり整数部が先頭の5ビットに収まる.

11'番地ではアキュムレータを残したまま積を戻し, 12'番地の A F で6'番地で入れた0番地の-1/16を足す, つまり整数部が0だったら-1/16になり, 26'番地へ進む. 26'番地で先程引いた1/16を足し, O 31 θ で空白を出力, 20'へ戻る.

整数部が0でなければ, 13'番地のジャンプは起きず, 14'番地の T F でアキュムレータをクリアし, 15'番地の T F で0番地を0にする.

16'番地の O 5 F は, 4番地の長語の先頭の5ビットが数字のコードなので, それを出力する.

17'番地で4番地の長語を取り出し, F 4 F では O 5 F で出力レジスタにいれた数字のコードを4番地の先頭の5ビットに取り出す. 残りのビットは0. それを19'番地で引き, 整数部をなくしてから, L 4 F で左へ4ビットシフト, 4番地へ戻す.

EDSACのシフトは難しい. 命令の下位から見て最初の1のビットが現れるまでシフトする. だから L D は1ビット, L 1 F は2ビット, L 2 F は3ビット, L 4 F は4ビットシフトになる. それより上に1のビットがあっても影響しない.

22'番地からは5桁出力したかのチェック. 最初に入れた1番地の-5/16から3'番地の-1/16を引く. 初めの4回は負なので9'番地へ戻る. 5回済むと25'から主ルーチンに戻る.

思わず長い説明になったが, EDSACの頃のプログラムはこういう感じであった.

このプログラムの問題点は0を出力すると, 5個の空白になることだ. 当時はメモリーが少く, プログラムを短かくするのが重要であった.

EDSACについては, シミュレータのTutorial Guideを参照されたし. その23ページにP6の記述がある.

2012年9月9日日曜日

多面体描画道楽

このタイトルの前回のブログ(2012年7月21日), 大二十面体の絵はもっともらしいと思ったが, 実は違っていた. 最後の方の図の三角形ABCの面は, 芯だと考えた正十二面体の面であって, 星形の面が20あるのはいいのだが, 十二面体の12は余計であった. というより, その図でBはもっと深い位置にあるべきであった.

間違いに気づいたのは, 面を色分けにしようとして, 面の数が多すぎると分かったからだが, 正しいB点の位置の計算に手間取っていた. その計算もやっと計算が終わったら, 次は隠面消去のアルゴリズムに手間取った. 秋風も吹き始めた今頃になって, どうやら正しい図を描くことが出来た.

まず前回の失敗の図だ. 頂点0の正5/2角錐が立つ正5角形は存在しない.



それに対して, 正しい図は次のとおり.



20枚の面のうち, 頂点0に関わって, しかも見える6枚の面に色をつけたのが次の図だ.



先ほどの正5/2角錐は, 底面が平面ではなくなったが, 一応角錐ということにすると, その底は逆に窪んだ5角錐になり, 見ている範囲でも3色別々の3枚の平面で出来ていることが分かる.

結構面倒な図形であった.