復活祭公式

前回のブログの最後の方の式, つまり補正した日数から対応する月を求める式

(define (month d) (quotient (+ (* 12 d) 373) 367))

の係数についてである.

前回のブログでは, 各月mの最初の日dを決める式を, 閏月を平均的に挿入し,

d=floor((7*(m+10))/12-6)+30*(m-1)
=floor((7m+70-72)/12)+(360m-360)/12
=floor((367m-362)/12)

のようにして求めた. これを使うと

m = 1 2  3  4  5   6   7   8   9   10  11  12
d = 0 31 61 92 122 153 183 214 245 275 306 336

が得られ, gregorianからfixへの変換に使えた. 今回は逆に与えられたdからmを計算したい. 前回の最後の表

d = 0-30 31-60 61-91 92-121 122-152 153-182 183-213 214-244
m = 1    2     3     4      5       6       7       8

d = 245-274 275-305 306-335 336-366
m = 9       10      11      12

を見ると,

0 ≤ d < 31ならm=1,
31 ≤ d < 61ならm=2,
...

としたい. mからdを決める式を代入すると,

floor((367*1-362)/12) ≤ d < floor((367*2-362)/12) ならm=1,
floor((367*2-362)/12) ≤ d < floor((367*3-362)/12) ならm=2,
...

なので,

floor((367m-362)/12) ≤ d < floor((367*(m+1)-362)/12)

の関係が分かる.

等号を右に移動して ≤d< から <d<≤ の形にするには, 両端の式から1を引くか中央の式に1を足す. また下の図で考えると



floor(f(m)) ≤ d < floor(f(m+1))はdのような, 下が閉区間, 上が開区間の範囲である. 従って, d+1は f(m) < d+1 ≤ f(m+1) となり, floorがいらなくなる.

(367m-362)/12 < d+1 ≤ (367*(m+1)-362)/12

367m-362 < (d+1)*12 ≤ 367*(m+1)-362

367m < (d+1)*12+362 ≤ 367*(m+1)

(m+1)=ceiling((12*(d+1)+362)/367)

m=ceiling((12*(d+1)+362)/367)-1

(12*(d+1)+362)/367-1

=(12*d+12+362-367)/367

=(12*d+7)/367

ceilingをfloorにしたいから, 分子に366を足すと

m=floor((12*d+7+366)/367)
=floor((12*d+373)/367)

で前回の式になった.

復活祭公式

復活祭(Easter)の日取りは, 春分かその後の最初の満月の後の日曜で決るから, このアルゴリズムに関心を持つ人は少なくない.

