時の不定時法

落語の「時そば」でお馴染の九つ, 八つ, ..., 四つという時の表現には, 現代人には理解しかねるものがある. なぜ減るのか, なぜ一つ, 二つはないのか, なぜ九つと八つの間が八つ半ではなく九つ半なのか.

それはまず置くとして, 九つが真昼(12時)と真夜中(0時)もよいとして, 六つが6時と18時ではないのも解せない. 明六つは文字通り夜明けであって, 日の出ではなく, 暮六つは夕暮れであって, 日の入りではない. 六つから六つまでは等間隔らしいがそれも定かではない. お江戸日本橋七つ立ちは確実に暗いうちの出立であった.

これを不定時法というらしい. 英語ではtemporal hourというらしい. まぁとにかく時計がないのだから, 明るくなった, 暗くなったに頼るしかなかったのも頷ける. 国立科学博物館には, この不定時法で動く和時計があったと記憶する. でも季節毎に人手で振り子の錘の位置を調整するらしい. メトロノームの原理と同じだ. セイコー時計資料館のHP参照.

ところで岩波国語辞典の裏見返しに, 気になる図がある. 不定時法の時と現代の時との関係を示す図である. こういう図を見ると, その描き方を考えるのが私の癖で, 今回もそういう話題である.



不定時法の図の説明はこうだ. 季節によって変動する不定時を, 約1ヶ月ごとの点で示す. この図は日出・日没基準ではなく, 一般に行われた薄明(明六つ・暮六つ)基準によって作った. 円の中心点から各点を通過する直線を引くことによって, 今の時刻と対比できる.

なるほどというので, 天文年鑑をとりだす. 東京における日出没時と東京における天文薄明継続時間の表がある. 下の表は左から太陽の黄経30°おきの日付, 薄明継続時間, 日出時刻, 日没時刻である.

[ 320 125 546 1753]
[ 421 131 502 1819]
[ 521 142 432 1844]
[ 622 149 426 1901]
[ 722 143 441 1854]
[ 823 132 506 1822]
[ 924 126 530 1736]
[1024 125 555 1655]
[1123 129 624 1630]
[1221 131 647 1631]
[ 121 129 648 1657]
[ 218 126 625 1725]

これを参考に, 簡単だとばかりに書いてみた.

まったく違う図になった. これはしたり. 昔の人は, 標準時とか平均太陽時とか均時差などの概念はないのだ. 視太陽時が基準である. すると自分で計算しなければならない.



ところで不思議なのは, 天文薄明継続時間である. 天文年鑑には「太陽の高度が-18度になる瞬時と日出時または日没時との間の時間」と書いてある. しかし例えば北緯35°の地点の夏至では, 太陽がもっとも低くなる(北中か?)高度は, -35+23.5=-11.5だから, -18°にならない. この継続時間は1時間25分から1時間40分くらいである. してみると, 高度というのは, 赤緯にそって18°かと考えたが, 1時間が15°なので, 18°は1時間12分にしかならず, 不可思議のままである.

とりあえず, 夏至の昼間(あるいは夜間)の長さを計算してみることにした. 図の上の円は地球を西側から見たもので, 左右の直径が北緯35°の地平線. 傾きの線は夏至, 二分(春分, 秋分), 冬至の太陽の軌跡である. 右上方に南中時に太陽がある. 色は太陽が地上にある部分を示す. 緑は夏至, 赤は二分, 青は冬至である.



右下の円は, 赤道を北から見た図で, これで昼間の時間を計算しようと思った. その結果, θは17.73°であり, それから夏至の昼間の時間, 14時間22分が得られた. 日出, 日没は, 太陽の上の辺が地平線にかかる時刻なのだが, 中心として概算すると, 夏至の昼/夜は, 14時間20分, 9時間40分なので, OKであろう.

一般のδ(太陽の赤緯)に対する真北からの偏角θは

(define (foo delta) (define lat 35)
 (let* ((a (sin (d2r delta))) (b (* a (tan (d2r lat)))))
  (r2d (asin (/ b (cos (d2r delta)))))))

で計算できる.

太陽の黄経が30°,60°のときのθも知りたい. これに対応する赤緯は, 11.5° 20.0°くらいなので, 8.19°, 14.76°であった. 薄明時間を赤緯方向に20°として図を描くと以下のようになった.



最初の図によく似ているではないか.

