ユリウス日

国立天文台の暦計算室のウェブページには興味ある記事が満載だ. その中に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 と得られた.

暦の使い方

私は暦が好きなので, Shakespeareの劇を読んでいても「暦」の文字があると立ち止ってしまう.

例えば, 誰でも思い出すであろう「夏の夜の夢」3幕1場で, 6人の職人が森の中で芝居の稽古をする場面:
スナウト おれたちが芝居をやる晩に月は出るのか?
ボトム 暦だ, 暦だ! 今年一年の暦を見て, 月が出るか調べるんだ, 月は出るか.
クインス うん, その晩月は出る.
(ここでの引用は, 私と同時代人である小田島さんの訳による.)

旧暦の国の人なら, 日付けで月のありなしは分る.
1日の前後は月出は夜明, 月入は夕方
8日の前後は月出は夜中, 月入は昼間
15日の前後は月出は夕方, 月入は夜明
22日の前後は月出は昼間, 月入は夜中
である. だから明治5年12月2日, 旧暦の最後の日は, 地球照で鎌のように細く光る三日月が 西の空にあったに違いない.

英国は, というかヨーロッパの殆んどは太陽暦だったから, 月の様子を知るには 暦が必要だったらしい.
しかしまだ疑問がある. 今我々の手元にある普通の暦は, 日付と曜日の対応を示 すのが主目的で, 次に祝日, 二十四節気, 雑節, さらに詳しいのには大安や友引の六曜の日付 も記載がある程度. 日出, 日入, 月出, 月入の時刻や月相は殆んど書いて ない. 天体の情報があるのは天文暦almanacだが, 英文「夏の夜の夢」の 当該部分はcalendarになっている.

東京の毎日の日出入, 月出入は, 理科年表の暦部や 暦象年表あるいは新聞の隅をよく見なければならない.

ところが1幕1場では:
シーシュース 楽しい日々をあと四日すごせば
新月の宵となる.
ヒポリタ そうすれば新月が, 天高く引きしぼられた
銀の弓のようにかかり, 私たちの婚儀の夜を
見まもってくれるでしょう.

だから, いつが新月か知る人は知っていた. 婚儀の日は新月なのに 「うん, その晩月は出る」というShakespeareの無頓着ぶりも楽しい.

ユリウス暦を採用した張本人「ジュリアス・シーザー」2幕1場では:

ブルータス そうだ, 明日は三月十五日であったかな?
ルーシャス さあどうでしたか.
ブルータス 暦を調べて, すぐ知らせてくれぬか.
ルーシャス かしこまりました.
ルーシャス 旦那様, 三月は十四日すんだところです.

正確に管理された日めくりカレンダーでなければ, 今日が何日かを知るの も普通には難しい. 最近ならスマホが手元にあれば分るが... さらにローマの暦の 日付は非常にややこしいから, 何日経過したかを彼等はどう計算したのか 興味がある.
ローマ暦のややこしい日付けに関しては この解説を見て欲しい.
その3幕1場で:
シーザー 三月十五日がきたぞ.
占い師 きましたが, シーザー, まだすぎさってはいません.

この「まだすぎさってはいません」がずっしり重い.

Shakespeareの史劇に登場する薔薇戦争時代で最も有名で残忍で話題も多いのは 「リチャード三世」であろう. 最近でも遺骨が発見されて改葬されたり, エドワード四世 の2人の王子はロンドン塔で殺害されたのではないというニュースがあったりした.

「リチャード三世」3幕8場:
リチャード おい, 暦をとってくれ.
だれか今日太陽を見たものはおるか?

これは日出の時刻を知ろうとしたものらしい.

「冬物語」の4幕4場, ごろつきのオートリカスが道化に, お前は酷い目に会うぞと驚かす場面 では:
オートリカス 一年中で一番暑い日と暦に出ている日に煉瓦塀にもたせかけられ, 太陽の カッと燃える目に見つめられながら, ...

7月22日に大暑と書いてあるようなことをオートリカスはいったのだろう か. 理科年表の「気温の最高および最低記録」のようなものを指すのか.
ついでだが, 昔の人は時刻を何で知ったか? ヨーロッパの町にいると, 教会の 鐘の音で「ヨーロッパに来たなぁ」と思う. でも教会から遠い戦場では, 太陽の 方向以外に時刻を知るのは難しかったであろう.

ロミオとジュリエット 1幕1場:
ベンヴォーリオ おはよう, ロミオ.
ロミオ 早いのか, まだ?
ベンヴォーリオ 九時だよ.

