Christopher StracheyのGPM

12 Days of Christmas

昨年12月12日のブログにクリスマスの歌を書いたのは, GPMのこのマクロを書いた直後であった.

それまでいろいろなプログラム言語で書いてみたので, GPMならどう書けるかと考えるのはまぁ自然だ.

こういう詩のようなものを出力するには, 空白や改行の制御が必要になる. GPMはマクロ評価部分以外は入力をそのまま出力するから簡単かもしれない. 一方, Tarai関数で述べたように, 空白や改行は邪魔なこともあるので, 私のシミュレータは空白改行を無視する/しないの切替え機能を持っている.

またこのマクロは, 本質的なところだけ再現出来ればよいという趣旨で, 元の詩の一部だけを使った.

とはいえ, 難題がある. このGPMの環境では, 通常は0から9までしか数えられないのに, 12日をどう扱うかだ.

とりあえずマクロはこうなっている.

$def,1-,<$-1,0,1,2,3,4,5,6,7,8,9,:,$def,-1,<~>~1;;>;
$def,getn,<$11th,10th,9th,8th,7th,6th,5th,4th,3rd,2nd,1st,
$def,11th,<~>~1;;>;
$def,gets,<$11 Pipers
,10 Lords
,9 Ladies
,8 Maids
,7 Swans
,6 Geese
,5 Rings
,4 Birds
,3 Hens
,2 Turtle Doves
,a Partridge
,$def,11 Pipers
,<~>~1;;>;
$def,song,<$~1,
$def,~1,<On the $getn,>~1<; day
$gets,>~1<;>~2<
$song,$1-,>~1<;,$gets,>~1<;>~2<;>;,
$def,-1,<On the 12th day
12 Drummers
>~2;;>;
$song,:,;

最初のマクロ1-はこれまでの1-とそっくりで, 引数が0なら-1, 1なら0, ... が返る. これは~nを読むとnのアスキーコードから48を引き何番目かを決める. ところでアスキーの表で9の次が:(コロン)なのを利用して:を10として扱うことにした. では;(セミコロン)を11として使えるかというと;はマクロ呼出しの終りを示すからこれは使えない. だが0から9の世界を0から10の世界に拡張できた. それでもまだ12日には問題が残る.

getnは引数が0なら11th, 1なら10th, ..., 10なら1stが返る. 次のgetsは同様で, 0なら11 Pipers改行, 1なら10 Lords改行, ... ,10ならa Partridge改行が返る. カンマが改行の次にあるからだ. 何行か飛ばして最後の行のsongの呼出しを見よう.

$song,:,;とマクロ呼出しになっている. つまり第1引数は:(コロン, 上述のように10), 第2引数は空文字列になっている.

そこでsongの定義をみると, 第1引数のifになっており, -1なら$def,-1,にあるようにOn the 12th day改行12 Drummers改行にして第2引数をつなげる. else部はOn the (getnの返した1st day改行みたいなもの), 続いて$gets, ~1;続いて第2引数. さらにsongを呼ぶようになっている.

という次第で

On the 1st day
a Partridge

On the 2nd day
2 Turtle Doves
a Partridge
...(3rdから11thまで省略)
On the 12th day
12 Drummers
11 Pipers
10 Lords
9 Ladies
8 Maids
7 Swans
6 Geese
5 Rings
4 Birds
3 Hens
2 Turtle Doves
a Partridge

の出力が得られた.

Color RingのRGB

つぎのような虹の両端を繋げた絵が描きたかった. 各色のRGBは円周外に併記してあるようになっている. 赤字がRの値, 緑字がGの値, 青字がBの値である. それぞれの値は0から0xffまでである.



円周外の文字は見難いいので, 変化の様子を図にすると



上の円形の図に右方向の赤が横軸の0で, RGBがff0000, そこから60度の間に緑が00からffへ増える(横軸30の地点). 次は赤が減り, 青が増え, 緑が減り, 赤が増え, 青が減るような変化である.

これは円周360度を36分割してあるが, さらに連続的な絵も描けるわけで, 2度置きに色を変えるとすると

for(i=0; i<180; i++) {
var ang0=-i*Math.PI/90-Math.PI/180;
var ang1=-i*Math.PI/90+Math.PI/180;
context.beginPath();
context.moveTo(r*Math.cos(ang1),r*Math.sin(ang1));
context.lineTo(R*Math.cos(ang1),R*Math.sin(ang1));
context.arc(0,0,R,ang1,ang0,true);
context.lineTo(r*Math.cos(ang0),r*Math.sin(ang0));
context.closePath();
context.fillStyle=cs[i];
context.fill();}