Knuth先生のThe Art of Computer Programming(TAOCP)にも演習問題として使われている. (TAOCP Vol.1,3rd Ed.,ex1.3.2-14, Fasc1.1,ex1.3.2'-32). 私も興味を持っていて, 数セミの1982年9月号に書いたりしたので, それを知った杉本敏夫さんや新井正夫さんが資料を送ってくださったことがある.

春分とか満月とか天体の運行に関連して日が決るように思えるが, 実は春分は3月21日に固定するとか, 年初の月齢を使うとか, アルゴリズムで決っているのである. とはいえ, TAOCPの演習問題でも分かるように, 結構面倒である.

先日, 以前頂いた資料をのんびりと見ていたら, Gaussの公式というのがあった. これは意外と簡単であった.

(define (easter20 y)
 (let* ((m 24) (n 5) (a (modulo y 19)) (b (modulo y 4))
   (c (modulo y 7)) (d (modulo (+ (* 19 a) m) 30))
   (e (modulo (+ (* 2 b) (* 4 c) (* 6 d) n) 7)))
  (if (> (+ d e) 9) (list 'april (+ d e -9))
   (list 'march (+ 22 d e)))))

という具合いに剰余しかとらない. let*の始めの方のmとnは世紀に依存する定数で, 2000年から2099年までは24と5である. そのため関数名に20をつけた. 杉本さんの資料に詳しい説明があるが, その話はまたいつかにして, 今回はReingoldとDershowitzのCalendrical Calculationsにある復活祭公式の話だ.

簡単にいえば, 計算法はKnuthのアルゴリズムと同様である. ただ計算の記述が本書では一風変っているので, その辺の話をしたい.

この本では, 日の計算に, fixed day number (fixed dateとかR.D.ともいう) の通日を使う. ユリウス日のようなものだが, Gregorian暦が過去まで使われたとして, その1年1月1日を第1日とするのである. その日は月曜なので, fixed dateを7で除した剰余が日曜を0とする曜日と重なって都合が良い,

Gregorian暦は400年毎に繰り返すから, 2001年1月のカレンダーをcal 1 2001で出力してみると

    January 2001
 S  M Tu  W Th  F  S
    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

でGregorianの1年1月1日が月曜であることが確認出来る.

Calendrical Calculationのeasterのアルゴリズムはこうだ. Schemeに書き直してある.

(define (easter g-year)
 (let* ((century (+ (quotient g-year 100) 1))
   (shifted-epact (modulo (+ 14 (* 11 (modulo g-year 19))
     (- (quotient (* 3 century) 4))
     (quotient (+ 5 (* 8 century)) 25)) 30))
   (adjusted-epact
    (if
     (or (= shifted-epact 0)
      (and (= shifted-epact 1) (< 10 (modulo g-year 19))))
       (+ shifted-epact 1) shifted-epact))
   (paschal-moon (- (fixed-from-gregorian g-year april 19)
    adjusted-epact)))
  (kday-after paschal-moon sunday)))

最後の行から見ると, easterはpaschal-moon(春分満月)の後の日曜を探すようになっている. kdayというのは, kがsunやmonや...やsaturで, それらをまとめてkdayという. (k曜日ということだ.) つまり春分満月の後の日曜を求めるのだ.

その前のlet*の中は, century(世紀), epact(歳首月齢), paschal-moonなどを計算する. paschal-moonは(fixed-from-gregorian year month day)で, 引数の(Gregorian暦の)年月日のR.D.を計算し, adjusted-epactを引いたものが春分満月であり, その直後の日曜(のfixed date)を計算する.

2012年の復活祭を計算すると, 734601と得られる.

Gregorian暦からR.D.への変換は以下の通り. まず定数を定義する. (ここには今必要なものしか書かない.) 次にその年が閏年であるかみる関数gregorian-leap-year?. そして本体だ.

(define sunday 0) (define april 4)
(define january 1) (define march 3)
(define gregorian-epoch 1)

(define (gregorian-leap-year? g-year)
 (and (= (modulo g-year 4) 0)
  (not (member (modulo g-year 400) '(100 200 300)))))

(define (fixed-from-gregorian year month day)
 (+ gregorian-epoch -1 (* 365 (- year 1))
 (quotient (- year 1) 4)
 (- (quotient (- year 1) 100)) (quotient (- year 1) 400)
 (quotient (- (* 367 month) 362) 12)
 (cond ((<= month 2) 0)
       ((and (> month 2) (gregorian-leap-year? year)) -1)
       (else -2))
 day))

結構オーソドックスに書いてある. ここで面白いのは, 前月末日までの年内の総日数の計算である. 一応2月が30日あるとして計算する式を用意し, 3月以降を平年か閏年かで補正する.

(map (lambda (m) (quotient (- (* 367 m) 362) 12)) (a2b 1 13))

をやってみると次の上の段が得られる. 次の段は平年の正しい値で, 1,2月はよいが, 3月以降は平年は2, 閏年は1多いので, 補正項を計算する. この結構きわどい係数をどのように見つけたのか.

 (0 31 61 92 122 153 183 214 245 275 306 336)
 (0 31 59 90 120 151 181 212 243 373 304 334)
月1  2  3  4   5   6   7   8   9  10  11  12

Calendrical Calculationsにはごたごた書いてあるが, 要するにこう導出したらしい.

前述のように2月が30日あるとすると, 1年は367日. 12ヶ月に31日の閏月(大の月)を7回平均的に置くとすると, 次のような図が描ける.



これは勾配7/12の線のfloorをとったもので, 1段上がるときを閏月だと思って1年の月を当てはめると, 図のように大の月と小の月が対応する. すこし位相をずらして見ると

(map (lambda (m) (floor (- (* (/ 7. 12) (+ m 11)) 6))) (a2b 0 12))
=> (0. 1. 1. 2. 2. 3. 3. 4. 5. 5. 6. 6.)

が得られ, 上のリストの1桁目と一致する. これに毎月30日ずつ足せば, 前月までの日数dになるわけだ.

従って上の式のm=0のところを1月に対応させると,

d=floor((7*(m+10))/12-6)+30*(m-1)
=floor((7m+70-72)/12)+(360m-360)/12
=floor((367m-362)/12)

大の月が7, 8月, 12, 1月と並んでいるのも, 結構平均的であったと改めて気づく.

dateの後のk曜の計算は, kday-on-or-beforeを使う.

(define (kday-after date k)
 (kday-on-or-before (+ date 7) k))

(define (kday-on-or-before date k)
 (- date (day-of-week-from-fixed (- date k))))

(define (day-of-week-from-fixed date) (modulo date 7))

たとえば3月20日の前の水曜(=3)は20日から3を引いた日(17日)の曜日(木=4)を知り, 20日からその4を引いた16日が水曜という計算をする.

さて, 復活祭の日がR.D.で教えられても有難くないから, Gregorian暦に変換する. まずR.D.のある年を計算する.

(define (gregorian-year-from-fixed date)
 (let* ((d0 (- date gregorian-epoch))
  (n400 (quotient d0 146097)) (d1 (modulo d0 146097))
  (n100 (quotient d1 36524)) (d2 (modulo d1 36524))
  (n4 (quotient d2 1461)) (d3 (modulo d2 1461))
  (n1 (quotient d3 365))
  (year (+ (* 400 n400) (* 100 n100) (* 4 n4) n1)))
 (if (or (= n100 4) (= n1 4)) year (+ year 1))))

n400は400年のあった回数. d1は400年内の日数. n100は100年のあった回数. d2は100年内の日数. n4は4年の回数. d3は4年内での日数. n1は年の回数. するとdateの年は, yearのようになるが, 閏の調整を最後に行う.

年が判明すると, R.D.の年月日が計算できる.

(define (gregorian-from-fixed date)
 (let* ((year (gregorian-year-from-fixed date))
  (prior-days
   (- date (fixed-from-gregorian year january 1)))
  (correction
   (cond ((< date (fixed-from-gregorian year march 1)) 0)
    ((and (>= date (fixed-from-gregorian year march 1))
     (gregorian-leap-year? year)) 1)
    (else 2)))
  (month (quotient
   (+ (* 12 (+ prior-days correction)) 373) 367))
  (day (- date (fixed-from-gregorian year month 1) -1)))
 (list year month day)))

fixed-from-gregorianを何回も使うが, なんとかしたい.

ここでも年内の日数から月を見つける方法がこの本の独特な点である. correctionを別として, 年内日数dからmonthが得られる関数である. 係数が前の式と微妙に違うが, 似たように計算出来よう.

(define (month d) (quotient (+ (* 12 d) 373) 367))

でクリティカルなところを調べるとちゃんと計算出来ている.

0-30 31-60 61-91 92-121 122-152 153-182 183-213 214-244
1    2     3     4      5       6       7       8

245-274 275-305 306-335 336-366
9       10      11      12

やはり3月以降は正しい値と違うので, 補正を先にしてから計算に取り掛かる. 先ほどのR.D.の復活祭は

(gregorian-from-fixed 734601)=>(2012 4 8)

で4月8日と判明した.

再帰曲線

3次元のHilbert曲線. 前回で2次までは描けることが分かった. まぁしかし, 3次を描くとどうなりそうかにも興味があってやってみた. ご用とお急ぎの方のために, 先ず結果をお目にかけるとこのようになる.



中央少し下の赤い点が出発点, 一番上の緑の点が到着点である. 今回はこの絵の描き方を説明したい.

天下り的だが,

(define (z+y^ d x y z)
 (x+z^ (/ d 2) x y z) (y+z^ (/ d 2) (+ x d) y z)
 (y+z^ (/ d 2) (+ x d) (+ y d) z)
 (z+xv (/ d 2) x (+ y d) z) (z+xv (/ d 2) x (+ y d) (+ z d))
 (y-zv (/ d 2) (+ x d) (+ y d) (+ z d))
 (y-zv (/ d 2) (+ x d) y (+ z d)) (x-zv (/ d 2) x y (+ z d)))

は1辺の長さdのz+y^をx,y,zを起点にして描く起動関数である. 描くといっても, 下請け関数が節点の座標を出力するだけで, 実際に描くのはPostScriptにまかせる.

この関数は下請けを8回呼ぶ. 最初のx+z^の場合は, 下請けの起点は貰ってきたx,y,zである. 次のy+z^の場合は, xを(+ x d), つまりxの位置をdだけずらす. z+y^のGrayコードが000, 001, 011,...と始まったのを反映している. だからその次はyも(+ y d)として呼ぶ. このようにして8回呼ぶと終わる.

その下請けは, やはりそれぞれのGaryコードに基づき,

(define (x+z^ d x y z)
 (p x y z) (p x (+ y d) z) (p x (+ y d) (+ z d)) 
 (p x y (+ z d)) (p (+ x d) y (+ z d))
 (p (+ x d) (+ y d) (+ z d)) (p (+ x d) (+ y d) z)
 (p (+ x d) y z))

(define (y+z^ d x y z)
 (p x y z) (p (+ x d) y z) (p (+ x d) y (+ z d))
 (p x y (+ z d)) (p x (+ y d) (+ z d))
 (p (+ x d) (+ y d) (+ z d))
 (p (+ x d) (+ y d) z) (p x (+ y d) z))
...

のようにして, x+z^, y+z^, z+xv, y-zv, x-zvの5個用意する. (p x y z)は(display (list x y z))である. 座標の8個出力後の改行もpの担当である. (z+y^ 2 0 0 0)と起動関数を呼ぶと, 64個の座標

(0 0 0)(0 1 0)(0 1 1)(0 0 1)(1 0 1)(1 1 1)(1 1 0)(1 0 0)
(2 0 0)(3 0 0)(3 0 1)(2 0 1)(2 1 1)(3 1 1)(3 1 0)(2 1 0)
(2 2 0)(3 2 0)(3 2 1)(2 2 1)(2 3 1)(3 3 1)(3 3 0)(2 3 0)
(1 3 0)(1 2 0)(0 2 0)(0 3 0)(0 3 1)(0 2 1)(1 2 1)(1 3 1)
(1 3 2)(1 2 2)(0 2 2)(0 3 2)(0 3 3)(0 2 3)(1 2 3)(1 3 3)
(2 3 3)(3 3 3)(3 3 2)(2 3 2)(2 2 2)(3 2 2)(3 2 3)(2 2 3)
(2 1 3)(3 1 3)(3 1 2)(2 1 2)(2 0 2)(3 0 2)(3 0 3)(2 0 3)
(1 0 3)(1 1 3)(1 1 2)(1 0 2)(0 0 2)(0 1 2)(0 1 3)(0 0 3)

が得られ, この座標に適当な倍率を掛けて, PostScriptで描くと


この絵は2次だが, さて, 3次まで描くには, 上にあった起動の関数と下請け関数を統合し,

x+y^ x+yv x+z^ x+zv x-y^ x-yv x-z^ x-zv
y+z^ y+zv y+x^ y+xv y-z^ y-zv y-x^ y-xv
z+x^ z+xv z+y^ z+yv z-x^ z-xv z-y^ z-yv

の24個について関数を書く必要がある. もちろん1個ずつ手書きするわけにはいかない. 当然関数発生関数を書くことになる. それを使い,

(3dgen 2 0 1 0) =>
(define (z+y^ d x y z) (if (= d 1) (begin (p x y z)
(p (+ x d) y z) (p (+ x d) (+ y d) z) (p x (+ y d) z)
(p x (+ y d) (+ z d)) (p (+ x d) (+ y d) (+ z d))
(p (+ x d) y (+ z d)) (p x y (+ z d)))
(begin (x+z^ (/ d 2) x y z) (y+z^ (/ d 2) (+ x d) y z)
(y+z^ (/ d 2) (+ x d) (+ y d) z) (z+xv (/ d 2) x (+ y d) z)
(z+xv (/ d 2) x (+ y d) (+ z d))
(y-zv (/ d 2) (+ x d) (+ y d) (+ z d))
(y-zv (/ d 2) (+ x d) y (+ z d)) (x-zv (/ d 2) x y (+ z d)))))

のような関数を24個揃える. そして(z+y^ 4 0 0 0)と起動すると, 512個の座標

(0 0 0)(0 0 1)(1 0 1)(1 0 0)(1 1 0)(1 1 1)(0 1 1)(0 1 0)
(0 2 0)(0 3 0)(1 3 0)(1 2 0)(1 2 1)(1 3 1)(0 3 1)(0 2 1)
(0 2 2)(0 3 2)(1 3 2)(1 2 2)(1 2 3)(1 3 3)(0 3 3)(0 2 3)
(0 1 3)(0 1 2)(0 0 2)(0 0 3)(1 0 3)(1 0 2)(1 1 2)(1 1 3)
(2 1 3)(2 1 2)(2 0 2)(2 0 3)(3 0 3)(3 0 2)(3 1 2)(3 1 3)
(3 2 3)(3 3 3)(2 3 3)(2 2 3)(2 2 2)(2 3 2)(3 3 2)(3 2 2)
(3 2 1)(3 3 1)(2 3 1)(2 2 1)(2 2 0)(2 3 0)(3 3 0)(3 2 0)
(3 1 0)(3 1 1)(2 1 1)(2 1 0)(2 0 0)(2 0 1)(3 0 1)(3 0 0)

(4 0 0)(4 0 1)(4 1 1)(4 1 0)(5 1 0)(5 1 1)(5 0 1)(5 0 0)
(6 0 0)(7 0 0)(7 1 0)(6 1 0)(6 1 1)(7 1 1)(7 0 1)(6 0 1)
(6 0 2)(7 0 2)(7 1 2)(6 1 2)(6 1 3)(7 1 3)(7 0 3)(6 0 3)
(5 0 3)(5 0 2)(4 0 2)(4 0 3)(4 1 3)(4 1 2)(5 1 2)(5 1 3)
(5 2 3)(5 2 2)(4 2 2)(4 2 3)(4 3 3)(4 3 2)(5 3 2)(5 3 3)
(6 3 3)(7 3 3)(7 2 3)(6 2 3)(6 2 2)(7 2 2)(7 3 2)(6 3 2)
(6 3 1)(7 3 1)(7 2 1)(6 2 1)(6 2 0)(7 2 0)(7 3 0)(6 3 0)
(5 3 0)(5 3 1)(5 2 1)(5 2 0)(4 2 0)(4 2 1)(4 3 1)(4 3 0)

(4 4 0)(4 4 1)(4 5 1)(4 5 0)(5 5 0)(5 5 1)(5 4 1)(5 4 0)
(6 4 0)(7 4 0)(7 5 0)(6 5 0)(6 5 1)(7 5 1)(7 4 1)(6 4 1)
(6 4 2)(7 4 2)(7 5 2)(6 5 2)(6 5 3)(7 5 3)(7 4 3)(6 4 3)
(5 4 3)(5 4 2)(4 4 2)(4 4 3)(4 5 3)(4 5 2)(5 5 2)(5 5 3)
(5 6 3)(5 6 2)(4 6 2)(4 6 3)(4 7 3)(4 7 2)(5 7 2)(5 7 3)
(6 7 3)(7 7 3)(7 6 3)(6 6 3)(6 6 2)(7 6 2)(7 7 2)(6 7 2)
(6 7 1)(7 7 1)(7 6 1)(6 6 1)(6 6 0)(7 6 0)(7 7 0)(6 7 0)
(5 7 0)(5 7 1)(5 6 1)(5 6 0)(4 6 0)(4 6 1)(4 7 1)(4 7 0)

(3 7 0)(2 7 0)(2 7 1)(3 7 1)(3 6 1)(2 6 1)(2 6 0)(3 6 0)
(3 5 0)(3 4 0)(3 4 1)(3 5 1)(2 5 1)(2 4 1)(2 4 0)(2 5 0)
(1 5 0)(1 4 0)(1 4 1)(1 5 1)(0 5 1)(0 4 1)(0 4 0)(0 5 0)
(0 6 0)(1 6 0)(1 7 0)(0 7 0)(0 7 1)(1 7 1)(1 6 1)(0 6 1)
(0 6 2)(1 6 2)(1 7 2)(0 7 2)(0 7 3)(1 7 3)(1 6 3)(0 6 3)
(0 5 3)(0 4 3)(0 4 2)(0 5 2)(1 5 2)(1 4 2)(1 4 3)(1 5 3)
(2 5 3)(2 4 3)(2 4 2)(2 5 2)(3 5 2)(3 4 2)(3 4 3)(3 5 3)
(3 6 3)(2 6 3)(2 6 2)(3 6 2)(3 7 2)(2 7 2)(2 7 3)(3 7 3)

(3 7 4)(2 7 4)(2 7 5)(3 7 5)(3 6 5)(2 6 5)(2 6 4)(3 6 4)
(3 5 4)(3 4 4)(3 4 5)(3 5 5)(2 5 5)(2 4 5)(2 4 4)(2 5 4)
(1 5 4)(1 4 4)(1 4 5)(1 5 5)(0 5 5)(0 4 5)(0 4 4)(0 5 4)
(0 6 4)(1 6 4)(1 7 4)(0 7 4)(0 7 5)(1 7 5)(1 6 5)(0 6 5)
(0 6 6)(1 6 6)(1 7 6)(0 7 6)(0 7 7)(1 7 7)(1 6 7)(0 6 7)
(0 5 7)(0 4 7)(0 4 6)(0 5 6)(1 5 6)(1 4 6)(1 4 7)(1 5 7)
(2 5 7)(2 4 7)(2 4 6)(2 5 6)(3 5 6)(3 4 6)(3 4 7)(3 5 7)
(3 6 7)(2 6 7)(2 6 6)(3 6 6)(3 7 6)(2 7 6)(2 7 7)(3 7 7)

(4 7 7)(4 7 6)(4 6 6)(4 6 7)(5 6 7)(5 6 6)(5 7 6)(5 7 7)
(6 7 7)(7 7 7)(7 6 7)(6 6 7)(6 6 6)(7 6 6)(7 7 6)(6 7 6)
(6 7 5)(7 7 5)(7 6 5)(6 6 5)(6 6 4)(7 6 4)(7 7 4)(6 7 4)
(5 7 4)(5 7 5)(4 7 5)(4 7 4)(4 6 4)(4 6 5)(5 6 5)(5 6 4)
(5 5 4)(5 5 5)(4 5 5)(4 5 4)(4 4 4)(4 4 5)(5 4 5)(5 4 4)
(6 4 4)(7 4 4)(7 5 4)(6 5 4)(6 5 5)(7 5 5)(7 4 5)(6 4 5)
(6 4 6)(7 4 6)(7 5 6)(6 5 6)(6 5 7)(7 5 7)(7 4 7)(6 4 7)
(5 4 7)(5 4 6)(5 5 6)(5 5 7)(4 5 7)(4 5 6)(4 4 6)(4 4 7)

(4 3 7)(4 3 6)(4 2 6)(4 2 7)(5 2 7)(5 2 6)(5 3 6)(5 3 7)
(6 3 7)(7 3 7)(7 2 7)(6 2 7)(6 2 6)(7 2 6)(7 3 6)(6 3 6)
(6 3 5)(7 3 5)(7 2 5)(6 2 5)(6 2 4)(7 2 4)(7 3 4)(6 3 4)
(5 3 4)(5 3 5)(4 3 5)(4 3 4)(4 2 4)(4 2 5)(5 2 5)(5 2 4)
(5 1 4)(5 1 5)(4 1 5)(4 1 4)(4 0 4)(4 0 5)(5 0 5)(5 0 4)
(6 0 4)(7 0 4)(7 1 4)(6 1 4)(6 1 5)(7 1 5)(7 0 5)(6 0 5)
(6 0 6)(7 0 6)(7 1 6)(6 1 6)(6 1 7)(7 1 7)(7 0 7)(6 0 7)
(5 0 7)(5 0 6)(5 1 6)(5 1 7)(4 1 7)(4 1 6)(4 0 6)(4 0 7)

(3 0 7)(3 0 6)(2 0 6)(2 0 7)(2 1 7)(2 1 6)(3 1 6)(3 1 7)
(3 2 7)(3 3 7)(2 3 7)(2 2 7)(2 2 6)(2 3 6)(3 3 6)(3 2 6)
(3 2 5)(3 3 5)(2 3 5)(2 2 5)(2 2 4)(2 3 4)(3 3 4)(3 2 4)
(3 1 4)(3 1 5)(3 0 5)(3 0 4)(2 0 4)(2 0 5)(2 1 5)(2 1 4)
(1 1 4)(1 1 5)(1 0 5)(1 0 4)(0 0 4)(0 0 5)(0 1 5)(0 1 4)
(0 2 4)(0 3 4)(1 3 4)(1 2 4)(1 2 5)(1 3 5)(0 3 5)(0 2 5)
(0 2 6)(0 3 6)(1 3 6)(1 2 6)(1 2 7)(1 3 7)(0 3 7)(0 2 7)
(0 1 7)(0 1 6)(1 1 6)(1 1 7)(1 0 7)(1 0 6)(0 0 6)(0 0 7)

が得られ, 始めの図が描けることになった. 折角描いてはみたが, 期待したほどには美しくなかったのが残念だ.

私がHilbert曲線の話を書いたのは2009年6月19日のブログであった. そこの図で分かるように, 2次元の曲線では, 1次の曲線が縦向きに出発するなら, 2次は横向き, 3次はまた縦向き,...のように交互に出発する.

今回, 3次元のHilbert曲線を3次まで描いて分かったのは, 前回のブログの最初の図で, 赤線の1次の曲線が0からx軸方向へ出発し, 2次が0からy軸方向へ出発したのに対し, 今回の図では, 赤点から予想通りz方向へ出発していることだ. そういうものだったのだ.

上の64行にわたる512の座標は, 8行ずつ8つのブロックにして書いてある. ブロックを上から順に0,1,2,...,7ということにすると, ブロック0の座標はxもyもzも0から3である.つまり原点に近い1/8の空間を占めている. ブロック1はxが4から7なので, x方向だけが4ずれた空間にある. 最後のブロック7はzが4から7で, z方向へずれた空間にあるわけだ.

ブロック0をもう一度書いてみる.

(0 0 0)(0 0 1)(1 0 1)(1 0 0)(1 1 0)(1 1 1)(0 1 1)(0 1 0)
(0 2 0)(0 3 0)(1 3 0)(1 2 0)(1 2 1)(1 3 1)(0 3 1)(0 2 1)
(0 2 2)(0 3 2)(1 3 2)(1 2 2)(1 2 3)(1 3 3)(0 3 3)(0 2 3)
(0 1 3)(0 1 2)(0 0 2)(0 0 3)(1 0 3)(1 0 2)(1 1 2)(1 1 3)
(2 1 3)(2 1 2)(2 0 2)(2 0 3)(3 0 3)(3 0 2)(3 1 2)(3 1 3)
(3 2 3)(3 3 3)(2 3 3)(2 2 3)(2 2 2)(2 3 2)(3 3 2)(3 2 2)
(3 2 1)(3 3 1)(2 3 1)(2 2 1)(2 2 0)(2 3 0)(3 3 0)(3 2 0)
(3 1 0)(3 1 1)(2 1 1)(2 1 0)(2 0 0)(2 0 1)(3 0 1)(3 0 0)

赤く書いてあるのは, すべての座標が2で割り切れるもので, 各行の1つずる現れる. 各行が1つのGrayコードになっていて, その原点に近いものがそれである. 下の図の8つの赤点に対応する.

それを2で割ると

(0 0 0)(0 1 0)(0 1 1)(0 0 1)(1 0 1)(1 1 1)(1 1 0)(1 0 0)

となり, 2次の64個の座標の1行目と同じになるが, 当たり前だ.

再帰曲線

今回は3次元のHilbert曲線が話題である. この2月は思わぬ風邪引きで, それも歳をとったせいか, 全治に時間がかかった.

それで家に引き籠もっている折, 例の「ハッカーのたのしみ」を眺めていたら, 3次元Hilbert曲線の基本部品という図を見つけた. 確かに3次元でも出来そうであるが, さてどうやったら描けるのか. しばらくそれが気になっていた.

2次元だと, 右折, 左折などを繰り返すと何となくアルゴリズムが見つかるが, 3次元の場合, これをどういえばよいか. 試行錯誤の結果, 基本部品に名前をつけ, それを組合せると描けそうだというところまで来た. もっともまだ最後まで出来たわけではなく, 今後挫折するかもしれないのだが.

とりあえず, 1次と2次の図を描いてみたのが下だ.



赤の太い線が1次, 黒の細いのが2次である. 通過した順が分かるように, 節点に0から番号を付けてある. 1次のは向こうの下の方0から始まる. この座標を000としよう. それから図で左下へ進む. この方向をxとする. つまり次の点1の座標はzyxの順に書けば001である. 今度は右下へ進む. この方向をyとする. 次の点2の座標は011だ. この伝でいくと, 010, 110, 111, 101, 100のように進むことになる.

もう1度揃えて書くと

zyx
000
001
011
010
110
111
101
100

これは誰が見てもGrayコードである.

ところで, 2次を見ると, 1次の0の回りを0→1→...→7と同様にGrayコードで進む8点で囲む. 1の回りを8→9→...→15の8点で囲む. ... 7の回りを56→57→...→63の8点で囲む. この8点の形はよく見ると5種類ある. つまり1と2(B), 3と4(C), 5と6(D)の回りの8点はそれぞれ同じで, その他に0(A)と7(E)で5種類である. Grayコードの部品がたくさんあるが, 全体がGrayコードになっているのではない.

A    B    C    D    E
zyx  zyx  zyx  zyx  zyx
000  000  011  110  101
010  001  001  111  111
110  101  000  011  011
100  100  010  010  001
101  110  110  000  000
111  111  100  001  010
011  011  101  101  110
001  010  111  100  100
x+z^ y+z^ z+xv y-zv x-zv

さてこれに分かりやすい名前を付けたい. Aを見ると最上段と最下段ではxが0から1と増えている. zは0から1になり0に戻る. yはふらふらして元に戻るから無視すると, xは増え(x+)zは010(^)の凸型ということで, x+z^と呼ぶことにする. 他のBからEも同じ要領で命名してある. 上にあった1次の形は, この方式でいうとz+y^である.

そのそれぞれを描くと次のようになる. 左上は1次の形である.



この1次の形z+y^は2次で描くとA, B, B, C, C, D, D, Eで出来ていたのだ. つまり

z+y^ = x+z^ 2(y+z^) 2(z+xv) 2(y-zv) x-zv.
        A      B       C       D     E

x座標はAで増え, Cで2回減増し, Eで減るから変化なし.
y座標はBで増え, Dで減るからy^である.
z座標はAで増減し, Cで増え, Eで減増するからz+である. 全体でz+y^となるわけだ.

2次のA, B, C, D, Eがこのように分解出来るなら, 3次も描けることになる.



分解を上の図を見ながらもう一度考えてみよう.

最初の+か-を伴う座標を主軸といおう. 次の^かvの座標を副軸といおう. 最後に名称に出てこないのを従軸といおう.

• 3と4は主軸の+, -により+か-になる. +の場合, 0,1,2では主軸は^, 5,6,7ではvになり, -なら逆になる.

• 副軸は^なら1と2で+, 5と6で-, vなら逆になる.

• 従軸の初期値は, 主軸と副軸の初期値と偶数パリティになるようにとる. それが0なら0で+, 3と4でv, 7で+. 1なら逆になる.

さっそく挑戦してみる.

x+z^ = y+x^ 2(z+x^) 2(x+yv) 2(z-xv) y-xv.

そう分かるとプログラムが書ける.

(define (xyz a) (list-ref '(x y z) a))
(define (+- b) (list-ref '(+ -) b))
(define (ud b) (list-ref '(^ v) b))

(define (bar a0 b0 a1 b1)
 (let ((a2 (- 3 a0 a1)) (b2 (modulo (+ b0 b1) 2)))
  (list
   (list (xyz a0) (+- b0) (xyz a1) (ud b1))
   (list (xyz a2) (+- b2) (xyz a0) (ud b0));0
   (list (xyz a1) (+- b1) (xyz a0) (ud b0));1,2
   (list (xyz a0) (+- b0) (xyz a2) (ud (- 1 b2)));3,4
   (list (xyz a1) (+- (- 1 b1)) (xyz a0) (ud (- 1 b0)));5,6
   (list (xyz a2) (+- (- 1 b2)) (xyz a0) (ud (- 1 b0))))));7

(z + y ^)は(bar 2 0 1 0)と入力する. すると

(bar 2 0 1 0) =>
((z + y ^) (x + z ^) (y + z ^) (z + x v) (y - z v) (x - z v))

と得られる. リストの先頭は入力パターンである. 従って上の図と同じになる.

Aを分解すると,

(bar 0 0 2 0) =>
((x + z ^) (y + x ^) (z + x ^) (x + y v) (z - x v) (y - x v))

一方, Grayコードを得るには次のプログラムでよい.

(define (foo a0 b0 a1 b1)
 (let ((b '()) (a (list 0 0 0)) (a2 (- 3 a0 a1))
  (b2 (modulo (+ b0 b1) 2)))
  (list-set! a a0 b0) (list-set! a a1 b1) (list-set! a a2 b2)
  (set! b (cons (reverse a) b))
  (list-set! a a2 (- 1 (list-ref a a2)))
  (set! b (cons (reverse a) b))
  (list-set! a a1 (- 1 (list-ref a a1)))
  (set! b (cons (reverse a) b))
  (list-set! a a2 (- 1 (list-ref a a2)))
  (set! b (cons (reverse a) b))
  (list-set! a a0 (- 1 (list-ref a a0)))
  (set! b (cons (reverse a) b))
  (list-set! a a2 (- 1 (list-ref a a2)))
  (set! b (cons (reverse a) b))
  (list-set! a a1 (- 1 (list-ref a a1)))
  (set! b (cons (reverse a) b))
  (list-set! a a2 (- 1 (list-ref a a2)))
  (set! b (cons (reverse a) b))
  (reverse b)))

初期値を設定してから, 従軸, 副軸, 従軸, 主軸, 従軸, 副軸, 従軸の順に1と0を交換する. "z+y^"は(foo 2 0 1 0)と入力する.

(foo 2 0 1 0) =>
((0 0 0)(0 0 1)(0 1 1)(0 1 0)(1 1 0)(1 1 1)(1 0 1)(1 0 0))

と得られる.

Aを１次とみて, その2次の絵を描いてみると,

なんだかうまく行っているらしい. 春から縁起がいいわ.

素因数探し

前回の素因数探しのブログ(1月23日)に「Edlerの作ったカードのセットはどこかに残っていないかしら」と書いた後, 探したら, factor stencilの画像はインターネットに存在した. 駄目元で探してみるものだ.

これ, 紙のカード?と思うが, その画像をコピーさせて頂く.





2004年にWashington D.C.のSmithsonian National Museum of American Historyの展示で撮影されたらしい.

1枚目の写真は, カードというよりは鳩サブレーのような厚みが感じられる. 2枚目のには文字がプリントしてあるが, それが読めないのが残念.


さて一方, 田中君と試作しようとしたstencil cardには, 2月3日に鎌倉のFabLabで再度挑戦した.

今回はこういう図を用意した.

これはR=2のプリントパターンである.



つぎはR=-1のカット用パターンだ.



今回試作するR=-1と2について, それぞれプリントパターンをカット用パターンを作った. とはいってもプリントパターンは, カード左上の見出しの数値だけが違う.

それぞれのパターンの四隅にある鉤の手は, 位置合わせのトンボだ.

まず多少厚手の白紙の中央に, プリントパターンを印刷する. 見出しの数値が2つあるのは, 内側のが完成したカードに残るもの, 外側がカット用パターンとの対を示す, 操作に安全を期すためのものである.

パターンを印刷した紙を, カッターマットに両面テープで張り付ける. しっかり止めないと, 精密な作業は出来ない. この結果, craftroboのカッターの刃先はマットにまで食い込むのである. また, 紙の端からトンボまで距離などを一応計測し, カットのソフトウェアに教える. そしてマシンをスタート.

craftroboはセンサーでトンボを位置を読み取り, 続いてカードの孔を切り始める. 孔が小さいから目にも止らぬ速さだ. ただ, どういう順にカットしているのか不思議である. 端から順に行くかと思うとそうでもない. そうこうしているうちに孔を切り終わり, 周囲を一周して終了する. この間約10分であった.

やってみると, カット圧力の調整が微妙なことが判明した. すなわち, 圧力が強いと孔から浮いたチップがカード上に散らばり, カッターの移動を妨害するのである. 従って, カットはするが, 紙片がはずれぬ程度の強さにしなければならぬ. そして, カッターが停止してから, ピンセットで孔を押し, 孔を開けることになる.

また反対に圧力が弱いと, 孔を手で押し開ける時に, きれいに切れないことになる. という次第で, 微妙であった.

また, トンボがあるにもかかわらず, カットした孔が, プリントした数字にきちんと重ならず, すこしずれ気味であったり, この辺の調整は難しいという経験をした.

カッターマットから外した紙からカードの形を抜き取り, カードの裏側(数字の印刷してない側)から, 孔を押し開ける. カードに約400個所の孔があるから, これも一仕事であった. (もちろん自宅へ戻ってからの仕事)

やはり, 計算機に接続されているIBMカードパンチを探さなければ.

完成したものは, 以下のとおり.





2月6日からローマで開催されたIFIP WG2.1の会合で, こんなもの作ったと見せたら, やんやの喝采を受けた. この作業では, SFCの田中浩也君に大変お世話になった. ありがとう.

除算サブルーチン

パラメトロン計算機 PC-1 が完成し稼働し始めたのは, 1958年3月26日ということになっているが, この時は乗算も除算も出来なかった. その後, それらの回路を設計し, やがて乗算が付き, 夏頃には除算も出来るようになった.

イニシアルオーダーは十進二進変換もするが, その際の10倍に乗算を使わず, シフトしたり足したりしているのは, 当時, 乗算命令が未実装だったからである.

従って, 始めの頃は, 乗算や除算のサブルーチンが存在した. 私は数値計算には興味がないので, 除算は必要なく, 高橋先生が作られたはずの除算のサブルーチンがどうなっていたかは知らない. 先生のライブラリのファイルが手元に来たので, 除算サブルーチンを探したが, 見つからなかった.

ところで, PC-1がお手本にした, Cambridge大学のEDSACには, 最後まで除算はなく, サブルーチンを使っていた. そういえば, EDSACでは除算はどうやったか知りたくなり, 例のEDSACの本を見ると, 思いの他簡単であった.

a₀←被除数,
c₀←除数-1
と初期化し, 漸化式
a_k+1←-a_kc_k+a_k,
c_k+1←-c_k².
で繰り返すと, c_k=0の時, a_kが商になる.

これをそのままSchemeでプログラムにする.

(define (div n d) ;n 被除数 d 除数
 (define (dv a c)
  (display (list a c)) ;途中の様子を見る
  (if (< (abs c) 0.0001) a
   (dv (- a (* a c)) (- (* c c)))))
 (dv n (- d 1)))

PC-1もEDSACも当時の計算機は, -1≤n<1の固定小数点数しか扱わない. 利用者は, 計算がこの範囲になるよう, 絶えず位取りに注意する必要があった.

そこで, 上のサブルーチンを(div 0.4 0.8), (div 0.4 0.5)でやってみると,

> (div 0.4 0.8)
(.4 -.19999999999999996)
(.48 -.03999999999999998)
(.4992 -1.5999999999999983e-3)
(.49999871999999995 -2.5599999999999945e-6)
=> .49999871999999995

> (div 0.4 0.5)
(.4 -.5)
(.6000000000000001 -.25)
(.7500000000000001 -.0625)
(.7968750000000001 -.00390625)
(.7999877929687501 -.0000152587890625)
=> .7999877929687501

正しい商は, 0.5と0.8である. 目がちらちらするから, 途中のa_k, c_kを整数でシミュレートしてから同様な表にすると,

0.4/0.8
k a          c           delta
0 0.4        -0.2        0.1
1 0.48       -0.04       0.02
2 0.4992     -0.0016     0.0008
3 0.4999872  -0.00000256 0.0000128
4 0.49999999 -0.00000000 0.00000000
    99967232    00065536   00032768

0.4/0.5
k a          c           delta
0 0.4        -0.5        0.4
1 0.60       -0.25       0.2
2 0.7500     -0.0625     0.05
3 0.796875   -0.00390625 0.003125
4 0.79998779 -0.00001525 0.00001220
    296875      87890625   703125

左端のkはaとcの添字, 右端のdeltaは正確なn/dとa_kとの差である. 最後のk=4の行は長いので, 少数点以下8桁目で折り返してある.

こういうのの性質をしらべるには, 各回におけるinvariantを発見するのが常道である. ためつすがめつするまでもなく, この例ではcとdeltaには比例関係があることが分かる. つまり, 上の計算ではdelta/cは-1/2, 下のではdelta/cは-4/5, つまり除算の商の符号を変えたものであった.

すなわち, n/d=a_k-(n/d)*c_kということだ.

ここまで分かると, 次は帰納法で証明するだけである.

k=0の時

a₀-(n/d)*c₀=n-(n/d)*(d-1)=n/d.

k=kで

a_k-(n/d)*c_k=n/d

が成立しているとする. その時k+1の式は

a_k+1-(a/d)*c_k+1
=-a_kc_k+a_k-(n/d)*(-c_k²)
=-a_kc_k+a_k+(n/d)*c_k²
=c_k(a_k-(n/d)*c_k)+a_k
=a_k-(n/d)*c_k)=n/d.

やはり, n/dになる. 一方 c_kはどんどん小さくなっていくから, c_k=0の時, a_k→n/d.

そういう風にEDSACの除算サブルーチンは出来ていたのであった.

世界一周プログラム

1月15日は私の恩師, 高橋秀俊先生の誕生日で, OBが先生のお宅に集まった. 先生が他界されたのは, 1985年だから, それから27年が経った. 成人の日は1月の第2月曜に変ったから, 15日ではなくなったが, 成人の日には今でもOBが某所に集まっている.

今回, その会で思わぬものを渡された. 高橋先生が東大にいらした頃の, パラメトロン計算機 PC-1 のライブラリのファイルである.



先生は1975年に東大を定年退官され, 慶応義塾大学へ移られた. 東大の教授室を後片付けした筑波大学の長谷部君が, そのファイルを見つけ, 自分の部屋で保管していたが, その長谷部君も定年になり, 筑波大学の教授室を片付けた. そしてそのファイルを私に委ねてきたのであった.

PC-1のライブラリを整備していた頃, 紙の質はまだ悪く, 日焼けで茶色になり, 縁もぼろぼろになってきた. 私が東大工学部にいた時, このままでは持たないと思って全部コピーしなおし, 製本した. その折, もう1部を作って高橋先生に差し上げたら大変喜ばれ, 「これを見ると, PC-1をやっていた頃が昨日のように思い出されます」という, 葉書の礼状を下さった.

先生がPC-1のライブラリファイルを, 東大に残していかれたとは, この正月まで知らなかった.

こうして手元にめぐってきたライブラリを見ると, 私の書いたイニシアルオーダーの前に, 先生が書かれたイニシアルオーダーがあった. 確かに1958年3月26日稼働開始の直後は先生のイニシアルオーダーが使われていたが, 長いし遅いので, 連休中に私が書き直したのが, その後1964年まで, 長く使われた.

先生のファイルには, "self reproducing program"もあった. これは計算機のメモリーをテストするため, 短いプログラムを走らせると, そのプログラムが自分を, アドレス部を修正しながら, 少し先のメモリーにコピーし, コピーが終わると, 新しい部分が走るようなものである.

Illinois大学のIlliac Iで使っていたと聞いていて, 高橋先生も書いてみられたのだろう.

Illiacのは, アドレス部がそれぞれその場所になっている, 2つのプログラムがあり, その差から新しいコピーを作っていると聞いたが, PC-1のは違っていた. (Illiacではこのプログラムをleapfrog(カエル跳び)と呼んでいたらしい.)

先生のスケッチは途中らしく, 正しくなさそうなので, 実際に使われていたプログラムを紹介しよう.



やはり, PC-1の命令から説明しなければなるまい. 1語は18ビットで, そのうち左端の6ビットが命令コード, 上のプログラムで, 各命令の先頭の"o"や"p"や"j"が入る. その右の1ビットが, l/sビットで, オペランドが長語(36ビット)か, 短語(18ビット)を区別し, 長語なら, 命令に"l"をつけ, そのビットが1になる. 残りの11ビットがアドレス部だが, PC-1は512語しかなかったので, 9ビットだけが必要, 上の2ビットは無視する.



0番地は"o"命令, outputの意味で, アキュムレータの左端の6ビットをテレタイプコードとして送り出す. アドレス部は, 命令の置いてある番地とすることになっているから, "0"である.

1番地は"p 16". "p"はpositive loadで, 16番地の短語をアキュムレータに取り出す. 従ってアキュムレータは, 16番地の命令"jl 17"になる.

2番地"jl"は, jump命令で, "j"の場合, lは長短ではなく, この場合は無条件ジャンプである. 従って, 3番地へ飛ぶ.

3番地"x 8"は, アキュムレータのアドレス部だけ, 8番地のアドレス部へ入れる命令. したがって, 8番地は "p 17"になる. 8番地のアドレスにかっこがついているのは, このアドレスは変更されるという意味である.

4番地の"a"はadd, 5番地の"s"はsubtractで, アキュムレータの"p 17"に, 16番地の"jl 17"を足し, 0番地の"o 0"を引く. 命令コード部も変るが今は関係なく, 17に17を足し0を引くから, アドレス部は17増えて34になり, 6番地, 7番地の"x"命令で, 9番地と12番地のアドレス部に入れられる.

このように, 定数を足すのにも, アドレス部は絶えず変るので, 定数を持っているわけにはいかず, 差を使う必要があった.

8番地に来て, "p 17", つまり"jl 3"をアキュムレータに持ってきて, それを9番地の"t", transfer命令で34番地にいれる. 10, 11番地は4, 5番地と同様で, アキュムレータのアドレス部を17増やし, 20にしてそれを34番地のアドレス部に入れる. 最終的に34番地は"jl 20"になる.

13番地は"p 8"だから, アキュムレータは"p 17"になる. 14, 15番地は, 14を足し, 15を引くから, 1引くことになり, アドレス部は16になり, 17番地へ飛び, さらに3番地へ飛ぶと, 今度は8番地, 9番地, 12番地のアドレス部がそれぞれ1減り, 16番地の"jl 17"が"jl 34"になって33番地に入る.

これを繰り返していると, やがて1番地の"p 16"が18番地に"p 33"として入り, 次のサイクルでは, "o 0"が"o 17"になって7番地の"jl 3"に上書きされる. コピーサイクルが17番地に到達した時点からコピーされたプログラムが走り出し, 同じことを繰り返すのである.

制御が新しいコピーに移ると, そこに出力命令があるので, プリンターが1文字打出す. つまりプリンターが文字を打ち続けている限り, メモリーは正常であるわけだ. 512/17=30だから, 30回のコピーでメモリーの始めに戻る. PC-1では, これを世界とみて, このプログラムを世界一周と呼んだ. 80日世界一周という映画があったから, ○○秒世界一周ともいったが, それが何秒だったかは, 残念ながら思い出せない.