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でも書くためである.

それはまた次回に.



2017年11月14日火曜日

Mercedes Euklid

前回のこのブログでは, Mercedes Euklid Model 1の使い方の説明をしたが, 今回は演算機構の中心とでもいうべき比例梃子の話しをしたい.

この下がインターネットから拾ってきた比例梃子の原理図である. 中央に横向きに並んだ0番から9番の10本のラックレールが見える. この図では手前へ1,2,3,...,9が順に右へ移動しているが,定常状態では, 10本全部の左右が同じ位置に揃っている.





左の方, ラックの下に丸に2番と丸に5番とある斜めに描かれた板が, 比例梃子で, これも定常状態では下の端が丸に11番に位置にある.

図の上部に6個の10歯の結果レジスタの歯車があることから判るように, この図は6桁の計算機になっている. その各桁に足すべき数が, 歯車の上にある8,0,5,9,6,2で, その各桁の下の方にラックと噛んでいる歯車がある. 8ならその歯車は8番のラックに乘っている.

最左端の桁に8を足すとき, 右に見えるモーターが半回転し, それに連結した棒が右端に移動することにより, 比例梃子が11の位置から5の位置へ移動する. それに従って, それぞれのラックは, その番号の歯数だけ右へ移動する. つまり8番のラックは8単位だけ移動し, その上の歯車は8だけ回転し, その桁に8が足されるのである.

図では比例梃子が右へ移動しきった時を示すので, すべてのラックの歯が縱方向に見ると揃っているのがわかる. 左端のずれも隣とは1単位ずれている.

モーターが図の位置まで半回転して足し終わると, 右上のカップリングギアが開放され, さらにモーターが半回転し, ラックに乘る歯車が逆回転するが, 結果レジスタには影響はない.


比例梃子方式は引き算も簡単である. 足し算の時は, 比例梃子は0番のラックを中心にして回転したが, 引き算の時は, 9番のラックを中心にして回転する. モーターからの棒は, 4と5のラックの間で比例梃子に接続しているから, 今度は0番のラックが9歯移動し, n番のラックは9-nだけ移動することになり, つまり9の補数を足すことで引き算を実現している.

しかし, 補数の計算に委しい人は, 1だけ違う筈だとすぐに気付くに違いない. たしかにその通りで, mercedes euklidではその補正も組み込まれている. それにはもう1枚の図を見る必要がある.

下の図は比例梃子方式の全体図で, 横線の一番上のZ0が0番のラック, 一番下のZ9が9番のラックである. 縱線も何本がある. 中央に比例梃子がごちゃごちゃ描いてあるので分り難いが, 等間隔に縱線があるとすると, 縱の線は17本になる.




その内左の7本がAlで, 次の9本が添字なしのA, 右端の1本がArである. Aの線にはRと添えられた太い線分が見えるが, これがラックに乘る歯車である. つまり0から9の任意のラックに乗せることができ, その値だけ足すことができる.

さてAlとArにある太線が9の補数対策である. 足し算用の値数レジスタの範囲がAの9桁であるが, 値数レジスタの範囲を越えたAlの部分にも0番のラック上に歯車が用意してある. これは他のラックの上へ移動することはなく0番に固定である.

足し算の時は0のラックにあるから, これらの桁は何の影響も与えない. しかし9の補数による引き算の時は, 結果レジスタの範囲では値数レジスタの範囲を越えて繰上げを左端まで伝搬しなければならず, そのための9が並んでいるのである.

一方, Arにも0番のラックの上に太線が見える. よく見るとこれは0番のラックより上にずらせて描いてある. 実はここには0番のラックに貼り付いた特別なラックと歯車があり, 0番のラックが9歯だけ動くと, Arの軸は9/10回転ではなく, 1回転するようになっている.


その桁が1回転して繰り上がりが生じ, これが1の補正を行っている. なんとも巧妙な設計ではないか.

2017年11月9日木曜日

Mercedes Euklid

1905年にドイツで発売された, Mercedes Euklid Model 1という素敵な名前の機械式計算機があった. 自動車のメルセデスとは関係ないらしい.

