unix time

前回のunix timeのブログに書いたように, Calendricalの方法は, 確かに簡単ではあるが, 1から勘定を始めるプログラムに, 多少の違和感を持っていた私は, やはり0から始めるプログラムが書いてみたく, そのため, 366, 365, 365, 365のように, おなじものが並ぶ後続の前に, それらより1だけ多いか少ない要素が先頭にあるときの処理について考えていた.

つまり例えば1600年1月1日から1999年12月31日の400年間, 146097日には, 最初の100年は1600年がうるう年なので, 36525日あり, 後の3回の100年は, 最初の00年が平年なので, 36524日であって, 36525, 36524, 36524, 36524であり, またその36524日ある100年を4年ごと25組に分けると, 第2組以降は最初の4年が1日少なく, 1460, 1461, 1461, ..., 1461となる.

このような問題提起について, d日目を入力し, その該当する年yと, その年内の日数d'を計算する方法を探した.

何万日を扱うのも憂鬱なので, とりあえず, 先頭が1多い4日, 3日, 3日, 3日のパターンと, 先頭が1少ない, 3日, 4日, 4日, 4日のパターンを考える.

4,3,3,3のパターンでは, 総日数は13日で, 0≤d<13に対し,

d         y  d'
0,1,2,3   0  0,1,2,3
4,5,6     1  0,1,2
7,8,9     2  0,1,2
10,11,12  3  0,1,2

また, 3,4,4,4のパターンでは, 総日数は15日で, 0≤d<15に対し,

d             y  d'
0,1,2         0  0,1,2
3,4,5,6       1  0,1,2,3
7,8,9,10      2  0,1,2,3
11,12,13,14   3  0,1,2,3

となるようにしたい.

前者については

(define (foo d)
(let ((y (floor (/ d 3.25))))
(list y (+ d -13 (floor (* 3.25 (- 4 y)))))))

とする. 13は総日数, 4は年の幅, 怪しげな3.25は13/4である.

後者については

(define (bar d)
(let ((y (- 3 (floor (/ (- 14 d) 3.75)))))
(list y (- d (floor (* y 3.75))))))

とする. 14は総日数-1, 3.75は総日数の15を4で割ったものだ.

