2014年12月15日月曜日

微分解析機

このブログの8月25日のにフロントラッシュの話を書き, 遊星歯車の話にもなった. また遊星歯車が話題だ.

微分解析機の加算器は差分歯車や遊星歯車を使うといわれているが, 理科大の微分解析機の加算器は見慣れない形をしていた. どうなっているのか.

下の図も一種の遊星歯車である. 左が軸方向から見た図. 同径, 同歯数の歯車6枚があり, Bの2枚とCの2枚は遊星キャリアに乗っている.

中央にあるAとDは, キャリアと同じ軸に乗っているが, A,D,キャリアは軸には固定されていない. 歯車による束縛条件はあるが, その下で自由に回転できる.

この図を左方から見たのが右の図で, 遊星歯車のBとCは見ての通り, 深さ方向にずれて配置されている. つまりAはBと左図の緑の点で接し, BとCは黒い点で接し, CはDと橙色の点で接している.


まずキャリアを固定し, Aを反時計方向に回転する. するとBは時計方向に, Cは反時計方向に回転する. ゆえにDは時計方向に回転するわけだ. 歯車の歯数が同じだから, 回転の角度も同じである.

隣同士の歯車は逆に回る. 中間歯車が1個なら同じ方向に回る. AとDは間に歯車が2個あるから, 逆に回る.

回転角を記述するのに, 反時計回りを正として, 歯車Aの回転角をφ, Dのをψ, キャリアのをθと表すことにする. そこでこれを

θ=0の時φ=-ψ

と書いておく.

キャリアがAと一緒に回転すると, この図がこのまま回転するから, DもAと同様に回転する.

φ=θ=αの時ψ=α

と書いておく.

すこし様子が分って來たが, もう少し図を示すと, 下図の0は上の図のA, 右上のB, 右下のCとDを書いたものだ. AとDが重さなっているから, A,B,CとB,C,Dと図を分割してある. 緑色の点はA,Bの接点, 黒の点はB,Cの接点, 橙色の点はC,Dの接点で, これも上の図と同様. 歯車やキャリアが回転すると, これらの点は付いて一緒に回る.


一段下って1の図はAを60度回転したものである(φ=60°.) Dは固定してある(ψ=0°.) すなわちDの橙色の点は0の図と同じ位置にある. するとキャリアは30度回転する(θ=30°.) キャリアの向きはAやDから見たBとCの接点の方向である.

図2はφ=120°, ψ=0° θ=60°の時だ.

3では今度はDを逆回転してみた. φ=120°, ψ=-60° θ=30°

4ではDを更に逆回転した. φ=120°, ψ=-120° θ=0°

つまり, この遊星歯車では φ+ψ=2θ となっている.

私の2010年1月4日のブログでは, 太陽歯車の半径r, 回転角φ, 内歯車の半径R, 回転角ψ, キャリアの回転角θとすると,

θ=φ*r/(R+r)+ψ*R/(R+r)

とした. 今回の歯車は遊星歯車を2段にし, 内歯車を太陽歯車と同様にした. だからr=Rとなり, それから上の式が導ける.

さて下は東京理科大の微分解析機の加算器の写真である. ベイとベイを繋ぐバスボックス内にある3本(手前からX,Y,Zとする)のスタブシャフトに設置してある.



中央のYに乗っているのが, 上述の遊星歯車である. 手前のXの歯車と噛んでいる平歯車と一緒になっている真鍮のがAである. その左, Aと噛んでいる手前と向こうのがB, Aからは離れているが, Bと半分の幅で噛んでいるのがC, Cの奥にDがあるのだが, それは見えていない. DはY軸に固定されている.

真鍮の歯車群の左にあるのがキャリアで, 2:1で向こう側のZ軸の平歯車と噛んでいる. φ+ψ=2θだったので, θを2θにしているわけだ.

これを図にしたのが次だ. 円内下端の数は歯数である. 軸が破線で描kてあるのは, 固定していないことを示す.



これまで何回ものブログで微分解析機の殆どの仕掛けが分ったように思う.

2014年12月12日金曜日

クリスマスの歌

今年もクリスマスが近付いた. 最近の私にそういう機会はないが, 以前はクリスマスキャロルに誘われて聞きにいったりした.

そこで歌われる歌の中に, On the first day of Christmas,... で始まる歌がある. 出典によって文句や順番が多少違ったりするが, 大体はこうだ.