タイガー計算機の修理や保守をされていた, 渡邊祐三さんの「美 機械式計算機の世界」にMercedes Euklid構造というページ(pp 68,69)があり, そこにこれは往復ラック方式であると書いていあった. 通常の説明書ではproportional leverだから比例梃子とでもいおうか. Mercedes Euklidの計算機はあまり知られていないが, 城憲三先生の「計算機械」の8,9ページにも簡單な記述がある.

なかなか面白そうな機械なので, この程そのシミュレータを作り, 感触を確かめることにした.

まずインターネットで見付たこの計算機の図である. 結構大きいものだ. これとよく似た計算機が, 東京理科大学の近代科学資料館にもあるが, それはタイガー計算機製となっている. それは実は電動化されたもので, Mercedes Euklid Model 7らしい.

この計算機(Modle 1), 巾が37cm, 奥行18cm, 高さ8cm, 重さ12Kg. 楔形の台に載っていて, 手前に傾いている. 手前の1/3くらいの上部がいわゆるキャリージで, 左右へ移動出来る.

インターネットを探すと, この計算機の シミュレータもあった. それが下だ.



上の方にOBENとあるのはoverのことで, 上から見た図である. 下のVORNEはfrontで, 正面図である. 灰色の薄い部分がキャリージに対応する. これはやってみるといちおう使えるけれども, 嬉しくないのはレジスタの数字が小さいことである. また数値の入力も面倒至極であった. というわけで自己流のシミュレータ(下図)を作った次第だ.



このシミュレータでは, 上の3/5くらいが本体, 下の2/5くらいがキャリージで, 右の方に余白があるが, そこはキャリージが動いていく場所である.

上の本体部分には, 9桁の置数レジスタMがある. 下のキャリージには, 16桁の結果レジスタP(Produktenzählwerkes)と, 8桁の回転数レジスQ(Quotientenaählwerkes)がある.

本物の計算機には, 本体側に置数レジスタに数をいれるスライダがあるが, 私のシミュレータではそういう入力法は面倒なので, 数値はキーボードから入れる.

上の本体部分の右上の円Kは, 本物の図にあるクランクハンドルのつもりである. Kと書いてある中央あたりのクリックが, ハンドルの右1回転に相当し, 加減算1回が出来る.

本体部分の左上には2個のスイッチ, ASとNCがある. ASは上に倒すとAdd Mulになり(つまりMの値がPに足され), 下に倒すとSub Divになる(Mの値がPから引かれる.) もう一方のNCのNはNormalstellung(標準状態), つまりKの1回転でQが1増え, Cは(Controlの略らしいが), Kの1回転でQが1減る.

キャリージ部分には, 左下にAEスイッチがあるが, この使用法は後述する.

このシミュレータのキャリージ下部の両端にあるLとRの矢印は, ここのクリックで, キャリージが左や右へ移動する.

Pの結果レジスタは通常の機械式計算機では, 直接数値の入力は出来ないが, Mercedes Euklidでは, 上の写真に見えるつまみ状のもので直接に入力出来る. 従って私のシミュレータでもキーボードから数値が入れられる.

M, P, Qの矢印は, それをクリックするとレジスタが零にリセット出来る. さらにMとPとでは, 矢印が黒い間はキーボードで打った数字が右から送り込まれる. MとPでは, レジスタの箱をクリックすると, 箱が黒くなり, キーボードからその桁に数字が入力できる. 負数は入力出来ない(負数は存在しない).

以下が基本的な使い方である.

加減算

加減算には置数レジスタMと結果レジスタPを使う. 例えば123+456なら, 左上のASスイッチをA(Addition)にし(クリックする度にAとSが切り替わる.) P矢印を2回クリックしてPレジスタをリセットする. (クリックが1回だと, Pが入力モードになったままである.)

次にM矢印をクリックして入力モードにし, キーボードから123を入れ, Mを再びクリックして入力モードから脱出し, Kのハンドルをクリック, 加算を1回行う.

もう1度Mに456を入れ, Kをクリックすると, 結果レジスタに和579が出来る.

