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でプログラマー人生を始めた頃はまったくこの通りであった.

2014年11月30日日曜日

EDSACのプログラム技法

英国の計算機歴史学者Martin Campbell-KellyがProgramming the EDSAC: Early Programming Activity at the University of Cambridgeという長い長い論文をIEEE Annals of the History of Computing Vol.2 No.1 January 1980に書いている.

この最後の方にInfluence on programming for other machinesという章があり, 東大のTACのことも述べているが, E. Wada, who was responsible for programming, designed a new set of initial orders based on those of Wheeler, but incorporating floating addresses and some other small refinements.と書かれているのを知った.

今回のブログの目的はその論文にあるEDSACのプログラムの例を説明することである. MartinはEDSACの例にTPKプログラムを使った. なんの変哲もない問題だが, アセンブリ械語のプログラムとしては長さ的に適当かも知れない.

11個のデータを読み込み, 配列aの0から10に置き, t =a[t ]として,

y =√|t|+5×t 3 (*)

を計算, 印刷するものである.

t が<1なら, 最初の項は平方根だからtより小さくなるから無視し, 次の項は3乗なのでこれが結果の大きさに利く. t が±10位ならyは±5000程度のt3の曲線を考えればよい.

EDSACの頃の計算機は, 扱う値を±1以下としていたから, tが10程度までなら1/16にし, yも1/213=1/8192として計算する.

つまり *の式にy=213y', t=24t'を入れ, y' =2-11t' +5×2^-1t'3 を計算し, 後で213倍する. t=1とすると, t'=2^-4. その平方根は2^-2. 2^-11を掛けると2^-13. またtの3乗は2^-12. 5.2^-1を掛けて5.2^-13. 両方を足して2^13を掛けると6になる.

EDSACに用意してあるライブラリのサブルーチンもみなこの範囲と思って働く.
さてMartinの書いた主プログラムはこうだ.



上から数行目の 0 |A θ と書いてある行を0番地としたアドレスにはθを最後に付ける. 下の方 M 0|P 4 D の行を0番地とするときは, Mを付ける. そのMは38θであること, set M-parameterの行で設定する.

0と1行目はWheelerリンケージで, 56絶対番地のR1へサブルーチンジャンプする. 2行目 T20Dはサブルーチンに20Dからデータを格納せよとの指定である.

データはテープに15+8+6-95+23+99+21+21-6+0001+0002-Fとパンチしてあり, Fが來るまで0.15, 0.8, -0.6, ...のように20番地, 22番地, 24番地, ...,40番地に読み込まれる. そして3番地に戻ってくる.

O 10M, O 11Mは10MにCarriage return(復帰), 11MにLine feed(改行)のコードがあって, 出力装置を改行する. こういう機能鍵にはアルファベットが書いてないので, パンチの都合上ギリシア文字が割り当てられている.

4,5番地にT Dが2度あるが, 最初のはアキュムレータを0にするもの, 次はその0を0番地の長語に格納するものだ. そして7M番地にある10を0番地の短語に設定し, 9,10番地でP7へサブルーチンジャンプする. EDSACの短語は, EDSACのプログラム技法のブログの2012年9月22日の命令語の図にあるように, L/Sビットのところが整数の1の桁であり, アドレスはその上に1ビットずれているから, P5Fと書くと10になり, 8MのようにPD, つまりアドレス部が0でL/Sビットがあると, 整数の1になる. (ファンクション部のPは0である.)

11番地へ戻るとO 12Mでスペースを2文字出力する. 次はtを16で割りt'にするのだが, もともと読み込むデータが1/10になっているとしたので, (15は0.15と読まれたが, 1.5のつもりであった.) 10/16を掛けることになる.

EDSACの乗算は一方の数を乗算レジスタへ置き, ある番地にある数と掛けて積をアキュムレータに置くので, H 4Mで10/16を乗算レジスタに置く. 4Mを見ると J F である. F はL/Sビットが0のことだが, Jはなにかというと, 前のブログでEDSACの文字コードを見るとJは01010である. 左端の0の次に小数点があるから, これが10/16なのである.

14番地のV40Dで40番地にある長語のデータと乗算が出来た. それを8Dに置いてサブルーチンへ行き, y' を計算する. y' も8Dに置いてある.

まだまだ続くが, 今回はここまでにしよう.

2014年9月24日水曜日

割込み

私が2年ほど前にツィートした割込みの話がいまでも時々リツィートされているようだ.

今日はその割込み(英語ではinterrupt)の昔話しを書いてみたい.

私がいた研究室, つまり東大物理学科の高橋研究室では, 院生の後藤さんが発明したパラメトロンを使って計算機を自作した. 1958年のことである. Parametron Computer No.1ということで, 通称はPC-1という.

当時計算機とかプログラミングでは, 英国ケンブリッジ大学のEDSACの本があったので, 計算機を作ろうとすると, EDSACが手本であった. PC-1ではEDSACのアーキテクチャを十分検討し, ほとんどそのまま採用したものもあれば, 我々の哲学に従って変更したものもある.

入出力命令も変更した部類であった.

EDSACの入力命令 I n は, 5単位紙テープの1列を, n番地の右端に読み込む; 出力命令 O n は, n番地の左端の5ビットを書き出す. しかしこれは不便ではないかと我々は考えた. 読み込んだ文字コードを調べるには, これを演算装置, アキュムレータに取り出さなければならい. また文字コードをアキュムレータで作ってから書き出すこともあろう. そういうわけで, PC-1で は, 入出力はアキュムレータを経由することとした. さらに入力はアキュムレータの左端に入れることにした. そうするとサインビットの判定も楽だからだ.

もちろん, 入力命令コードは i, 出力は o である. ではこれらの命令の番地部は何に使うか. 計算機の処理速度に対して, 入出力は遲いから, 計算機は次々と入力や出力を続ける際, 待たされることになるだろう. そこで入出力命令を実行しようとした時, 相手の機器の準備が出来ていなければ, 別の計算に移れるように, 命令の番地部へとぶようにした. つまり条件ジャンプ命令にしたのである.

入出力をしながら, 同時に計算も出来るということで, 聖徳太子は十人の訴えを同時に聞いたと伝えられるから, この機能には聖徳太子というあだ名がついた.

しかし, 実際問題として, こういうプログラムを書くのは結構面倒である. 従って普通には命令は自分の番地にジャンプする, つまりダイナミックループをするように書かれていた.

下がPC-1が最初に実行したプログラム(neglect erase)である. 当然だが二進法で紙テープに穿孔してある.

0    i 0;
1    b 6;
2    zl 0;
3    b 6;
4    o 4;
5    jl 0;
6    111111000...0;

0番地のi命令の番地部は0, 4番地のo命令の番地部は4になっている. このプログラムは紙テープを読みそのまま出力テープに複写するが, 抹消のコード(6単位がすべて1, つまり全穿孔)は複写しない. 6番地には左端6ビットに1が入っており, 1番地のb命令でビット毎のxorをとり, 結果が0ならzl 0;で0番地に戻る. そうでなければ3番地で元に戻し, 4番地で出力, 5番地で0にジャンプバックする.

入出力と演算の同時実行のプログラムがなぜ面倒かというと, 入出力機器の準備ができていなくて演算をしていると, 今度は準備が出来たかどうかを絶えず見にいかなければならず, それはまた煩わしいのである.

例えていえば, なるべく早く食事したいが, 食堂がまだ準備中なので本を読んでいるとする. しかし食堂が開いたかどうか, 絶えず注意していると, 本を読むのも途切れ途切れで能率が悪い. 能率よく本を読みしかもすばやく食事にありつく対処法は食堂の準備が出来た時, 食堂の方から知らせてくれることである.

というわけで, 入出力機器, といっても出力のテレタイプの方だけだが, 準備ができたらプログラムに知らせる機能をつけることにした.

どうやってプログラムに知らせるか, プログラムの実行中の情報をどう保存するかを次に考えることになった.

PC-1の記憶装置は512語であったので, 510番地に次に実行する命令の番地をいれ, 511番地にジャンプすることになった. 割込む機器は1台であったが, 複数台になったとき, 割込みが続けておきると, 帰り番地の情報が失われる恐れがあるので, 割込みと同時に, その後の割込みを禁止する. この禁止はプログラムで解除する.

PC-1にはプログラムでセット, リセットできるフリップフロップが2個あった. 2014年3月13日の私のブログ「計算機による音楽演奏」で使ったのもこのフリップフロップである. その一つを割込み禁止に使った.

では割込み処理の説明に移ろう.

なにか長い計算をしていて, 途中時々数値を出力しているとする. 出力サブルーチンはその数値を十進法に変換し, そのそれぞれの数字や記号をテレタイプのコードに変換してテレタイプに送りだすのが通常の処理である.

それに対し, 割込みではコードに変換するまでは同じだが, それを1文字ずつリングバッファに入れるのである.

下の図が当時の割込み処理の流れ図である.



出力サブルーチンは, コード変換がすむと, 流れ図の左上からこの処理に入ってくる. まず割込みを禁止し, 帰り番地を510番地に設定する. 次に出力すべきテレタイプのコードをリングバッファに入れ, 入り口の指標indexを1増やす.

それから中央の線へ行き, 出口の指標exdexの指す一番古いコードをアキュムレータに取り出し, o 命令を実行する. 出力機器の準備が出来ていれば出力でき, yesの口から出る. 反対に準備がまだならnoへ出, バッファが満杯かどうかを見る. 満杯なら仕方ないから, 先程の古いコードを取り出すところへ戻り, 機器の準備を待つ. 満杯でなければ, 次のコードを入れられるから, 下へ抜けてもとのプログラムへ戻る.

1文字出力出来た場合は, 中央の線の一番上へ戻り, 出口の指標exdexを1増やす. バッファが空なら戻る. 空でなければ一応出力してみようと 試みる.

割込み処理から戻るときには, 割込み禁止を解除しておく.

一方割込みが掛かって処理ルーチンに来るときは右上から入る. この時は帰り番地の設定は済み, 割込み禁止も掛かっているが, なにかの実行中で飛んできたので, アキュムレータや他のレジスタ類を退避してから 出力のループに入る.

簡単だがこれはいちおうマルチプログラミングである. この実験をやったのは1959年の夏ころと思う. ちょうどその頃, パリでUNESCO主催のコンピュータの会議があり, Christopher Stracheyが割込み, マルチプログラミング, TSSの提案をしていた. 我々はその割込みを先取りしているので大いに気持ちがよかった.

