Life Game

Life Gameはuniversalであるといわれる. つまり計算機を構築する部品がいろいろ考案され, それを組合せればよいという議論らしい. ただ, von Neumannの自己増殖機械みたいに, 縦何セル, 横何セル, 遺伝子のセル何個で出来るというような具体的な数値がないので, これで時間さえかければ解けるよね, といわれている程度にしか納得できない.

さて, そういう部品の一つに希薄銃(thin gun)というのがある. Life Gameで作る計算機の情報(信号)は, 2次元世界を飛び回わるグライダーだが, グライダー銃は30クロックごとにグライダー1個を発射し, これでは濃すぎる(多すぎる)ので, もっと間引いたグライダー列を作りたい. その機構が希薄銃である. 本当に出来るかなと思いやってみた.

次の図が希薄銃の様子を示す. これはProcessingのシミュレータで動いているもののある時点でのスナップショット(をPostScriptで書き直したもの)である.



赤い丸で数字を囲ったのが3個所あり, それぞれグライダー銃がグライダーを発射する場所だ. 左上のグライダー銃0は, 右下へグライダーを出し続ける. 図に見えるグライダーに先頭からa, b, c,...と標識を付ける.

(グライダー銃はもっと複雑な形だが, こんな大きなものをシミュレータに書き込むと, 広い面積が必要になるので, このシミュレーションでは, シミュレータがどのタイミングにどの位置でどの方向へグライダーを発射するかのシナリオに従ってグライダーを発生する.)

右端の中ほどに1番のグライダー銃があり, これは左上へグライダーを送る. こちらのグライダーをf, g, h,...とする.

これらグライダーの2列の間で, eというグライダーが右上に向っている. これはConwayのいうキックバックで, eは2つのグライダー列の間を右上と左下で往復する. eは間もなくfと接触し, eとfは消滅し, eと逆向きのグライダーが出現する.

そうして戻ってくるキックバックのグライダーは, タイミング的には0番のグライダー列のdと接触し, dを消し, また右上へキックバックする.

eが前回右上に向う前に, 左下へ向っていたグライダーはaの1機前のグライダーに衝突し, 右上へキックバックしていたのである. 0番のグライダー列は, したがって, a, b, cとキックバック地点を通過した後, dは消滅し, その後また3個通過し, 1個消滅する,...  を繰り返す.

左上へ進んでいる1番のグライダー列も, 4個に1個はキックバックで消えるが, 3個はキックバックされずに通過する. しかしキックバック後のグライダーは要らないから, Eと書いたイーターで吸い込む.

上の図には中央あたりに2番のグライダー銃があり, グライダーを左下へ発射している. そこへキックバックを回避してきた0番のグライダー列が近づく. この衝突は, 両者が消滅する位相であって右下へ通過するグライダーはないが, すでにキックバックで消えた0番グライダーに対応するホールでは, 2番のグライダーは衝突を免れ, 左下へ進む. 図の左下のグライダーiは, そのように抜けてきたものである. その後3機の2番のグライダーは, キックバックを免れた0番のグライダーのために消滅してホールになる. jはaの前のグライダーがホールだったため, 抜けてきたのである.

このようにして, キックバックで3/4になった0番のグライダー列により, 1/4に希薄になったグライダー列が発生出来る.

次は, キックバックの起こし方だ. グライダーの位置決めのため, ある方向に進むグライダーの4位相パターンのそれぞれに, 注目点を決めよう.



この右下へ移動するグライダーの図で, 左上, 右上, 左下, 右下が赤字のようにそれぞれ位相0,1,2,3で, 一番進行方向に近い場所を注目点とする. 位相0はいわゆるハッカーエンブレムの形で, 4クロック後には右へも下へも1セル分移動する.

Conwayの示すキックバック機構の図は以下の通り.



白黒を問わず, 丸が生きているセルである. 下の数字0から8は相対クロックだ. クロック0では上に左下へ進むグライダー, 下に右下へ進むグライダーが見える. 黒丸は次のクロックでも生き残るセル. 白丸は次には消えるセル, 点は今は死んでいるが, 次のクロックで生まれるセルだ.

この反応を経て, クロック6,7,8では右上へキックバックされるグライダーが形成され, 8が丁度0と180度逆転している. つまり, ターンで8クロックかかるわけだ.