近くの教会の鐘が鳴ったか.

再び暦に戻って:
「マクベス」 4幕1場, マクベスといい合った魔女や幻影が消えた後:
マクベス どうした? 行ってしまったか? 畜生, この忌ましい時め, 暦の つづくかぎり呪われるがいい!

自然の「天地の続く限り」ではなく, 人工の「暦の続く限り」というのが 意外な発想である. Shakespeareの暦についての考えが少しは分ったように思う.

月齢カレンダー

私がこの数年愛用している予定表の話はこのブログに何回か書いた. ところで, 今年の予定表の表紙にある週日定数も月齢定数も共に0である. つまり下のようなのだ.

そうするとこの次にこうなるのは何年後かも気になる.

y年の週日定数, week(y)と月齢定数, moon(y)は次の式で計算する.

def week(y): return (y-2000+(y-2000)//4 
    +4)%7
def moon(y): return (((y-11)%19)*11)%29

そこで次に(0,0)になる年を計算してみたら, なんと2272年; 247サロス周期の後であった. 今年は昭和100年, 戦後80年, 阪神淡路大震災30年といわれるが, 暦屋にとっても誠に貴重な年なのである.

この表を見た序でに, 今年の復活祭を調べよう. まず3月21日の月齢は21. 次の満月は月齢45として4月14日. その日の曜日は1で月曜だから, 日曜はその6日後で4月20日になる.

もうひとつ序でに. Winning WaysのDoomsdayは3月0日の曜日なので, 今年は5, つまり金曜日である.

Heronの式

前回のブログ「四面体の体積」に名称だけ出てきたHeronの式というのがある. 三角形の3辺の長さをそれぞれa, b, cとし, (a+b+c)/2をSとすると, その三角形の 面積は√S(S-a)(S-b)(S-c)というのである.

たとえば, 辺の長さが3,4,5のPythagoras三角形ではS=6だから, 平方根号の中は 6×3×2×1=36で, 面積は6だ. 右の図, 1辺が2の正三角形では, S=3で, 根号内は3×1×1×1=3で, 面積は√3である.

Heronの式を覚えた時, これが正しいという証明を見たかどうかは分らない. 勿論, Sの次元は長さ, S-a,S-b,S-cも長さなので, 根号内の次元は長さの4乗; 開平すると 2乗になり, つまり面積だ. そういえば, 前回のブログにあった, 四面体の体積を6本の 辺の長さから求める式も体積の次元になっている. あたりまえだ.

でも定量的に合っているかを調べるのが今回の頭の体操である. 方針としては, 根号の あるのがウンシャンだから, 面積の二乗の値を考えることにする. 例題の三角形は こういう頂点名と辺名を持つ. 計算の途中で時々数値的にチェックするのにもこれを 使う. この三角形は∠Aが直角なので, 面積は√8×√18/2=6だ.

式を使わざるを得ないから, Texの出力のpdfにも助けてもらおう.

式(0)はHeronの式で, 2乗すると(1). Sを定義で戻し, 展開したのが式(2)である. この程度 の計算はWolframAlphaにやって貰う.

三角形ABCの各頂点の対辺をa,b,cとすると, その各々の2乗もpdfにある通り. これを 式(2)の分子に入れて展開するわけだが, これは中々面倒である. WolframAlphaも 時間がかかるとか文句をいうし, 我が家のRapsberry PiのMathematicaは今は ネットに繋っていないので, 使い勝手が悪い. 結局は多項式の乗算のプログラムを Schemeでちょこっと作り, それで展開した. 時々例題の値を入れて験算する.

そして展開した結果は項数21のpdfにある式である. これも目がくらくらするから, x_ix_jy_ky_lについて, ijと klの対の項の係数を表にしてみた. さらに例題の座標で験算すると576. 16で 割ると36だ. しめしめ.

さてこの式は本当に三角形の面積であろうか. 頂点の座標(x₀y₀), (x₁y₁),(x₂y₂)の三角形の 面積(頂点の回り方によってはその負数)の式は, 幸いにも私は高校生の頃から知っている. 次のpdfにある式だ. 行列式を計算式にしてその分子の2乗を展開すると, またもや項数21の多項式が得られた. 今回もその係数の表を作るとブラボー! 先の表の1/4の値である. という訳けで Heronの式は面積になるのであった.

だが待てよ. この行列式はなぜ正しいのか. こちらも気になったので, 計算してみた.