2014年8月25日月曜日

微分解析機

このブログでフロントラッシュの話を書いたのはもう1年も前だ.

今回の再生プロジェクトでは, 精度は二の次なので, フロントラッシュを取り付けるまでには至っていないが, Crankの本でその辺を読み直してみると, 出力軸は入力軸より10パーセント高速に回転すると書いてあった. どうしてそうなるか考えてみた.

前回 ブログに掲載したCrankの本の分解図は



だが, どうもきれいとはいえない図である. そこで描きなおしたのが下だ.



写真の右端のピンの部品は省いてある. その左のドラム, 遊星歯車キャリア, 出力軸に固定されている内歯車が, その順に描いてある. 歯車の歯は描いてないが, 歯車のつもりの円板の外径は相方の歯車の外径に合うようになっている.

さて入力軸の逆転時には, ピンがドラムのペグに当るまでの時間は, ドラムは入力軸と無関係に停止している.

しかしその左のキャリアは, 入力軸に固定されて回転しているから, キャリアの入力側の遊星は, 太陽歯車で自転させられながら, 公転する.

入力側の遊星歯車の自転はそのまま出力側の遊星歯車の自転になっている.

出力側の遊星歯車の公転と自転により, 内歯車が回転し, 出力軸も回転する仕掛けである.

上の図の太陽歯車, 入力側の遊星歯車, 出力側の遊星歯車, 内歯車の半径, つまり歯数は, 最初の写真の大きさに大体合せてある.

そこで, 入力軸に対して出力軸がどう回転するか見るための図を描いてみる.



太陽歯車 S
太陽歯車半径 SA r0
遊星歯車 P
入力側遊星歯車半径 PB r1
出力側遊星歯車半径 P'C'' r2
内歯車半径 SD'' r3

入力軸が反時計回りにθ0だけ回転すると, Pにあった遊星歯車はP'へ移動する. 入力側の歯車は始めABで合っていたが, A''B''で合うようになり, 先程のBはB'へ移動している. 弧AA''は弧B'B''と同じ長さだから, 遊星歯車の回転角θ1はr0θ0/r1である.

内歯車はキャリアに従って回転する他, 出力側遊星歯車の回転でさらにθ2だけ回転する. その角は青で示した弧の関係からr2θ1/r3になる.

従って内歯車=出力軸は, キャリア=入力軸がθ0回転するのに対してθ02回転する. θ2がθ0の10パーセントなら, Crankのいう通りだ.

上の図ではr0:r1:r2:r3が3:4:2:9になっているから, 1/6だけ速く回転していることになる. 実際のフロントラッシュの比率は, 理科大へ行って歯車の歯数を数えなければならないが, 近代科学資料館はアナコンの企画展の後, 夏休みになっているので, 再開したら見に行こう.


以下9月5日記

資料館が再開されたので, フロントラッシュの歯数を数えてきた.

太陽歯車 歯数 22
入力側遊星歯車 歯数 28
出力側遊星歯車 歯数 12
内歯車 歯数 62

従って
θ20r0r2/r1r3=0.15θ0

つまり 10パーセントではなく, 15パーセントの歯数比であった. (1/6とあまり違わないね.)

2014年8月2日土曜日

曜日の計算

数学セミナー1978年7月号に島内剛一先生が寄稿した「万年七曜表」に登場する計算尺については, すでに2回のブログで説明した.

今回は円形計算尺によるものを話題としたい. 私はこれが一番好きだ.

前回の円筒式のものを見れば, そのまま円形にもなりそうに思うが, 島内方式の円形計算尺は, また趣向が違って, 次のような形をしている. 島内流にいうと上からA,B,Cである. (私の流儀で書き直しているから数セミの図とは多少違う.)







これらを重ねるのでA,B,Cをそれぞれ黒, 青, 赤の各色で示す. A, B, Cは紙に書き込まれ, 共通の中心の回りに回転出来る. Aには時計でいえば1時方向, 3時方向, ..., 11時方向には, Cにある月名を見る穴が開いている. (島内式では5時と7時方向にしかない. 紙方式ではAとCの間にBがあるので, Bにも同じ半径で5時方向から7時方向までの円弧状の穴がある.)

Aの黒文字は西暦の上2桁で, 外周の0から15がユリウス暦, 内周の15から22がグレゴリオ暦だ. Bの青文字は外の一周が日付, 内側の3周は西暦の下2桁. Cの赤文字は, 外が週日, 内側が月名である. 週日のRは木曜のこと.

前回と同様, まず2014年7月について, これを使ってみるには, まず11時方向の黒の上2桁(20)と7時方向の青の下2桁(14)を合わせる. 下の図は黒はそのままで青を右方向に4時間分ほど回転している.



黒と青をそのように合せてから, 外側の赤の紙を, 週日名と日付が合わさるように回転し, 黒の月名用のどれかの穴から目的の月名が見えるようにする. この図では9時方向の穴に赤字の7が見える. 週日名は円周に沿って6回繰り返すので, 月名は6個の穴のどれで見えてもかまわない. このように回転すると7月1日が火曜なのが分る.

ところで, その穴で7が端の方のあるのは, ここで丸めているからで, 赤字の紙は青字の日付と赤字の曜日が重さなるように, 離散的に動かすから, このようになる.

5時方向の穴に4が見えるのは, 4月と7月は曜日の関係が同じだからだ.

再び前回同様に2000年1月について試みると



のよう黒の20と青の00が重さなり, 11時方向の枠に1が入り, 1日の相方に土曜の方のSが來るので, 1月1日は土曜であった.

例によってこの仕掛の説明である. 2014年7月の計算法について, mod 7の図を示すと次のようだ.



2014の上2桁は20なので, 黒の目盛のf(20)=5.875にマークする. 下2桁14ではg(14)=17なので, mod7をとり, 青の3にマークし, これらを合せる.

青は1.875分右に移動する. 今は赤も青と一緒に移動する.

このとき赤のb(7)=13.25は黒の目盛上では赤点の1.125に対応し, これを丸めると1になる. 従って7月0日は月曜だ.

5.875+17+13.25=36.125=1.125 mod 7

2000年1月で試みると



5.875+(-0.25)+6.75=12.375=5.375 これが黒目盛上の赤点の位置で, これを丸めて1月0日は金曜である.

調べようとする年月の0日の曜日までは分ったが, これをカレンダーに対応させるにはまたひと工夫がいる.

とりあえず青に下の図のように, 1,2,...,7を書き入れる. これが日付だ. また赤にも図のように曜日を書き込む. そして赤目盛を動かす.

また黒目盛の0(mod 7だから7も)の両側に0.5の幅をとり(つまりこの範囲を丸めると0になる), 赤目盛上の月の値の黒マークがその範囲に入るように, 青と目盛の位置を合わせながら移動するのである.

2014年7月で使ったふたつ上の図の赤目盛は一目盛分左にずれて, マークが7付近に近づき, 青の1日が赤の火曜に対応する.



2000年1月での図もこのようになる.



要するに日付と曜日はなんとなく書き入れておき, 赤目盛の黒点を黒目盛のどの値の範囲に置くと曜日と日付が対応するか試み, その範囲が0付近になるように日付を再調整したまでである.

なかなか楽しい島内方式の検討であった.

私は例によってPostScriptを活用したが, 島内さんは「円周を42等分するには, まずコンパスで6等分し, それをさらにコンパスまたはディバイダで7等分するのがよいだろう.」と書いているから, そのようにして作図していたのかもしれない.

いまやこの計算尺を作るにも文明の利器が出現してきているが, 一方計算尺がなくても曜日が簡単に分かる手段も増えていて, 計算尺の出番もない.

2014年7月23日水曜日

曜日の計算

前回のブログで島内先生の計算尺方式の万年七曜表を話題とした.

余談ながら現在東京理科大学の近代科学資料館で開催中の企画展「計算する器械たち---アナログコンピュータ展」には 円筒形の計算尺も展示されている. 今回はそういう円筒方式の万年七曜表の話だ.

下の図を見てほしい. ちょっと分り難いがこれが円筒式のもの(展開図ともいうべきもの)だ.



前回のブログを見た人なら思い出すだろうが, 最上段の右から0,1,...,15 少しあけて15,16,...,20とあるのは西暦の上2桁である. 右側の15まではユリウス暦, 左側の15からはグレゴリオ暦用のものだ.

その下, 5段位に左下がりに並ぶのは西暦の下2桁である. 段々と下るのは円筒の卷いたとき, 21に22を自然に繋げるためである. なかなか凝っている(次の図を見ると分る). 44は丁度切れ目になってしまった.

次は2, 8, 5,...と現れて, 4と7が重さなっているから, これは月名である. 展開図の範囲に4回出てくる.

最後は, 今度は右下がりだが, これは一箇月のカレンダーである.

もうひとつ青色の枠の図がある. 週日名が書いてあるほか. 枠の所が3カ所くり抜いてある.

島内さんの記事では, 西暦上2桁と月名が1枚の紙でCといい, 下2桁とカレンダーがもう1枚の紙でBといい, 青色の枠の紙はAである. 同じ幅に描いてあるが, A,B,Cの順の幅がすこしずつ短くなり, Aが一番外側, Bが真中. Cが内側でそれぞれ円筒状になって外中内の順で差し込まれている.

勿論底の場所は同じで, 展開図では左右に移動するように, 円筒は回転できる.

前回の例と同様に2014年7月の曜日を見るには,



のようにCの20の真下にBの14を合せる. そしてCの7月がAの月名を入れる青色の小さい枠に入るように合せると, 図のように青色の大きい枠がカレンダー部分に重さなり, 7月1日が火曜と判明するのである.

この図では, Bの紙が円筒になっているのが見えるであろう.

難しいのは例の丸めのところである. この図でも青色枠内の4と7が枠の左右の中央にないのが見てとれる. 4と7の中央は1.125という位置にあり, 枠の中央が1なのである. どこか中央かは, 下のカレンダー部分の枠の左右がカレンダーの数字の丁度真中に來るように離散的に合せることで分るのだ.

もうひとつの例



この例は2000年1月で, 閏年であり, 1月, 4月, 7月が同じ曜日の配置になる(Aの同じ枠に入る)ことが見てとれる.

さてこの計算尺の仕掛けの説明だが, mod 7の範囲だけの図にしたのが次だ.



最上段の線が西暦の上, その下の線が西暦の下で, 2014年になるようにその20と14を対応させた図になっている. つまりBの紙が右に移動している. (輪だから左かもしれない.) 移動していることを示すため赤で描いてある.

