Christopher StracheyのGPM

GPMの続きだ.

二進化

例えば3を11, 6を110, 9を1001へのように十進数dを二進数bに変換したいとする.

d % 2の左に⌊d/2⌋の二進化したもを置けばよいが, 除算も剰余もなければ2を繰り返し引くしかない. とりあえずSchemeで書くと

(define (b x y)
 (if (< x 2)
  (if (= y 0) (list x)
   (append (b y 0) (list x)))
  (b (- x 2) (+ y 1))))

というわけで, GPMにすると

$def,b,<$$lt,~1,2;,
 $def,t,<$>~2<,
  $def,>~2<,<$b,>>~2<<,0;>>~1<;,
  $def,0,>~1<;;>;,
 $def,f,
  <$b,$1-,$1-,>~1<;;,$1+,>~2<;;>;;>;
$b,0,0;,$b,1,0;,$b,2,0;,$b,3,0;,$b,4,0;, => 0,1,10,11,100,
$b,5,0;,$b,6,0;,$b,7,0;,$b,8,0;,$b,9,0; 
=> 101,110,111,1000,1001

二項係数

Pascal三角形を思うえば, C_n,mは両端にある時, つまりm=0かm=nの時は1, それ以外は一段上の左(C_n-1,m)と右(C_n-1,m-1)の和 にすればよい.

$def,b,<$~2,
 $def,~2,<$+,$b,$1-,>~1<;,>~2<;,
  $b,$1-,>~1<;,$1-,>~2<;;;>;,
 $def,0,1;,$def,~1,1;;>;
$def,binom,<$bb,0,
 $def,bb,<$~1,
  $def,~1,<$b,>>~1<<,>~1<;,
  $bb,$1+,>~1<;;>;,
  $def,>~1<,1;;>;;>;
$b,4,0;,$b,4,1;,$b,4,2;,$b,4,3;,$b,4,4; => 1,4,6,4,1
$binom,3; => 1,3,3,1

上のbが二項係数で, $~2でmの値を見る. 下の方$def,0,1;はm=0の時, $def,~1,1;はm=nの時の値を返す. $def,~2, がその他の場合を計算する.

binomはPascal三角形のn段目を計算するもので, bbでmを0からnまで回している.

素数テスト

基本演算で定義した剰余を利用する. nの素数性はxを2から順に増やしながらnをxで割り剰余が0なら偽, x=nになったら真とする.

