2018年8月22日水曜日

e1000桁の計算


60年くらい前によく行われていた, Napier's constant eの値を1000桁計算するパラメトロン計算機PC-1のプログラムについてのファイルがあったので, 私のブログに載せておくことにした.

自然対数の底 e の値を小数点以下1000桁まで計算することは,1960年代の計算機で流行した. Taylor展開

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

で計算するのだが, n = 460くらいで1000桁の精度がある.
   
実際には上の式で計算するのではない.

e = 1 + 1/1(1 + 1/2(1 + 1/3(1 + ... +1/459(1 + 1/460)...)))

を使う. さらに例えば4項ずつまとめて計算することを考える.

Qa-4 = 1 + 1/(a - 3)(1 + 1/(a - 2)(1 + 1/(a - 1)(1 + 1/a・Qa)))

a = 460, 456, 452, ... 0

Q460 = 1

e = Q0

これはさらに簡単になって

Qa-4 = {a(1 + (a - 1)(1 + (a - 2)(1 + (a - 3)))) + Qa}

/{a(a - 1)(a - 2)(a - 3)}

しかし 460×459×458×457 = 44192868840 となり, 36ビットのPC-1のレジスタの容量235 = 34359738368には収まらない. 連続する4個の整数の積は4で割れるから, 2ビット右シフトして計算すれば, この問題は解決する.以下のプログラムはそのようにしてeの値を計算している.

次は10進数での出力である. 小数を10進法で出力するには, 全体を10倍し,整数部分にでた数を印字する. これでは乗算に時間をかけて1桁しか印字できないので,以下のプログラムでは, 小数部分に一度に1010をかけ,10桁分の整数部を作り出す. また多倍長の乗算を10区間に分割し,各区間毎に10桁の整数を1桁ずつ印字するようになっている.

作業場所は10進で1000桁なので, 3300ビット必要. PC-1の36ビットの長語を100語使う
(300〜498 300番地の方が小数点側).

作業場所は10進で1000桁なので, 3300ビット必要. PC-1の36ビットの長語を100語使う (300~498 300番地の方が小数点側).

0h=80,0i=100,0p=180,0n=272,0m=300,

70:
 0   it
 1   jl   h   作業場所クリアへ
 2   it
 3   jl   i   計算プログラムへ
 4   it
 5   jl   p   印字プログラムへ
 6   0

0h:
 0   a   8n ┐
 1   x  13r ┘ plant link
 2   p   2n
 3   x   7r
 4   s   3n
 5   z  11r
 6   l  35    アキュムレータを0に
 7   tl(  )   作業場所をクリア
 8   p   7r
 9   a   8n
10  jl  3r
11   pl 16n   Q460を1に
12   tl   m
13   jl(  )

0i:
 0   a   8n ┐ e1000桁計算ルーチン
 1   x   3r ┘ plant link
 2   pl 18n ┐
 3   zl(  ) ┘ a = 0なら帰る
 4   sl 16n   a - 1
 5   tl 20n  
 6   sl 16n   a - 2
 7   tl 22n  
 8   sl 16n   a - 3
 9   tl 24n  
10   vl 22n   ×(a - 2)
11   ll 34   
12   vl 20n   ×(a - 1)
13   ll 35
14   vl 18n   ×a
15   ll 34    a(a-1)(a-2)(a-3)
16   tl 26n   を26nLへ
17   p   2n
18   x  22r
19   x  25r
20   s   3n
21   z  29r
22   ql(  ) ┐
23   r  33  │ Qを4で割る
24   ll 33  │
25   tl(  ) ┘
26   p  22r
27   a   8n
28   jl 18r
29   pl 24n   a - 3
30   al 16n   + 1
31   vl 22n   x(a - 2)
32   ll 35
33   al 16n   + 1
34   vl 20n   x(a - 1)
35   ll 35
36   al 16n   + 1
37   vl 18n   x a, aは4の倍数なので
38   ll 33    最後の2ビットは0
39   al   m   Rレジスタは0
40   tl   m
41   p   2n
42   x  46r
43   x  48r
44   s   3n
45   z  52r
46   ql(  ) ┐ 多倍精度/単精度
47   dl 26n │ a(a-1)..(a-3)で割る
48   tl(  ) ┘
49   p  46r
50   a   8n
51   jl 42r
52   pl 24n
53   sl 16n
54   tl 18n
55   jl  3r