後記
青木信仰「時と暦」によると, 明け六つは日出前定時法の二刻半, 暮れ六つは日入後二刻半(一刻は一昼夜の1/100. 二刻半は36分). 不定時法から定時法になったのは明治5年11月末, 太陽暦の採用と同時. 懐中時計の輸入普及が原因と書いてある.
Edward M. ReingoldとNachum Dershaowitzの"Calendrical Calculations"にはThe ancient Egyptians --- as well as the Greeks and Romans in classical times --- divided the day and night separately into 12 equal "hours" each. Because, except the equator, the length daylight and nighttime varies with the seasons, the length of such daytime and nighttime hours also vary with the season. These seasonally varying temporal (or seasonal) hours are still used for ritual purposes among Jews. の記述あり.

京王線調布駅

京王線の調布駅の列車運用は実に見事である. 普段は下車すると直ぐに改札を出るし, あまり何も考えずにアナウンスに従って乗車するが, 時刻表を調べてみた.

朝晩はややこしいが, 平日の昼間の20分間のダイヤは以下のようであった. つまり20分のネットダイヤである.

京王線下り        上り  相模原線下り        上り
     着 発       着 発         着 発       着 発
準特 04 05  各停 03 03    快速 03 06  快速 07 09
各停 07 08  準特 07 08    急行 13 15  各停 13
特急 14 15  各停 13 13    各停    19  急行 17 19
各停 17 18  特急 17 18

図に示すとこうなる. 特急, 準特急, 急行, 快速は時刻表の表示と同じである. また破線は相模原線を示す. 調布の上り方は京王線だが, 途中で実線に変えるのは面倒なので, 破線のままだ.



これで分かるように, 相模原線の上り列車は7分, 13分, 17分に京王線の上り列車の4番線への到着と同時に, 平面交差して3番線に到着する. 4番線から8分, 18分に準特急, 特急が出発した1分後に3番線から快速, 急行が出発する. 13分には各停が京王線, 相模原線から同時着し, 京王線が出発した後, 相模原線から来た各停は, 本線を布田方へ引き上げ折り返しの準備をする. 下りに注目すると, 京王線は17分に各停が到着, 18分に出発するが, その直前に折り返しの車が1番線に戻ってきていて, 19分に出発する. ここが最も神業のところだ. 急いで折り返せるように, 運転士がもう1人乗っている. 後は3分に快速が到着, 4分に準特急が到着, 5分に準特急が出発した後, 6分に快速が相模原線へ出発. 13分に急行が到着, 14分に特急が到着, 15分に同時発する. これが何事もないように20分毎に繰り返される. なお, 相模原線内で区間運転の各停には, 都営の車も使われている.

周知のように, 調布付近は連続立体交差の工事が進行中で, 2012年度に工事が完成すると, この光景は見られなくなる.

参考: 調布駅構内の配線は以下の通り.

ビットスワップ

Knuth先生のTAOCPは多くのことを盛り込みたいが故に, 説明は最小限であって, 手を動かし, プログラムを書き, 考えないとなんのこっちゃである. 7.1.3項演習問題52, 53もその類いである.

まずδswap. 下の図のxはビット列で, そのブロックAとC, BとDを交換し, x''としたい. 交換するブロックの長さはそれぞれ等しく, ブロック間の距離δはすべて等しい. その方法は, まずxを右にδ桁シフトしxと排他和をとり, 右側のブロックの形111...1でマスクし, yとする. ABA⊕Bのこと. xとyの排他和をx'とする. yを左にδ桁シフトしx'と排他和をとりx''とする. これをδswapという.



TAOCPではこの話題の前に, ビット列のi番目とj番目を交換したい. 本書の解を見る前に自分の方法を考えよとある. 私の考えはこれと同じであった.

この方法で64ビットのビット列を反転するには, 左右32ビットずつをそれぞれ反転し, 32ビットのブロックを交換する. それには0³²1³²のマスクがいる(0が32個並び1が32個並ぶパターンをこう書く). となり同士のビット交換には, ...010101のマスクがいる.
2^kビットの0と2^kビットの1で出来た(つまり...0^2k1^2k)長さ2^dビットのマスクをμ_d,kと書く. 上の32ビットごとのマスクはμ_6,5であり, ...010101はμ_6,0である. ...010101は...111111を3で割ればよい. ...00110011は...11111111を5で割る. つまりμ_d,k=2^2d-1/2^2k+1である. これはパラメトロン計算機の頃からの常識である. この応用で1/3は二進法では0.010101..., 十進法の0.1は二進法では0.0001100110011...なのが分かる.

