ユリウス日

前回のこのブログでは, 暦象年表のトピックスにあったグレゴリオ暦か らユリウス日(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であった.