0p:
 0   a   8n ┐ 出力ルーチン
 1   x  68r ┘ plant link
 2   q  16n
 3   p   4n
 4   x   8r
 5   x  12r
 6   s   5n
 7   z  16r
 8   pl(  ) ┐  
 9   wl 12n │
10   tl 28  │ 1回目の1010
11   ql 28  │
12   tl(  ) ┘
13   p   8r
14   s   8n
15   jl  4r
16   p   4n ┐
17   tllt   │ 2. と印字
18   s  70r │
19   tllt   ┘
20   p   6n   20行用カウンタ設定
21   x   9n
22   p   7n   5ブロック用カウンタ設定
23   x  10n
24   z  65r
25   pl 28
26   tl 68
27   p    n   1010倍10回カウンタ
28   x  11n   初回は0に設定
29   q  16n
30   p   4n
31   x  39r
32   x  43r
33   s   5n
34   z  62r
35   p  11n
36   z  47r   10回カウンタを調べる
37   s  17n
38   x  11n
39   pl(  ) ┐
40   wl 12n │
41   tl 28  │ 1010
42   ql 28  │
43   tl(  ) ┘
44   p  39r
45   s   8n
46   jl 31r
47   ql 68    1010倍10回終り
48   tl 28    1文字 印字へ
49   p    n
50   dl 14n   109で割る
51   l  18
52   a  56r
53   x  54r
54   p (  )   文字コードを取る 
55   tllt     印字
56   q   2n   アドレス部分は52で使う
57   vl   n
58   ql 28
59   tl 68
60   p   1n
61   jl 37r
62   p  10n
63   s  17n
64   jl 23r
65   p   9n
66   s  17n
67   z (  )   link
68   x   9n
69   p   8n
70   l  12    CR
71   tllt
72   l   2    LF
73   tllt
74   r   1    SP
75   tllt
76   tllt
77   jl 22r

0n:
 0   0       ┐
 1   10      ┘ 10(長語)
 2   229676   (0) 300 以下文字コード
 3   233972   (1) 500 と定数
 4   197106   (2) 498
 5   172332   (3) 300
 6   135188   (4)  20
 7   217093   (5)   5
 8   245762   (6)   2
 9   180224   (7) (  ) カウンタ
10   143360   (8) (  ) カウンタ
11   184320   (9) (  ) カウンタ
12   38146   ┐
13   254976  ┘ 1010 
14   3814    ┐
15   182784  ┘ 109
16   0       ┐
17   1       ┘ 1 (長語)
18   0       ┐
19   460     ┘ a (長語)
                この後27nまで作業場所
jl70.

2018年4月8日日曜日

PDP8のプログラム技法

前回のブログで基本ループの動くのは分ったが, calのプログラムを書くにはまだmon0, mon1のデータを用意する仕事がある. そのプログラムの話だ.

mon0は各月の1日の曜日が決ったとき, その週の最左端の枠, つまり日曜の日付けである. 1日の曜日の計算は月の定数にその年の定数と(日付けの)1を足す.

その年の定数は, 年の下2桁をyとして, (y+floor(y/4)+世紀の定数)%7である. 世紀の定数は, 年の上2桁を4で除した剰余の0,1,2,3について, それぞれ+4, +2, 0, -2である.

この辺までを計算するのが次のプログラムの000行から034行である.



000行はプログラムの200番地からの格納の指定. 001行はアキュムレータとリンクのリセット. 002,3行で1文字(1数字)を読みaへ置く. 009行までで4文字(4数字)を読み, 1000,100,10,1の各桁がa,b,c,dへ入る. 次のgetchは改行を読み, eへ入れるがこれはそれを捨てアキュムレータをクリアする.
 
012行からが上のy+floor(y/4)の計算. 016行まででcを10倍し, dを足す. それをeへ格納し, 20行までで2桁右シフト. それにeを足して計算完了.

024行からは上2桁の4で除した剰余をとる. それにはa+a+bを3でマスクする. 一旦fに入れ(028行), 031行までで剰余を5倍し, 世紀内の定数eを足しさらに033行で5を足す.

0,1,2,3から7を法として4,2,0,-2を作るのに, 5倍して4を足した. つまり4,9,14,19は7を法として4,2,0,-2だからである. 4を足すのではなく5を足すのは, 1日の分を足すからである. これをyconstに置く(034行.)

035行から051行は閏年の判定である.