のように扇型のような領域を下にあるcs[i]のRGBで塗り分けることになる.

ここでのマクロはcs[i]の値を生成するものである.

var cs = new Array (
"#ff0000", "#ff0900", "#ff1100", "#ff1a00", "#ff2200", "#ff2b00", 
"#ff3300", "#ff3c00", "#ff4400", "#ff4d00", "#ff5500", "#ff5e00", 
"#ff6600", "#ff6f00", "#ff7700", "#ff8000", "#ff8900", "#ff9100", 
"#ff9a00", "#ffa200", "#ffab00", "#ffb300", "#ffbc00", "#ffc400", 
...
"#ff0033", "#ff002b", "#ff0022", "#ff001a", "#ff0011", "#ff0009",
0);

60度の間に0から255まで変化する. 2度ずつだと255を30で割るから1回の増減量は8.5に相当する. 緑の増え方をみると, 最初が0, 次が9(8.5を四捨五入した), 次が 十六進で11つまり17(8.5*2), ... となっていることが分る.

しかしこの計算はせず, この昇順と降順の値をマクロで定義した.

$def,00^,09;$def,09^,11;$def,11^,1a;$def,1a^,22;$def,22^,2b;
$def,2b^,33;$def,33^,3c;$def,3c^,44;$def,44^,4d;$def,4d^,55;
$def,55^,5e;$def,5e^,66;$def,66^,6f;$def,6f^,77;$def,77^,80;
$def,80^,89;$def,89^,91;$def,91^,9a;$def,9a^,a2;$def,a2^,ab;
$def,ab^,b3;$def,b3^,bc;$def,bc^,c4;$def,c4^,cd;$def,cd^,d5;
$def,d5^,de;$def,de^,e6;$def,e6^,ef;$def,ef^,f7;$def,f7^,ff;
$def,09_,00;$def,11_,09;$def,1a_,11;$def,22_,1a;$def,2b_,22;
$def,33_,2b;$def,3c_,33;$def,44_,3c;$def,4d_,44;$def,55_,4d;
$def,5e_,55;$def,66_,5e;$def,6f_,66;$def,77_,6f;$def,80_,77;
$def,89_,80;$def,91_,89;$def,9a_,91;$def,a2_,9a;$def,ab_,a2;
$def,b3_,ab;$def,bc_,b3;$def,c4_,bc;$def,cd_,c4;$def,d5_,cd;
$def,de_,d5;$def,e6_,de;$def,ef_,e6;$def,f7_,ef;$def,ff_,f7;

Hilbert曲線の時と同様, ^は上向き, _は下向きとし, $def,00^,09;は上向きで00の次は09という定義である.

$def,next,<$~1^;>;
$def,last,<$~1_;>;

と定義すると

$next,00; => 09
$next,$next,00;; => 11

のようになる. 次にfとbを定義する.

$def,f,<$~2, 
 $def,~2,<"#>~1~2~3<",$f,>~1<,$next,>~2<;,>~3<;>;,
 $def,ff,;;>; 
$def,b,<$~2, 
 $def,~2,<"#>~1~2~3<",$b,>~1<,$last,>~2<;,>~3<;>;,
 $def,00,;;>; 

そして

$f,ff,00,00;

と呼ぶと, 第1引数はff, 第2引数は00, 第3引数も00である. fの定義は第2引数=ffなら何もしない; そうでないなら"#の次に3つの引数をつなげ",を出力しfを続けて呼ぶ. 但し第2引数をnextにする となっている. 従って

"#ff0000","#ff0900","#ff1100",...

が得られる. 続いてbを$b,,ff,ff00;と呼ぶと, 第1引数は空文字列, 第2引数はff, 第3引数はff00だから, RGBのRの値を次々と減らした列が得られるのである.

こうして出来た配列で描いたColor Ringがこれだ.

HP-16Cのプログラム技法

私の個人用電卓で他の電卓に見当たらない便利な機能はユリウス日の計算である.

-4712年1月1日を0とした通日である. 例えば2015年3月2日は20150302と入力し, JDのキーを押すと2457084が得られる. これに1を足し7で割った剰余を計算すると1になり月曜であるのが分る.

これをHP-16Cのプログラムにしたいが, グレゴリオ改暦の前を使うことはまずないので(-4712年1月1日が0になるか確認するくらいだ), 今回はFixed Day Numberを計算することにした.