グライダーの列は30クロック間隔なので, 向こうでもう一度キックバックされ, 30の倍数で戻ってくれば, 後続のグライダーとうまく衝突出来る. こちらのターンで8, 向こうで8かかると, 16クロック. 1歩進むのに4クロックかかるから, その距離往路と復路で8クロック. 16足す何倍かの8クロックの和で30の倍数になるのは,

30=16+8*1.75
60=16+8*5.5
90=16+8*9.25
120=16+8*13

だから, 最後のが解で13ステップの距離を往復し, 片側で120クロック毎にキックバックを起すことになる. 120はグライダー間隔30クロックも4倍なので, 4機に1機が犠牲になる.

この計算から分かるように, 240クロックでもグライダー間隔を28にすればキックバック出来る. Conwayの記述には, 1/Nに希釈出来るが, Nは4の倍数でなければならないとある.

先ほどの図で, グライダー0とキックバックしたとして, 次のグライダー1とのキックバックはどこだろうか. タイミング0の下のグライダーの注目点(下の3x3の正方形の右下)をx=0,y=0とすると, タイミング8の右上に出発しようとするグライダーの注目点はx=0,y=5である(yは上向きに測る.)

次の衝突のタイミング0は13セル間隔離れているから, x=13, y=18, つまりタイミング0の上のグライダーの注目点がそこである. 従って, その時点で1番のグライダーの注目点は, 0の下のグライダーと180度回転した形で, xが1少なく、yが4多い. x=12, y=22にグライダーがあれば良い. そのタイミングでグライダーを発射するように1番のグライダー銃を配置する. または, xとyを1ずつ遠ざけると, 発射の位相を4早める.

キックバックの図を描くプログラムを利用して, 2個のグライダー消滅の様子も描いてみよう.



この図のように, 左上から来るグライダー(下)と右上から来るの(上)とが, この位置で出会うと4クロック後に消滅する. 右上の2番から出続けるグライダーは, 左上のグライダーがあると消滅し, ないと飛び続けるから, これをNot回路という.

シミュレータを使うノウハウも結構たまったので, またいろいろ遊べそうだ.

PC-1シミュレータ

パラメトロン計算機 PC-1のシミュレータは何回も書いたことがある. 数年前にC言語とxlibで実装したのは, 計算機のインディケータの点滅も再現していて, 見ていると楽しい. それを公開したいと思っていたが, うまい方法が見つからず, 延び延びになっていた.

しかし, 情報処理学会コンピュータ博物館のウェブページで, 「PC-1 パラメトロン式計算機」を見ると, 「現在PC-1の現物は存在しないが, 和田英一によりノートPC上でPC-1シミュレータが実現され, 当時のプログラムが動作している.」と書いてあるので, 早く公開しなければと少しずつ公開用の実装を進めていた.

今回, やっとProcessingで実装した同様のシミュレータが完成したので, お目にかけたい. ただ, 私の新しい4CPUのMacBook Airでは順調に走るのだが, 古い方のMacBookAirでは, インディケータが高速には点滅せず, 臨場感がいまいちだ. また私のウェブページパラメトロン計算機 PC-1も参考にされたい.

さて, 今回のシミュレータである.

http://playground.iijlab.net/~ew/pc1sim/pc1sim.html

にアクセスすると, 下のような画面が現れる. これがシミュレータのインターフェースだ.



中央左よりの緑の小箱の列は, インディケータで, 最上段の左18ビットが命令レジスタ, その右11ビットがSCC(シーケンスコントロールカウンタ), その下で36個並ぶのは, 上からメモリーレジスタ, アキュムレータ, Rレジスタである. PC-1は1語が18ビットだが, 2語つなげて36ビットを長語として演算が出来た.

各レジスタの内容の0(白で)と1(オレンジ色で)が表示される.

情報処理学会のコンピュータ博物館にパラメトロン計算機PC-1のウェブページがあり, 後藤さんと高橋先生がPC-1の前にいる写真がある. そのPC-1の上の方にネオン管が並んでいるが, それがこのインディケータである.

