2015年3月28日土曜日

復活祭公式

今年の復活祭が間もなくだ. 今年はそれもその前の晩に月食があるという豪華版だ. とはいえ私は宗教には無関心で, 復活祭はこのブログに何度か書いたようにその日取りについてだけ興味がある.

復活祭が何かも知らぬ小学校の上級のころ, 大叔父が「今日は復活祭だよ」といい, 手元のSOD(Standard Oxford Dictionary)を開きながら「春分の後の満月の後の日曜」と不思議なことをいったのが興味を持ったきっかけであった.

長ずるに及び, 折がある度に復活祭の日程計算の資料を取り込んできた. ずいぶんいろいろなことが分かってきたが, 今回話題とするのは, 我が愛読書「Calendrical Calculations」の8章The Ecclesiastical Calendarsの扉のページにあるFinger calculation for the date of Easterについてである.


この図の周囲には1582年から1699年, 1700年から1799年, 1800年から1899年の黄金数に対する歳首月齢(epact)の表とドミノ文字(dominical letter)の表があり, もちろん何語か知れぬが(恐らくスペイン語)使い方の記述もある.

最初の表は下のようだ.


[1582]つまりグレゴリオ改暦の年は, 黄金数が6, ドミノ文字がGであったことが分かる. 他にもこういう年の表示が5カ所あるが, 解像度が悪く判読できない.

ドミノ文字は毎年と閏日に1文字ずつ前へ戻る. 従って表を右端から眺めるとまずD, その前年が閏年で3月以降がE, 1,2月がF, 更にその前年がGのようになる.

最後がBからGになるのは改暦のせいだ. 改暦で10日抜けたから, 普通なら1582年のGから1583年はFになるのだが, G,A,Bと3(10 mod 7)日分進み, Bになっている. ここをCに書き直すと, この表は28個の場所があるから, ドミノ文字の一周分になることが分かる.

さて, 1600年はこの表の右端で, 黄金数は(modulo 1600 19)→4だからそれに1足して5. よって歳首月齢は, 5の下の段を見て15. 1月1日の曜日は(modulo (+ 1 (jd 1600 1 1)) 7)→6 だから土曜で, B. 表では右から10個めのABと上下に並ぶのの下3月は下のAになる.

1600年の復活祭は別の資料を見ると4月2日であった.

ここで手の掌の図を見る. 指の内側にあるのが歳首月齢なので, 15を探すと人差し指の左縁のDとある内側に15が見える. そこから右に指の縁を辿り, 最初のドミノ文字Aを探し, その位置を覚えておく.

次に親指の左縁のDのところに22とあるが, これを3月22日として右へ23日, 24日,...と先ほど覚えた位置まで進むと, そこが4月2日になるのである.

指の形はどうでもいいらしいから, 指の縁を直線に延ばして描いたのが次である. 上の線の右端から下の線の左端へ繋がる.

これを見ると歳首月齢が一通り書いてあり(*や24,25だったり,XXVは今は無視), その右隣からドミノ文字が一週分書いてあるのが分かる.




歳首月齢が分かると春分満月も一意に決まる. それがこの線図に赤字で書いたものである. 1600年の例でいえば, 歳首月齢が15だと春分満月が3月29日であり, その直後の日曜は右にドミノ文字を探せばよいわけである.

という次第で, finger calculationとはいうものの, 要はこの対応表だけであった.

西暦年から黄金数は簡単に分かるが, ある年のドミノ文字は西暦年の下2桁を28で割った剰余で表を引くのであろうか. その辺は何語かで書いてあるし, 鮮明さも欠くので, 想像するしかないようだ.

2015年3月17日火曜日

Christopher StracheyのGPM

cal

calというのはある年(ある月)のカレンダーを出力してくれるものである.

cal 3 2015

と入力すると
     March 2015
Su Mo Tu We Th Fr Sa
 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
が出力される. これは文字ベースだが, pcalはPostScriptで出力が得られるものである. 私は自分で書いたpcalも持っている.

2007年の夏のプログラミング・シンポジウムで私は「いろいろなプログラミング言語によるcalのプログラミング」という題でMMIX, Haskell, PostScript, Scheme, Teco, Texでcalを書いた.

GPMを使っていみていると, GPMでもcalを書いてみたくなるのは人情である. で早速やってみた.

まずはSchemeでロジックを確認しよう. ある月のカレンダーを次のように出力するとして, パラメータa, b, cを次のものとする.