Fixed Day Number(RD, Rata Die)は改暦よりも前までグレゴリオ暦が使われていたとした, 1年1月1日 (月曜)を1にするものである. 2015年3月2日RDは735659である. ユリウス日より1721425日少ない.

Calendrical Calculationのアルゴリズムは次のようだ(同書ではCommon Lispを使う)

(defun fixed-from-gregorian(year month day)
(+ (* 365 (1- year))
   (quotient (1- year) 4)
   (- (quotient (1- year) 100))
   (quotient (1- year) 400)
   (quotient (- (* 367 month) 362) 12)
   (if (<= month 2) 0
    (if (gregorian-leap-year? year) -1 -2))
   day))

要するに年数に365を掛け, 閏年の調整をし, 月始めまでの日数を求め, 日付を足す. HP-16Cのプログラムはこうだ.

001 43.22.b  LBL B ;20150302 GSB B
002 44.00    STO 0 ;ymd -> 0
003 43.05.00 CF 0  ;flag 0 閏年
004 43.05.01 CF 1  ;flag 1 m<=2
005 04       4
006 00       0
007 00       0
008 44.04    STO 4 ;400->4
009 42.b     SR
010 42.b     SR
011 44.03    STO 3 ;100->3
012 42.09    RMD   ;ymd%100
013 44.02    STO 2 ;d -> 2
014 02       2     ;m<=2の比較用
015 45.00    RCL 0 ;ymd
016 45.03    RCL 3 ;100
017 10       /
018 44.00    STO 0 ;ym
019 45.03    RCL 3 ;100
020 42.09    RMD   ;ym%100
021 44.01    STO 1 ;m -> 1
022 43.01    x<=y
023 43.04.01 SF 1  ;m<=2ならF1をセット
024 45.00    RCL 0 ;ym
025 45.03    RCL 3 ;100
026 10       /
027 44.00    STO 0 ;y -> 0
028 45.03    RCL 3
029 42.09    RMD
030 43.40    x=0   ;y%100=0?
031 22.00    GTO 0
032 45.00    RCL 0 ;y
033 04       4     ;4
034 22.01    GTO 1
035 43.22.00 LBL 0
036 45.00    RCL 0 ;y
037 45.04    RCL 4 ;400
038 43.22.01 LBL 1
039 42.09    RMD
040 43.30    x=0
041 43.04.00 SF 0 ;y%(4or400)=0
042 03       3
043 06       6
044 07       7
045 44.05    STO 5 ;367 -> 5
046 05       5
047 30       -
048 44.06    STO 6 ;362 -> 6
049 03       3
050 40       +     ;365
051 45.00    RCL 0 ;y
052 01       1
053 30       -     ;-1
054 44.00    STO 0 ;y-1 -> 0
055 20       *     ;365*(y-1)
056 45.00    RCL 0
057 04       4
058 10       /
059 40       +     ;+(y-1)/4
060 45.00    RCL 0
061 45.03    RCL 3
062 10       /
063 30       -     ;-(y-1)/100
064 45.00    RCL 0
065 45.04    RCL 4
066 10       /
067 40       +     ;+(y-1)/400
068 45.01    RCL 1 
069 45.05    RCL 5
070 20       *     ;m*367
071 45.06    RCL 6
072 30       -     ;-362
073 01       1
074 02       2
075 10       /     ;/12
076 40       +
077 45.02    RCL 2
078 40       +     ; +day
079 43.06.01 F? 1
080 43.21    RTN   ;m<=2なら帰る
081 02       2     ;平年なら2
082 43.06.01 F? 0
083 42.b     SR    ;閏年なら1
084 30       -     ;引く
085 43.21    RTN   ;帰る

テストはクリティカルな日付で行う

20150228 => 735657
20150301 => 735658

20120229 => 734562
20120301 => 734563

20141231 => 735598
20150101 => 735599

00010101 => 1

RDは7での剰余がそのまま曜日になる(1年1月1日月曜を1にしたから.)

HP-16Cのプログラム技法

2月21日のブログで疑似素数発生のプログラムが出来たから, 次は素因数分解のプログラムを書こう. 私がSchemeの環境で使っているfactorizeはこういうものだ.

 (factorize 32) => (2 2 2 2 2)
 (factorize 12) => (3 2 2)
 (factorize 108) => (3 3 3 2 2)
 (factorize 97) => (97)