その下の線が月名で, これは最上段と同じ紙, 同じ位置にある.

さて赤の線の基準に対して7月の位置を計算すると, g(14)とf(20)が合っているから, 黒の基準位置はg(14)-f(20). 従って7月の位置はg(14)-f(20)+b(7). ところがfのスケールは反対に取ってあるから, g(14)+f(20)+b(7)になって, 実際の値を入れると
17+5.875+6.25=29.125=1.125 mod 7
従ってこれを丸め, 7月0日の値が1. すなわち7月1日は2(火曜)となる.

7月のすぐ右に9月がある. だから9月0日の値は0, 9月1日は1(月曜)だ.

従って青色の枠は下のように作ればよい.



上の小さい枠に7月があると, その1日は火曜. 上の枠がひとつ右に移ると, 曜日の名前が右にひとつずれて, 9月1日が月曜になるという仕掛である.

2014年7月14日月曜日

曜日の計算

私が2011年5月25日や6月5日のブログに書いた故島内剛一(しまうちたかかず)先生の曜日の計算は興味深いものであった. 先生の数学セミナーの記事(1978年7月号)には驚くようなアルゴリズムが沢山登場する.

その中で, しばらく私が読み飛ばしてはいたものの, 気になっていたのは, 「計算尺方式の万年七曜表」であった.

先日来, ちょっと読み直してみたが, もとの図を島内さんがどのようにして書いたか不明である. 雑誌に掲載された図の正確さもよく分らないので, 使い方が分かった時に, 自分で図を書き直してみた. (元の記事には「日のひと目盛の1/8=0.125までの精度はどうしても必要」と書いてある.)

ざっといえば, 下のようなものである.



上下2枚の図があり, どちらも半分が上下逆になっているのは, 中央で折って表と裏にするのである. 島内さんは上をB, 下をCといった.

Bの上半分は横に2倍にしたカレンダーである. 左側のどの列から始めても一週間分があるようになっている. 下半分は西暦の下2桁の表だ. こちらもかなり冗長である.

Cの下半分は月名である. 左の方, 3と11が縦に書いてあるのは, 3月と11月の曜日が同じことをしめす.

上半分は西暦の上2桁になっている. 0から15までの部分と15から22までの部分があるが, 前者はユリウス暦, 後者はグレゴリオ暦に対応する.

私は島内さんの記事のこれらの図をコピーし, 切り抜き, BとCを計算尺のように滑らせながら操作してみて, 使い方が分ってきた.

逆さまの図は見にくいから, 正立させた図にしよう. 上下の余白も省いた.



この上の図では上と下が表面と裏面に相当する. それぞれの面の途中の横線が計算尺の滑る位置だ. 上のカレンダーの中央付近に緑色の枠が見えるのは, 通常の計算尺にあるカーソルのようなもので, その内部がその年その月のカレンダーになる.

前の説明の繰り返しだが, 上から月, 日, 年の下2桁, 年の上2桁である.

さて, 2014年7月のカレンダーを見るには下の図のようにする.



下の西暦上2桁の20と, 上の西暦下2桁の14が合うように下の目盛を右にずらす. すると上の月名の目盛は(上下反転したので)同じ長さだけ左にずれる.

日の表の上端に曜日の区間の境界を示す青い線があり, 7の下の対応する区間の中央にM(monthのつもり)がくるように緑の枠を滑らすと, 7月1日が火曜になり, 7月のカレンダーが得られる仕掛けである.

目盛合わせが正確にできるよう, この図には多くの縦線が引いてある. ユリウス暦とグレゴリオ暦の間には微妙な差があるので, グレゴリオ暦の部分は青線で表示する.

この計算尺でカレンダーが出来る魔法のような理由は次のとおり.

式が導出された経過は省略するが, 島内方式で西暦y年m月d日の曜日は

round(f(y idiv 100) + g (y % 100) + b (m) + d) % 7

で計算する.

y idiv 100は西暦の上2桁. y % 100は下2桁である.

f(y)= - floor(3(y+1)/4) - y + 6 + (y % 4 == 0 ? -0.125 : +0.125)
(ユリウス暦に対するfの式は省略)
g(y)= y==0? -0.25 : (y + floor (a/4) - y % 4 ==0 ? -0.5 : 0)
bは1月から12月について
6.75, 2.75, 3.25, 6.25, 1.25, 4.25, 6.25, 2.25, 5.25, 0.25, 3.25, 5.25
である.

y=15から22までのfの値は
0.125, 5.875, 4.125, 2.125, 0.125, 5.875, 4.125, 2.125

y=0から15までのgの値は
-0.25, 1, 2, 3, 4.5, 6, 7, 8, 9.5, 11, 12, 13, 14.5, 16, 17, 18

アルゴリズムを解明すべく, これを次のような目盛に記入する. gは7の法をとってある. 目盛を下からf, g, m, dということにする. dの線の上の1〜7はカレンダーの一番上の数である.



いま2014年の計算をしようとして, 下の図のようにfの目盛の20とgの目盛の14を合わせると(赤い線), gの0はfの1.875に合う.

これはf(20)=5.875, g(14)=3の和 5.875+3=8.875=1.875(mod 7) である.



通常の計算尺でA尺とB尺を使い, log x + log y を計算するには, A尺のlog xのところにB尺の1(原点)を合わせ, B尺のlog yに対応するA尺を見るとlog x + log yになっている.

しかし今の計算では, gの目盛が逆向きにとってあるので, 通常の計算尺の減算の要領で加算が出来るようになっている. 結局fが右に1.875ずれているわけで, 同時に一番上のmも左1.875(=f+g)ずれている. mでの7月のb=6.25に対応するdの位置(緑の線)は, f+g+bの値f(20)=5.875, g(14)=3, b(7)=6.25の7を法とした和1.125)になっており, 丸めると1になる.

7月はこの値が1. ということは1日はd=1を足すから2になって火曜になる. 従ってカレンダーの中央の水曜を2日の上に置くことになる. dの値の1の上に2とある理由だ. 丸めた値が2だと1日は3だから水曜で, カレンダーの中央を1日に置く. dの目盛のカレンダーの値はそのように決めてある.

この説明も回りくどいが, 実は私も理解に到達するのに結構苦労した. 島内さんの発想に驚くばかりである.

2014年6月30日月曜日

微分解析機

東京理科大学近代科学資料館で企画展「計算する器械たち」を開催中だ.

そこで毎日微分解析機を実演している. デモは例のサークルテストだが, やはり一周すると2パーセントくらいのエラーが出る.

このデモの様子を, 昔微分解析機を使っていらした, 渡辺勝先生にお見せしたら, 「サークルテストでの反転時のバックラッシュが小さいようにお見受けしております」というメイルを頂いた.

ところで手回し計算機や電動計算機が出回ると, Runge Kuttaという計算法で連立m元常微分方程式を解くようになった. 我々がプログラミングを勉強したケンブリッジ大学のEDSACの本の例題にRunge Kutta Gillがあったので, 当時のコンピュータ屋にはお馴染の方法である. この方法については「伊理正夫,松谷泰行,Runge-Kutta-Gill法について,情報処理,Vol.8,No.2,pp.103-107」に解説がある. 今回はそれを使ってサークルテストを計算してみた. Schemeのプログラムは次のようだ.
(define (runge-kutta-gill m h n)
 ;m 変数の個数, h 分点間隔, n 分点数
(let ((y (make-vector m)) (f (make-vector m))
  (q (make-vector m)) (k 0) (qi 0) (s 0) (r 0)
  (c0 (- 1 (/ 1 (sqrt 2)))) (c1 (+ 1 (/ 1 (sqrt 2)))))
 (define (inity)    ;yの初期値
   (vector-set! y 0 0)
   (vector-set! y 1 0)
   (vector-set! y 2 1))
 (define (calcf)    ;fを計算
   (vector-set! f 0 1)
   (vector-set! f 1 (vector-ref y 2))
   (vector-set! f 2 (- (vector-ref y 1))))
  (define (printy)  ;yの出力
   (display (list 
     (vector-ref y 0)
     (vector-ref y 1)
     (vector-ref y 2))) (newline))
  (newline) (inity)
  (do ((i 0 (+ i 1))) ((= i m)) (vector-set! q i 0))
   ;qiの初期化
  (calcf) (printy)
  (do ((j 0 (+ j 1))) ((= j n))
  (for-each (lambda (e c) 
   (do ((i 0 (+ i 1))) ((= i m))
    (set! k (* h (vector-ref f i)))
    (set! qi (vector-ref q i))
    (set! r (e))
    (set! s (vector-ref y i))
    (vector-set! y i (+ s r))
    (set! r (- (vector-ref y i) s))
    (vector-set! q i (+ qi (* 3.0 r) (* c k))))
   (calcf))
   (list (lambda () (- (* 0.5 k) qi))
     (lambda () (* c0 (- k qi)))
     (lambda () (* c1 (- k qi)))
     (lambda () (/ (- k (* 2.0 qi)) 6.0)))
   (list -0.5 (- c0) (- c1) -0.5))
  (printy))))

(runge-kutta-gill 3 (/ (atan 1) 2.5) 20)
普通とちょっと違うのはy[0]がxであることだ. 従ってf(というのはdy/dxのことだが)のf[0]=1にしてある.

yの初期値では, y[0](xのこと)=0, y[1](yのこと)=0, y[2](dy/dxのこと)=1. fの計算では f[0]=1, f[1]=y[2], f[2]=-y[1]とする.

Runge Kuttaでは分点を1個進めるのに, x, x+h/2, x+h/2, x+hについて4回fの値を計算するが, そのループはfor-eachで回している. ループごとに変る値については, (lambda (e c) ..)の引数で渡す. 特にeの値は, 呼ばれた時に計算する必要があるので, 引数の方は(lambda () ..)の形にし, 使う時に呼出す.

initiy, calcf, printfなどのサブルーチンが, 関数runge-kutta-gillの中で定義してあるのは, 問題ありだが, yやf読み書きするので, これで我慢している.

最後の行で呼び出しているが, 0から2πまでを20分割で計算する. つまりh=2π/20; πを4*(atan 1)で計算したいからh=(atan 1)/2.5としてある.

たったの20分割なのだが, 一周したときのy[2]が.999868...だから1万分の2の程度の精度だ.

その出力を貰ってサークルを描くPostScriptのプログラムは次の通り.