さて64ビットのビット列の任意の置換は, 適切なマスクθ_k, θ'_kを用い,
k=0,1,...,5について x←θ_kによるxの2^kswap.
k=4,3,...,0について x←θ'_kによるxの2^kswap.
で出来ると書いてある. 64ビットの反転はこの前半だけでよい.

この応用が演習問題52である. 次の各場合の置換のためのθ_k, θ'_kを求めるのだ. もとの各ビット位置をj=(j₅j₄...j₀)₂とする.

a)完全シャッフル
jπ=(j₀j₅j₄...j₁)₂
b)8×8ビット行列の転置
jπ=(j₂j₁j₀j₅j₄j₃)₂
c)4×16ビット行列の転置
jπ=(j₁j₀j₅j₄j₃j₂)₂
d)FFT(Fast Fourier Transform)に出てくるパターン
jπ=(j₀j₁j₂j₃j₄j₅)₂

例を示すのに64ビットに名前をつける．RFC1341のBase64 Alphabetを用い

とする．左端Aの位置が111111, 右端/の位置が000000, 中央fが100000, gが011111である. 従ってa)の場合はgのj₅の0が右端へ移動し111110に変わり，62になるからAの右に来る．そう考えると

にしたいのだ．

今は理解が先だから，まず解答を見てみる．

a) 0<=k<5についてθ_k=μ_6,k＆μ_6,5, θ'_k=μ_6,k＆(μ_6,k+1⊕μ_6,k-1);
θ₅=θ₄,θ_-1=0
と書いてある．

ではやってみよう．64文字はブログの画面では長すぎるので, 32文字のところで分割する. 各マスクの右はシフト数である．

一体どうなったか. まず32ビットシフトの直前までは, 右半分32ビットの反転である. これは簡単. そこで左, 右の32ビットの落ち着き先をみると, Aは63, fは1, /は0, gは62なので, 63-1(奇数のみ) 0-62(偶数のみ)である. 63-1 0-62 と書く. それぞれのブロックを中央で分けると, 落ち着き先は
63-33 31-1 0-30 32-62 [3 2 1 0]. 右の[と]の中はブロック番号. ブロック2と0を交換すると
63-31 32-62 0-30 31-1 になる. クイックソートと同じで, 今後中央を越えての移動はない.
またそれぞれを半分にする.
63-49 47-31 32-46 48-62 0-14 16-30 31-17 15-1 [7 6 5 4 3 2 1 0]
ブロック6と4, 3と1を交換. 左半分は右半分のほとんど鏡像だから, 右半分に注目.
31-17 16-30 0-14 15-1 になる. 半分にする.
31-25 23-17 16-22 24-30 0-6 8-14 15-9 7-1 [7 6 5 4 3 2 1 0]
6と4, 3と1を交換. 右半分に注目.
15-9 8-14 0-6 7-1
半分にする.
15-13 11-9 8-10 12-14 0-2 4-6 7-5 3-1 [7 6 5 4 3 2 1 0]
交換. 右半分を見ると
7-5 4-6 0-2 3-1
で, 文字では cd98/+ef. 次はdと8, /とeを交換 (2)
7 6 4 5 3 2 0 1 になる. c89de+/f
dと9, /とfを交換 (1) 7 6 5 4 3 2 1 0 になった.

とまぁこういうことをやったのである.

他の場合については別の機会に.

入れ子のかっこ

Gosperのハック

定数ステップで次を求めるかっこの足跡のアルゴリズムは, Gosperのハックにヒントがあるらしい. GosperのハックはMITのAIラボのHAKMEMのitem 175にあり, To get the next higher number with the same number of 1 bits のプログラムである. 元はPDP10の機械語だが, TAOCPのex7.1.3--20に通常の書き方がしてある. それをまたSchemeにすると

(define (gosper x)
 (let* ((u (iand x (- x)))
        (v (+ x u)))
  (+  v (irsh (quotient (ixor v x) u) 2))))

iand, ixor, irshは符号なし整数値を二進数とみてビット演算する. 途中結果も含めて計算の進行状況を見るとこうなる.

   x      u      v   xor x    / u     y
000111 000001 001000 001111 001111 001011
001011 000001 001100 000111 000111 001101
001101 000001 001110 000011 000011 001110
001110 000010 010000 011110 001111 010011
010011 000001 010100 000111 000111 010101
010101 000001 010110 000011 000011 010110
010110 000010 011000 001110 000111 011001
011001 000001 011010 000011 000011 011010
011010 000010 011100 000110 000011 011100
011100 000100 100000 111100 001111 100011
100011 000001 100100 000111 000111 100101
100101 000001 100110 000011 000011 100110
100110 000010 101000 001110 000111 101001
101001 000001 101010 000011 000011 101010
101010 000010 101100 000110 000011 101100
101100 000100 110000 011100 000111 110001