だから素数ならこのリストの長さが1になるわけで

 (define (prime? n)
  (= (length (factorize n)) 1))

 (map prime? (a2b 2 10)) => (#t #t #f #t #f #t #f #f)

だが, 問題はn=1の時だ.

 (factorize 1) => (1)

と1も素数に判定されてしまう. それはさておき...

個人用電卓Happy Hacking Calculatorでは, n>=2の素因数分解は, pをnの最小の因数としてn=p^c*q となるようなp, c, qをスタックに積むようにした. つまり

 n fact -> p, c, q
 2 fact -> 2, 1, 1
 4 fact -> 2, 2, 1
 6 fact -> 2, 1, 3
 8 fact -> 2, 3, 1
12 fact -> 2, 2, 3

のようになる. 計算が終了した時, 最小素因数pが見える. それが最初に入れた数nなら素数だった. そうでないなら, その最小素因数が何回あったかcがスタックのすぐ次で分かる. さらに次に最小素因数で割り切った残りqがあるから, 続けてその素因数分解を始める. だから, 108(=2 * 2 * 3 * 3 * 3)は

108 fact -> 2, 2, 27 popを2回
27 fact -> 3, 3, 1

で108=2²*3³ということが分るのである.

これが便利だったので, HP-16Cのプログラムでも同様なインターフェースとした.

Happy Hacking Calculatorの素因数分解のアルゴリズムをSchemeで示すと

(define (hhcfactorize n)
 (let ((fl 疑似素数差分リスト) (d 2) (c 0) (i 0))
  (define (floop)
   (define (gloop) (set! n (/ n d)) (set! c (+ c 1))
    (if (= (modulo n d) 0) (gloop)))
   (if (< (quotient n d) d) (list n 1 1)
    (if (= (modulo n d) 0) (begin (gloop) (list d c n))
     (begin (if (= n 485) (set! i 5))
       (set! d (+ d (list-ref fl i))) (set! i (+ i 1)) (floop)))))
(floop)))

のように書いたので, 今回はこれをHP-16C用に書き換える. 上のプログラムでは配列flが疑似素数表だが, それに前回のブログのものを使う.

001 43.22.A LBL A
002 44.0    STO 0 ;n->0
003 6       6     ;差分を設定
004 44.01   STO 1
005 44.03   STO 3
006 4       4
007 44.04   STO 4
008 44.06   STO 6
009 44.08   STO 8
010 2       2
011 44.05   STO 5
012 44.07   STO 7
013 44.09   STO 9
014 44.A    STO A
015 44.d    STO D ;2->d
016    1    1
017 44.b    STO B
018 43.35   CLx
019 44.c    STO C ;0->c
020 1       1
021 1       1
022 44.32   STO I ;11->i
023 43.22.0 LBL 0 ;ループバック地点
024 45.0    RCL 0 ;n
025 45.d    RCL D ;d
026 10      /
027 45.d    RCL D ;d in x n/d in y
028 43.01   x<=y ;d<=n/d?
029 22.01   GTO 1
030 45.b    RCL B ;商が除数より小さくなった
031 45.b    RCL B ;1,1,nを積んでもどる
032 45.00   RCL 0 ;n
033 43.21   RTN
034 43.22.1 LBL 1 ;n/d>=d
035 43.06.4 F? 4  ;剰余!=0フラッグ
036 22.03   GTO 3
037 43.22.2 LBL 2 ;divisible
038 45.00   RCL 0
039 45.d    RCL D
040 10      /
041 44.00   STO 0 ;n/d->n
042 45.c    RCL C
043 1       1
044 40      +
045 44.c    STO C ;c+1->c
046 45.0    RCL 0 ;n
047 45.d    RCL D ;d
048 42.09   RMD ;remainder
049 43.40   x=0
050 22.02   GTO 2
051 45.00   RCL 0 ;n
052 45.c    RCL C ;c
053 45.d    RCL D ;d
054 43.21   RTN
055 43.22.3 LBL 3 ;割り切れなかったので次の除数を作る
056 45.d    RCL D
057 45.31   RCL (I)
058 40      +
059 44.d    STO D
060 43.23   DSZ
061 22.0    GTO 0
062 8       8
063 44.32   STO I
064 22.00   GTO 0

017までは1からBのレジスタに差分を入れる. 2の時には最初の除数もDへ入れる(015). 019でカウンタCをクリアする. 022でIレジスタの初期値を11に設定. 023のLBL 0がループバックの場所.

n/d<dを調べ, 終わりなら030から1,1,nをスタックに積んでRTN.

035は026の除算で剰余が0でない, 割り切れなかった時はフラッグ4が立っているから, F? 4で調べ, 立っていればLBL 3へ飛ぶ.

050までは同じ除数で割れるだけ割り(cも増やしつつ), 割れなくなったら051からn, c, dをスタックに積んでRTN.

そこから下は前回の手法による除数の更新である.

このプログラムが消費したメモリーをしらべると, プログラムは64バイトだが, 7の倍数ずつ領域を確保するから, 70バイト使う. レジスタは0からDまで32ビット= 4バイトを14個つかったら56バイト, 合わせて126バイトだったから, 全体で203バイトのメモリーの62パーセントも消費した.

このプログラムをHP-16Cの実機で実行すると, 不思議なことに分解する数の大きさにあまり関係なく30秒ちょっと程度で計算できる. DM-16はさすがに現代の電卓だけあって, 2,3秒で計算する. ただ誰もがいうようにこの電卓のキーの押し心地は最低だ.

前回紹介したシミュレータは大きい数はかなりの時間がかかり, 実用にならぬ. しかし除数が徐々に増えていくのが窓に見えるので待つのも苦にならない.

やはり私の作った個人用電卓の組込み機能が最高だなぁ.

Christopher StracheyのGPM

Ackermann関数

再帰呼出しというとすぐに出てくるのがAckermann関数である.

(define (a x y)
 (cond ((= x 0) (+ y 1))
       ((= y 0) (a (- x 1) 1))
       (else (a (- x 1)
                (a x (- y 1))))))

SICPの定義は少し違う. (演習問題1.10)

(define (a x y)
 (cond ((= y 0) 0)
       ((= x 0) (* 2 y))
       ((= y 1) 2)
       (else (a (- x 1)
                (a x (- y 1))))))

SICPの著者のSussman君はもちろんこの方が正しいというが, ここでは巷間に伝わる前の定義に従う.

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

1桁の世界ではこの辺までしか計算できないから, あっという間に終るのは有難い. 特にコメントすることもあるまい.

Catalan三角形

TAOCPの式7.2.1.6-(22)にCatalan数とかCatalan三角形というのがある. このブログに「入れ子のかっこ」という話を書いたが, 4組の(,)の並べ方が14通りあり, それがCatalan数Cを使いC(4)=14と表現出来るという話題だ. 三角形の方は



のような形で, Catalan数C_nは図のC_nnである.

C₀₀=1,
C_p,q=C_p(q-1)+C_(p-1)q 0≤p≤q≠0の時
      =0 p<0 または p>qの時

だからSchemeでは

(define (c p q)
 (cond ((and (= p 0) (= q 0)) 1)
       ((or (< p 0) (> p q)) 0)
       (else (+ (c p (- q 1)) (c (- p 1) q)))))

GPMでは

$def,c,<$~1~2,
 $def,~1~2,<$$|,$lt,>~1<,0;,$lt,>~2<,>~1<;;,
  $def,t,0;,
  $def,f,<$+,$c,>>~1<<,$1-,>>~2<<;;,$c,$1-,>>~1<<;,>>~2<<;;>;;>;,
 $def,00,1;;>;
$c,0,0;,$c,1,1;,$c,2,2;, => 1,1,2,
$c,0,3;,$c,1,3;,$c,2,3;,$c,3,3;, => 1,3,5,5,
$c,0,4;,$c,1,4;,$c,2,4; => 1,4,9

と書ける.

論理関数

andやorはどうなるだろうか. 四則演算も自分で定義したくらいだから, 論理演算も出来るであろう. and, or, notはSchemeで書くと

(define (and x y)
 (if x y #f))
(define (or x y)
  (if x #t y))
(define (not x)
 (if x #f #t))

つまり(and x y)はxが真なら結果はyの真偽による; xが偽ならもちろん偽という定義でよい. (or x y)はxが真なら即真, 偽ならyに依存する. (not x)も見ての通り. これをGPMのif then elseで実装する. 真偽値はltと時と同様, tとfとする.

$def,&,<$~1,$def,~1,f;,$def,t,~2;;>;
$def,|,<$~1,$def,~1,t;,$def,f,~2;;>;
$def,\,<$~1,$def,f,t;,$def,t,f;;>;

マクロ名には大体の文字が使えるから, andは&, orは|, notは\を使った.

$&,f,f;,$&,f,t;,$&,t,f;,$&,t,t; => f,f,f,t
$|,f,f;,$|,f,t;,$|,t,f;,$|,t,t; => f,t,t,t
$\,f;,$\,t; => t,f

となる. これを使って多数決maj(majority)と排他的論理和xorを定義する.

$def,maj,<$|,$|,$&,~1,~2;,
  $&,~2,~3;;,$&,~3,~1;;>;
$def,xor,<$|,$&,~1,$\,~2;;,
  $&,$\,~1;,~2;;>;
$maj,f,f,f;,$maj,f,f,t;, => f,f,
$maj,f,t,t;,$maj,t,t,t; => t,t
$xor,f,f;,$xor,f,t;,$xor,t,f;,$xor,t,t; =>f,t,t,f

という次第だ. 二進の加算器が得られた気分だ.

Hilbert曲線

GPMはもともとプログラムの前処理系だったから, これを使ってプログラムの一部を発生させてみたい.

PostScriptでHilbert曲線を描いてみよう. もちろんPostScriptには再帰呼出しの機能があるから, PostScriptだけで簡潔にプログラムできるわけだが, GPMを 使うのも一興である.

まずはHilbert曲線の描き方の復習から.



左下のpは左から進んで来て, +が示すように左折する(+は反時計方向に曲り, - は時計方向に曲る), 左中央のfを左回転したf+分だけ進み, -で右折, fと進み, -で右折, 右回転したf-と進み +で左折して右へ抜ける.

左上のqの描き方も同様に出来ている.

これらp, qから中央のp'を描くのも図の通りで, +で左折, 左回転したq+, 左回転したf+, 右折, p, f, pと進み, 右折, f-, q-と進んで左折する.

q'も同様.

次に大きいp'', q''は今のp, qをp', q' に変更するだけである.

これだけ分ればSchemeで実装できる.

(define (^ x)
 (modulo (+ x 1) 4))
(define (_ x)
 (modulo (+ x 3) 4))
(define (f d)
 (display (list-ref '(
 "l 0 rlineto "
 "0 l rlineto "
 "l neg 0 rlineto " 
 "0 l neg rlineto ") d)))
(define (p n d)
 (if (> n 0)
  (begin
   (q (- n 1) (^ d)) (f (^ d))
   (p (- n 1) d) (f d)
   (p (- n 1) d) (f (_ d))
   (q (- n 1) (_ d)))))
(define (q n d)
 (if (> n 0)
  (begin
   (p (- n 1) (_ d)) (f (_ d))
   (q (- n 1) d) (f d)
   (q (- n 1) d) (f (^ d))
   (p (- n 1) (^ d)))))

関数 ^ は4を法として1を足す, _ は4を法として1を引く.

(f d)はd向きのfを引くPostScriptの命令を返す. つまりd=0なら l 0 rlineto (x 方向へl ,y 方向へ0だけ移動,) d=1なら 0 l rlineto(x 方向へ0, y 方向へl だけ移動)のようになっている. negはその前の値の符号を変える.

(p n d)は大きさnのpを描く. q+, f+ ...だったから(q (- n 1) (^ d)), 1小さいqを(^ d)向きに描き, (f (^ d)), fを(^ d)向きに描き, ...と定義する.

(q n d)も同じだ.

これをGPMにしたのが次だ.

$def,1-,<$-1,0,1,2,3,4,5,6,7,8,
  $def,-1,<~>~1;;>;
$def,^,<$1,2,3,0,$def,1,<~>~1;;>;
$def,_,<$3,0,1,2,$def,3,<~>~1;;>;
$def,f,<$~1,
 $def,0,<l 0 rlineto >;,
 $def,1,<0 l rlineto >;,
 $def,2,<l neg 0 rlineto >;,
 $def,3,<0 l neg rlineto >;;>;
$def,p,<$~1,
 $def,~1,
 <$q,$1-,>~1<;,$^,>~2<;;$f,$^,>~2<;;
  $p,$1-,>~1<;,>~2<;$f,>~2<;
  $p,$1-,>~1<;,>~2<;$f,$_,>~2<;;
  $q,$1-,>~1<;,$_,>~2<;;>;,
 $def,0,;;>;
$def,q,<$~1,
 $def,~1,
 <$p,$1-,>~1<;,$_,>~2<;;$f,$_,>~2<;;
  $q,$1-,>~1<;,>~2<;$f,>~2<;
  $q,$1-,>~1<;,>~2<;$f,$^,>~2<;;
  $p,$1-,>~1<;,$^,>~2<;;>;,
 $def,0,;;>;
/l 20 def 100 100 moveto 
$p,3,0;
stroke

大部分はマクロの定義で, 最後の方 /l 20 def 100 100 moveto はこのまま出力にコピーされる. これはPostScriptのプログラムで, 変数lの値を20に設定する; 座標100 100 に移動する; という指令だ.

次の$p,3,0;でサイズ3, 進行方向0のHilbert曲線を描く命令

0 l rlineto l 0 rlineto 
0 l neg rlineto l 0 rlineto 
l 0 rlineto 0 l rlineto 
(この後同様な28行省略)
0 l neg rlineto 

が出力され, 最後に strokeの命令がコピーされる.

出来上がった図は



上に2つ並ぶ凹のような形に足のついたものの辺の数は15で, それが全部で4個あり, 左と上と右のそれらを繋ぐ辺が3こあるから, この図全体の辺は63個.

上の出力は28行の省略があるから全体で32行. 最後の行以外はrlinetoが2個ずつあるからrlineto, つまり辺の数は63個あるわけだ.

GPMの例題はまだ続く.

HP-16Cのプログラム技法

1980年代にHewlett Packardが出荷していたHP-16Cという電卓の名器があった. 下のような面構えである.



私は30有余年愛用していて, 今も毎日持ち歩いているが, 十進法の計算はiPhoneを使うことが多くなった. HP-16Cの最大の特徴はRPN, すなわちReverse Polish Notationで, (1+2)*(3+4)を1,2+3,4+*の順に入力する. これは大変理にかなっていて, 私が以前実装したHappy Hacking Calculatorも勿論そうなっている. ただHP-16Cはスタックが4段で, 手前から(HPのマニュアルでは下から)x,y,z,tしかない. つまり上の例ではまず1がxに入り, カンマで示したENTERキーでその1がyへ移動. 2がxに入り, +で和の3がxに出来る. 3を入力すると, 和の3がyへ移動. xが新入力の3になり, ENTERでyの3がzに移動, xの3がyに移動, 4の入力でxに4が入り, +でxが7. yが3になり, *でxに21が出来る.

もう1つの特徴が, 数値の表示が二進, 八進, 十進, 十六進に切り替わることだ. 内部表現はいうまでもなく二進である.

iPhoneで真似が出来ないのがHP-16Cのプログラム機能だ. 現今の標準から見れば計算速度も相当遲いし記憶容量も小さいので大したことは出来ないが, こういうプログラムもあるということは, ソフトウェア歴史上無視できない.

今回のブログでは, HP-16Cのプログラムのスタイルを話したい.

HP-16Cには203バイトのメモリーが塔載されている. 0番地の方からデータの記憶場所, 反対側からプログラムの記憶場所として使う. プログラムの方は1ステップが1バイトだが, データの方は1ビットから64ビットまで長さが任意で, プログラムに取られた場所以外を様々なサイズでデータ用の番地を取っていく. これをレジスタという. 私には計算に使うx, yなどがレジスタという気分なのだが, そちらはスタックというのだろう. この可変長のレジスタがどういう実装なのか私には分らない.

HP-16CのプログラミングはKeystroke Programmingといって, 通常手で計算するキーの押し方の順を記憶するもでのある. この他にラベルを置くこと, ラベルへジャンプすること, 条件判断で1ステップ飛ばすことなどが出来る.

またIで表記するインデックスレジスタが1個ある. (I)と表記するとIの内容という意味である.

以下で扱うHP-16Cのプログラムは, 素因数分解のプログラムで使う予定の疑似素数発生噐である. つまり2,3,5,7,11,13,...を次々を作るのだが, 素数だけを作るには大量に記憶するか, 毎回素数性を検査するとか, 実用にならないので, 途中から先は, たまに合成数があっても苦しうないというものだ.

2,3,5で割れない数は30個で循環して現れる. 0から29までの数を書き, 2,3,5で割れるものは1, 割れないものは0として, 縦に0が揃うものを*で表わし, その間隔を 調べると以下のようになる.

2 101010101010101010101010101010
3 100100100100100100100100100100
5 100001000010000100001000010000
&  *     *   * *   * *   *     *
   2     6   4 2   4 2   4     6

また, 始めの方の素数の間隔は

2,3,5,7,11,13,17,19,23,29,31,37,41,43,
 1 2 2 4  2  4  2  4  6  2  6  4  2

だから, 2に1,2,2と足して3,5,7が得られた後は, 上の循環数列の最初の4からを順に足す. こうするとそのうち49が現れるが驚かない.

ではプログラムに移ろう. 循環の最後を調べるのを容易にするには, インデックスレジスタIを使う. 1足したり引いたりした結果が0になると, 次のステップを 飛ばす機能があるから, レジスタには次のように逆順に差分を入れておく. インデックスは最初11に設定し, DSZ(decrement skip if zero)で1ずつ減らし, 0になったら8に再設定する.

0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11
   6, 2, 6, 4, 2, 4, 2, 4, 2,  2,  1

プログラムは不思議なことに001番地から始まる.

001 43.22.00 LBL 0
002 43.34    PSE
003 54.31    RCL (I)
004 40       +
005 43.23    DSZ
006 22.00    GTO 0
007 34       x<>y
008 44.32    STO I
009 34       x<>y
010 22.00    GTO 0

001行目はLBL 0(ラベル 0)だ. つまりここに0のラベルをつける. 002行目はPSE(pause)で一旦小休止してxの内容をしめす. すぐに走りだすから, 目を凝らす必要がある. 003行目 RCL(I) インデックスレジスタの内容をxにもってくる, xにあった疑似素数はyに送られている. 004行目 +. これで前の疑似素数に差分を足す. 結果はxに出来ている. 005行目 DSZ. インデックスレジスタから1を引き, 結果が0なら007行目へ飛ぶ. 006行目GTO 0. ラベル0へジャンプ. 007行目. xに出来ている疑似素数とyにある数を交換する. 008行目 STO I. xの値をインデックスレジスタへ入れる. 009行目. x<->y. 010行目GTO 0. 0へ飛ぶ. という風に走る.

従って,
11をインデックスレジスタに入れる.
8をxに入れる. ENTERしてyに送る.
最初の疑似素数2をxにいれる.
GTO 0でスタートへ.
R/Sで起動する.

HP-16Cが手元に無ければシミュレータが使える. 例えばここにある.

これらのキーにポインタを合せてクリックすると入力できる. たとえば最初の例. 1,2+3,4+*だが, 1を押し, 中央のENTERを押し, 2を押し, +を押し, 3を押し, ENTERを押し, +, *を押す.

HEXを押せば15に, OCTを押せば25に, BINを押せば10101に表示が変る. その時の基数を窓の右端にh,d,o,bで表示する.

各キーには上にオレンジ色, 下に青色で機能が示されている. その方を使うには, その前に左下のfかgを押しておく.

というわけで疑似素数発生プログラムをやってみよう.

まず差分を入れる. レジスタの長さは8ビットでよさそうだからDECとしてから8を押しfを押し, オレンジ色のWSIZEを押す.

HP-16Cではキーの位置を2桁(行と列)の十進数で表わす. 行番号は上から1,2,3,4. 列番号は左から1,2,3,...,9,0である. fの位置は42, WSIZEの位置は44だ. (2行にわたるENTERは36)

レジスタ1に6を入れるには6を押してxに置き, STO 1と押す. こうしてレジスタ9に2まで入れたら, 次に2と1を10と11に入れるためにHEXにする. 2 STO A 1 STO Bで差分の入力は終わる

プログラムを入力するには, gと31のキー, P/Rを押す. そうすると窓に000が表示される. いかにも000に入りそうだが, 表示の次の位置の入るので要注意だ. プログラム入力モードになっているので, 窓にPRGMの表示がある.

LBL 0をいれる. 窓には 001 43.22.00 と出る. 001番地にキー位置43 (つまりg), 22 (つまりLBL), 00 (これはキー位置ではなく0)が格納されたことが分かる.

このようにして上のプログラムリストの右端の部分を次々と入力する. それが済んだら, 入力モードのままで, GTO . 001と入力すると001番地の命令をみることができる. 以下続く番地の内容を見るにはSSTを押す. こうしてプログラムを確認できる.

プログラムを走らせるには, P/Rを押して入力モードから抜ける. つまり計算モードにする. DECであるのを確認し, 11 STO I としてインデックスレジスタを11に初期設定する. 次に8をENTERし, 2をxに入れて, GTO 0としR/Sを押すとプログラムが走りだす.

各々方はうまく走らせることが出来たであろうか. 私の場合, HP-16Cの実機ではうまく行くが, 件のシミュレータでは, R/Sで停止しない. ONを押して電源を切ると止まってくれた.

HP-16Cのプログラムはかようなものである.

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

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