割込み

私が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の提案をしていた. 我々はその割込みを先取りしているので大いに気持ちがよかった.

微分解析機

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

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

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



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



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

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

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

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

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

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

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



太陽歯車 S
太陽歯車半径 SA r₀
遊星歯車 P
入力側遊星歯車半径 PB r₁
出力側遊星歯車半径 P'C'' r₂
内歯車半径 SD'' r₃

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

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

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

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


以下9月5日記

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

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

従って
θ₂=θ₀r₀r₂/r₁r₃=0.15θ₀

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

曜日の計算

数学セミナー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等分するのがよいだろう.」と書いているから, そのようにして作図していたのかもしれない.

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

曜日の計算

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

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

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



前回のブログを見た人なら思い出すだろうが, 最上段の右から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日が月曜になるという仕掛である.

曜日の計算

私が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の目盛のカレンダーの値はそのように決めてある.

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

微分解析機

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

そこで毎日微分解析機を実演している. デモは例のサークルテストだが, やはり一周すると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

図はこのようだ.

鶴亀鴉算

小学生のころ鶴亀算を習った. 例えば鶴亀合せて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と判明した.

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