ユリウス日

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

直感的にもこちらの方が面倒なのは分る. 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行で年と月が正しく得られるという仕掛けであった.

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