図の三角形の, 座標軸に平行な線での外接長方形を作る. それから長方形の辺と 三角形の辺で作る3個の三角形の面積を引く. そういう式を作ってみる.

WoflramAlpha様に

Expand[(x2-x1)(y0-y2)-(x0-x1)(y0-y1)/2-(x2-x1)(y1-y2)/2-(x2-x0)(y0-y2)/2]

と伺いを立てたら

(x0y1)/2-(x0y2)/2-(x1y0)/2+(x1y2)/2+(x2 y0)/2-(x2y1)/2

のお告げが出た. 整理すると

(x0y1+x1y2+x2y0-x0y2-x1y0-x2y1)/2

となり, 行列式の値と一致している. もう一度ブラボー! ちゃんとした証明でもないが, 私にとっては十分納得出来る頭の体操であった.

四面体の体積

筑波大の三谷さんのツィッターに四面体の体積を計算する話があった.

角封筒の上の辺の中央と下の辺の両端を結ぶ折目をつけ, それを山折りにすると, 合同な 二等辺三角形4枚による四面体が出来る. 元の封筒の横と縦が1と√2の時の 体積を計算せよという問題である. ちょっとやってみたところ, 存外簡単であった.

下の図で, 赤い線は元の封筒を展開したもの. 下の辺と上の辺は繋っている. この斜めの 線を谷折りにすると, ABCを底面, ACD'を手前の面, ABDを向うの面, 辺DEとD'E'は 繋っていて斜面だが上の面になる.

錐体の体積を計算するには, 底面積と高さが必要だが, 底面積は簡単だ. 辺BCが1, 辺OAは辺ABが√2だから, √7/2である.

高さについては, 三角形ABDは, 辺ABを折目として頂点Dが起き上がるから, Dに対応する 平面上の点はABと直交する破線で示す線に沿って下がり, AO上のGまで来ると下から上って来たD’と 出会って頂点になる. GAの距離はADの距離にtan αを掛けたものだが, tan αはOB/OAだしADは1だから, GA=1/√7である. よって 高さは√6/√7.

求める体積は (√7/2)*(√6)/((√7)*2*3)=1/2√6

三谷さんのツィッターにはヒントと称して, 四面体の6辺の長さからその体積を 求める式があった. 平面三角形のHeronの公式の三次元版である. そこにある式の 添字の付け方を, 私流に修正したのが次の図である.

このように添字を付けると, 四面体の体積Vの2乗は

の一番上の式で得られる. しかしこれも目がくらくらするので, 始めの3項と終りの4項を それぞれ6変数のtfと3変数のtgに分け, 2行目以下のように定義する. 最初のV²を計算する式を書いたLatexは次の様だ. 計算式と同様に分解してある.