On the 1st day of Christmas,
my true love sent to me:
a Partridge in a Pear Tree.

On the 2nd day of Christmas,
my true love sent to me:
2 Turtle Doves
and a Partridge in a Pear Tree.

クリスマスの1日目, 私の恋人は
ナシの木の1羽のヤマウズラをくれた.

クリスマスの2日目, 私の恋人は
2羽のヤマバトと
ナシの木の1羽のヤマウズラをくれた.
つまり
クリスマスのn日目, 私の恋人は
n個のなにか, 
n-1個のなにか,
...
2羽のヤマバトと
ナシの木の1羽のヤマウズラをくれた.
と12日目まで続く.

結構規則的なので, 記号処理プログラムで生成するのが簡単であり, その昔Snobol 4の手引き書に出ていたのを覚えている. さらにこの歌詞をAlgol 68での階乗のプログラムにしたという話もあるが(The Most Contrived Factorial Program, Allgol Bulletin 42 (1978)), それはまたいつか書くこともあろう.

私自身はエディタのTecoで書いたことがあった. (bit Vo. 14, No. 10, 1982) 最近のことだが, GPMでも書いてみた. いろいろ苦労したが, 超絶技巧な技を利用. それについてもまたいつか書きたい.

TecoとかGPMのようなプリミティブなプログラム言語とは別に, Schemeで書くとどうなるか. やはり超簡単の単であった. 以下の通り.
(define a '
("12 Drummers Drumming\n"
 "11 Pipers Piping\n"
 "10 Lords a Leaping\n"
 "9 Ladies Dancing\n"
 "8 Maids a Milking\n"
 "7 Swans a Swimming\n"
 "6 Geese a Laying\n"
 "5 Golden Rings\n"
 "4 Colly Birds\n"
 "3 French Hens\n"
 "2 Turtle Doves\nand "
 "a Partridge in a Pear Tree.\n"))
(define b '
("12th" "11th" "10th" "9th" "8th" "7th" 
 "6th" "5th" "4th" "3rd" "2nd" "1st"))
(do ((i 11 (- i 1))) ((< i 0))
(display (string-append
"\nOn the " (list-ref b i) 
" day of Christmas,\nmy true love sent to me:\n"
(apply string-append (list-tail a i)))))
ところで私は何となく贈りものの数を1+2+...+12=78だと思っていたが, ヤマウズラは12日間贈り, 2羽のヤマバトは11日間贈るので, 総計は1*12+2*11+3*10+...+12*1だったらしい. 丁度途中の6*7+7*6のところで折り返すから,1*12+2*11+3*10+4*9+5*8+6*7まで計算し2倍する.

4*9+5*8+6*7は簡単かもしれない.


赤の枠が5*8, 緑の枠が4*9. 青の枠が6*7である. 青の斜線部分を左に1だけ移動し, 赤枠内に移動すると, 緑の斜線部分と合せて赤枠の面積マイナス1になり, 緑の残り部分と青の残り部分も赤枠-1になるから, 5*8*3-2=118だ.

1*12+2*11+3*10も同様に2*11*3-2=64. 全体は(118+64)*2=364だ. 1年の日数に1日足りないのか.

2014年12月2日火曜日

EDSACのプログラム技法

前回のブログではy' を計算するサブルーチンにジャンプするところまで述べた.

今回はそのサブルーチンの中身の説明から始めよう.



0行目と1行目はlinkageで主ルーチンにもどる命令を作り, 22行目に入れる.

2行目は8Dにあるt' をアキュムレータに取り出す. 3行目のEはアキュムレータが正ならジャンプする条件ジャンプ命令なので, t' が正なら6行目(6番地)へ行く. 負の場合は8Dを2回引いて正にし, 6行目で合流した時には絶対値になっている.

平方根をとるサブルーチンS2は, 被開平数を4Dで渡すのでT4D. そしてS2へジャンプ. 結果は4Dにある.

9行目からはt' の3乗を計算する. まずH8Dで8Dにあるt' を乗数レジスタに置く. V8Dで2乗をつくり, TDで0D番地へ置く. VDで3乗が出来る.

13行目のRDはアキュムレータの1ビット右シフトである. EDSACのシフト命令のシフト数は, 命令の一番右のビットの位置で決まる. だからRDなら1ビット右シフト. R1Fなら2ビット, R2Fなら3ビット, R4Fなら4ビットである. 番地部が2nならn+2ビットシフトする. 命令語の上の方のビットは影響しないから, 1ビットシフトは番地部が奇数ならなんでもよい.

