ユリウス日

前回のこのブログでは, 暦象年表のトピックスにあったグレゴリオ暦か らユリウス日(jdn)への変換, またその逆の変換のプログラムの仕掛けを 解明した. かなり巧妙に出来ていて, ビューティフルコードといってもよ さそうである.

しかし私は, ビューティフルコードとしては, 見る人が「このプログラム は確かにこの問題を解く」と確信出来るようでありたいと願う. そういう こともあって, また最近 グレゴリオ暦の年月日からユリウス日への変換とそ の逆変換のプログラムを書いたのでお見せしたい. 先頃のブログ, Fixed Day Numberで紹介した式も活用している.

グレゴリオ暦の年月日からユリウス日への変換

def greg2jdn(y,m,d):
    def leap(y):
        return y%400==0 if y%100==0 else y%4==0
    return(
        y*365+y//4-y//100+y//400+d+1721425  
        -((4777-367*m)//12
          -(0 if m>2 else(1 if leap(y)else 2))))

まぁこれは見た通りであろう. 局所関数leap(y)は, y年に閏日があるとい う真理値を返す.

本体はy*365がy年末まで平年としての日数. y//4がその間の閏年の 数. -y//100が, 100の倍数の年の数で, グレゴリオ暦で閏にしない年の数. しかしこれでは, 1600年や2000年なども引くので, y//400で修正する.

その後+dは, 前月末までの日数にその月内の日数を足す. +1721425は -4712年からの値にする為の補正. 次の行の長い式が, 前月末までの 日数を計算する. このプログラムは一目瞭然だと思う.

例題をやってみる.

print(greg2jdn(1582,10,15)) => 2299161
    (グレゴリオ暦の初日)
print(greg2jdn(1,1,1)) => 1721426
    (Fixed Day Numberの初日) 

逆問題の方は, 局所関数としてminx(x,xs)を使う. これはxからxsの数を 次々と引き, この項を引くと値が負になる直前で止め, 引けた項数と, 最後の剰余を返す.

def minx(x,xs):
    i=0;
    while(i<len(xs)and x>=xs[i]):x=x-xs[i];i=i+1
    return [i,x]

例えば, 8と[1,3,5,7,9,11,13,15,17,19]が与えられた時, 先ず1を引いて 7が残り, 3を引いて4が残り, 5を引くと-1になるので, これは引かず, 1 と3の2項をを引き, 4が残る場合なので, [2,4]が返る. つまり前の項2が整数 平方根; 後の4が剰余. 8=2²+4を示す. xを0から19まで 変えて計算してみる.

list(map(lambda x: 
         minx(x,[1,3,5,7,9,11,13,15,17,19]),
         range(0,20)))

とすると

[[0, 0], [1, 0], [1, 1], [1, 2], [2, 0], [2, 1],
 [2, 2], [2, 3], [2, 4], [3, 0], [3, 1], [3, 2],
 [3, 3], [3, 4], [3, 5], [3, 6], [4, 0], [4, 1],
 [4, 2], [4, 3]]

となる. この逆変換の基本方針では, 残りの年の先頭が平年なら, ２月28日の後に 無理やりに2月29日に相当する日を挿入する. つまり残りのユリウス日が >=59なら, それに1を足すのである. 一見面倒なようだが, これが結構う まくいく.

ユリウス日からグレゴリオ暦の年月日への変換

def jdn2greg(jd):
    def minx(x,xs):
        i=0;
        while(i<len(xs)and x>=xs[i]):
            x=x-xs[i];i=i+1
        return [i,x]
    jd=jd-1721060
    a,jd=divmod(jd,146097)
    b,jd=minx(jd,[36525,36524,36524,36524])
    if(b>0 and jd>=59):jd=jd+1
    c,jd=divmod(jd,1461)
    d,jd=minx(jd,[366,365,365,365])
    if(d>0 and jd>=59):jd=jd+1
    m,jd=minx(
        jd,[31,29,31,30,31,30,31,31,30,31,30,31])
    return([a*400+b*100+c*4+d,m+1,jd+1])

局所関数minxは前述のもの.

jd=jd-1721060の引く値は, 0年1月1日からにユリウス暦だったして, そこを 起源にする変換である.

1721060=greg2jdn(1600,1,1)*5-greg2jdn(2000,1,1)*4

a,jd=divmod(jd,146097) は, 146097=365*400+97が4世紀の日数で, aが 4世紀が何回あったの数, jdは現在の4世紀になってからの日数.