\newcommand{\sq}[1]{l_#1^2}
\newcommand{\tf}[6]{\sq#1\sq#2(\sq#3+
\sq#4+\sq#5+\sq#6-\sq#1-\sq#2)}
\newcommand{\tg}[3]{-\sq#1\sq#2\sq#3}
$V^2=\frac{1}{144}[\tf012345+\tf234501\\
+\tf450123\tg025\tg034\tg124\tg135]$

これを使ってpythonで計算したのが次だ.

import math
def sq(x):
  return x*x
def tf(a,b,c,d,e,f):
  return sq(a)*sq(b)
    *(sq(c)+sq(d)+sq(e)+sq(f)-sq(a)-sq(b))
def tg(a,b,c):
  return -sq(a)*sq(b)*sq(c)
def v2(l0,l1,l2,l3,l4,l5):
  return (tf(l0,l1,l2,l3,l4,l5)+tf(l2,l3,l4,l5,l0,l1)
    +tf(l4,l5,l0,l1,l2,l3)
    +tg(l0,l2,l5)+tg(l0,l3,l4)+tg(l1,l2,l4)
    +tg(l1,l3,l5))/144
l0=1;l1=1;l2=math.sqrt(2);l3=math.sqrt(2)
l4=math.sqrt(2);l5=math.sqrt(2)
print(math.sqrt(v2(l0,l1,l2,l3,l4,l5)))

=> 0.20412414523193143

この値は確かに1/2√6 である.

ところで, 最初の封筒の四面体の体積は, 形状が特殊であったせいか, 簡単であった. 底面の3点の座標と, その各点から頂点までの距離が与えられている時にも体積は 得られるだろうか. 以下はその計算である. 式が沢山あるので, Latexでpdfに したもので示す. 底辺の3点, A, B, Cの座標をそれぞれ(x₀,y₀,0), (x₁,y₁,0),(x₂,y₂,0); 頂点Dの座標を(x,y,z)とする. A,B,Cから頂点までの距離をそれぞれ l₀,l₁,l₂とする.

今回もpythonで計算すると

x0=0;y0=-1/2;x1=math.sqrt(7)/2;y1=0;x2=0;y2=1/2
l0=math.sqrt(2);l1=1;l2=math.sqrt(2)
a0=2*(x0-x1); b0=2*(y0-y1)
c0=sq(l1)-sq(l0)+sq(x0)-sq(x1)+sq(y0)-sq(y1)
a1=2*(x1-x2); b1=2*(y1-y2)
c1=sq(l2)-sq(l1)+sq(x1)-sq(x2)+sq(y1)-sq(y2)

x=(c0*b1-c1*b0)/(a0*b1-a1*b0)
y=(a0*c1-a1*c0)/(a0*b1-a1*b0)

z=math.sqrt(sq(l0)-sq(x0-x)-sq(y0-y))
print(x,y,z)

#=> 0.9449111825230684 0.0 0.9258200997725514

期待通りの値が得られた. でもこんな実数の値より, 1/2√6 の方が嬉しい.

Fixed Day Number

前回のブログのアルゴリズムに, 年初から前月末までの総日数を計算する式:

(367 * month - 362) // 12

があった. monthを1から12まで変えて値を計算すると.

[0, 31, 61, 92, 122, 153, 183, 214, 245, 275, 306, 336]

が得られる. 演算子「//」で小数点以下を落しているが, その落された値を調べると, ( (4777 - 367 * month) // 12 の図も右に示す.)

のような図になる. どちらも2月は30日あるとする.

まず左の図に注目する. 1月から2月は下がり, 2月から3月は上がり, 4月へは下がり, ... 8月と9月へは下がり, ... つまり大の月の後は下がり, 小の月の後は上がるらしい.

計算して見ると下がる量は0.4167, 上がる量は0.5833である.

もう一度式を見ると, monthが1増えると, 367 / 12 つまり 30.5833 増える. これと30及び31の差であった. これで納得できるが, Calendrical Calculation にはこういう説明があった.

(7 * m - 2) // 12 + 30 * (m - 1)

この第2項 30 * (m - 1) はすべて小の月として前月末までの総日数.

第1項はm = 1, 2, ..., 12について,

0, 1, 1, 2, 2, 3, 3, 4, 5, 5, 6, 6
  31 30 31 30 31 30 31 31 30 31 30

である(上の行), 下の行は1月から11月までの日数で, こうして見ると, 大の月の 上で上の行の値は1増える. つまり小の月とした時の誤差を得ている.

そこで上の第1項と第2項を足すと, (367 * m - 362) // 12になるのであった.

2月を30日にしたという発想がよかったのだ. 一方, 年末から計算する方の式はこう導く.

欲しい値は1月, 2月, ..., 12月で, 2月も30日とするから

[367, 336, 306, 275, 245, 214, 184, 153, 122, 92, 61, 31]

である. 30の倍数の列

[360, 330, 300, 270, 240, 210, 180, 150, 120, 90, 60, 30]

に誤差の列

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

を足したものだ. これは7から年初からの計算の誤差を引いたのもだ. つまり,
7-(7 * m - 2) // 12
でこれを(13 - m) * 30 に足す. これを一纏めにしたい.

まず下の計算の1行目のように, 床関数を天井関数にする. そして1を足して床関数に戻す. しかしmの値によっては, 整数になることがあり, 1を足すと失敗する. 従って, 2行目のように, 1の代りに 11/12を足してみる. 3行目はそれに30日の倍数を足したところ. 最後に分数の外に ある部分を分子に足す. すると(4777-367m)//12が得られる.

どうだろうか.

Fixed Day Number

私の好きな本に, Edward M. ReingoldとNachum DershowitzのCalendrical Caluculation がある. 随分昔に購入したもので, 手元にあるのは2004年のThird Printing.

要するに暦に関するアルゴリズムが満載で面白い. この本では暦の変換の基準は, 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月とか１３月が与え られても計算してしまう恐れがある. 表なら範囲外としてエラーに出来る. それを 心配したとしても, 日の値に範囲外が与えられる可能性はあるのだから, まぁ我慢することに しよう.