まずcとd(下2桁)を足し, 和が0(100の倍数)ならhへ. そうでないならc+c+dを持ってjへ. hではa+a+bを持ってjで合流. これらが4の倍数なら閏年である. 3でマスクし, 0でなければ048行で-1を置き, 閏年なら049へスキップしてきて1増やすから, 閏年なら1が出来, 平年なら-1が1増えるから0が出来る.

050行は上の結果をマイナスにし, 051行でそれをmleapへ入れる.
 
これまででyconstとmleapが出来たので, mon0とmon1の表の修正を開始する.



修正はまず月の定数のうち, 1月と2月から1を引くのと, 月の日数のうち, 2月を1増やすことである. 052でmon0のアドレスをeに入れ, 054で1月の値を取り出し, mleapを足してもとへ戻す. iszでアドレスを増やし, 2月についても同様にする.

061行からは2月の日数の調整である. mon1の表は負数になっていたから, 2月の値から1引くわけで, それが065行だ.

mon0はyconstを足した(072行)のが1日の曜日であった. これからその月の最初の枠の日付けの計算に移る. 067行のm14は十進では-12のことで, 12回ループの準備をする.

この計算はmon0の値にyconstを足し, 7の法をとり(073行,) その後, 各月の1日の値がnなら, -n+1で置き換える. -nの操作が074行. 1足すのが075行である. これを実行すると, 前回のブログにあったmon0とmon1の表が得られる.

次の081行から097行は定数で, m7770, m7774はand命令で使うマスクビットである.


 
085行からのmon1は各月の日数だが八進法なのでちょっと異様にも見える.

098行からは1文字読込み(getch)と7の剰余をとるサブルーチン(mod7)である.



getchでは西暦の年号を読むだけなので, 数字のコードの八進260を引く.
 
mod7は法をとる披除数をアキュムレータに置き, 上9ビットをa, 下3ビットをbとし, aが0でなくなるまでaを3ビット右シフトしてbに足す操作を繰り返す. aが0になったらbからを引き, 結果が正, つまり0ならそのまま, 負なら7を足して戻る.

これでcalのブログラムの主要な部分は大体完成した.

2018年4月1日日曜日

PDP8のプログラム技法

私がいろいろなプログラム言語で書いたcalのプログラムは, unixのcalのコマンドで印刷されるのと同じものだ.

例えば

cal 2000

とすると, 2000年の12ヶ月のカレンダーが印字される.


こういうプログラムをpdp-8用に作るための工夫を書いてみたい.

かなり前のことだが, アメリカにJerry Weinbergという心理学者が活躍していて, 私も何回か会ったことがある. 当時プログラムを書くのにトップダウンがいいかボトムアップがいいかという無益な議論があった. その時Jerryが「難しいところから書く」を主張していて, 私と同じ考えの人がいると思った.

というわけで, 一番中心部分のプログラミングから始めよう. 見出しの類いはあとからゆっくり追加すればよい.

まずは各月の第1日の曜日と各月の長さが必要である.

このブログでも何度か述べたように, 1月から12月について次の表が重要になる. (ある月の定数は前の月の定数に前の月の日数を足し, 7による法をとったもの(と等価なもの)だ.)

 
 
1234567819101112
2558361417250


この表とその年の定数から各月の0日の曜日が判る. 2000年の定数は4なので, 例えば2000年8月15日は(4+4+15)mod 7=2で火曜だ. 注意すべきは2000年が閏年だったので, 1月と2月はこの表から1引くことである. 従ってこの表は補正後は, 1,4,5,8,3,... のように始まる.

それに年の定数を足して7の法をとると, 5,1,2,5,0,... になる. すると1月1日は4+1+1=6で土曜である. 1日が土曜だと, 金曜が0, 木曜が-1, 水曜が-2, 火曜が-3, 月曜が-4, 日曜が-5に相当し, 1月の第1週は-5日から始まる. そこから1日ずつ増えて, 1になるのが土曜である. 従って次の図のmon0と書いた表の最上段(1月)に7773, つまり10000-5を入れる. mon0はこのように出来ている. 年の定数が変るとか, 年の定数が同じでも閏年であったりするから, この表の内容は年により違う.





右側にあるmon1の表は, 上から順に各月の日数の負の値である. 日付が順に増えてゆき, 日数を過ぎると表示を空白にしなければならず, それには日数を減算するが, PDP-8には減算がないので, 最初から負数にしてある.

最初のカレンダーの構造からわかるように, このプログラムは大局的には3ヶ月のカレンダーの4回の繰返しである. その1回分は, 1ヶ月の3回の繰返し; 1ヶ月は週の6回の繰返し, 週は日の7回の繰返しなので, 以下のプログラムは4重のループで出来ている.