これをみて理由を考える. まずxがα01^a0^bであったとする. ただしa≥1, b≥0である. (1^aは1がa個並んでいること.) 上からの4行についていえば
x=000111, α=00, a=3, b=0
x=001011, α=001, a=2, b=0
x=001101, α=0011, a=1, b=0
x=001110, α=0, a=3, b=1
である. uは右端の1だから u=10^b. vはxの右端の1に1を足すから, a桁連続している1が0になり, 繰上げが出る. v=α10^a+b. xとvをxorするとαがキャンセルされて, 1^a+10^bになる. uで割ると0^bがなくなり, 2ビット右シフトすると, 1^a+1が1^a-1になり, 結局 α10^b+11^a-1が得られる.
xのαの右にあった0は1に変り, bが1増え, aが1減り, 場所が交代した.
結果的に1の数は変らない.

Gosparハックの逆関数もある.(TAOCPex7.1.3--21)

(define (igosper y)
 (let* ((t (+ y 1)) (u (ixor t y)) (v (iand t y)))
  (- v (quotient (iand v (- v)) (+ u 1)))))

次にように計算が進行する.

   y      t      u      v   and -v   u+1    x
111000 111001 000001 111000 001000 000010 110100
110100 110101 000001 110100 000100 000010 110010
110010 110011 000001 110010 000010 000010 110001
110001 110010 000011 110000 010000 000100 101100
101100 101101 000001 101100 000100 000010 101010
101010 101011 000001 101010 000010 000010 101001
101001 101010 000011 101000 001000 000100 100110
100110 100111 000001 100110 000010 000010 100101
100101 100110 000011 100100 000100 000100 100011
100011 100100 000111 100000 100000 001000 011100
011100 011101 000001 011100 000100 000010 011010
011010 011011 000001 011010 000010 000010 011001
011001 011010 000011 011000 001000 000100 010110
010110 010111 000001 010110 000010 000010 010101
010101 010110 000011 010100 000100 000100 010011
010011 010100 000111 010000 010000 001000 001110

方法としては, yの右端が0の時, 最も右の`10'を'01'にする. yの右端が1の塊の時, その左の0を挟んだ左の1から右を01と1の塊にし, その右に0を詰める.
これをifで分岐せずに実行するのが味噌であるが, 説明は場合を分ける.

yの右端が0の時, y=α10^b(b≥1), t=α10^b-11, u=1, v=y, v and -v=10^b, quotient (v and -v) (+ u 1)=10^b-1; x=α010^b-1.

yの右端が1の時, 1が全部右に寄ったら終りなので, 右の1の塊の左に0があり, その左に1があるとする. つまりy=α10^a1^b(a,b≥1)とする. t=α10^a-110^b, u=1^b+1, v=α10^a+b, v&-v=10^a+b, u+1=10^b+1 quotient (v&-v)(u+1)=10^a-1, これをvから引くとα10^a+b-10^a-1=α01^b+10^a-1.

入れ子のかっこ

Lisp屋には違和感が全くない, 入れ子のかっこの話題だ. TAOCPにかっこの足跡(parenthesis trace)という話がある. (-sesと複数ではなく, -sisと単数なのはなぜか.)
例えば, 左かっこ4個, 右かっこ4個で, 正当なかっこの組み合せを作ると, ()()()(), ()()(()), ..., (((())))の14通りが出来る. n=4で14通りになるというのは, 2009年8月のブログ, 「投票数」に書いた通り.
₈C₄-₈C₃=70-56=14だ.

さて, かっこの足跡は左かっこを0, 右かっこを1で表した二進数である. 従って, ()()()()は01010101, (((())))は00001111となる. この方法で辞書式順で書くと

色の線で囲ったのは, 同じパターンが見えるところだ.

かっこ構造を保ったまま, 辞書式順で次の構造に移るにはどうするか.　二進法の表記では, 右から0と1を数えながら左へ`01'を探す. (0と1の数を#0, #1と書こう.)　ただしこれまで通過した#0<#1でなければならない. そういう`01'を見つけたらそれを`10'と交換し, その右に#0だけ0を並べ, さらにその右に#1だけ1を並べる. これは辞書式順で最小の数を作るためである.

テストするには, 左から探したいので, 左右逆転のリストを使い,

(define x '(1 1 1 1 0 0 0 0))

(define (next x)
 (let ((z 0) (o 0))
  (define (dup n a) (if (= n 0) '() 
   (cons a (dup (- n 1) a))))
  (define (nx x)
   (if x (if (and (= (car x) 1) (= (cadr x) 0) (< z o))
           (append (dup o 1) (dup z 0) '(0 1) (cddr x))
           (begin (if (= (car x) 0) 
              (set! z (+ z 1)) (set! o (+ o 1)))
            (nx (cdr x)))) 'ok))
  (nx x)))

で(set! x (next x))を, x=okになるまで次々と実行する.

実はこれを一気に, 定数ステップで計算するアルゴリズムがあるらしい. これこそプログラムハックである.



μ₀は...010101というパターンの二進数である.

その計算の進行の様子を下に示す.

上の段の左端は説明用の行数である. 左上のxで, 赤い`01'は交換する場所である. 次のozは, それより右の1と0の数. xのビット位置を右から0,1,2,...と数える. #0=#1となる可能性があるのは, 奇数番目の位置にある0なので, 3行目の右端の`01'の0, 8行目の`0101'の0が要注意だ. 奇数番目を取り出すべく, μ₀を利用する.

tはxの右から見て最初の奇数桁目の1の右を0にする. uは反対に最初の奇数桁目の1の左を0に, 右を1にする. つまり右の1は要注意の0を隠す. この1の数は偶数だが, (#0+1)*2であることにも注意. vはuとxの∨でxの交換すべき`01'の右がすべて1になった. 交換すべき`01'の左は現状のまま. これに1を加えれば, `01'が`10'になり, 右は#0+#1個の0になる. 次のv∧¬wは右から#0+#1+1個の1なので, #0+1桁右シフトすれば, 右端に#1個の1が得られる. それが√(u+1)で割る理由である.