14行目からは5倍の作業だ. まずUDで0D番地に入れる. UはTと同様な格納命令だが, アキュムレータをクリアしない. L1Fで左へ2ビットシフトで4倍する. それに0Dにあった元の数を足して5倍を実現する. それを一旦0Dへ退避し, 4Dにあった平方根をとりだし, R512Fで(512は2^9なので)右へ11ビットシフトし, 2-11倍する. それに5.2-1t' 3を加え, y' を8Dに置いて主ルーチンに戻る. T8Dをやったので, アキュムレータは0になり, 0は正なので, 正ジャンプのE命令で戻ってきた.

以後の作業には前回のブログのプログラムを見てほしい.

18行目. y' を乗算レジスタに入れる. 19行目. アキュムレータに取り出す.

次はyが400を超えるときはtoo largeと出力する仕様なので, y' を400.2^-13と比較する.

アドレスのところが2^-15だから, P400*2^2Fと比較すればよく, 5Mのような定数になる. ただ5MにはP1600Dとあるから, 400.125*2^-13と比較している. ちょうど400の時はtoo largeにしたくないというつもりだろうか.

21行目のGは負ジャンプ. 400を超えた場合は22番地で999.2^-13を乗算レジスタに入れる. これも定数はアドレス部に999の4倍が書いてある. (too largeと出力する代わりに999で代用する.)

ところで, これは正の方向に400を超ているを見ているだけで, -400より小さいときは別になにもしないようだ. 不思議な仕様である. 他の例を見ても正の場合しか考慮していない.

続いて213=8192を掛けるのだが, 全体を小数にする必要があるので, 十進の小数の4桁目に小数点が來るように, 213/10000を掛ける.

乗算VnDは, 乗算レジスタにある数とnDにある数を掛けてアキュムレータに足すから, TDでアミュムレータをクリアしておく. そして2MDにある数を掛ける. しかし2MDのようにコードレターを2つ書くことは, EDSACのイニシアルオーダーでは出来なかったので, Dの代わりにπを書いていた.

2Mと3Mにあるのが, 213/10000である. なぜそうかというと

0.819210= 0.110100011011011100010111010110001110001000011001011012

これをEDSACの短語に対応させてみると, EDSACでは偶数語が右, それ+1の奇数語が左なので



のようになる. この図で見ると2Mの語の最後はDだが, プログラムではFなのはこの辺で切り捨てたのであろうか.

掛け算が済んだので, 0Dに入れ, 小数の印刷ルーチンP14を呼ぶ.

このP14は28行目にあるプログラムパラメタで, レイアウトできる.

例えば小数0.123456789を
のように出力するには, スペースの位置の下の数と, 最後の桁の上の数を足したものがパラメタの値であり, 28行目は3072+32=3104となっている.

この後はt を取り出すアドレスと, カウンタi の変更が残っている.

29行目は14行目の命令を40D, 38D, ...と減らすもので, A14θ, S 9M, T 14θ で変更する. 32行目からは, A 7M, S 8M, U 7Mとカウンタを減らし, 35行目で, アキュムレータに正の数が残っていれば, 3番地へ戻る. 負になったらジャンプせず, 36行目のZ命令で停止する.

論文にあった出力はこうだ.
10  0000 04472
 9  0000 03162
 8  0999 00000
 7 -0044 85586
 6  0047 75414
 5  0999 00000
 4  0062 35157
 3  0999 00000
 2 -1077 55051
 1  0999 00000
    0018 09974
最下行 0 がプリントされていない. 論文には印刷ルーチンの虫だろうと書いてある. そうだろうか.

私がSchemeで計算したのは次だ.
(define (tpk t)
(+ (sqrt (abs t)) (* 5 (expt t 3))))

(define data 
'(1.5 8 -6 9.5 2.3 9.9 2.1 -2.1 6 0.001 -0.002))

(map tpk data)
=>
(18.09974487139159 2562.828427124746 -1077.5505102572167
 4289.957207001485 62.35157508881029 4854.641426544511
 47.75413767461895 -44.85586232538106 1082.4494897427833
 .03162278160168379 .044721319549995794)
という次第でこのプログラムは解明出来た.

いまから思うと大変な手間であったが, 私がパラメトロン計算機PC-1でプログラマー人生を始めた頃はまったくこの通りであった.