00行はこのプログラムが八進の200番地から格納することを示す. 02行からmon0pの値をaへ, mon1pの値をbへ置く. 06行から4,6,3,7のカウンタを初期化する. 15行 aの指す内容(今は1月の先頭の日付け)の1未満を調べるため, アキュムレータに-1を置きmon0の値いを足す. それが負なら空白を出力するためl1へ. (0も-1に  なるので, 飛ぶ) 次にこの月の日数を引き(18行), 正になればこれも空白出力へ進む.

日付を出力できる条件を満すと, 21行からそれを出力し26行へ行き日付を進める. 日付の値は負から始まり, 正へ進むが, isz命令は0になった場合, 次の命令をスキップするから, isz命令の次にjmp .+1かnop(7000)を置く.

この後1週間分, つまり7回出力したかを見, そうなら次の月の1週目に進むから, aとbのポインタを進める(31,32行目.) 33行目はそれを3ヶ月やったかの確認で, 3ヶ月済むとまた左端の月に戻るので, 35行目からaとbを3ずつ減らす. これで1行出力したから, 改行する(41行目.)

42行めは6週間のテスト. 市販のカレンダーでは1日が金曜か土曜だと, 30日や31日は23日や24日と一枠を共有するが, calではそうせず, 次の行にするので, 1ヶ月は最大6行になる. その6行が終ったかを調べる.

それが終わると続く3ヶ月に移るので, 44行目からは, aとbを3増やす段取りだ.

 

53行目からは定数や作業場所, 基本的なルーチン類になる. 54行, 57行が上の図のmon0とmon1である. 59行からはカウンタの初期値, 64行からのa, b, c, d, e, fが作業場所になる.

このプログラムの出力が以下だ. 横7文字目, 14文字目の次に縦線を引き, 6行, 12行, 18行の下に横線を引くと, 最初のカレンダーと同じ情報が得られる. 1日がA, 2日がB, 26日がZだから, アルファベット各文字が何番目かを熟知しているひとには, カレンダーとして十分役立つと思われる(^^).

@@@@@@A@@ABCDE@@@ABCD
BCDEFGHFGHIJKLEFGHIJK
IJKLMNOMNOPQRSLMNOPQR
PQRSTUVTUVWXYZSTUVWXY
WXYZ[\][\]@@@@Z[\]^_@
^_@@@@@@@@@@@@@@@@@@@
@@@@@@A@ABCDEF@@@@ABC
BCDEFGHGHIJKLMDEFGHIJ
IJKLMNONOPQRSTKLMNOPQ
PQRSTUVUVWXYZ[RSTUVWX
WXYZ[\]\]^_@@@YZ[\]^@
^@@@@@@@@@@@@@@@@@@@@
@@@@@@A@@ABCDE@@@@@AB
BCDEFGHFGHIJKLCDEFGHI
IJKLMNOMNOPQRSJKLMNOP
PQRSTUVTUVWXYZQRSTUVW
WXYZ[\][\]^_@@XYZ[\]^
^_@@@@@@@@@@@@@@@@@@@
ABCDEFG@@@ABCD@@@@@AB
HIJKLMNEFGHIJKCDEFGHI
OPQRSTULMNOPQRJKLMNOP
VWXYZ[\STUVWXYQRSTUVW
]^_@@@@Z[\]^@@XYZ[\]^
@@@@@@@@@@@@@@_@@@@@@

2018年3月24日土曜日

PDP8のプログラム技法

今から50年程前にDEC(Digital Equipment Corporation)が発売していたPDP8というミニコンがあった. 思うにこれはMITのTX-0を源とする似たようなアーキテクチャの計算機である.

簡単に概要を示すと, 1語は12ビット, 記憶装置は12ビットの4096語. 1ページが128語の32ページで構成される.

演算装置の中心は12ビットのアキュムレータで, 2の補数(4096の補数)で演算される. アキュムレータの左に1ビットのリンクがあり, 加算の結果繰上がりがあると, リンクのビットは反転される. アキュムレータのシフトは回転シフトであり, リンクを経由する. アキュムレータのクリアはリンクに影響しない.


図の2段目は記憶装置参照命令で12ビット. 下の図のような構成である.

左端の3ビットが命令コード. 右端の7ビットがページ内アドレスである. ビット4が0の時は第0ページの番地を示し, 1の時はその命令のあるページ内でのアドレスを示す.