/ps[[0 0 1]
[.3141592653589793 .30899155257892935 .9510578492071948]
[.6283185307179586 .5877376828378169 .8090352529734782]
[.9424777960769379 .8089575954451165 .5878333484965556]
[1.2566370614359172 .9510010098334781 .30910245672629366]
[1.5707963267948966 .9999670230159169 1.2303914619280843e-4]
[1.8849555921538759 .9510645042444495 -.30886434562367165]
[2.199114857512855 .8090808881734999 -.5876187580148391]
[2.5132741228718345 .5879134969774285 -.8088585919500648]
[2.827433388230814 .3092092738117766 -.9509316169859521]
[3.141592653589793 2.4607017746574233e-4 -.9999340319806838]
[3.4557519189487724 -.30873714204448566 -.951071143410806]
[3.7699111843077517 -.5874998314987273 -.8091265072362017]
[4.084070449666731 -.8087595818584492 -.5879936306339407]
[4.39822971502571 -.9508622132841098 -.30931607883671747]
[4.71238898038469 -.9999010268957623 -3.6909309235748844e-4]
[5.026548245743669 -.951077766707203 .30860994184321267]
[5.340707511102648 -.8091721101619074 .5873809032915227]
[5.654866776461628 -.5880737494657694 .8086605651723112]
[5.969026041820607 -.3094228718001786 .9507927987297933]
[6.283185307179586 -4.921078894070952e-4 .9998680077626148]
] def
280 360 translate
/mv {200 mul exch 200 exch mul moveto} def
/ln {200 mul exch 200 exch mul lineto} def
ps 0 get dup 1 get exch 2 get mv
1 1 20 {ps exch get dup 1 get exch 2 get ln} for stroke
図はこのようだ.



2014年6月8日日曜日

鶴亀鴉算

小学生のころ鶴亀算を習った. 例えば鶴亀合せて10匹. 足は合せて30本. 鶴と亀はそれぞれ何匹か. 全部鶴とすると足は20本. 残りの10本は鶴より2本足が多い亀によるから, 亀は5匹. 従って鶴も5匹. 足の本数を検算すると2*5+4*5=30.

矢野健太郎先生がアメリカの子供にこのような話をしたら, 鶴と亀とは見ただけで分り, 足だってまったく違うという反応だったらしい. そこで矢野先生は2ドルのケーキと4ドルのケーキを合せて10個買い, 合計を30ドルにするには2ドルと4ドルのケーキを何個ずつ買うかと話したそうだ.

小学6年の私は近くに住む中学1年の先輩と話をしていた. 先輩は「鶴をx, 亀をyとすると, x+y=10, 2x+4y=30. これを解くとx=5, y=5が得られる」と得意そうに話した. へぇーすごいなぁと思ったが, 1年後には私もその先輩と同じ中学(詳しくいえば7年制高等学校尋常科)に入学し, 同じように解けるようになった.

さて東京理科大学近代科学資料館で, 6月19日から8月8日まで「計算する器械たち —アナログコンピュータ展—」が開催される.

それに門脇廉君(現 九州大学)が作った機械式アナログ計算機の「三元連立方程式求解機」も展示される. 私は展示準備中のその機械を見たが, 非常に美しく作られていて感動した.

企画展開催中に説明を担当する理科大の学生さんたちと三元連立方程式の例題を作って解いてみたりしていたが, 見学に来た小学生に連立方程式の意味を教えるにはどうするかということになった.

そこで私が考えたのが, このブログの始めにあった鶴亀算であった. でもそれは二元であるが, もう1種類の変わった動物を登場させようと思った.

それがサッカーの3本足の鴉である. ボールを掴んでいるのが3本目の足だ.



これが三元連立方程式求解機用問題.

t, 亀k, 鴉cとする. ただし鴉はサッカーJFAのシンボル 3本足のものである. 頭の数が6, 羽の数が10, 足の数が16のとき, 鶴, 亀, 鴉それぞれ何匹か.

つまり連立方程式
t+k+c=6 (0)
2t+0k+2c=10 (1)
2t+4k+3c=16 (2)
を解くわけだ.

早速三元連立方程式求解機のシミュレータでやってみる.



4枚のプレートがあり, 左からt, k, c, 定数に対応している. また制約するテープは内側から順に上の式(0), (1), (2)に対応している.

プレートの傾斜は右端の定数のを10°にしたので, 左からの傾斜は31.395°, 10.0°, 20.322°,10°である. つまり(t, k, c)=(3, 1, 2).

この図で分かるのは定数項の値が係数に較べて大きいので, 右端の定数のプーリーが分散しているのに対し, 係数の方のプーリーがごちゃごちゃと固まっていることである.

これを解決するには上の連立方程式の変数のそれぞれから1, または2を引くのである.

1を引くと
t+k+c=3 (0)
2t+0k+2c=6 (1)
2t+4k+3c=7 (2)

2を引くと
t+k+c=0 (0)
2t+0k+2c=2 (1)
2t+4k+3c=-2 (2)

2を引たので解いたのが下の図である. 各プーリーが分散している.


傾斜が左から10°, -10°, 0°, 10°になり, (t, k, c)=(1, -1, 0)であるから, 鶴が3, 亀が1, 鴉が2と判明した.

却って難しくなったかしら.

2014年5月31日土曜日

横方向の和の和

昨日のブログではさぼったnusum (2r-1) から nusum 2rを計算してみた.

前回のような表を書くと

nnusumνnes
2r-1nusum 2r-1r(r-1 r-2 r-3 ... 0)
2rnusum 2r1(r)


羃の列が(r-1 r-2 r-3 ... 0)のnusum (2r-1) の式は

nusum (2r-1) = (r-1)2r-2 + (r-2+2)2r-3 + (r-3+4)2r-4 + ... + (0+2r-2)2-1

で, これにrを足したものが

nusum 2r = r.2r-1

になればよい.

nusum (2r-1)
={各項の係数の中を計算する}
(r-1)2r-2+r2r-3+(r+1)2r-4+...+(r+(r-2))2-1
={降順の羃の列で同じ係数を持つものをならべる}
(r-1)2r-2+(r-1)2r-3+(r-1)2r-4+...+(r-1)2-1  {1行目}
             +2r-3+2r-4+...+2-1                    {2行目}
                    +2r-4+...+2-1                    {3行目}
                          +2r-5+...+2-1              {4行目}
                          +...
                          +2-1                          {r行目}
={各行を2p+2p-1+... +2q=2p+1-2q で置き換える}
(r-1)(2r-1-2-1)
+2r-2-2-1
+2r-3-2-1
+...
+20-2-1
={1行目を展開し2行目以下の第1項を足し合せる 2行目以下のr-1個の2-1を足す}
(r-1)2r-1-(r-1)2-1+(2r-1-20) -(r-1)2-1
={いろいろ消えて}
r2r-1-2r-1-(r-1)+2r-1-1
=r2r-1-r.

rを足せばnusum 2rになる.

存外簡単であった.

2014年5月30日金曜日

横方向の和の和

TAOCPの演習問題713-42にこういうのがある.

「問題42. [M21] e1 > ... > er ≥ 0について, n = 2e1 + ... + 2er なら, 和 ∑k=0n-1νkを, 冪e1, ..., erを使って表せ.」

ここでνkkを二進法で表わした時の1の数であり, 横方向の和といったりする.

ν0=0; ν1=1; ν2=1; ν3=2; ν4=1; ν5=2; ν6=2; ν7=3; ν8=1; ν9=2; ...

従って問題はnが与えれた時, n-1までの各整数にあった1の総和nusumを計算するものである.

n=10とすれば, 上のνの値を使いnusum 10 = 0+1+1+2+1+2+2+3+1+2=15

しかし途中のνを全部計算して求めるのはしんどいから, 10を二進法で1010と表わし, 10 = 8 + 2 = 23 + 21. そのe1=3, e2=1からnusumを得たいという問題である.

私はこの解答にある J.-P. Allouche, J. Shallitの本 Automatic Sequences (2003)にあるというフラクタルの図が見たかったが, その本が簡単に見られそうもなかったので, PostScriptで描いてみることにした. そこで解答は次のようだ.

「解答. rについての帰納法により,

e12e1-1+ (e2+2)2e2-1 + ... + (er+2r-2)2er-1.

D. E. Knuth, Proc. IFIP Congress (1971), 1, 19--27. この和のフラクタル模様はAlloucheとShallitの本の図3.1と3.2に示してある.] またSn'(1)も考察せよ. ただし Sn(z) = ∑k=0n-1zνk = (1+z)e1+ z(1+z)e2+... +zr-1(1+z)er.」

つまりe1=3, e2=1だったから 3*23-1 + (1+2)*21-1 = 3*22 + 3*20 = 12 + 3 = 15なわけだ.

もう少し計算してみる. 下の表は左端がn. 次がnusum. 各行のnusumは 1行上のnusumに1行上のνnを足したものになっている. その次のかっこ 内はe1...er. さらにその右は上の 式による計算とその和を示す.

nnusumνnes
101(0)0.2-10
211(1)1.201
322(1 0)1.20 + 2.2-12
441(2)2.214
552(2 0)2.21 + 2.2-15
672(2 1)2.21 + 3.207
793(2 1 0)2.21 + 3.20 + 4.2-19
8121(3)3.2212
9132(3 0)3.22 + 2.2-113
1015(3 1)3.22 + 3.2015

さてこれでrによる帰納法が出来るであろうか.

偶数のnからそれより1大きいn'に進むとき, eの列は 最後にr+1項目として0が追加される. だから

nusum n+1 = e12e1-1+ (e2+2)2e2-1 + ... + (er+2r-2)2er-1 +(er+1+2(r+1)-2)2er +1-1 .

er+1=0とすると

nusum n+1 = e12e1-1+ (e2+2)2e2-1 + ... + (er+2r-2)2er-1 +2r.2-1
=nusum n+r = nusum nn = nusum n'.

一方n=7から8へのように繰上がりがあるときはどうするか.

31から32の計算をしてみると, eは(4 3 2 1 0) から(5)になる.

このときのnusumは

4.23+5.22+6.21+7.20+8.2-1

これを次のように計算してみる.

3(23+22+21+20+2-1)
+(23+22+21+20+2-1)
+(22+21+20+2-1)
+(21+20+2-1)
+(20+2-1)
+(2-1)
=
3(24-2-1)
+(24-2-1)
+(23-2-1)
+(22-2-1)
+(21-2-1)
+(20-2-1)
=
3(24-2-1)
+(25-20)
-5(2-1)
=
3*15.5+31-2.5
=75
nusum 32は(5)だから, 5.24 = 80. さっきの75にν31=5を足して 80になる.