インディケータの右で3行3列に並ぶ枠がスイッチである. 上の3個は左からクリアスタート(Clear Start), イニシアルスタート(Init Start), リスタート(Restart).  リスタートは中断したプログラムを再開する. PC-1では, 停止したとき, 次の命令が命令レジスタに読み出されているから, それを棄て, SCCの示す場所の命令からスタートするのがイニシアルスタートである. さらにSCCを0にした上でイニシアルスタートするのがクリアスタートである.

このシミュレータは先読みしないので, 中と右のスタートは同じである.

中段の3個は左からイニシアルロード(Init Load), フリーラン(Free Run), ストップ(Stop)で, フリーランは連続実行モードをオンにし, ストップはオフにする. オフの時のスタートスイッチは, 命令のステップ実行になる. 左端のイニシアルロードは, イニシアルオーダーのテープを読込む.

下段は左からイニシアルオーダーテープを用意する(Place R0); e1000桁のテープを用意する(Place e1000)スイッチである. 右端は, 他のデモテープを選択する(Place Tape)スイッチで, 現れたメニューから選ぶ.

テープ選択スイッチを押すと,



のように右上にメニューが出るので, 使いたいデモプログラムをクリックする.

デモテープは次の4つが用意してある.

a. factorize32 2^31-1, 2^31-3, ..., 2^31-19 を素因数分解する (5分55秒)
b. eratosthenes 篩を使い256から11327までの素数の表を作る (7分)
c. lucaslehmer 3から199までの素数pについてMpの素数性を調べる
d. e1000 自然対数の底eの値を1000桁計算する (6分15秒)

lucas lehmer以外のそれぞれの最後は私が最近ウェブでテストしたの終了までの時間である.

シミュレータを使うには, まず左下のLoad R0スイッチで, イニシアルオーダーR0 のテープを用意する. 次にその上のInit Loadスイッチで, R0を読込む. テープの6列分がメモリーレジスタに読込まれ, 長語(36ビット)として, SCCの表示する0番地, 2番地, ..., に読込まれ, 同時にアキュムレータの各ビットとXORをとってアキュムレータに戻す. 読み終わったとき, アキュムレータが0になるようにパリティビットがテープに入っている.

e1000桁は当時ポピュラーだったデモプログラムである. フリーランにして, クリアスタートでテープ読込みと計算が始まる.

このプログラムは460の階乗分の1までで計算するが, 4回分を一度にまとめて計算するので, インディケータの規則的な繰り返しが115回で計算は終了. その後は結果の二進十進変換になる. このプログラムはR0を壊さないから, 次のデモはR0の再読込みなしで始められる.

factorize31は 2³¹-1, 2³¹-3, 2³¹-5,...2³¹-19を素因数分解するデモプログラムである. このテープは入力サブルーチンのテープを読込んだ後, 停止するから, フリーランをクリックし, リスタートしなければならない. 結果は以下のようになる.

2147483647.=2147483647.
2147483645.=5.19.22605091.
2147483643.=3.715827881.
2147483641.=2699.795659.
2147483639.=7.17.18046081.
2147483637.=3.3.3..13.6118187
2147483635.=5.11.337.115861.
2147483633.=5843.367531.
2147483631.=3.137.263.19867.
2147483629.=2147483629.

最初と最後が素数であり, 特に最初はMersenne素数である. 2³¹の平方根 46340までの疑似素数で割ってみる. 割る値がRレジスタに見えている. 46340は八進法では132404なので, Rレジスタの中央より右が132000くらいになると素数でも計算は終了する.

このデモプログラムはR0を破壊するので, 次のデモの前にはR0を再読込みする必要がある.

eratosthenesの篩 256から11327までの素数を篩で探す. 篩には36ビット長語の32ビットを使い, 残りの4ビットには疑似素数の差を保持する.

最初3で篩い, 次に5で篩い, ..., 11327の平方根106で篩うまで篩えばよいのだが, このプログラムは, 疑似素数の差の循環時に終了テストをするから, 231(八進で347)になるまで篩う. 篩う数はRレジスタに見え, 篩う数が大きくなると, 篩うのに要する時間短くなるのが分かる.

篩終わると, 出力ルーチンが読込まれ, 結果の1315個の素数を持つ素数表257. 263. ... . 11321.  が印字される. このデモプログラムはR0を破壊しない.