(define (isprime? n)
 (define (p x)
  (cond ((= x n) #t)
        ((= (modulo n x) 0) #f)
        (else (p (1+ x)))))
 (p 2))

従って, isprime?をp?と書くと

$def,p?,<$p,2,
 $def,p,<$~1,
  $def,~1,<$$r,>>~1<<,>~1<;,
   $def,$r,>>~1<<,>~1<;,
    <$p,$1+,>>~1<<;;>;,
   $def,0,f;;>;,
  $def,>~1<,t;;>;;>;
$p?,2;,$p?,3;,$p?,4;,$p?,5;, => t,t,f,t,
$p?,6;,$p?,7;,$p?,8;,$p?,9; => f,t,f,f

p?の定義はまず$p,2,と(p 2)を実行し, すぐにpの定義が続く. $def,p,<$~1, の~1はSchemeのプログラムのxである. その次の行の$def,~1,はelseの部分. 一番下の 行の$def,$gt;~1<<,の~1はn, すなわちx=nならtという定義がこの行だ.

else部分に戻ると$$r,>>~1<<,>~1<,;,とあるが, ここがnをxで割った剰余を計算するところで, >>~1<<がn, >~1<がxである. その剰余でマクロ呼出しし, 下の方の$def,0,f;が割り切れた場合は偽と定義する. elseの定義は$def,の後でもう一度剰余を計算する. やることは$p,$1+,x;である.

tarai関数

竹内郁雄君の発案したtarai関数はGPMでやってみるには都合がよい. tarai関数の定義は

(define (tarai x y z)
  (if (<= x y) y
  (tarai (tarai (- x 1) y z)
         (tarai (- y 1) z x)
         (tarai (- z 1) x y))))

で, (trai 4 2 0) とかやってみると, 盥まわしの様子が分かる. しかしz=0なので早速-1が現れて問題となる. $1-,0;は実行出来て-1になる. 基本関数のltの 引数の下が-1まで使えるようになっているのは, ここで使いたかったからである.

次は(<= x y). 基本関数にあったのはleではなく, ltであったが, もちろん(< y x)の形で使い, then部分とelse部分を交換して書いておく.

従ってtarai関数は

$def,tarai,<$$lt,~2,~1;,
 $def,f,~2;,
 $def,t,<$tarai,
  $tarai,$1-,>~1<;,>~2<,>~3<;,
  $tarai,$1-,>~2<;,>~3<,>~1<;,
  $tarai,$1-,>~3<;,>~1<,>~2<;;>;;>;
$tarai,4,3,2; => 4
$tarai,4,2,0; => 4

この辺でGPMの空白改行問題を説明しなければならない. StracheyのGPMの論文には, マクロ呼出しは評価の文字列に置き換わるがそれ以外は入力がそのまま出力されると書いてある. アセンブリ言語の前処理用としてはその通りであるが, tarai関数のマクロ定義をこのように整形しておくと実はうまく走らないのである.

(tarai 4 3 2)の実行され方を見ると, (<= 4 3) は #fなので(tarai (- 4 1) 3 2), (tarai (- 3 1) 2 4), (tarai (- 2 1) 4 3) つまり(tarai 3 3 2), (tarai 2 2 4), (tarai 1 4 3)をまず計算する. この3個はどれも(<= x y)なので, それぞれ3, 2, 4であり, 次に(tarai 3 2 4)を計算しなければならない.

これも(<= 3 2)ではないので, (tarai 2 2 4), (tarai 1 4 3), (tarai 3 3 2)を計算し, これらは直接終わるので(tarai 2 4 3)の計算に移り, (<= 2 4)だから4となるわけだ.

問題は

$tarai,
 3,
 2,
 4;

になった時に改行や空白が邪魔になることである. このtaraiの第1引数は`改行空白3', 第2引数は`改行空白2', 第3引数は`改行空白4'であり, これらが$1-に渡されてしまう.

そういう次第で, GPMの処理系では改行や空白は無視するようにしてあるが, この後で出て来る例題では改行や空白をそのまま使いたいものもあって, その辺は どう対処するのがよいか疑問である.

BCPLを設計したMartin RichardsはBGPMというGPM処理系を使っているそうだが, そのBGPMではバッククォート(`)をエスケープに使い, バッククォートから後その行の最後までと次の行から空白を無視する仕様になっているとそうだ.

Christopher StracheyのGPM

前回はHanoiの塔まで説明した. また続きのマクロを示す.

Fibonacci数

GPMで遊ぶのに適している例題の一つがFibonacci数である.

(define (fib n)
 (cond ((= n 0) 0)
       ((= n 1) 1)
       (else (+
        (fib (- n 1))
        (fib (- n 2))))))

とりあえずGPM風にすると

$def,fib,<$n,
 $def,n,<$+,
  $fib,$1-,n;;,
  $fib,$1-,$1-,n;;;;>;
 $def,1,1;,
 $def,0,0;;>;

+のマクロは基本演算で定義した.

nはもちろん~1にする. <,>が2重に使われているが, 内側のクォートの内部では ~n は>~n< にする.

従って

$def,fib,<$~1,
 $def,~1,<$+,
  $fib,$1-,>~1<;;,
  $fib,$1-,$1-,>~1<;;;;>;
 $def,1,1;,
 $def,0,0;;>;

やってみると

$fib,0; => 0
$fib,1; => 1
$fib,2; => 1
$fib,6; => 8

階乗

次は階乗. マクロ名に!が使えて嬉しい. *は基本演算参照

$def,!,<$~1,
 $def,~1,<$*,>~1<,$!,$1-,>~1<;;;>;,
 $def,0,1;;>;

と簡単だ. 実行例は

$!,0; => 1
$!,1; => 1
$!,2; => 2
$!,3; => 6

中央値

英語ではmedian. TAOCPの7.1.1に登場する. 奇数個の値をソートしてa₀,a₁,...,a_2nが得られた時のa_nを値とする.

<1,0,4,2,3>=2だ. 奇数個の値がfalseとtrueだけとし, false<trueとした時, 中央値は多数決になる. <false,false,true>=false,<true,false,true>=true. 中央値は多数決を一般化したものである.

ここでは3個の値の中央値を見つける.

(define (med a b c)
 (if (< b a)
  (if (< c a)
   (if (< c b) b c)
   a)
  (if (< c b)
   (if (< c a) a c)
   b)))

このGPM版は

$def,med,<$$lt,~2,~1;,
  $def,t,$$lt,~3,~1;,$def,t,$$lt,~3,~2;,$def,t,~3~2~1;,
                                        $def,f,~2~3~1;;;,
                     $def,f,~2~1~3;;;,
  $def,f,$$lt,~3,~2;,$def,t,$$lt,~3,~1;,$def,t,~3~1~2;,
                                        $def,f,~1~3~2;;;,
                     $def,f,~1~2~3;;;;>;

実行すると

$med,0,0,0;,$med,0,0,1;,$med,0,0,2;,$med,0,1,0;,$med,0,1,1;,
$med,0,1,2;,$med,0,2,0;,$med,0,2,1;,$med,0,2,2;,$med,1,0,0;,
$med,1,0,1;,$med,1,0,2;,$med,1,1,0;,$med,1,1,1;,$med,1,1,2;,
$med,1,2,0;,$med,1,2,1;,$med,1,2,2;,$med,2,0,0;,$med,2,0,1;,
$med,2,0,2;,$med,2,1,0;,$med,2,1,1;,$med,2,1,2;,$med,2,2,0;,
$med,2,2,1;,$med,2,2,2;
=>
0,0,0,
0,1,1,
0,1,2,
0,1,1,
1,1,1,
1,1,2,
0,1,2,
1,1,2,
2,2,2

GCD

Euclidの互除法が有名だが, 9までの数を扱うこのGPMの世界では引き算で計算できる.

(define (gcd a b)
 (cond ((= a b) a)
       ((< a b)
        (gcd a (- b a)))
       ((> a b)
        (gcd (- a b) b))))

GPMに書き直す.

$def,gcd,<$~2,
 $def,~2,
  <$$lt,>~1<,>~2<;,
  $def,f,
   <$gcd,$-,>>~1<<,>>~2<<;,>>~2<<;>;,
  $def,t,
   <$gcd,>>~1<<,$-,>>~2<<,>>~1<<;;>;;>;,
 $def,~1,~1;;>;
$gcd,2,4;,$gcd,5,3;,$gcd,6,3; => 2,1,3

この辺までは簡単の単だ.

Christopher StracheyのGPM

このブログでGPMのことを書いたのは2013年12月だから1年以上前になる. 英国の計算機保存団体の発行するComputer Resurrectionの昨年の夏号にStrachey's General Purpose Macrogeneratorという記事を見付けたのを切っ掛けに元の論文にあった例題以外のマクロを書こうと思いたった.

かなりたくさんのマクロを書いたので, GPMの楽しさを紹介しようとして今年のプログラミング・シンポジウムでそれに関する話をした.

これとそれに続く何回かのブログでは, それらのマクロについて書きたい.

まず簡単におさらいをしておく.

"$def,a,<b~1d>;" はマクロ "a" が "<b~1d>" であると定義する. "$a,c;" でマクロ "a" を実引数 "c" をもって呼び出すと, マクロ本体の第1引数 "~1" に "c" が代入されて "bcd" が返る. "~0" は第0引数を表わすから, マクロ "a'" を "$def,a',<~0b~1d>;" と定義し "$a',c;" と呼び出すと "a'bcd" が返る.

マクロ呼出しは "$マクロ名, 実引数1, 実引数2, ... ;" の形で, マクロ名も実引数もまず評価される, つまりマクロ呼出しがあれば実行し, 引数"~n" があれば置き換える. それ以外はそのまま. 評価を避けたい時は"<" と ">" で囲むが評価が済むと両端の "<" と ">" がなくなる.

上の "a" の定義があった時, "$a,$a,c;;" と呼び出すと, 第1引数 "$a,c;" がまず評価されて "bcd" になり, それで "$a,bcd;" と呼び出すから "bbcdd" が返る.

マクロ呼出しの実引数の並びの中に局所定義を書くことが出来る. 定義されてから呼び出すから, "$a,c,$def,a,<b~1d>;;" のように書くと, "bcd" が返る. 局所定義は呼出しが終わると消える.

同じマクロ名の定義はスタック状に記憶され, 最新のものを使う. 定義が消えるともとの定義が見えるようになる.

"$def,a,b;" と定義しておき, "$a,$def,a,c;,$def,a,d;;" と呼び出すと, "a" の最後の定義が "$def,a,d;" だから "d" が返り, それが済むと最初の定義がスタックの上に来るから, ここで "$a;" と呼び出すと "b" が返る.

この機能を使うと条件式が書ける.

"(if (eq? α β) γ δ)" は"$α,$def,α,δ;,$def,β,γ;;" と書く. αとβが同じ時は, 2つの局所定義の後のもαになるからそちらを使い, γになり, 違う時は, 前の局所定義を使い, δになる.

これだけ知っていれば, GPMのマクロを書くことが出来る.

基本演算

GPMには算術演算もないから, その辺から始めなければならない. 最初は引数に1を足すマクロ "1+". 定義と使用例は次のとおり.

$def,1+,<$1,2,3,4,5,6,7,8,9,10,$def,1,<~>~1;;>;
$1+,0; => 1
$1+,4; => 5
$1+,9; => 10

"$1+,4;" の場合, マクロ "1+" の本体 "$1,2,3,4,5,6,7,8,9,10,$def,1,<~>~1;;" が評価される. つまりマクロ "1" を実引数 "2,3,4,5,6,7,8,9,10" で呼び出す. その"1"の局所定義は後にあり, "~1" に "4" が代入されて "$def,1,~4;" になっている. 従って第4引数 "5" が返る.

$def,1-,<$-1,0,1,2,3,4,5,6,7,8,$def,-1,<~>~1;;>;
$1-,0; => -1
$1-,9; => 8

"1-"も上のマクロ "1+" と殆ど同じ. 引数が"~"の後の1文字なので, ここで使う整数は0から9までということになる.

$def,+,<$~1,
 $def,~1,<$1+,$+,$1-,>~1<;,>~2<;;>;,
 $def,0,~2;;>;
$+,0,3; => 3
$+,3,5; => 8

マクロ "+" は条件式になっている. 最後の行 "def,0,~2;" は, 第1引数が "0" なら第2引数が答, そうでないなら, 第1引数から1を引いたものと第2引数と足してそれに1を足すと再帰的に定義する.

Scheme流に書けば
(λ (x y) (if (= x 0) y (1+ (+ (1- x) y))))

$def,-,<$~2,
 $def,~2,<$-,$1-,>~1<;,$1-,>~2<;;>;,
 $def,0,~1;;>;
$-,3,0; => 3
$-,6,4; => 2

第2引数による条件式で, 0なら第1引数が答. そうでないなら第1引数から1を引いたものから第2引数から1を引くと再帰的に定義する.

(λ (x y) (if (= y 0) x (- (1- x) (1- y))))

$def,*,<$~1,
  $def,~1,<$+,$*,$1-,>~1<;,>~2<;,>~2<;>;
  $def,0,0;;>;
$*,0,5; => 0
$*,1,4; => 4
$*,2,3; => 6

(λ (x y) (if (= x 0) 0 (+ (* (1- x) y) y)))

$def,lt,<$~1,
  $def,~1,<$p,>~1<,>~2<,$def,p,$lt,$1-,>~1<;,>~2<;;;>;
  $def,-1,t;$def,~2,f;;>;
$lt,-1,-1; => f
$lt,-1,0;  => t
$lt,0,-1;  => f

ltつまり<は, この世界の下が-1までであることを利用している.

(λ(x y) (cond ((= x y) f) ((= x -1) t) (else (lt (1- x) y))))

これを使って剰余を書く.

$def,r,<$$lt,~1,~2;,
 $def,t,~1;,$def,f,<$r,$-,>~1<,>~2<;,>~2<;>;;>;
$r,9,5; => 4

(λ (x y) (if (< x y) x (r (- x y) y)))

基本演算ばかりでは面白くないので, 応用として情報科学標準問題のHanoiの塔をやってみよう.

Hanoiの塔

$def,1-,<$-1,0,1,2,3,4,5,6,7,8,
 $def,-1,<~>~1;;>;
$def,hanoi,<$~4,
 $def,~4,<$hanoi,>~1<,>~3<,>~2<,$1-,>~4<;;
  +>~1<->~3<+
  $hanoi,>~2<,>~1<,>~3<,$1-,>~4<;;>;
 $def,0,<>~1<->~3<>;;>;
$hanoi,a,b,c,0; => a-c
$hanoi,a,b,c,1; => a-b+a-c+b-c
$hanoi,a,b,c,2; 
=> a-c+a-b+c-b+a-c+b-a+b-c+a-c

$hanoi,a,b,c,n;はn枚の円板をaからcへbを中継点として移動する手続きで, nが0なら$def,0,にあるように ~1-~3とする. -は引き算ではなく, 引数の間に-を書くことである. そうでないなら, $hanoi,~1,~3,~2,$1-,~4;; をやり +~1-~3+とし, $hanoi,~2,~1,~3,$1-,~4;; をやると定義してある.

結果は上のようだ. 次のブログへ続く.

Piの近似値

Martin GardnerのPiの話を読んでいたら, 下のような図があり, おやっと思った.


円の半径を1とすると, FGの長さが16/113になるというのだ. 本当かと思い, 計算してみた.

図はこのように描いてある. ACGは半径1の円の1/4である.

BはCAのCから7/8の点, DはBGのGから1/2の点である.

DからACと平行にDEを引き, CG上の点をEとする. DからBEと平行にDFを引き, CG上の点をFとする.

BGの長さを計算すると BG²=(7/8)²+1²=(49+64)/8² =113/8²

BCGとDEGは相似だから EG/CG=DG/BG ∴ EG=CG*DG/BG

BEGとDFGは相似だから FG/EG=DG/BG ∴ FG=EG*DG/BG

この2つから FG=CG*DG/BG*DG/BG

CG=1だから FG=DG²/BG²=(1/4)/113/8²=16/113

もう何だか分った人も多いだろうが, この値に3を足すと(3*113+16)/113=355/113で円周率の近似値になるのである. つまりFGは円周率の小数部分であった.

微分解析機

このブログの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てあるのは, 固定していないことを示す.



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

クリスマスの歌

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

そこで歌われる歌の中に, 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日足りないのか.

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ビットである. 番地部が2ⁿなら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' ³を加え, 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より小さいときは別になにもしないようだ. 不思議な仕様である. 他の例を見ても正の場合しか考慮していない.

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

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

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

0.8192₁₀= 0.11010001101101110001011101011000111000100001100101101₂

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