一般的には上の(4 3 2 1 0)が(r-1 r-2 ... 0)なのでr について同様の計算をしなければならない. しかしもう出来そうになったから フラクタル図の方へ急ごう.

nusumの式をPostScriptで書くことになる.

/es {4 dict begin /s exch def /m exch def /n exch def
/l s length def
n 0 eq{s}{n 2 mod 0 eq{n 2 idiv m 1 add s es}
{n 2 idiv m 1 add m s aload pop l 1 add array astore es}
ifelse}ifelse end} def

/nusum {4 dict begin /n exch def
/ee n 0 [] es def /sum 0 def
0 1 ee length 1 sub{/i exch def /e ee i get def
/sum e i 2 mul add 2 e 1 sub exp mul sum add def} for
sum end} def

0.33 setlinewidth
0 1 560{/x exch def
x 0 moveto 0 x nusum 4 div rlineto stroke} for
最初の手続きesはnの冪の列を返す. n 0 [] esのように呼ぶ. 10 0 [] es => [3 1]

手続きnusumが総和の計算である. 32 nusum => 80.0, 560 nusum => 2480.0

n=0から560までを描いたのが下の図だ. フラクタルと言われればそう見えなくもないが.



横軸の目盛は100ごと, 縦軸の目盛は500ごとである.

2014年4月14日月曜日

ビットごとの秘法と技法 から

TAOCPの演習問題7.1.3-20にGosperのハックというのがある.

正の整数xについて

ux&-x;
vx + u;
yv + (((xv)) »2);
としてyを返す.

とにかくやってみよう.
(define (gosper x)
 (let* ((u (band x (- (expt 2 16) x)))
        (v (+ x u))
        (y (+ v (>> (/ (bxor v x) u) 2))))
 y))

(gosper 1) => 2
(gosper 2) => 4
(gosper 4) => 8
(gosper 8) => 16
(gosper 3) => 5
(gosper 5) => 6
(gosper 6) => 9
(gosper 7) => 11
(gosper 11) => 13
という次第で, 1のビットの数が同じの次に大きい整数が得られる.

その理由を考えてみた.

x=α01a0b

とする. つまりxは左に適当な0と1の列αがあり, その右に0が1個, その右に1がa >0, その右に0がb ≥0個あるとする.

uの式は有名な最も右の1を取り出すものだから

u=10b

である.

v=x + uは1の並びの右端に1を足すから

v=α10a0b.

xv は左のαがキャンセルされるから

xv=11a0b

になり, これをuで割るから右の0の列が消えて

(xv)/u=11a

つまり1が右端にa +1個ならぶ. これを2ビット右シフトするから右端に1がa -1個ならぶ. これをvに足すから

y = α10b01a-1     (a >0に注意)

で次のyが得られたのであった. たしかにハックだね. これはMITのHAKMEM175番にある.

2014年4月6日日曜日

ビットごとの秘法と技法 から

最も左のビットの位置を知る法

TACPにはもっと凄いλ関数の計算法がある. 今回はそれを説明しよう.

xのサイズは知らなければならないが, それがどんなに大きくてもbignumのような2adic integer(2個進数)での処理だ.

まず準備としてあるアルゴリズムを考える. 具体的に4ビットとしよう. 4ビットだから対象は0≤x<16 である. c≤8と最も左のビットだけが1のh=8が登場する.

h|x-cは引かれる方が8≤ <16だから, 0≤x<16について, 差は8-c≤ <16-cの2回繰り返しになる.

たとえばc=3とすれば, 5,6,...,12,5,6,...,12になる. つまりそれぞれの最初のc個だけが<8になる. これにxをビットごとorすれば, 右半分はxの最も左のビットが1なので≥8になる. 従って全体では左端のc個だけが<8, それ以外は≥8になって, hでビットごとandをとると左端のc個が0, それ以外は8になる. c=3の例では0,0,0,8,8,...,8だ. このアルゴリズムで重要なのは, 4ビットの範囲の左から繰下げがないことである.

以下のアルゴリズムにはこれを多用する.

xが如何に大きくても, gを2の羃, ng2とし, x<2nの⌊lg (x)⌋を計算する.



以下の説明ではg=4, n=16とする. そうすると最初の行にあるlは216-1, つまり1111111111111111を24-1, つまり1111で割るから0001000100010001になる. 従ってhはそれを3ビット左へシフトした1000100010001000になる. 最初の行は16ビットを4ビットの区間4つのそれぞれについて, 4ビットがオール0かそうでないかを上のアルゴリズムで調べる. c=1の場合だ. 各区間はその区間がオール0か否かにより, 0か8になる. その結果がt1だ. だからt1は各々の区間の結果をp,q,r,sのビットで表すと, p000q000r000s000になる. p,q,r,sの少なくても1つは1である. 次の行はまずaを計算する. 212-1つまり111111111111を23-1つまり111で割るからa=001001001001である. これを先ほどのt1に掛けるから, 積は
         p000q000r000s000
      p000q000r000s000
   p000q000r000s000
p000q000r000s000
の和,
p00pq0pqrpqrsqrs0rs00s000
222221111111111
4321098765432109876543210←ビット位置

になり, この下から12,13,14,15ビットをとるからyはpqrsにlを掛けたもの, pqrspqrspqrspqrsになる. 言い換えれば区間の情報を4ビットに詰め込んで4(=g)回コピーしたものが出来た.

今度はt2の計算だ. bは220-1を25-1で割るからb=1000010000100001となる. だから式中のbは左端の区間から順に1000, 0100, 0010, 0001で, これらをcとしてまたあのアルゴリズムを使う.

pqrs≥1000>100>10>1ならすべての区間で8になり,
1000>pqrs≥100>10>1なら0888になり,
1000>100>pqrs≥10>1なら0088になり,
1000>100>10>pqrs≥1なら0008になる.
従ってt2は二進法では
1000100010001000
0000100010001000
0000000010001000
0000000000001000
のいずれかになり, このそれぞれからm
1111111111111111
0000111111111111
0000000011111111
0000000000001111
になり, 4ビット右シフトしたものとxorするとmとして
1111000000000000
0000111100000000
0000000011110000
0000000000001111
が出来る.

次の行は最も左の1を持つ区間をmでマスクしてとりだし, lを掛けて左へ4個ならべ16ビットに区切ってから12ビット右シフトして右端の4ビットに置き, さらにlを掛けて4区間全体に展開する. 先ほどのpqrsから最も左の区間を選んだように, 今度は最も左のビットの位置がt3に出来る.

整理すると最も左の1のある区間が右端にあればt2は0008, その1つ左にあれば0088, 左端の区間なら8888である.

その区間で最も左の1が右端にあればt3は0008, 左端にあれば8888である.

これからビットの位置λを計算するにはこうする.

t2を5ビット右シフトするから, 0,4,44,444のいずれかになる.
t3を7ビット右シフトするから, 0,1,11,111のいずれかになる.

これらを上下から取って足し, lを掛けて12,13,14,15ビット目をとるとそれがλである. たとえば444と1だとすると和は445になり,
   445
  445
 445
445
この4+4+5で13が得られる. 下からの繰り上げはないか心配だが, この計算は十六進であり, 最大の444と111でも5+5+5になって繰り上げはない.

このアルゴリズムは繰り返しがないのでO(1)と書いてあるが, 乗算が5回もあるから実用にはならないとTAOCPはいう.

次はDoループはあるがステップ間のジャンプはないアルゴリズムである.



d=4とするとn=d・2d=64になる. この値を使ってアルゴリズムを調べてみよう.

B1. λ=0とした後, ⌈lg n⌉=6だから, k=5,4についてxを32ビット, 16ビットの幅に1があるかどうかでxを右シフトし, λを増やす. このステップが終わるとxは最も左のビットを含む16ビットの区間になり, その区間がどこにあるかに従ってλは0, 16, 32, 48のいずれかになっている.

B2. 16ビットの区間に縮められたxを4つの区間にコピーする.

次にμd,kというのがでてくるが, 右から2kごとに1と0を繰り返す長さ2dのビット列である.

μd,k=(22d-1)/(22k+1)

今はd=4だから
μ4,3=0000000011111111
μ4,2=0000111100001111
μ4,1=0011001100110011
μ4,0=0101010101010101

μ4,3の0の部分を左半分, 1の部分を右半分といおう. 同様に
μ4,2の0の部分を左四半分, 1の部分を右四半分といい,
μ4,1の0の部分を左八半分, 1の部分を右八半分といい,
μ4,0の0の部分を左十六半分, 1の部分を右十六半分といおう.

すると下の図で灰色の部分が左なんとか半分, 白い部分が右なんとか半分である.


B3. これらのμを繋げてビットごとにnotをとったものとxのandをyとする. つまりxの左なんとか半分の1をとりだす.

B4. xyは右なんとか半分の1になるわけだ. そしてh=8000800080008000としておなじみのアルゴリズムを使う.

そうすると左部分をアルゴリズムのx, 右部分をcとすることになるので, x<cの時に0, xcの時に8だったように, 左≥右なら0, そうでなければ8000になる.

上の図で最上段の0〜fはビット番号で, そのビット6に最も左の1があるとすると上から順に
左<右
左>右
左>右
左<右
だから, B4のtはそれぞれの区画で0, 8000, 8000, 0になる.

B5. ⌈lgd⌉=2なので, 0,1のkについて
24-20=15ビット,
25-21=30ビット
ずつ左シフトして加える. 8の立つビットを左からp,q,r,sとするとそれらのビット位置はそれぞれp=63, q=47, r=31, s=15である.

まず15ビット左へ移るから, qは62, sは30になりpq, rsは並ぶ. それをさらに30ビット左シフトすると31,30にあったrsは61,60になり, pqrsは並ぶ.

264の剰余をとるから64ビットレジスタの左端に出来る.

先ほどの例ではその4ビットに0110, つまり最も左のビット位置の6ができる. それをn-d=60ビット右シフトし, すでに得られていたλと足すと最も左の1のビット位置が得られるのである.

ところで理解して頂けたであろうか. 私はまずSchemeでプログラムを書いて実行したりして, 様子を見ながらアルゴリズムがどうなっているかを考えた. 分かってしまえば簡単なのだが, 最初に述べたアルゴリズムが何をしているか分からないうちは, なんとも不思議に思うアルゴリズムであった.