b,jd=minx(jd,[36525,36524,36524,36524]) は, 現在の4世紀で今はどの 世紀かを見る. 0<=b<4. jdは世紀内の日数. ここでminxを使う.

bが1,2,3の時は, 先頭の年が100年で割れて閏年なのに閏日なしなので, if(b>0 and jd>=59):jd=jd+1 で閏があったように補正する.

c,jd=divmod(jd,1461) は世紀をまた4年組に分ける. 1461=365*4+1; 0<=c<25. d,jd=minx(jd,[366,365,365,365])

if(d>0 and jd>=59):jd=jd+1 は閏を入れる補正である.

こちらも例題.

jdn2greg(2299161) => [1582, 10, 15]
jdn2greg(1721426) => [1, 1, 1]

初めの例題でのa,b,c,dの値はa=3, b=3, c=20, d=2であった.

ユリウス日

今回のブログも国立天文台, 暦象年表のトピックスにあったユリウス日のアルゴリズムの続きだ. 今回は反対に, ある日のユリウス日からそのグレゴリオ暦の年月日を計算するものである.

直感的にもこちらの方が面倒なのは分る. 3つの数を掛けた積から, 元の数を 探すようなものだから.

その方法はこうだ.

def jdn2greg(jdn):
    l = jdn + 68569
    n = (4 * l) // 146097
    l = l - (146097 * n + 3) // 4
    i = (4000 * (l + 1)) // 1461001
    l = l - (1461 * i) // 4 + 31
    j = 80 *l // 2447
    d = l - 2447 * j // 80 
    l = j // 11
    m = j + 2 - 12 * l
    y = 100 * (n - 49) + i + l
    return [y,m,d]

前回のプログラムもそうだったが, 今回もif文がなく, 切り捨て除算ばかり なのが目立つ. O'BeirneのPuzzles & Paradoxesにある復活祭公式を 見る気分だ. (そういえば, 今年は明日4月20日が復活祭.)

まず分りそうなところから調べると, 下の方に m = j + 2 - 12 * lというのがあり, 前のプログラムのように, 1,2月は前年の11, 12月にし, 3〜12月を1〜10月に していることが分る.

また, y = 100 * (n - 49) + i + lから, nは, -4900年を基準にした世紀数 らしい. そこで前回同様の絵を描く. この図の左下, 1500年3月1日の(まだユリウス暦の 時代だがレゴリオ暦とした)jdnは青文字の2268983. この箱内の日数は赤文字の 2337552. 従って, -4900年3月1日のjdnは-68569. プログラムの先頭の 定数はこれであった. つまりlは-4900年からの日数であった.

次の行 n = (4 * l) // 146097 の146097は, 1460が365の4倍で, その 100倍にグレゴリオ暦の400年間の閏日を足したものである. lを100年の日数で割って 世紀数を出す代りに, lを4倍して400年の日数で割る. 理由は分らず. l = l - (146097 * n + 3) // 4 は世紀内の日数の計算.
i = (4000 * (l + 1)) // 1461001 は世紀内の年数の計算. ここも誤差処理の 為か, lを365では割らずに, 4000倍してから365の4000倍で割る. このiを 最後の年の計算で n+i のように使う.

l = l - (1461 * i) // 4 + 31 は年内の日数を計算する. 1461は閏日を 最後にした4年の日数で, 世紀内数に4年分を掛けてから4で割る. 最後に31を 足すのは, 1月を0ではなく, 1と出す為と思われる.

j = 80 *l // 2447 は 367 * m // 12 と同じように, 30と31が, 7,8月と 1,12月で連続する他は, 交互に現れるような式である.

f(l)= 80*l // 2447を, f(l)が1増えるようなlがどこからかを調べた表を 下に示す.

l = j // 11 で11, 12月(元々の1,2月)ならl=1, そうでなければ l=0に する.

最後の2行で年と月が正しく得られるという仕掛けであった.

誤差対策までは詳しく追わなかったが, このプログラムがどういう計算を しているかが判明したわけである. 中々面白い.

ユリウス日

国立天文台の暦計算室のウェブページには興味ある記事が満載だ. その中に2023年の暦象年表のトピックスとして「ユリウス日について」>というのがあった.

その最後の方に, グレゴリオ暦の年月日からその日のユリウス日(Julian Day Number)を 計算する式があり, その計算の仕方が面白かったので, 説明したい.