Lucas Lehmerテスト ある素数pについて, M_p=2^p-1が素数の時, M_pをMersenne素数という. この素数性を調べるLucas Lehmerテストというのがある. s₀=0, s_n=s_n-1²-2(mod M_p)を計算し, s_p-2=0(mod M_p)ならM_pは素数である.

デモプログラムlucas lehmerはp=2,3,...,199の素数について, M_pの素数性をしらべ, 素数ならpを, 合成数ならcを印字する.

最後まで実行すると相当時間がかかるので, 適当なところで中断しよう. PC-1が活動していた頃, 分かっていたMersenne素数はM₃₂₁₇くらいまでであった.

デモプログラムはplayground.iijlab.net/~ew/pc1sim/{e1000, factorize31, eratosthenes, llt}で見られる.

これらのシミュレーションは50年前の実機より格段に速い. しかしネオン管の明滅は当時を彷彿とさせる. 最近のCPUはブラックボックスになって, なにがどうなっているか分からないが, 当時はこのようなインディケータを眺めて計算の進行状況を知り, 一喜一憂していたものだ. 出力も印刷電信機で, 大きな音をたて, 机をゆさゆさ揺らしながら, 文字を出力していた. あれから半世紀が過ぎた.

再帰曲線

東京大学新聞に2012年度後期日程入試問題が掲載されていた. その総合科目IIにフラクタルの絵があった. 問題はフラクタル次元を計算するものだが, 私としてはその絵を描くプログラムに関心があった.

まず図はこういうものだ.


左からF₁, F₂, F₃

左からG₁, G₂, G₃
問題に記述はこうだ.

はじめに, 1辺の長さ1の正方形をF₀とする. F₀を, 1辺の長さ1/3の9個の小正方形に分割し, 中央の小正方形を取り去ったものをF₁とする.  また, それぞれの小正方形を1辺の長さ1/3のユニットと呼ぶ. 次に, F₁ を構成する8個のユニットのそれぞれについて, これを9個の小正方形に分割し, そこから中央の小正方形を取り去ったものをF₂とする.  また, この分割で得られたそれぞれの小正方形を1辺の長さ1/9のユニットと呼ぶ. 以下, 同様の操作を繰り返して得られる図形をF₃, F₄, ...とする.

PostScriptのこのプログラムは次のようだ.

/n 5 def
/l 540 def
/xs [0 1 2 0 2 0 1 2] def
/ys [0 0 0 1 1 2 2 2] def
/a {n 0 eq {0 0 moveto l 0 rlineto 0 l rlineto l neg 0 rlineto
 closepath fill}
 {/n n 1 sub def
 gsave 1 3 div dup scale
 0 1 7 {/i exch def
 xs i get l mul ys i get l mul gsave translate a grestore} for
 grestore
 /n n 1 add def} ifelse} def
 40 40 translate a

まず全体の次数nを決める. この例では5. 1辺の長さlも決める. 次に8個の正方形を描くので, その順に左下のxとy座標のリストxs, ysを書く. そしてユニットを描く手続きaの定義だ.



n=0なら, 長さlの正方形を書く. そうでないならnを1減らし, スケールを1/3にする. dupはx座標とy座標の両方を1/3にするためのもの. gsaveとそれに見合うgrestoreは, スケールや原点を変えるとき, 以前の環境をスタックするものである. xsとysから各小正方形の原点の座標をとり, 図を各手続きaを呼ぶ.

こうしてn=5を描いたのがこれだ.



次のは少し手ごわい.

はじめに, 1辺の長さ1の正方形G₀を, 1辺の長さ1/2の正方形ユニット4個に分け, 右上のユニットについては, 左下の1辺1/4の正方形を取り去る. こうして得られた図形をG₁とする. 次に, 残った3個の, 欠陥のないユニットのそれぞれについて, 同様の操作を行なう. すなわち, それぞれを4個の1辺の長さ1/4のユニットに分け, 右上にユニットについては左下の1辺の長さ1/8の正方形を取り去る. 得られた図形G₂とする. 続いて, G₂に含まれる1辺の長さ1/4の, 欠陥のない各正方形ユニットについて, 同様の操作を行なう. 得られた図形をG₃とする. これを繰り返して得られる図形の列をG₁, G₂, G₃, ...とする.

