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