要するに暦に関するアルゴリズムが満載で面白い. この本では暦の変換の基準は, Fixed Day Numberといい, Gregorian暦が過去に遡って使われたとして, その1年1月1日(の正子)を 第1日とするものである. いちいちfixed day numberといわず, R.D.と いうらしい. (Rata Die)
似たようのものに, -4712年1月1日正午から起算するユリウス日があるが, これは途中でJulian暦 からGregolian暦になるので変換が面倒である.
Calendrical Caluculationにある, Gregorian暦の年月日からそのR.Dを計算する アルゴリズムを見てみよう. 本書のアルゴリズムはCommon Lispで記述されているが, 今日はPythonに翻訳する.
def fixeddaynumber(year, month, day):
return (365 * (year - 1) + (year - 1) // 4
- ((year - 1) // 100) + (year - 1) // 400
+ (367 * month - 362) // 12
+ (0 if (month <= 2) else
(-1 if leap_year (year) else -2))
+ day)
最初の365 * (year - 1)は前年までがすべて平年であったとしての, 前年の
最後の日までの総日数である. しかし, 4年に一度閏年があるから, その分を足す.
+ (year - 1) // 4. Gregorian暦では, 100で整除される年は平年に
するのであったからそれを引く. - ((year - 1) // 100). しかし400で
整除されるなら閏年にするからそれを足す. + (year - 1) // 400.
前年末までの総日数はこれで計算出来た. 次は前月末までの日数を計算する. 目標としては, 1月なら0. 2月なら31, 3月は 平年なら59, 閏年なら60, ... 閏年の補正は後でやることにし, 本書にあるアルゴリズム は(367 * month - 362) // 12というものだ. monthを 1から12まで変えてこの値を計算すると:
print(list(map (lambda m:(367*m-362)//12, range(1,13))))
[0, 31, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336]
3月からは, 平年なら2, 閏年なら1多いことがわかる. 従ってこの関数の
month番目を使い, monthが1, 2ならこの値そのまま, 3以降はこの値から
平年は2, 閏年は1を引いて使う.これで前年末までと前月末までの総日数が分ったから, 後はdayを足せばよい.
閏年の補正に本書ではこういう式を使う.
def leap_year (year):
return (year % 4) == 0 & ((year % 400)
in [100, 200, 300])
これで1年1月1日を計算すると:
print (fixeddaynumber (1, 1, 1)) 1なるほど. また, おととし, 開業150年と話題になった鉄道開通の日, 1872年10月14日(明治5年には日本は まだ旧暦を使っていたが, これは太陽暦に換算してある)のR.D.を計算してみる.
各々の項の値は682915, 467, -18, 4, 275, -1, 14で, R.D.は683656 である.
ところで私は何度も暦の計算のプログラムを書いたが, その年の終までの日数を計算し, その月の前までの日数を引いて日を足すのが好きだ. 上のアルゴリズムに year - 1が4回もあるのが気に食わないのである.
それで普通は年末から今月0日までの日数を引くことになるが, 大体はそこは定数の 表を利用する. しかし, Calendricalの本にあるような式が使えないかと 考えた. Calendricalの式と同様, 1,2月は特別扱いにしてもよいとして, 定数など を変えてテストした. 欲しい値は:
[[365, 366], [334, 335], 306, 275, 245, 214,
184, 153, 122, 92, 61, 31]
である. 左端が1月と2月で, かっこ内は左が平年, 右が閏年. 右端が12月.
その結果, (4777 - 367 * m) // 12がよいらしかった. つまり
print(list(map (lambda m:(4777 - 367 * m) // 12,
range(1,13))))
[367, 336, 306, 275, 245, 214, 184, 153, 122, 92, 61, 31]
1月, 2月は平年は2を引き, 閏年は1を引く. ここで使う定数4777は非常にクリティカルで, 4776や 4778に変えてみると, 前者は7月が小さく, 後者は2月が大きい.
print(list(map (lambda m:(4776 - 367 * m) // 12,
range(1,13))))
[367, 336, 306, 275, 245, 214, 183, 153, 122, 92, 61, 31]
print(list(map (lambda m:(4778 - 367 * m) // 12,
range(1,13))))
[367, 337, 306, 275, 245, 214, 184, 153, 122, 92, 61, 31]
これを利用して書いたfixed day number関数は次の通り:
def fixeddaynumber2(year, month, day):
return (365 * year + year // 4 - (year // 100)
+ year // 400 - (4777 - 367 * month) // 12
+ (0 if (month > 2) else
(-1 (year % 400 == 0)
if (year % 100 == 0) else
(year % 4 == 0) -2))
+ day)
こちらで鉄道記念日を計算したのと, 本年6月30日の値の7で割った剰余とを見ると:
print (fixeddaynumber2 (1872, 10, 14)) 683656 print ((fixeddaynumber2 (2024, 6, 30)) % 7) 0だから, R.D.を7で割った剰余が曜日になるわけだ. つまり, R.D.1の日, 1年1月1日は月曜であった.
年末から月始めまでの日数を今回のように式で計算すると, 実引数に0月とか１３月が与え られても計算してしまう恐れがある. 表なら範囲外としてエラーに出来る. それを 心配したとしても, 日の値に範囲外が与えられる可能性はあるのだから, まぁ我慢することに しよう.
0 件のコメント:
コメントを投稿