目で見ると簡単だが, 欠陥のない正方形をどう判定するか.

私の考え方はこうだ.

各正方形を4区画に分け, 左下, 右下, 左上, 右上の順に下請けに渡すとする. G₀を描く手続きをaとすると, 1/2のサイズで, a,a,aと3回呼び, 最後に左下を取り除いた正方形を描く手続きbを呼ぶ.

aはnを1減らして, a,a,a,bと呼べば良い. bはnを1減らし, 左下を除いて, a,a,aと呼べば良い. そう考えて書いたのが, このプログラムだ. 上のプログラムが理解出来れば, こちらも分かるだろう.

/n 7 def
/l 540 def /l2 l 2 div def
/a {n 0 eq{0 0 moveto l 0 rlineto 0 l rlineto l neg 0 rlineto
 closepath fill}
 {/n n 1 sub def
 gsave
 0.5 dup scale
 gsave 0 0 translate a grestore
 gsave l 0 translate a grestore
 gsave 0 l translate a grestore
 gsave l l translate b grestore
 grestore
 /n n 1 add def} ifelse} def
/b {n 0 eq{
 l2 0 moveto l2 0 rlineto 0 l rlineto l neg 0 rlineto
 0 l2 neg rlineto l2 0 rlineto closepath fill}
 {/n n 1 sub def
 gsave 0.5 dup scale
 gsave l 0 translate a grestore
 gsave 0 l translate a grestore
 gsave l l translate a grestore
 grestore
 /n n 1 add def} ifelse} def
 40 40 translate a

この結果, n=7で描いたのが, 下の図である.



PostScripでこういう図を描くのは思ったより易しい.

太陽太陰暦

英語ではlunisolar calendarというらしいから, 日本語とは順が反対だが, 太陰太陽暦でもいいようだ.

要するに月齢に合わせて月を決めるが, ときどき閏月をおいて, 季節を合わせる. その規則はどこにでも見つけることができる.

• 新月の起きる日を, 新しい月の朔日(ついたち, 1日)とする.

• 24節気のうち, 太陽の黄経が30度の倍数になる中気に各月を対応させる. 中気の対応しない月を閏月とする.

昔の本を読むと, 春一月, 夏四月, 秋七月, 冬十月と書いてあるが, 旧暦では1, 2, 3月が春, 4, 5, 6月が夏, 7, 8, 9月が秋, 10, 11, 12月が冬である. その各月に中気が対応し, 中気のない月が閏と書いてあるから, 手元の天文年鑑を見ながら計算してみた.

新月と節気の日付と時刻を順に書いてみると,

新月 2月22日  7h35m       ---   a月1日(2月)
啓蟄 3月 5日 13h21m 345度
春分 3月20日 14h14m   0度
新月 3月22日 23h37m       ---   b月1日
清明 4月 4日 18h 6m  15度
穀雨 4月20日  1h12m  30度
新月 4月21日 16h18m       ---   c月1日
立夏 5月 5日 11h20m  45度
小満 5月21日  0h16m  60度
新月 5月21日  8h47m       ---   d月1日
芒種 6月 5日 15h26m  75度
新月 6月20日  0h 2m       ---   e月1日(5月)
夏至 6月21日  8h 9m  90度

aの月は春分があるから2月である. eの月は夏至があるから5月だ.




b,c,dが3月と4月とその閏月になる. bの月には黄経30度の穀雨があるから3月になる. cの月は小満があるから4月で, dの月には30度の倍数になる節気がないから, 閏4月だと思った.

ところが旧暦に関するウェブページを見てみると違うのである. 閏は4月21日からの3月の方であって, dの月が4月であった.

なぜかというと, 小満と5月の新月は, 小満の方が7時間半ほど早いが, 同じ日である. 新月が何時でも, その日の午前0時から5月になるのである. 従って小満はd月1日にあることになっているらしい.

天体の運行に基づくカレンダーはこういうところが微妙だ. 平朔とかカレンダー用のアルゴリズムをCalendrical Calculationsでもっと勉強する必要がある.

復活祭公式

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

(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行目と同じになるが, 当たり前だ.