(define (pcal a b c)  ;0<=a,b<7, 3<=c<6
(define (ss) (display "   "))  ;空白3文字出力
(define (dd x y) (display x) (display y) (display " "));xy空白
(define (nn) (newline)) ;改行
(define (r n x y)  ;月末処理 最後の行 xyからn回出力
 (if (> n 0) (begin (dd x y)
  (if (= y 9) (r (- n 1) (+ x 1) 0) (r (- n 1) x (+ y 1))))))
(define (q m n x y)  ;m行目n曜日にxyを出力
 (if (= m 0) (r (- 7 b) x y) ;月末処理へ
  (begin (dd x y)
   (cond 
    ((and (= n 0) (= y 9)) (nn) (q (- m 1) 6 (+ x 1) 0))
    ((and (= n 0) (< y 9)) (nn) (q (- m 1) 6 x (+ y 1)))
    ((and (> n 0) (= y 9)) (q m (- n 1) (+ x 1) 0))
    ((and (> n 0) (< y 9)) (q m (- n 1) x (+ y 1)))))))
(define (p n)     ;月初処理 始めの空白日を出力
 (if (= n 0) (q c (- 6 a) 0 1)
  (begin (ss) (p (- n 1)))))
(p a)
)
(newline)
(pcal 3 2 4)@
これで上手くいったのでgpmへ書き換える. マクロ中の空白や改行に注意.
$def,pcal,<
$def,1+,<$~1,
 $def,~1,<$1,2,3,4,5,6,7,8,9,10,$def,1,<~>>~1<;;>;,
 $def, ,1;;>;
$def,1-,<$-1,0,1,2,3,4,5,6,7,8,$def,-1,<~>~1;;>;
$def,-,<$~2,$def,~2,<$-,$1-,>~1<;,$1-,>~2<;;>;,$def,0,~1;;>;
$def,nn,
;
$def,ss,   ;
$def,dd,<~1~2 >;
$def,r,<$~1,
 $def,~1,<$dd,>~2<,>~3<;$>~3<,
  $def,>~3<,<$r,$1-,>>~1<<;,>>~2<<,$1+,>>~3<<;;>;,
  $def,9,<$r,$1-,>>~1<<;,$1+,>>~2<<;,0;>;;>;,
 $def,0,;;>;
$def,q,<$~1,
 $def,~1,<$dd,>~3<,>~4<;$>~2~4<,
  $def,>~2~4<,<$q,>>~1<<,$1-,>>~2<<;,>>~3<<,$1+,>>~4<<;;>;,
  $def,>~2<9,<$q,>>~1<<,$1-,>>~2<<;,$1+,>>~3<<;,0;>;,
  $def,0>~4<,<$nn;$q,$1-,>>~1<<;,6,>>~3<<,$1+,>>~4<<;;>;,
  $def,09,<$nn;$q,$1-,>>~1<<;,6,$1+,>>~3<<;,0;>;;>;,
 $def,0,<$r,$-,7,>>~2<<;,>~3<,>~4<;>;;>;
$def,p,<$~1,
 $def,~1,<$ss;$p,$1-,>~1<;;>;,
 $def,0,<$q,>>~3<<,$-,6,>>~1<<;, ,1;>;;>;