ハンドルを回す度にQレジスタの内容も変化するが, いまは関係ないので無視しよう.

減算も同様だが, ASスイッチをS側にしておく. こういう計算機には負数というものはない. 引いた結果が負になると, 結果レジスタが16桁なので, 1016の補数表示になる. 例えば結果レジスタを0にしておき, 1を引くと, 結果レジスタは9999999999999999になる. 最上位まで繰下げが伝わるのは驚きである. タイガー計算機やBrunsvigaなどOdhner型の計算機では, 繰上げ, 繰下げは途中で止まってしまう.

乗除算

123を3倍したければ, 置数レジスタに123を置き, 結果レジスタをリセットし, ハンドルkを3回クリックする. 23倍したければクリックを23回やってもよいが, キャリージを右に1桁移動し, 置数レジスタが結果レジスタの10の桁に足せるようにして, 10の桁で2回足す. 乗算の積は結果レジスタに出る. 結果レジスタがPという名前なのは, (Product)だからであろう.

何回足したかわかるように, 回転数レジスタがある. キャリージがどの位置にあっても, 置数レジスタの1の位置に対応する回転数レジスタの桁で, NCスイッチがNの時は1が足され, Cの時は1が引かれる. その際, 回転数レジスタ内でも繰上げや繰下げは最上位まで正しく処理される.

従って通常の乗算の時は, ASスイッチはA, NCスイッチはNにしておく.