とりあえず, これでdに対して(y d')を求めてみると

(map foo (a2b 0 13))
=>((0 0.) (0. 1.) (0. 2.) (0. 3.) (1. 0.) (1. 1.) (1. 2.)
(2. 0.) (2. 1.) (2. 2.) (3. 0.) (3. 1.) (3. 2.))

(map bar (a2b 0 15))
=>((0. 0.) (0. 1.) (0. 2.) (1. 0.) (1. 1.) (1. 2.) (1. 3.)
(2. 0.) (2. 1.) (2. 2.) (2. 3.) (3. 0.) (3. 1.) (3. 2.) (3 3.))

とうまく行きそうである.

fooとbarの図を描いてみるとこうなる. 上がfoo, 下がbar.



図でははっきりしないが, dに対するfloorを取る前の値は
(0 .31 .62 .92 1.23 1.54 1.85 2.15 2.46 2.77 3.08 3.38 3.69)
と
(3.73 3.47 3.2 2.93 2.67 2.4 2.13 1.87 1.6 1.33 1.07 .8
.53 .27 0)
である.

そこで, この3とか4とかを実際の値に変えて関数を書き, クリティカルな値に対してテストすると次の通りだ.

(define (foo d)
(let ((y (floor (/ d 36524.25))))
(list y (+ d -146097 (floor (* 36524.25 (- 4 y)))))))
(map foo '(0 36524 36525 73048 73049 109572 109573 146096))
=>((0 0.) (0. 36524.) (1. 0.) (1. 36523.) (2. 0.) (2. 36523.)
  (3. 0.) (3. 36523.))

(define (bar d)
(let ((y (- 24 (floor (/ (- 36523 d) 1460.96)))))
(list y (- d (floor (* y 1460.96))))))
(map bar '(0 1459 1460 2920 2921 4381 4382 5842 33602 35062
  35063 36523))
=>((0. 0.) (0. 1459.) (1. 0.) (1. 1460.) (2. 0.) (2. 1460.)
  (3. 0.) (3. 1460.) (23. 0.) (23. 1460.) (24. 0.) (24 1460.))

(define (foo d)
(let ((y (floor (/ d 365.25))))
(list y (+ d -1461 (floor (* 365.25 (- 4 y)))))))
(map foo '(0 365 366 730 731 1095 1096 1460))
=>((0 0.) (0. 365.) (1. 0.) (1. 364.) (2. 0.) (2. 364.)
  (3. 0.) (3. 364.))

あとはこれを繋げるだけ. 注意すべきは, Calendricalの方は1年1月1日から始まるのに対し, こちらは0年から399年までの400年から始まるので, 0年分の366日を足す; Calendricalは1月1日のfixed dateが1なのに対し, こちらは, 0年1月1日が0なので1日分の修正が必要だ. この修正の結果がd0である. 400年のところは, 常に146097日なので, 計算は簡単だ. 残りは上の関数を利用している.

(define (gr rd)
(let* ((d0 (+ rd 366 -1))
       (y400 (quotient d0 146097))
       (d1 (modulo d0 146097))
       (y100 (floor (/ d1 36524.25)))
       (d2 (+ d1 -146097 (floor (* 36524.25 (- 4 y100)))))
       (y4 (if (= y100 0) (quotient d2 1461)
           (- 24 (floor (/ (- 36523 d2) 1460.96)))))
       (d3 (if (= y100 0) (modulo d2 1461)
           (- d2 (floor (* y4 1460.96)))))
       (y1 (if (and (> y100 0) (= y4 0)) (quotient d3 365)
           (floor (/ d3 365.25))))
       (d4 (if (and (> y100 0) (= y4 0)) (modulo d3 365)
           (+ d3 -1461 (floor (* 365.25 (- 4 y1))))))
       (y (+ (* y400 400) (* y100 100) (* y4 4) y1))
       (leap (if (= (modulo y 100) 0) (= (modulo y 400) 0)
           (= (modulo y 4) 0)))
       (e (if leap
        '(0 31 60 91 121 152 182 213 244 274 305 335)
        '(0 31 59 90 120 151 181 212 243 273 304 334)))
       (mon (apply + (map (lambda (f) (if (<= f d4) 1 0)) e)))
       (d (- d4 (list-ref e (- mon 1)) -1)))
    (list y mon d)))

テスト

(gr 1) => (1. 1 1.)
(gr 734274) => (2011. 5 16.)

こんな具合いだが, まぁいいか.

unix time

前回のblog以後, ユリウス日を年月日に変換するプログラムも書きたくなっている.

私の愛読書の1つ, E.M.Reingold, N. DershowitzのCalendrical Calculationsでその辺の計算法を見てみよう. ユリウス日は, Julian暦-4712年1月1日が起点だが, 本書のカレンダーの起点は, Gregorian暦を昔の方へ外挿した1年1月1日を1とし, それからのFixed Day Numberで表わす. この日数をR.D.(Rata Die, fixed dateのラテン語)という.

本書の関数は, 本文では不思議な構文で書いてあり, 付録にはCommon Lispによる実装があるが, Scheme風にすると,

(define gregorian-epoch 1)
(define (gregorian-leap? y)
 (and (= (modulo y 4) 0)
   (not (member (modulo y 400) '(100 200 300)))))

(define (fixed-from-gregorian y m d)
 (+ gregorian-epoch -1 (* 365 (- y 1)) (quotient (- y 1) 4)
    (- (quotient (- y 1) 100)) (quotient (- y 1) 400)
    (quotient (- (* 367 m) 362) 12)
    (cond ((<= m 2) 0)
      ((gregorian-leap? y) -1)
      (else -2)) d))

(fixed-from-gregorian 1 1 1) => 1
(fixed-from-gregorian 2011 5 16) => 734273

この本では, 年内の日数と月の変換は, 2月も30日あるとして計算し, 補正する

(map (lambda (m) (quotient (- (* 367 m) 362) 12)) (a2b 1 13))
=> (0 31 61 92 122 153 183 214 245 275 306 336)

floorをとらないと, 以下のような値である.

(.417 31. 61.583 92.167 122.75 153.333 183.917 214.5
245.083 275.667 306.25 336.833)

隣り同士の差を取ってみると,

(map (lambda (a b) (- a b))
 '(31 61 92 122 153 183 214 245 275 306 336)
 '(0 31 61 92 122 153 183 214 245 275 306))
=>
(31 30 31 30 31 30 31 31 30 31 30)

補正は, 1月2月はこのまま, 3月から後はうるう年なら1を引き, 平年なら2を引く.

さて, Gregorian暦のy, m, dからfixed dateを得る最初の関数fixed-from-gregorianの解説である.

前年までの経過日数が必要なので, (- y 1)が頻出する.

(* 365 (- y 1))      ;前年までの平日の日数
(quotient (- y 1) 4) ;Julian暦のうるう日の日数
(- (quotient (- y 1) -100)) ;100年の倍数の年はうるうをやめる
(quotient (- y 1) 400) ;しかし400年の倍数なら, やはりうる
う年にする
(quotient (- (* 367 m) 362) 12) ;2月を30日と仮定して, m月
の前月までの日数
(cond ((<= m 2) 0) ;補正
 ((gregorian-leap? y) -1)
 (else -2)) d))      ;dを足す.

最初のgregorian-epochは, 1年1月1日を1にするためである.
最後のテスト例のようにうまく行く.
予想通り, 逆は難しい. プログラムは以下のようだ.

(define (gregorian-year-from-fixed rd)
  (let* ((d0 (- rd 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))
         (y (+ (* 400 n400) (* 100 n100) (* 4 n4) n1)))
  (if (or (= n100 4) (= n1 4)) y (+ y 1))))

(define (gregorian-from-fixed rd)
  (let* ((y (gregorian-year-from-fixed rd))
   (prior-days (- rd (fixed-from-gregorian y 1 1))     
　 (correction (cond ((< rd (fixed-from-gregorian y 3 1)) 0)
                ((gregorian-leap? y) 1
                (else 2)))
   (m (quotient (+ (* 12 (+ prior-days correction)) 373)
             367))
   (d (+ (- rd (fixed-from-gregorian y m 1)) 1)))
  (list y m d)))

(gregorian-from-fixed 1) => (1 1 1)
(gregorian-from-fixed 734273) => (2011 5 16)

と定義しておき, (gregorian-year-from-fixed rd)
で rd に対する年 y を計算する.

(n400 (quotient d0 146097)) :d0に400年の日数が何回あるか見る. 146097は400年の日数.
(d1 (modulo d0 146097)) :その400年内の日数

(n100 (quotient d1 36524)) ;その日数に100年は何回あるか. しかし400年の100年の日数は, 最初の100年の先頭の00年はうるう年なので, 36525日あるから

(d2 (modulo d1 36524)) ;その100年内の日数
(n4 (quotient d2 1461)) ;その日数内の4年の数
(d3 (modulo d2 1461)) ;その4年内の日数
(n1 (quotient d3 365)) ;4年内の1年の数
(n1 (quotient d3 365)) ;4年内の1年の数

基数変換をやっているみたいに簡単なのに驚く. 400年の中の日数d1のいろいろな値から得られるyをみてみる.

(define (gregorian-year-from-fixed rd)
  (let* ((n100 (quotient rd 36524))
         (d2 (modulo d 36524))
         (n4 (quotient d2 1461))
         (d3 (modulo d2 1461))
         (n1 (quotient d3 365))
         (y (+ (* 100 n100) (* 4 n4) n1)))
  (if (or (= n100 4) (= n1 4)) y (+ y 1))))

すると, 以下のようになっていることが分かる.

       rd       y                 rd       y
     0-   364   1   365     35794- 36158  99   365
   365-   729   2   365     36159- 36523 100   365
   730-  1094   3   365     36524- 36888 101   365
  1095-  1460   4   366
  1461-  1825   5   365     72318- 72682 199   365
  1826-  2190   6   365     72683- 73047 200   365
  2191-  2555   7   365     73088- 73412 201   365
  2556-  2921   8   366
  2922-  3286   9   365    108842-109206 299   365
  3287-  3651  10   365    109207-109571 300   365
  3652-  4016  11   365    109572-109936 301   365
  4017-  4382  12   366
  4383   4747  13   365    145001-145365 398   365
  4748-  5112  14   365    145366-145730 399   365
  5113-  5477  15   365    145731-146096 400   366
  5478-  5843  16   366

0から勘定を始めるのが基本と思っている私なら, 下の図の左のような関数を書くところだが, 本書の流儀は違う. d1=0ならy=1が返る. d1=146096なら, y=400であった. なるほど! (図の太い横線は366日の年を示す. 横線の左の黒丸は閉区間, 右の白丸は開区間を示す.)




別の図を描くと下のようになる. つまり, 図で網掛けの例外を最後に置くので計算が簡単になっていたのだ.




それなら, 私の流儀でも, 400年紀末から逆に計算すればおなじわけだった.


ところで, これではまだfixedから, 年が得られただけである. 本番のプログラムはこれを下請けに使う(gregorian-from-fixed d)である. その解法はこうだ.

とりあえずこのrdの落ちるyを求める. 次のその年の1月1日のfixed dateを求め, prior-daysとする. その年の3月1日のfixed dateを求め, prior-daysがそれより小さければ, 補正は0, そうでなくて, うるう年なら補正は1, 平年なら補正は2である. 次に前に日数を計算した式で, 月を見つけ, その月の1日のfixed dateとの差から, 日が分かるのである.

fixed-dateを何回も使うが, それだけ分かりやすいアルゴリズムになっている.

こういうプログラムの解読もなかなか面白い.

unix time

unixには1970年1月1日正子(0時0分)からの延秒数を数えている32ビットの時計がある. 2038年1月19日にサインビットが立つといわれ, 2000年問題みたいになにか起きかもしれないが, 私は多分もうこの世にはいず, 状況を知ることはかなわぬ.

この時計の元は, MITのMulticsではないか. Multicsには, 1900年1月1日正子からのマイクロ秒を数える52ビット時計があった. マイクロは10^-6だから, 20ビット程度であり, unixの32ビットに対して52ビットなのは分かる. 私がMITに滞在したのは, 1973年9 月から74年7月までだが, その時計のサインビットが立ったのは, その少し前のたしか5月だったとある院生から聞いた.

まず脱線して, それがいつだったか計算してみよう. 例の個人用電卓が活躍する. 2⁵¹は2251799813685248. これを1日のマイクロ秒864000000000で割る.

商は26062, 剰余は43013685248. つまり1900年1月1日から26062日後を知りたい. それには1900年1月1日のユリウス日2415021に26062を足し, その2441083がユリウス日になる日を知ればよい.

ここから先は電卓から離れ, 理科年表のユリウス日の表による.

すると, 1971年5月11日がその日であることが判明. 日以下を計算すると, 11時 56分 53秒 685248マイクロ秒であった.

jdをユリウス日を計算する関数として, Schemeで検算すると,

(+ (* (- (jd 1971 5 11) (jd 1900 1 1)) 86400000000)
   (* 11 3600000000) (* 56 60000000) (* 53 1000000) 685248)
=> 2251799813685248

(factorize 2251799813685248)
=>
(2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2)

(length (factorize 2251799813685248)) => 51

本題へ戻り, 通常のy年mon月d日h時min分s秒からunix timeへの変換の方は, Multicsの場合と同様に簡単だ.

(JD(y,mon,d)-JD(1970,1,1))*86400+h*3600+min*60+s

となる.

一方, unix time tからy,mon,d,h,min,sへの変換は, ユリウス日からy,mon,dへの変換関数があればなんでもないが, まだそういうプログラムは書いたことがない.

カレンダーの変換で一番面倒なのは, Grogorian暦でのうるう年の計算である. しかし, ことunix timeに限れば, 2000年が通常のうるう年なのが幸いだ. そこで書いたのが次のSchemeのプログラムである.

(define (unixtime t)
 (let* ((s (modulo t 60)) (m (quotient (modulo t 3600) 60))
  (h (quotient (modulo t 86400) 3600)) (d (quotient t 86400))
  (y (- (floor (/ (+ d 731) 365.25)) 2))
  (c (- d (floor (+ (* y 365.25) 0.25))))
  (e (if (= (modulo y 4) 2)
        '(0 31 60 91 121 152 182 213 244 274 305 335)
        '(0 31 59 90 120 151 181 212 243 273 304 334)))
  (n (apply + (map (lambda (f) (if (<= f c) 1 0)) e))))
 (list (inexact->exact (+ 1970 y)) n
  (inexact->exact (- c (list-ref e (- n 1)) -1)) h m s)))

引数のtがunixtimeである. 最初にs(秒), m(分), h(時)を取り出す. dは通算の日数になる. その後の変数は, yが1970年以降の年数, cがその年内の日数, nが月だ. yはdを365で割ってfloorを取ればよいが, 閏年があるからそうは問屋が卸ろさない.

まず, y年について, 前の年の終りまでの日数は,



欲しい値は,



365の代りに365.25で割ればよさそうに見えるので, 2.25でテストしてみる.

(map (lambda (n) (floor (/ n 2.25))) (a2b 0 20))
=> (0 0. 0. 1. 1. 2. 2. 3. 3. 4. 4. 4. 5. 5. 6. 6. 7. 7.
   8. 8.)

4が3個並ぶのがうるう年に対応し, これを2年にしたいから, 3つの0と2つの1の5個をスキップするために, nの代りに(+ n 5)とし, 最後に2を引く.

(map (lambda (n) (- (floor (/ (+ n 5) 2.25)) 2)) (a2b 0 20))
=> (0. 0. 1. 1. 2. 2. 2. 3. 3. 4. 4. 5. 5. 6. 6. 6. 7. 7.
   8. 8.)

なるほどうまくいくので, 356.25に修正し, テストする.

(map (lambda (n) (- (floor (/ (+ n 366 365) 365.25)) 2))
'(0 364 365 729 730 1095 1096 1460 1461 1825 1826 2190
   2191 2555))
=> (0. 0. 1. 1. 2. 2. 3. 3. 4. 4. 5. 5. 6. 6.)

うまくいく. これを1970に足せばよい.

次に上の表にあった前年までの日数の和を計算するには,

(lambda (y) (+ (* y 365) (quotient (+ y 1) 4)))

か

(lambda (y) (floor (+ (* y 365.25) 0.25)))

とする. これをdから引くと, その年内の日数cが得られる. 前の年の終りまでの日数の和のように, 前の月の終りまでの日数の和のリストを, うるう年か否かで変数eに用意する.

月nは, このリストで, cが越えるものの数として得る. nが決れば, 月内の日数は, cから先ほどのリストの要素を引いて作る.

これで完成. テストしてみる.

(unixtime 0) => (1970 1 1 0 0 0)
(unixtime (expt 2 31)) => (2038 1 19 3 14 8)