ビット3が1の時は, その語の内容を番地とする.

記憶場所の八進法で0010から0017までの場所を間接番地で参照すると, この番地の内容を1増やしてから使われるから要注意である.
 
記憶場所参照命令は, 命令コードの0から5が次のようである.
0 AND
1 TAD  
2 ISZ
3 DCA
4 JMS
5 JMP

ANDは番地の指す記憶場所の内容とアキュムレータのビット毎のANDをとり, 結果をキュムレータに置く.
 
TADは番地の指す記憶場所の内容をアキュムレータに加える. 繰上がりはリンクを反転する.

ISZは番地の指す記憶場所の内容を1増やし, その結果が0になれば, 次の命令をスキップする.

DCAはアキュムレータの内容を番地の指す記憶場所に入れ, アキュムレータをクリアする.

JMSはこの命令の次ぎの番地を番地の指す記憶場所に入れ, その次の命令に進む. つまりサブルーチンジャンプする.

JMPは番地の指す記憶場所の命令に進む.

命令コード6は入出力命令で, テレタイプのキーボードから読んだり, プリンタに印字したり, 紙テープリーダから読み込んだりする. この話しは省略する.

命令コード7はマイクロプログラム命令で, 図に示すようにビット3の0, 1に従って, グループ1とグループ2とがある.

グループ1はビットの4,5,6,7が1だと, アキュムレータをクリア(CLA), リンクをクリア(CLL), アキュムレータの各ビットを反転(CMA), リンクを反転(CML)する.

ビット8, 9はアキュムレータとリンクを繋げた左(RAL)または右(RAR)回転で, ビット10の0/1に従って1ビットか2ビット(RTL, RTR)になる. ビット11が1だと, アキュムレータを1増やす(IAC).

マイクロプログラムでは実行順が重要である. 最初がクリア, 次が反転, 次がIAC, 最後が回転となる.

グループ2は, ビット4が1の時アキュムレータのクリア(CLA), ビット5,8が1,0の時アキュムレータの負(ビット0が1)でスキップ(SMA), 1,1の時アキュムレータの正でスキップ(SPA). ビット6,8が1,0の時アキュムレータの=0でスキップ(SZA), 1,1の時アキュムレーの≠0でスキップ(SNA). ビット7,8が1,0の時リンクの≠0でスキップ(SNL), 1,1の時リンの=0でスキップ(SNL)である.

ビット9が1の時, コンソールの12ビットのスイッチをアキュムレータとビットごとORする. ビット10が1の時計算機を停止する(HLT).

実行順はまずスキップ, クリア, スイッチ読込みと停止である.

プログラムでは, CMA IACのように書くだけである. こう書くとアキュムレータがビット反転され, 1足されるので, アキュムレータの値を負にする. PDP8には減算命令はなく, A-BはBをアキュムレータに置き, CMAで1の歩数を取ってからIACで1を足す.

ところでPDP8ではプログラムは原則八進法で書く. 格納開始番地は星印(*)を使い
*200
のように書く. 記号アドレスの定義は
start,
のようにコンマ(,)を使う. 命令のアドレス部には記号アドレスや, この命令の番地を示すピリオド(.)や, プラス, マイナス(+, -)に数を続ける.

斜線(/)はコメントの開始を示す. プログラムの最後には($)を書く.

例えばアキュムレータにある正を数を十進法の10で割り, 商をq, 剰余をrに置くサブルーチンdiv10を八進の100番地からに入れるには

*100
div10,0  /サブルーチンの入口. 帰り番地が入る
dca r    /披除数をrに入れる
dca q    /商の場所qをクリアする
tad r    /部分剰余をアキュムレータに置く
tad m10  /10を引く
spa      /剰余が正ならスキップ      
jmp .+3  /負になったので後始末へ
isz q    /剰余を1増やす
jmp .-4  /4命令前へ戻る
tad p10  /引き過ぎた剰余を戻す
dca r    /剰余をrに入れる
jmp i div10 /主プログラムへ戻る
m10,-12  /定数 -10
p10,12   /定数 10
r,0      /剰余の場所
q,0      /商の場所
のように書く.

PDP8はこの程度の計算機だから, シミュレータやアセンブラを書くのも簡単だ. 今回PDP8を説明したのは, 私がいろいろなプログラミング言語で書いている例題, 与えられた年のカレンダーを出力するcalをPDP8でも書くためである.

それはまた次回に.