除算では, 披除数を結果レジスタに, 除数を置数レジスタに置き, キャリージで引く位置を調整しながら, 引き算を繰り返す. (見た目は悪いが, 結果レジスタに披除数が直接入れられるのは, 非常に便利である. 引いた回数が商になるから, NCスイッチはCにする.

披除数の最上位と除数の最上位の位置が合うようにキャリージを合わせ, 筆算でやるように, 各桁で引けなくなるまで, つまり結果レジスタ≤置数レジスタになるまで引き算を繰り返す. 商は回転数レジスタに得られる. このレジスタがQなのは, ここに商(Quotient)が出るからであろう.

短絡乗算

ある桁の乗数が8とか9とか大きい時は, 8回や9回まわすのではなく, 上の桁を1回多く足しておき, この桁では2回か1回引くことで済ませられる.

例えば256×4096なら, Mに256を置き, 3桁右シフト. 切替スイッチをAとNにして4回(4000倍を足した), 1桁左シフト, 1回(100倍を足した), 2桁左シフト, 切替スイッチをSとCにして4回まわす(4倍を引いた). この時, Qレジスタは4096, Pレジスタは1048576(=2^20)になっている.

自動除算

Mercedes Euklidは除算を容易にするために考案されたといってもいいような計算機であった. 例として1048576を144で割ることを考えよう. この除算の商は7281, 剰余は112である.

Pレジスタに1048576, Mレジスタに144を置く.

1048の下に144が來るようにキャリージを右3桁シフト.
SNにして負になるまで回転. 8回で Pレジスタは999...896576で負になり, Qレジスタは8000になる. この様子を下に示す.
Mercedes Euklidでは, 加減算の後, 結果レジスタが正から負, または負から正へ変化すると, ハンドルが固定され回転できない仕掛けyになっている. ー計算器ではベルがなるだけであったが, Mercedes Euklidではハンドル動かなくなることで,いわゆる引放し除算を実現している.



Mercedes Euklidにはさらに工夫があり, キャリージを右シフトした時に内部でスプリングを卷きあげ, 左シフトの時はラッチを外すだけで, スプリングの力で移動する. またAEスイッチをEの側にしておくと, キャリージが動くと同時にA/S, N/Cが反転する. Youtubeに除算をしている動画がある. この速さ除算が出来るのは見事というしかない. ただ白状すれば, どこでラッチを外しているか, 私にはよく判っていない.

2017年4月26日水曜日

したたり算法

この話題の前回のブログでは, 円周率を表すmixed radixの式を示し, 右の項から 正規化しながら十進の小数点以下の数字を求める話を書いた.

通常の計算ならそうするのだが, したたり算法では左から計算したい. それも 可能であって, 次のように計算が進行する.

下のtex出力の図で, 左上の式を見てほしい. 左端のvとuにすべての情報を 記憶させながら, 左側(外側)からかっこをはずす. 話の発端となった円周率の式 は, 左側の2行目のものである. 2+で始まっていたが, 外から1を掛けて標準形に してある.


最初の仕事は左かっこの外にあるvを, かっこ内の2つの項に掛けることだ. 今は 1だから左の3行目のように, vが消えた以外はなにも変らない.

次は先頭に来たuに相当する2を, 1/3の右のかっこの中にいれることである. それには この2に1/3の逆数, 3/1を掛ける. そうして出来た2×(3/1)+2をまとめて8 にする. 下から3行目のように, 1/3(8+ ... が出来て, 標準形であるが, かっこが 外側から1組減った.

また同じように1/3をかっこ内の8と2/5に掛ける. そうして出来た8/3を 2/15の逆数を掛けてかっこ内に入れる. すると右側の上の行のようになる. 後は 一々書かぬが, 同様に進むのである.

これを手でやるのは面倒なので, 以後は計算機にやってもらう. かっこ内は左からu とvと, 標準形の式にあったiの値である. 最後の項が5に相当するところ(4行目)までが 上の出力に見られる.

(8 1/3 2)
(22 2/15 3)
(160/3 2/35 4)
(122 8/315 5)
(1352/5 8/693 6)
(8818/15 16/3003 7)
(8832/7 16/6435 8)
(18782/7 128/109395 9)
(356984/63 128/230945 10)
...
これを更に進めると
(864779799556983658/40970475 67108864/450883717216034179 31)
のようになり, このuとvを掛けると3.1415926533011596になる. 小数点以下9 桁目まで合っている.

これが3.1415...なのは計算機が十進法にしてくれたからで, したたり算法では 自分で計算しなければならない. その過程が次の出力である.


最初の行は, 先程の長々とした分数のuとv, それとiから計算した残りの項 (31/63(2+...))である. u*vの整数部分を計算してみると
(floor->exact
 (*  67108864/450883717216034179
   864779799556983658/40970475)) => 3
と3になるので, 次の2行にわたる式のように, 先頭に3を括りだす. そしてかっこの中から3にvの逆数を掛けたものを引いて辻褄を合せる. かっこ 内を計算したものが最後の行である.

この最後の行のvとuと十進にしたいので10を掛けると, 整数部分に1が得られる. 3.14...の1のことだ.
(floor->exact
(* 10 67108864/450883717216034179
  2615629766097082765349437/2749482034790400)) => 1
ここでもTexで書くのは止め, 計算機にやってもらう. 次々のuとvと出力される 数である.
(2615629766097082765349437/2749482034790400
   67108864/450883717216034179 3)
(1536675519372845944725869/5498964069580800
   671088640/450883717216034179 1)
(58841914244318110336667/5498964069580800
   6710886400/450883717216034179 4)
(437921482322098289538739/109979281391616000
   67108864000/450883717216034179 1)
(136926162079932661882877/219958562783232000
   671088640000/450883717216034179 5)
(196056880918257839392441/10997928139161600000
   6710886400000/450883717216034179 9)
(60341900506756319941901/13747410173952000000
   67108864000000/450883717216034179 2)
(196925612577461046092237/549896406958080000000
   671088640000000/450883717216034179 6)
(48785647745580267174347/2199585627832320000000
   6710886400000000/450883717216034179 5)
(222531979586221607133547/109979281391616000000000
   67108864000000000/450883717216034179 3)
(8569388169424319751667/1099792813916160000000000
   671088640000000000/450883717216034179 3)
今日の計算はまずかっこを取り払い, つまりかっこの次の2を消し, この処理をi=31に なるまで進め, 今度は十進の各桁と順次に括りだすものであった. 大体は 左から右へ処理できたが, 本来のしたたり算法では, 出力できる状態に なると同時に出力するわけで, 次回はその話をしたい.