$p,~1;
>;
$pcal,3,2,4;
本来のcal, pcalにするには年月をパラメータa,b,cに変換する必要がある. これには400など我々のGPM世界をはみ出す定数が現れるので, 取りあえずはSchemeで書いてある.
年月からpcalのパラメータを求めるプログラム
(define (ym2pcal y m)
 (let* ((p (lambda (n) (= (modulo y n) 0)))
  (q (lambda (n) (quotient y n)))
  (s `(0 3 ,(if (if (p 100) (p 400) (p 4))
   1 0) 3 2 3 2 3 3 2 3 2 3))
  (d (modulo (- (+ 1 y (q 4) (q -100)
   (q 400)) (apply + (list-tail s m))) 7))
  (z (+ (list-ref s m) 28)))
 (list d (modulo (- 42 (+ d z)) 7) 
  (quotient (- (+ d z) 1) 7))))
(map (lambda (m) (ym2pcal 2015 m)) (a2b 1 13)) =>
((4 0 4) (0 0 3) (0 4 4) (3 2 4) (5 6 5) (1 4 4)
 (3 1 4) (6 5 5) (2 3 4) (4 0 4) (0 5 4) (2 2 4))
$pcal,a,b,cがうまくいったので, 今回はこれで満足しよう.

2015年3月13日金曜日

O'Beirneのキューブ

O'Beirneがこれらのピースの形を得た道筋を考えてみた.

24個の単位キューブが2×3×4の空間に並んでいたとする. (下の図の一番上のI) それに図のようにaからxまでのラベルをつける. 左端の2×3の枠が最下段. 右端の枠が最上段で, 順に4段に積んである. それぞれの段には6個の単位キューブがあり, ラベルは図の通りとする.



Iの図の赤線のところで分割し, ずらして合体させると, 再下段では右上のbcfがadeの上に来るからJの図の左端のように変る. 他の段でも同様な図になる. こうしてJが得られる. JからKへの分割面は上から見た図ではうまく表現できないが, 下から2段目のhiが再下段の上へ降りてきて, Kの左のようになる.

周囲の方向の壁にも名前をつける. 上はtop, 下はbot, 水平方向は北, 東, 南, 西で, nor, eas, sou, wesとした.

そして外面と分割面に向う単位キューブを数え上げると次のようなリストが得られる.
(define wall '(
;I
(bot a b c d e f)
(top s t u v w x)
(nor a b c g h i m n o s t u e k q w)
(eas c f i l o r u x a e g k m q s w)
(sou d e f j k l p q r v w x b h n t)
(wes a d g j m p s v b f h l n r t x)

;J
(bot b c a f d e h i m r v w)
(top t u s x v w b c g l p q)
(nor b c h i n o t u p q g l)
(eas c f e i l k o r q u x w)
(sou d e j l p q v w m r h i)
(wex b a c h g j n m p t s v)

;K
(bot h i b c a f d e n m g j u x w q)
(top t u s x v w p q h b a d o r l k)
(nor h i n o t u)
(eas i c f e o r l k u x w q n m g j)
(sou d e j k p q)
(wes h b a d n m g j t s v p o r l k)

;L
(bot n h i m b c g a f j d e)
(top t u o s x r v w l p q k)
(nor n h i t u o m a e s w k)
(eas i c f e o r l k m a s w)
(sou j d e p q k n b f t x l)
(wes n m g j t s v p b f x l)

;M
(bot m n h i g a b c j d e f v w x r)
(top s t u o v w x r p q k l g a b c)
(nor m n h i s t u o g a b c p q k l)
(eas i c f o r l)
(sou j d e f p q k l m n h i v w x r)
(wes m g j s v p)

;N
(bot g a b c j d e f n q u x)
(top s t u o v w x r a d h k)
(nor g a b c m n h i s t u o)
(eas c f i l o r g j n q u x)
(sou j d e f p q k l v w x r)
(wes g j m p s v a d h k o r)
))
ちょっとしたプログラムで処理すると, それぞれのキューブで外面に向いていない方向が得られる.

((a sou) (b eas) (c sou wes) (d nor eas) (e top wes) (f top nor)
 (g sou) (h eas) (i top wes) (j top nor) (k bot) (l bot) (m top)
 (n top) (o bot sou) (p bot eas) (q wes) (r nor) (s bot sou)
 (t bot eas) (u sou wes) (v nor eas) (w wes) (x nor))
ここはキューブが離ればなれにならない向き, つまりキューブが繋がっている方向なので, Iの図に接続の様子を赤字で書き込んでみる. 矢印はそこで水平方向に 繋がっていることを示す. 丸の中の点とプラスは電流の記号と同じで, 丸に点は上向きに, 丸にプラスは下向きに繋がっていることを示す.



これを眺めながら, 前回のブログの6種のピースと対応させると
adek C+
bcfl A-
gjpq B+
hior B-
msvw A+
ntux C-
うまくいったなぁ. O'Beirneもこう考えたのかな.

2015年3月12日木曜日

O'Beirneのキューブ

ネットで何かを調べていたら, 思いがけずO'Beirneのキューブ(O'Beirne's Cube)というのに出くわした.

O'Beirneの正しい読み方はよくは分らぬが, とりあえず「オベアーン」ということにしよう. 上のウェブページで動き続けている形はスクリーンショットをとると下の6枚のようになっている. 上左(I) → 上右(J) → 中左(K) → 中右(L) → 下左(M) → 下右(N) → と移り変って左上のIへ戻る.



目をこらして図を眺め, それぞれのサイズを見ると(奥行き, 幅, 高さの単位長は6,4,3なので)
  奥行きx  幅y  高さz
I6×2=124×3=123×4=12
J6×3=184×2=83×4=12
K6×4=244×2=83×3=9
L6×4=244×3=123×2=6
M6×3=184×4=163×2=6
N6×2=124×4=163×3=9
で体積は当然だがすべて1728になる. 移り変るにつれて, 奥行き, 幅, 高さのうち2つだけが±1変化し, 残りが不変で, Gray Codeみたいである.

しっかり分りたかったので, まず模型を作ることにした.

6種類のピースで出来ているが基本的には3種類の立体が必要で, そのそれぞれの展開図を下の図のように描き, 30センチ×40センチの工作用紙(30円)を使って4個ずつこしらえた.





これらは元々6×4×3の直方体を6の方向へ2倍(12×4×3になる. 上の写真左, 以下X), 4の方向へ2倍(6×8×3になる. 写真中, Y), 3の方向へ2倍(6×4×6になる. 写真右, Z)したものである. この2個ずつからA+, B+, C+, C-, B-, A-の6個のピースを作る. 下の写真の前列が左からA+, B+, C+, 後列がC-, B-, A-である.



写真では分り難いから, 図で示すとこうなる.



これらを積み重ねた下の写真の左側は上からA+, B+, C+, 右側は上からC–, B–, A–である. これを左右から合体させると最初の写真のIが出来る.



右のを少し奥にずらしてから合体させたのがJである.

模型をいろいろこねくりまわして, 移り変わりを描くことができた. 黒線はピースの境界. ジグザグの赤線は分割し, 1段ずらして再度合体させると次の図形に進む線だ. 青線は逆方向に回転する分割線. 6個のピースのすべてが見えているのが意外であった. 黒いジグザグ線も見えるが, これはその面の中心を通っていないから, ずらして合体しても直方体にならぬ. そもそも断面が向う側まで届いていない.



とにかくどう動くかは分った. 実にうまく出来ていると思った. Tom O'Beirneがどうしてこれに辿りついたかが, 次なる疑問だ.

2015年3月6日金曜日

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

HP-16Cにあって私の個人用電卓にない機能のひとつが実数演算である. なにしろHP電卓設計チームのひとりがIEEE 754の代表のWillian Kahanなのだからあるのは当たり前か.

しかしHP-16Cの内部表現は56ビットということを別として, 私には分かっていることもまだ分からないこともある.

整数から実数に演算を切り換えると, y*2xなる値がxレジスタに入る. 反対に実数から整数に切り換えると, その時のxの値がy*2xなるようにxとyが設定される.

たとえば円周率3.141592...は

(* 1686629713. (expt 2 -29)) => 3.1415926534682512

だから, 1686629713 ENTER 29 CHS FLOAT とすると窓に 3.141592 00 が現れる.

ここでまたDECにもどすと, xが-30, (normalizeされた)yが 3373259425 となり

(* 3373259425. (expt 2 -30)) => 3.1415926525369287

実数演算でなにかやってみよう. 円周率の連分数展開などが面白そうだ.
3.141592 これから 整数部の3を引く.
1.415926 -01 1/xをとる
7.062513 00 7を引く
6.251333 -02 1/xをとる
1.599658 01 15を引く
9.965869 -01 1/xをとる
1.003424 00 1を引く
3.424719 -03 1/xをとる
2.919947 02 ...
という次第で連分数は3,7,15,1,291,...

したがってπの近似値はまず3, 次は3+1/7=22/7, 次は3+1/(7+1/15) = 333/106, 次は有名な113/355. 捨てたのが1/292という項だから精度は抜群である.

さてこれをプログラムに書き直そう. 1/xはそういう命令があるから簡単だ. 問題は逆数をとったあと, 整数部を引くことである. この電卓にはfloorとかroundの命令は見当たらない.

いろいろ考えた末, 最初の例にあるように, yにある整数に, 2のxにある整数乗を掛けたのが実数だから, yにある整数をxにある整数の絶対値だけ右シフトすれば整数部が得られることが分った.

ということで書いた連分数展開のプログラムが以下のである.
001 43.22.C  LBL C
002 44.0     STO 0 ;0に格納しておく
003 24       DEC   ;整数計算に
004 44.32    STO I ;羃数をIに於く
005 33       R↓  ;yを取り出す
006 43.22.0  LBL 0 ;ループバック地点
007 42.b     SR    ;右1ビットシフト
008 43.24    ISZ   ;インデックス-1 結果0?
009 22.0     GTO 0
010 31       R/S   ;一旦停止
011 0        0     ;羃を0にして
012 42.45.48 float.;実数演算に
013 45.0     RCL 0 ;元の数を取り出し
014 34       x<>y  ;x,yを交換して
015 30       -     ;引く
016 43.26    1/x   ;逆数をとる
017 22.C     GTO C ;次の桁の計算に戻る
FLOATにしてからレジスタに3.141592....を入力し, GSB Cとするとしばらくrunningした後3が現れて停止する. 値を読み取り, R/Sを押すとまた走りだす. 後は7, 1, 291と続く.

TAOCPの第2巻に連分数の例がある.
π=3+//7,15,1,292,1,1,1,2,....//
e=2+//1,2,1,1,4,1,1,6,...//
φ=1+//1,1,1,1,.....//

8/29=//3,1,1,1,2//
√8/29=//1,1,9,2,2,3,2,2,9,1,2,1,9...//
など
これらを計算してみると楽しいが, 始めの円周率の例で見たように, 292のところがそろそろ限界で, HP-16Cの桁数(56ビット)では破綻し始める.

2015年3月4日水曜日

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がこれだ.

2015年3月2日月曜日

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にしたから.)

2015年3月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=pc*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=22*33ということが分るのである.

これが便利だったので, 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秒で計算する. ただ誰もがいうようにこの電卓のキーの押し心地は最低だ.

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

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