2014年3月22日土曜日

ビットごとの秘法と技法 から

最も左のビットの位置を知る法

二進法で表したxの最も左の1のビットの位置λ(x)は, 2を底とする対数lgがあれば簡単だ.

λ(x)=⌊lg (x)⌋

たしかに
 
(lg 1) => 0
(lg 2) => 1.
(lg 3) => 1.5849625007211563
(lg 4) => 2.
(lg 5) => 2.321928094887362
(lg 6) => 2.584962500721156
(lg 7) => 2.807354922057604
(lg 8) => 3.
になっている.

TAOCPでの最初の方法は浮動小数点演算命令による.
 FLOTU y,ROUND_DOWN,x; SUB y,y,fone; SR lam,y,52
ここで fone=#3ff0000000000000.

これでうまく行く理由だが, MMIXの浮動小数点はIEEE/ANSI standard 754で, その形は



Sは1ビット, Eは11ビット, Fは52ビットの整数で, S=0なら正, =1なら負. 0<E<2047なら
±2E-1023(1+F/252)
を表す.

だからx=1ならFの部分は0で, Eの部分は1023; x=2ならFの部分 は0で, Eの部分は1024; のようになる. 従って上の命令のように, FOUTU(convert fixed to floating unsigned) で変換し, Eの部分から1023を引いて52ビット右シフトしたのがλになるわけだ.

一方, 浮動小数点を使わないMMIX流の方法は次の通り.
 SRU y,x,32; ZSNZ lam,y,32;
 ADD t,lam,16; SRU y,x,t; CSNZ lam,y,t;
 ADD t,lam,8; SRU y,x,t; CSNZ lam,y,t;
 SRU y,x,lam; LDB t,lamtab,y; ADD lam,lam,t;
まずxを32ビット右シフトしてyに置く.

ZSNZ lam,y,32はyが≠0なら(つまりxの左半分に1があれば)lamに32を, そうでないなら0を入れる. 次はそのlamに16を足してtに置いておく. lamは16か48である. xを16か48ビット右シフトしてyに置く. そのyが≠0ならlamにtを入れ, そうでないならlamはそのまま. 次は8ビットの範囲で同様なことをするので, lamは最も左の1を含む8ビットの範囲の右の境界の値になっている.

SRU y,x,lamでその範囲をレジスタの右端に移動し, 256バイトの表lamtabを引いてlamに足すのである. lamtabは0の場所は使わない. 1のところは0, 2,3のところは1, 4〜7は2, 8〜15は3, ..., 2k〜2k+1-1はk, 128〜255は7になっている.

最も右のビットを取り出すには x & -x とやったが, 左に対してはうまい方法はない. 「ハッカーのたのしみ」の著者 Warrenの考えた方法というのはこうだ.

yxとする. レジスタの長さを2dとして0 ≤ k <dについて yy | (y » 2k)とする. y - (y » 1) がxの最も左ビットである.

yにあるxの最も左のビットは次々とシフトされて右方向に2ビット, 4ビット, 8ビットと広がり, レジスタの右端まで1で埋める(塗りつぶす). 他の右の方にある1はこの塗りつぶしに埋没する. 左端の0の並びはそのままである. そこでyからyを1ビット右にシフトしたものを引くと, xの最も左のビットがとれるわけだ.

TAOCPには

λx = λy   if and only if    xyx & y

という式がある. これも証明は簡単だ.

x, y の最も左のビットの位置が同じとすると, xyのその位置は0になり, x & yのその位置は1になるから, ⊕の方が小さい. 最も左の1の位置が異なると, より左にある1の位置では, ⊕は1, &は0なので⊕の方が大きいことになる.

気になるのは等号の場合だ. x, yが共に0の時, λxもλyも不定になる. 不定同士を=にしたいとすると, ⊕も&も0だからこの場合は等号になるのであろう.

2014年3月21日金曜日

ビットごとの秘法と技法 から

最も右のビットの位置を知る法

TAOCPにある第二の方法はde Bruijnサイクルを利用するものだ.

4ビットのde Bruijnサイクルは例えば次のようなものである.

0000111101001011
サイクルというから円状に描いたのがこの図で, 11時の方向にある0から時計回りに上の数字が並んでいる. 各数字から右方向へ4文字で0から15までのいずれかの二進数になっている.


内側の色の付いた円弧は, その上のビットの二進数が下の同じ色の十進数であることを示す.

この作り方は次の通り.

まず0000と書く. つまり0だ. この4ビットの左端を削除し, 右端に0か1を挿入すると0000と0001が出来る. つまり0の次は0と1だ.

1000からも0と1が出来る. 8の次も0と1だ.

要するに0から15までの数について, 8以上なら8を引き, 8未満ならそのままで, それを2倍したものと2倍して1を足したものが次として出来るわけだ.

0 → 0,18 → 0,1
1 → 2,39 → 2,3
2 → 4,510 → 4,5
3 → 6,711 → 6,7
4 → 8,912 → 8,9
5 → 10,1113 → 10,11
6 → 12,1314 → 12,13
7 → 14,1515 → 14,15
この遷移を有向グラフにしたのが次の図である. それぞれのノードは入次数も出次数も2の正則グラフである.


このグラフですべてのノードを経由して元へ戻るHamilton経路の一例が赤い矢印で, それが上にあったde Bruijnサイクルである. この例は 0000 の右に 1111 を置き, 0100 の後に反対の 1011 を置いた形になっている.

4ビットのde Bruijnサイクルがどのくらいあるか, 手元に計算機があるから全解探索をしてみた. その結果が下の図で16通りあった. 上の目の子で探した解は最下段の左から2つ目である.


そのそれぞれに対するサイクルは次のようだ.

0000100110101111
0000100111101011
0000101001101111
0000101001111011
0000101100111101
0000101101001111
0000101111001101
0000101111010011
0000110010111101
0000110100101111
0000110101111001
0000110111100101
0000111100101101
0000111101001011
0000111101011001
0000111101100101

そこでρ関数の話になる. 最も右の1のビットにするのは常套手段x & -xで, それにde Bruijnサイクルの定数を掛け, 64ビットの左端から6ビットを取ると, 最初の1の位置により, 0から63のいずれかが得られる.

それを先程の4ビットの例でいうと, サイクルは0000111101001011だったから, 1を掛けた時は0000が, 2を掛けた時は1ビット左にシフトしたのの左端4ビットだから0001, 4を掛けた時は0011, ... のように得られる. その様子を次の図で示す.



横長の枠がde Bruijnサイクルを何ビットか左シフトした位置である. 16ビットのレジスタの左端の4ビットに相当する場所にだけ二進数が書いてある. 下の方, 213, 214, 215を掛ける辺りはサイクルの右端の4ビットがレジスタからはみ出しているが, サイクルの左端が0000だったのが幸いしてちょうど輪になっている.

さてこれから分るように

×1は20を掛けたのだから, 表の0番には0を置く.
×2は21を掛けたのだから, 表の1番には1を置く.
×4は22を掛けたのだから, 表の3番には2を置く.
×8は23を掛けたのだから, 表の7番には3を置く.
... とやって出来上がった表が

0, 1, 10, 2, 8, 11, 13, 3, 15, 9, 7, 12, 14, 6, 5, 4

であり, ρ関数は最も右端の1にde Bruijnサイクルの定数を掛け, 左端の4ビットの位置を表で見ると得られるのである.

TAOCPの記述はこうだ.