すばらしい!

細線化アルゴリズム

細線化アルゴリズムがどう働くか, 「和」の字を細くしただけではよく分からない. そこでまた実験をした. 横40ピクセル, 縦20ピクセルの黒い長方形を作り, このアルゴリズムで細線化した. その様子を下の図に示す.

玉ねぎの皮を剥くように「骨から肉が削がれる」様子が見える. 第1回目は1と書いた矢印の先の赤のピクセルたちが消える. 上から下へ各行を処理, 行内では左から右へ処理するとしよう.
f(x_NW,x_N,x_NE,x_W,x,x_E,x_SW,x_S,x_SE)=x∧¬g(x_NW,...,x_W,x_S,...,x_SE)
もともと0のピクセルは, fの式に x∧ があるから0のままである. この長方形の左上の隅で1のピクセルに最初に出会う. この時使うgのパターンは下の図に左端のものだ.


gの前に ¬ があるから, このパターンの時は要らないと読む. 従って左上隅は消える. 次に上の縁を処理するが, これはgのパターンの左から2番目で消える. パターンの下の数字は, 前回のブログの図のどこにあったかを, 左端を0として示す. このようにして, 上の縁, 右上隅, 右の縁, 右下隅の赤い地帯が消えていく.

2回目はgのパターンを180度廻転するから, 外周の青で示す左と下が無くなる. 3回目はオレンジ色の部分が消え, 4回目に緑が消える. こうして20行あった横の列は, 19回で骨だけになるのであった.

楕円でも事情は同様であった.


この楕円のデータは, TAOCPに楕円曲線の塗り潰しの例にあったものを使った. TAOCPによると, 楕円や円など, 楕円曲線の塗り潰しはことの他簡単という. そのうちやってみよう.

細線化アルゴリズム

TAOCP 7.1.3項の式159のGuo & Hallの細線化アルゴリズムがある. Life Gameのように, ピクセルのキング近隣(セルオートマトンでは,
Moore Neighbourhoodという. 自分と周りの8隣り)を見, 次の時刻の状態を
f(x_NW,x_N,x_NE,x_W,x,x_E,x_SW,x_S,x_SE)=x∧¬g(x_NW,...,x_W,x_S,...,x_SE)
で決める. ただし, gは, 近隣が

の時, 1とする.
上の操作を奇数回目とすると, 次の偶数回目は, gのパターンの向きを180度廻転して使う. それで2回連続で変化がなければ停止する.

本当かなぁ, と思い, 実装, テストした. 和田研フォントの「和」の字を200×200のグリッドにしてやってみた. その結果が次である.



まずまずはうまく行った. 演習問題に, M行N列の黒四角を細線化するとどうなるかというのがあり, それもやってみた. 30行40列だと, 中央あたりに11個の横1列が出来て終わりであった. 当たり前か.