Pythonで書いたその式はこうだ. //は切り捨て除算である.

  def jdn(y, m, d):
    k = (14 - m) // 12
    return (
        (- k + y + 4800) * 1461 // 4 
        + (k * 12 + m -2) * 367 // 12 
        - ((- k + y + 4900) // 100) * 3 // 4
        + d - 32075)

私は普段, Schemeを使うので, Schemeで書くとこうだ.

  (define (jdn y m d)
  (let ((k (quotient (- 14 m) 12)))
  (+ (quotient (* (+ (- k) y 4800) 1461) 4)
     (quotient (* (+ (* k 12) m -2) 367) 12)
     (- (quotient (* (quotient (+ (- k) y 4900) 100) 3) 4))
     d -32075)))

暦の計算でなんとも鬱陶しいのは2月の末にあったりなかったりする閏日だ. 今回話題にする プログラムでは, 1年は3月から始まり, 1月と2月は前年の13月と14月にする扱いである.

また紀元前(BC, BCE)はマイナス何年の表記と1年ずれるのも面倒で, 暦の計算ではBCなど は使わないのが普通である.

さて, 「ある日」dのユリウス日jdn(d)の計算には, どこか基準になる日d₀とjdnの分かって いる日d₁を決め, d₀からd₁までの通日n(d₁)と jdn(d₁)の差(Δ_d)を, dの通日n(d)に足す.

Δ_d=jdn(d₁)-n(d₁), jdn(d)=n(d)+Δ_d

このプログラムでは, d₀をユリウス日の開始日の-4712年より前で, 400で整除出来る年の -4800年3月1日として, その日からの通日を使う.

通日の計算は次の図を見て欲しい. いまjdnの分かっている日d₁を, グレゴリオ暦への改暦 (1572年)後の1600年3月1日にし, 「ある日」dを1600年10月21日にして, ユリウス日の計算の進行を見よう. この日は 日本の歴史では, 旧暦の慶長5年9月15日, 関ヶ原の合戦があった.

私の「個人用電卓」 で10月21日のjdnを予め計算すると2305742である.

まず-4800年3月1日から1599年14月28日までの, 閏日を無視した長方形を 描く. 閏日は, 4行ごとの最後に入る. 100年毎や400年毎の補正は, 100行ごとの最後, 400行ごとの最後に入る. 長方形の最下行の右端に 飛び出しているのが, 1600年2月(1599年14月)29日で, 1600年は100 で整除出来るが, 400年でも整除出来るから, 閏日は入れる.

365*6400=2336000
4年毎の閏日は 6400/4=1600,
100年毎の閏日は 6400/100=64,
400年毎の閏日は 6400/400=16.

従ってこの(多少の飛出しのある)長方形は, 233600+1600-64+16=2337552. 左下の赤字の値で, 3月1日の 通日である.

次は1600年内の計算で, 3月0日から「ある日」までの通日を計算する. それには3月0日から m月0日までの通日が必要になる. これに,

f(n)=367*n//12

が使えるかもしれない. nを0,1,2,...と替えながらこの値と, f(n+1)とf(n)との差Δ_fを 計算してみると, 次のようになる.

Δ_fは赤線のように30と31の12個の列が循環し, 途中31が連続する箇所が2回ある. これが7月, 8月および, 12月, 13月に相当するわけで, mの行に対応する月を記入した. f(n)が各月までの通日になる. 3月までの通日を 0にしたいところだが, それは最後の補正の時にやることにして, これを使う.

大体の様子が分ったので, もう一度プログラムを見る.

  def jdn(y, m, d):
    k = (14 - m) // 12
    return (
        (- k + y + 4800) * 1461 // 4 
        + (k * 12 + m -2) * 367 // 12 
        - ((- k + y + 4900) // 100) * 3 // 4
        + d - 32075)

最初のkの値は, mの1から12に対して, 1,1,0,0,0,0,0,0,0,0,0,0,0,0 となり, 1月, 2月を別扱いにする時に使う.
次の行, (- k + y + 4800) * 1461 // 4 は, y - k で, 1,2月を前年に繰上げ, -4800を足して基準を-4800年にし, 1461(=365+365+365+366)を掛けて4で 切り捨て除算し, 4年毎の閏日を処理する.

次の行, + (k * 12 + m -2) * 367 // 12 は. m=3, 4, ..., 12を. 1, 2, ..., 14に変換して, m月0日までの通日を得る. 但し3月の欄で見るように30日多い. これは Δに組み込むことにする.

その次の行は100年の補正である. 年を100で割り世紀ごとの補正にする. 0,1,2,3,...,7 について, * 3 // 4をやると, 0, 0, 1, 2, 3, 3, 4, 5,...が得られる. 1,2,3と 増えるのは, 100で割った値が4で割れないで, 閏にしない年, 2行上では4で割れる年はいつも閏にしていたのから これを引いて補正した. 0,0や3,3のように繰り返すのは閏にする年で, 補正にはならない. ただこれを見ると, 補正しないのは-4800年の次の世紀年で, これでは1世紀遅すぎる. その修正のため, ここでは4800ではなく, 4900を足すのである.

改で左下を見ると赤字は前年までの通日, オレンジ色は3月1日までの通日, 従ってn(d₁)は 2337583になる. 一方3月1日のユリウス日jdn(d₁)は2305508でΔは プログラムにあるように, 32075になる.

1600年10月21のjdnは 2337552+244+21-32075=2305742 と得られた.

ユリウス日

数日前に書いた個人用電卓にユリウス日のキーがあった. これは特に週日を求めるためのものである. ユリウス日に1を足し, 7で法をとると日曜を0, 月曜を1, ..., 土曜を6とする週日が得られるのである. またこの電卓は, 除算のとき, スタックトップには商が, そのすぐ下には剰余が出るので, 20081020 [JD] 1 [+] 7［/] [Pop] -> 1 となり, 本年10月20日は月曜と分かる. -4712年1月1日も月曜であったわけだ. ユリウス日は, -4712年1月1日 世界時正午を0.0とし, それからの通日である. 2008年10月20日のユリウス日が2454760というのは, 正確にはその日の世界時正午, JSTで午後9時のユリウス日が2454760.0であるということだ. それを省略してその日のユリウス日という. -4712年1月2日のユリウス日は1である. 理科年表や天文年鑑にも計算用の表がある.

A表              B表    *   C表      *
-1000 1355808     0  0,-1    1   0   0
 -900 1392333     1   365    2  31  31
 -800 1428858     2   730    3  59  60
 ...              3  1095    4  90  91
 1300 2195883     4  1460   ...
 1400 2232408     5  1826   11 304 305
 1500 2268933    ...        12 334 335
 1500 2268923    98 35794
 1600 2305448    99 36159
...
 2000 2451545

例えば2004年4月23日のユリウス日ならA表から2000年の2451545, B表から4年の1460, C表から4月の右の閏年用の91をとり, それに23を足して(+ 2451545 1460 91 23) -> 2453119が得られる. このA表は, 例えば2000年なら, -4712年1月1日から1999年12月31日までの日数である. 従ってこれは2001年1月1日のユリウス日のになりそうだが, B表の最初0年に, 閏年なら-1, 平年なら0とあるのが重要で, 2000年は閏年なので, 1引くから結局1月0日のユリウス日になる. 変換の式もサーチエンジンで探せるが, 私が情報処理学会誌のHaskell Programmingに書いた,

mon0, mon1 :: [Int]
mon0 = [0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334]
mon1 = [0, 31, 60, 91, 121, 152, 182, 213, 244, 274, 305, 335]

gleap, jleap :: Int -> Bool -- Gregorian, Julian暦のうるう年
gleap y = if y `mod` 100 == 0 then y `mod` 400 == 0
                              else y `mod` 4 == 0
jleap y = y `mod` 4 == 0

julianDate :: Int -> Int -> Int -> Int
julianDate y m d =                   -- aからgは途中の値
 let a = (y + 4712) * 365
     b = (y + 4712 + 3) `div` 4
     c = if y > 1601 then y' `div` 400 - y' `div` 100 else 0
           where y' = y - 1601
     e = if [y,m,d] >= [1582,10,15] then -10 else 0
     f = (if leap y then mon1 else mon0) !! (m - 1)
           where leap = if y > 1600 then gleap else jleap
     g = d - 1
  in a + b + c + e + f + g           -- 最後の結果

のようにも書ける. インターネットなら, https://www.onlineconversion.com/julian_date.htmにJulian Date Converterがあり, Year, Month, Day(, Hour, Minute, Second)を入れ, Computer Julian dateを押すと結果が返ってくる. 個人用電卓では結局上の式を使った. 電卓は実はやっとiPod touchに実装した. 電車の中でもテストが出来て楽しいが, 実用化にはまだ改良の余地は多々ある.