TAOCPのレジスタは64ビットだから6ビットのde Bruijnサイクルを使う. それは#03f79d71b4ca8b09である. やはり左端が000000で始まっていることに注意しよう. まず6ビットずつ取ると0から63が得られることを確認しよう.
(define as
(map (lambda (n) 
 (quotient (modulo (* #x03f79d71b4ca8b09 (expt 2 n)) 
  (expt 2 64)) (expt 2 58)))
 (a2b 0 64)))

as =>
(0 1 3 7 15 31 63 62 61 59 55 47 30 60 57 51 39 14 29 58 53
43 23 46 28 56 49 35 6 13 27 54 45 26 52 41 19 38 12 25 50 37
10 21 42 20 40 17 34 5 11 22 44 24 48 33 2 4 9 18 36 8 16 32)
ソートして0から63が1回ずつなことを確かめる.
(sort as <)
=>
(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
63)
decodeで使う表を作る.
(define bs 
(map (lambda (n) (- 64 (length (member n as)))) (a2b 0 64)))

bs
=>
(0 1 56 2 57 49 28 3 61 58 42 50 38 29 17 4 62 47 59 36 45 43
51 22 53 39 33 30 24 18 12 5 63 55 48 27 60 41 37 16 46 35 44
21 52 32 23 11 54 26 40 15 34 20 31 10 25 14 19 9 13 8 7 6)
もちろんTAOCPにある表と合っている.

2014年3月19日水曜日

ビットごとの秘法と技法 から

最も右のビットの位置を知る法

TAOCP 7.1.3はbitwise tricks and techniques(ビットごとの秘法と技法)で楽しい話題が満載だ. このブログで以前取り上げたビットスワップもそこにある.

今回の話は計算機にある語xの最も右にある1のビットの位置を計算するもので, x = 0なら1のビットはないから, エラーか不定とするか無限大とするかだが, その辺はどうでもいいのでx≠0 の場合を考える.

xに対してこの位置をρ(x)で表すので, 位置を求めるアルゴリズムをρ関数という. ρは右(right)のRに対応するギリシア語のアルファベット(小文字)である. これに対し最も左のビットはTAOCPではleftのLのギリシア語のλという. λはLisp屋にはちょっと困る命名だ.

二進法での計算に慣れている人は, 最も右の1のビットを取り出すのが簡単なことは知っている. x & -xでそのビットが得られる. x = 0ならどのビットも立たず0である.

例えばxが十進法の10なら, 二進法では1010, -10は二進法では ...0110(左の方は1が何桁も並んでいる)だからandをとると...0010となる. 二進法のビットを右からa0, a1,...と番号をつけると, ...0010の1はa1だから位置としては1とする. ρ(10)=1だ.

奇数は右端が1なのでρ(奇数)=0である.

ρ関数は原始的には右端のビットが1になるまでxを右シフトすればよい.

(define (rho x)
 (do ((i 0 (+ i 1))) ((odd? x) i)
  (set! x (quotient x 2))))

(map rho (a2b 1 17))
=> (0 1 0 2 0 1 0 3 0 1 0 2 0 1 0 4)
x & -xで最も右のビットだけが得られたら, その2を底とする対数をとるという方法も考えられる. TAOCPでは2を底とする対数をlgと書くから,
(define (lg x)
 (/ (log x) (log 2)))
2の63乗までのlgを取ってみると
(do ((i 0 (+ i 1))) ((= i 64))
 (display (lg (expt 2 i))) (newline))
0
1.
2.
3.
4.
5.
6.
7.
8.
9.
10.
11.
12.
13.
14.
15.
16.
17.
18.
19.
20.
21.
22.
23.
24.
25.
26.
27.
28.
29.000000000000004
30.
31.000000000000004
32.
33.
34.
35.
36.
37.
38.
39.00000000000001
40.
41.
42.
43.
44.
45.
46.
47.00000000000001
48.
49.
50.
51.00000000000001
52.
53.
54.
55.00000000000001
56.
57.
58.00000000000001
59.00000000000001
60.
61.
62.00000000000001
63.
ところどころに誤差が出るからroundしinexact->exactすると
(define (rho x)
 (inexact->exact (round (/ (log x) (log 2)))))

(map (lambda (i) (rho (expt 2 i))) (a2b 0 64))
=>
(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42
 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62
 63)
しかしこれらはアルゴリズムとしてはお粗末で, いろいろな工夫がされている. 例えば語の中の1を数える命令(sideway addition SADD)を計算機が持つているなら, x & -xで右端の1だけが取れた後, 1を引くと右端の1は0になり, その右の0は1になるから1の数を数えればよい.

MMIXのSADDはこのためにあるような命令で,
 SUBU t,x,1; SADD rho,t,x
SUBU でxから1を引き, SADDはtのビットが1で, xのビットが0のものの和をとるからrhoに答えが得られる. x = 0の時は, 1を引くとtがオール1になり, xが0だから64ビットの1が足されて64が得られる. TAOCPにはこれを無限大とみてよいと書いてあった.

SADDのような奇妙な命令がない場合のTAOCPの方法はこうだ.

MMIXというアセンブリ言語で書いてあるがそれを写すと
m0 GREG #5555555555555555 ;m1 GREG #3333333333333333;
m2 GREG #0f0f0f0f0f0f0f0f ;m3 GREG #00ff00ff00ff00ff;
m4 GREG #0000ffff0000ffff ;m5 GREG #0000ffff0000ffff;
 NEGU y,x; AND y,x,y; AND q,y,m5; ZSZ rho,q 32;
 AND q,y,m4; ADD t,rho,16; CSZ rho,q,t;
 AND q,y,m3; ADD t,rho,8; CSZ rho,q,t;
 AND q,y,m2; ADD t,rho,4; CSZ rho,q,t;
 AND q,y,m1; ADD t,rho,2; CSZ rho,q,t;
 AND q,y,m0; ADD t,rho,1; CSZ rho,q,t;
最初の3行はマスク用の定数である. GREGは大域レジスタを設定する命令でm0というラベルの場所を十六進(#で示す)の5555555555555555にするというような命令が並ぶ.

1行目の中程「 ;m1 」のセミコロンの次にすぐm1が続いているのは, m1をラベルにしたいからである. 4行目から命令の前に空白があるのは, ラベルがないことを示す. 行頭やセミコロンの直後に文字があるとラベルに扱われる.

さて命令の部分に入って, NEGU y,x; はxのunsigned negativeをとってyに入れる. AND y,x,y;はx & -xを作ってそれをyとした. yには最も右の1のビットだけが立っている.

AND q,y,m5; はyのビットが64ビットの左半分にあればqが0になり, 右半分にあれば, yのままの場所に1が残る.

ZSZ rho,q,32; はzero or set if zeroで, qが0ならrhoに32を入れ, そうでないなら0を入れる. 従ってyのビットが左半分にあればrhoは32になる.


次の行はm4のマスクで調べ, 以前のrhoに一応16を足したものをtに用意した上で, 上の図の16の範囲にあれば前のrhoにtを入れる. CSZは条件付きでqが0ならtをrhoに入れるが, そうでないなら何もしない.

そういう次第で, 上の図の上の2段目の区画のどこにyのビットがあるかでrhoに0, 16, 32, 48 のいずれかが入る.

同様にして8ビットごとの区間のどこにあるかでrhoに8を足す足さないとしてyのある区画を決める.

後も同様で, 変数rhoには答えの位置の番号が入るのである.

TAOCPのこのアルゴリズムの後にはまだ何か書いてある. 「最後の3行を
 SRU y,y,rho; LDB t,rhotab,y; ADD rho,rho,t;
としてもよい. rhotabは129バイトの表(その8個だけを使う)の原点である.」 これは何か.

最後の3行に来たときにはyのビットは上の図の最下段のどれかの8ビットの区画にあり, その右端のビット位置がrhoに入っている. 従ってyを右にrhoビットシフトすると, yのビットは右端の8ビットの区画に移動する. つまり

のどれかになっている. この一番上なら(1なら)0, その次なら(2なら)1, 4なら2, ..., 128なら7をrhoに足せばよい. 従って次のような表を用意すればよいわけである. たしかに129の場所があり, そのうち8カ所だけを使っている. 

MIT Schemeのbit-stringを使って実装してみた.
(define (bmake x) (signed-integer->bit-string 64 x))
(define (band x y) (bit-string-and x y))
(define (bzero? x) (bit-string-zero? x))

(define m0 (bmake #x5555555555555555))
(define m1 (bmake #x3333333333333333))
(define m2 (bmake #x0f0f0f0f0f0f0f0f))
(define m3 (bmake #x00ff00ff00ff00ff))
(define m4 (bmake #x0000ffff0000ffff))
(define m5 (bmake #x00000000ffffffff))

(define (rho x)
 (let ((bx (bmake x)) (by (bmake (- x))) (bq 0) (rho 0) (t 0))
  (set! by (band bx by))
  (set! bq (band by m5))
  (set! rho (if (bzero? bq) 32 0))
  (set! bq (band by m4))
  (set! t (+ rho 16))
  (if (bzero? bq) (set! rho t))
  (set! bq (band by m3))
  (set! t (+ rho 8))
  (if (bzero? bq) (set! rho t))
  (set! bq (band by m2))
  (set! t (+ rho 4))
  (if (bzero? bq) (set! rho t))
  (set! bq (band by m1))
  (set! t (+ rho 2))
  (if (bzero? bq) (set! rho t))
  (set! bq (band by m0))
  (set! t (+ rho 1))
  (if (bzero? bq) (set! rho t))
  rho))

(rho 1000)=>3
(rho 10000)=>4
(rho (+ (expt 2 62) (expt 2 32)))=>32
≥263ではbit-stringを作るところでエラーになったが, 他はうまくいっている.

2014年3月13日木曜日

計算機による音楽演奏

3月11日から13日まで, 東京電機大学東京千住キャンパスで情報処理学会第76回全国大会が開催された. 私は「〜コンピュータパイオニアが語る〜「私の詩と真実」」のセッションで昔話をさせられた.

その中で思い出深いパラメトロン計算機PC-1による音楽演奏にも触れた. この事は約半世紀前には一部の計算機屋によく知られていたが, 最近は知る人もほとんどいないであろうから, ブログに記録していおこうと思う.

話は56年前, つまり1958年にさかのぼる. 私は秋の半ばからニューヨークへ出張した. ニューヨーク滞在中にあちこちに出掛けたが, 12月だったかボストンへ行きMITを訪問した. ちょうどMITにいた高橋研究室の先輩がMITの名所を案内してくれた. 記憶に鮮明に残っているもののひとつは最先端のトランジスタ計算機TX-0であった. いろいろなデモを見せて貰ったなかに, コンソールのタイプライタのキーを押すとキー毎にいろいろな音が出て音楽が演奏できるのがあった. (TX-0は私が1973年から74年にMITでファカルティメンバーとして在籍した頃も学内のどこかにあったらしく, TX-0で出力したというレポートをもってきた学生がいた.)

それはとても楽しそうだったから, 東京に戻ったらなんとかしてPC-1でもやってみたいと思いつつ帰国した.

という次第で以下に述べる音楽演奏の実験をしたのは1959年の2月か3月頃であったろうか.

PC-1は理学部の研究者の科学計算にも使われていたが, もともとは高橋研の実験機だったから, 何かに使えるだろうということでフリップフロップが2個組み込まれていた. 30と31という名前である.

PC-1の命令はアルファベット1文字か, それにバリエーションを示すl(エル)がついているかである. 通常の演算には使わないアルファベットyがあったので, y30とy31がそれぞれ30と31のフリップフロップをセットする命令, yl30とyl31がリセットする命令である.

私はこのフリップフロップに目をつけた. PC-1に研究室に転がっていたマグネチックスピーカを固定し, フリップフロップのセット, リセットでスピーカのコーンを押したり引いたりしたら 音が出せるのではないか.

その話を同僚の相馬君に伝えると早速回路を作ってくれた. それが下の図である.



この辺でパラメトロン回路の見方を説明しなければならない. 図の左下がパラメトロン回路で丸がパラメトロンを表す. パラメトロンは左側が入力, 右側が出力だ. 丸の中のプラスやマイナスも入力だが, これらは定数で, プラスは常に1, マイナスは常に0である.

もうひとつ大事なことはパラメトロンは多数決素子であって, 入力は定数も含めて奇数になっている. 真を1, 偽を0で表すと, 定数がマイナスのパラメトロンはAndであり, 定数がプラスのはOrである.

入力の直前の線と直交する短かい線分はNotを示す.

またパラメトロンは3拍励振といって, それぞれのパラメトロンはI, II, IIIの3つの相のどれかであり, I相の出力がII相の入力, II相の出力がIII相の入力, III相の出力がI相の入力になる. 従って上の図でff31の下の3個のパラメトロンは, 左のプラスのあるパラメトロンをI相とすると, その右のマイナスのパラメトロンがII相, そのさらに右の白丸がIII相になる. また入出力の関係から, 左端のyやlのパラメトロンはII相になる.

他の命令を解読実行している時は左端のyのパラメトロンは0の状態で, その右隣りのAndの出力も0である. この0がフリップフロップの最左の素子に入るが, 定数の1と打ち消してIII相からの情報がそのまま出力になる, つまりフリップフロップの状態は保たれる.

フリップフロップのII相のマイナスのある素子は, 左下からの入力が0だとNotによって1になり, 定数と打ち消される.

yが1, lが0だと上のAndの出力は1になり, フリップフロップは1にセットされる.

yもlも1の時は下のAndが1になり, フリップフロップのII相の素子は2つの入力が0になるから, フリップフロップはリセットされる.

パラメトロン素子の0と1は, 実際には発振の位相が基準のパラメトロンの位相と合っているとき1, 反対位相のとき0という.

そこでIII相の素子の下のプラスだけのパラメトロンは, 1の位相で発振しており, フリップフロップも1の位相のとき, その右のトランスから出力波が得られる. 0だと上下2つの波が消し合い, 出力はない. こういう出力を整流増幅してスピーカに入れている.

プログラムとしては, 一定の時間ごとにy 31とyl 31の命令を繰り返えせばよいから簡単である.

それが下の図だ. 0番地と3番地にy 31とyl 31がある. それらの命令の次に時間を調整するための左シフト命令l nがある. この命令の実行時間はシフトする量nによって変るのでこういう時に便利である. 各命令の右端の{,}内は実行時間で, 基本的な命令は4τである. 但しτはクロック時間, 各相が出力を出す間隔であり, 100マイクロ秒くらいである.



プログラムの5番地には0番地に戻るためにo rという命令が置いてあり, これと時間を合せるために2番地に3番地へジャンプするjl 3rという命令を置く.

ここで今度は命令oを説明しなければならない. oはoutputのoで, 本来はアキュムレータの最上位6ビットをテレタイプに送る命令である. 情報の元はアキュムレータだからこの命令にはアドレス部はいらいない筈である. しかしテレタイプは計算に比べて遲いから, 次のビットを送ろうとしても, テレタイプはまだ前の文字の処理中である公算は大きい. そこでo命令は出力しようとしたとき, テレタイプが準備出来ていなければ, その番地部の示す命令にジャンプするという方式設計になっていた. (入出力と計算を同時に出来るということから, 十人の訴えを一度に聞いたといわれる「聖徳太子の機能」といった.)

仮にテレタイプが1秒に10文字印字するなら, o命令は一旦ビットを送信した後は100ミリ秒は番地の示す先にジャンプする. したがって先のプログラムでは100ミリ秒のループが作れる.

100ミリ秒経つとダミーの情報をテレタイプに送り, ループから下へ抜ける. そこでは次の音高に従い, nの値を設定してまたループに戻ることになる. (ビットが送られると テレタイプはがちゃがちゃいってうるさいから電源を切っておく.)

最後に決めるのはnの値である.



この表の左frequencyの欄はdoに対する各音の周波数を示す. しかし我々は上の図のように周期で考える方が分り易いので, 右の欄を使う. いろいろ試行したが, 下のdoの周期を120にし, reをその9/10倍, miをreの8/9倍, ...のようにすると, reのところだけ胡麻化しだが, なんとか収まるようである. それに対応するnの値を一番右の欄に示した.

doを例にすると, 1周期は120τ, 半分の波の時間は60τ, 命令の固定の時間は4τ + 6τ + 4τ=14τだから, 60から14を引いてn=46になる. 先程のプログラムのシフト命令にこれらのnの値を順々にいれて走らせるとスケールの音が聞こえるのである.

計算機で音を出すことには成功したが, 四分音符しか演奏できないから, 実用にはならなかった. でも手元に計算機があると面白いことが出来るという見本である.

2014年1月14日火曜日

開平法

CurtaのウェブページにDibble-Dabble法という平方根の アルゴリズムがあった.
 2. Dibble-dabble method

2.A Shift the carriage full CCW and zeroize.

2.B Enter "1" in the most significant input slider (IS).

2.C Add ONCE. If the value in the total register exceeds the 
target number, then subtract out the last number added and go 
to step F.

2.D Increase the value in the currently active IS by mentally 
adding "2". If adding "2" to "9" change the "9" to "1" and 
add "1" the the next left input slider (not relevant first 
time through).

2.E Repeat steps C & D.

2.F Reduce the setting of the current IS by "1" and shift the 
carriage CW one step. If the limit of carriage travel is 
reached, go to step I.

2.G Enter "1" in the next most significant IS (the one to the 
right of the one previously used).

2.H Repeat steps C through G.

2.I The number in the "turns counter" register is the true 
square root.
やってみよう. その前に日本語にすると

2.A キャリージを反時計回り(counter clock wise)に一杯に 回し, 零にする(zeroize).

2.B 置数レジスタ(input slider)の最上位に1を置く. 2.C 一度加算する. 結果レジスタ(total register)の値が, 目標の数値を越えるなら, 最後に足した数を引き, Fへ進む.

2.D 置数レジスタを2増やす. 2を9に足す時は, 9を1にし, 置数レジスタの左の 桁に1足す.

2.E CとDのステップを繰り返す.

2.F 置数レジスタの値を1減らし, キャリージを時計回りに1ステップ回す. キャリージ が回らないならIへ進む.

2.G 置数レジスタの次の最上位(前回の桁の右の桁)に1を置く.

2.H CからGのステップを繰り返す.

2.I 回転レジスタの数が平方根になる.

ここに Curtaの精巧なシミュレータがあるからそれを使って2の平方根を求めてみよう.

Curta計算機の部分の名称は下の図を参照してほしい.



Curtaの操作の基本は次の通り.

加減する数は左の絵の下の置数レジスタにつまみを上下して入れる.

加減算の結果は右の絵の結果レジスタに出る. レジスタへの加減する位置は キャリージを回して決める. それにはキャリージの左右の矢印をクリックする.

減算の時は左の絵の上の操作ハンドルを上に引き上げ, 加算の時は 押し下げ, 右の絵の操作ハンドルを回す.

結果レジスタや回転レジスタは右の絵のクリアレバーを回す.

そこで2の平方根を小数点以下5桁求めるために20000000000の平方根を 計算する. (下の説明の行末の番号は後のトレースの行番号.)

まずキャリージの右の矢印を5回押してキャリージを回す.

置数レジスタの6の桁を1にして加算する. 
結果レジスタは10000000000, 回転レジスタは100000になる. (1)
置数レジスタの6の桁を3にして加算する. 
結果レジスタは40000000000, 回転レジスタは200000になる. (2)
40000000000は20000000000より大きいから減算する. (3)
6の桁を1減らして2にし, キャリージを1桁戻す.

置数レジスタの5の桁を1にして(置数レジスタは210000になる)加算
する.(4)
結果レジスタは12100000000, 回転レジスタは110000になる.
置数レジスタの5の桁を3にして加算する. 5にして, 7にして, 9にし
て加算する.
結果レジスタは22500000000, 回転レジスタは150000になる.(8)
20000000000より大きいから減算する.(9)
5の桁を1減らして8にし, キャリージを1桁戻す.

置数レジスタの4の桁を1にして(置数レジスタは281000になる)加算す
る.(10)
3にして加算する.(11)
結果レジスタは20164000000, 回転レジスタは142000になる.
20000000000より大きいから減算する.(12)
4の桁を1減らして2にし, キャリージを1桁戻す.

置数レジスタの3の桁を1にして(置数レジスタは282100になる)加算す
る.
9まで加算すると(17)
結果レジスタは20022250000, 回転レジスタは141500になる.

20000000000より大きいから減算する.(18)
3の桁を1減らして8にし, キャリージを1桁戻す.
この辺までで大体分かったからプログラムを書いてみる.
(define (dibbledabble radicand)
 (let ((s0 1) (s1 1) (s 0) (a 0) (c 0))
  (define (loopp) (display "+") (display s) (display " ")
   (set! a (+ a s)) (set! c (+ c s1)) (display a)
   (display " ") (display c) (newline)
   (if (> a radicand)
    (begin (display "-") (display s) (display " ") 
     (set! a (- a s)) (set! c (- c s1)) (display a)
     (display " ") (display c) (newline)
     (set! s (/ (- s s0) 10))
     (if (> s0 1)
      (begin (set! s0 (/ s0 100)) (set! s1 (/ s1 10)) 
       (set! s (+ s s0)) (loopp)) c))
    (begin (set! s (+ s s0 s0)) (loopp))))
  (do ((n 100 (* n 100))) ((> n radicand))
   (set! s0 (* s0 100)) (set! s1 (* s1 10)))
  (display radicand) (display " ") (display s0)
  (display " ") (display s1) (newline) (set! s s0)
  (loopp)))

(dibbledabble 20000000000)
上のプログラムでletで用意したs0は置数レジスタを増減する数, s1は回転 カウンタを増減する数, sは結果レジスタに加減する数, aは結果レジスタ, cは回転カウンタ, radicandは平方根をとる数である.

radicandを20000000000とした時, まずs0とs1を用意する. それが本体の do loopである. s0は10000000000, s1は100000になる.

上の方のloopが1桁分の計算で, aにsを, cにs1を足す. (if (> a radicand) は結果レジスタが目標を超えた場合で, 置数レジスタと回転レジスタを 戻し, sからも1引いてキャリージを回すように1/10にする. またs0を1/100, s1を1/10にする.

超えない時はsにs0を2回足して加算を繰り返す.

このプログラムを実行した結果が次だ. 左端の斜体の行番号は後から追加した. 0行目はtarget, s0, s1を示す. 後の行はsを足したか引いたかを+, -で示し, s, a, cを順に示している.
0 20000000000 10000000000 100000
1+10000000000 10000000000 100000
2+30000000000 40000000000 200000
3-30000000000 10000000000 100000
4+2100000000 12100000000 110000
5+2300000000 14400000000 120000
6+2500000000 16900000000 130000
7+2700000000 19600000000 140000
8+2900000000 22500000000 150000
9-2900000000 19600000000 140000
10+281000000 19881000000 141000
11+283000000 20164000000 142000
12-283000000 19881000000 141000
13+28210000 19909210000 141100
14+28230000 19937440000 141200
15+28250000 19965690000 141300
16+28270000 19993960000 141400
17+28290000 20022250000 141500
18-28290000 19993960000 141400
19+2828100 19996788100 141410
20+2828300 19999616400 141420
21+2828500 20002444900 141430
22-2828500 19999616400 141420
23+282841 19999899241 141421
24+282843 20000182084 141422
25-282843 19999899241 141421
この1, 3, 5,...と引くのは要するに手回しの機械式計算機の平方根の計算法で あって, それがdibble-dabbleという名前だとは知らなかった. 以前のブログの 5倍してから5, 15, 25,...と引くのは電動計算機の方法というべきであったかも しれない.