満月の高さ

子供の頃から, 夏の月は低く, 冬の月は高いと知っていた. 太陽について いえば, 夏の太陽は高く, 冬の太陽は低い. なぜ反対なのか. 私は小学校 高学年のころ, 原田三夫「子供の天文学」が愛読書で, 天文学の基本的な ことは知っていたが, 月の高さについてはよく分らなかった. 直感的にか 出任せにかさかしらにかは知らないが, 理由は「太陽を見るのは昼で, 月 を見るのは夜でだから」と説明していた.

それから何十年の経った今から3年前の学士会会報(No.962, 2023-V)に, 国立天文台(当時, 現在は京都産業大学神山宇宙科学研究所長)の渡部潤一 さんが「中秋の名月」なる一文を寄稿されていて, それに下のような図が あった. (この図は学士会会報の図を見ながら, 私がPostScriptで描き直 した.) 図ではしっかりと冬至の満月は高く, 夏至の満月は低いと説明さ れている. そうなると, 私もこのことを確信したくなった. これまできち んと調べなかったのは, 月の軌道である白道が, 太陽の軌道である黄道と 交差していることから, 月の位置は捉え処がなく, どう動くか調べようも ないだろうと思っていたところ, 確信を以って「夏至の満月は低い」とあっ たので, 自分でも調べてみる気分になった.

手初めに, 天球に於ける月の軌道を描くことにした. 月の位置の赤経赤緯 の表は, 理科年表や天文年鑑にあるが, 紙に印刷してあるから取り出して データとして使うことは出来ない. しかしネット上には国立天文台の 「暦象年表」のページがあり, そこから2026年の毎日の太陽と月の赤経赤 緯を取り込むことが出来る. これを利用して, 太陽と月の軌道を図に描く ことにした.

暦象年表とは, 理科年表の暦部のことである. 私も詳しいことは覚えてい ないが, 終戦の頃やその後の数年は理科年表は刊行されず, 暦象年表が市 販された. 理科年表は暦部の他は毎年殆んど同じ内容なので, 必要なのは 暦象年表だけである. しかし数年経つと, 理科年表が刊行されるようにな り, 暦象年表は書店から姿を消した.

ところがある時, どこかの天文台で, 思いがけなく暦象年表を見付た. 一 般への市販はなくなったが, 天文屋さんのためには出版されていたらしい. そうこうするうちに, 暦象年表はインターネットで見られるようになっ た. ありがたいことである. (有料だが, 理科年表もインターネットで見 ることが出来る.)

この図は時間とは無関係に, 黄道(赤の曲線)と白道(青の曲線)の天球にお ける位置の動きを示す. 赤線の右端は, 黄道が赤経0hの時(春分点の時)の 赤緯は0°で, 右から1/4は, 黄道が赤経6hの時(夏至の時)の赤緯は +23.5°であるのように見る. 白道についても同じ. ただ白道の線は, 14本 くらいが殆んど重なっている.

曲線上にある点は, 暦の各月の満月の日に於ける赤経と赤緯の位置を示す. 5が2回あるのは, 5月に満月が２回あったからだ. これを見ると, 同じ暦 の月の満月は太陽と丁度12時間ずれていることも分る.

これで判明したのは, 交差角が5°とはいえ, 白道は殆んど黄道と一致して いることだった.

この図からの結論をいえば, 太陽と新月の月は, その赤経と赤緯は同じで, 満月の月の赤経は12hずれ, 赤緯は南北反対だが, 絶対値は同じというこ とである. 従って夏至の満月は, 太陽が天の赤道の上, 23.5°の高いとこ ろで南中するのに対し, 月は天の赤道の下, 23.5°低いところで南中する. 冬至はその反対である. 目出度し目出度し.

今回の検討は, 私にもうひとつの宿題の解答も与えた. それは1960年かそ の少し後だっと思うが, 物理学科の何年か後輩の佐々木不可止君がスェー デンの大学の留学を終えて帰国し, 私のところへ来て雑談していた時のこ とだ. 佐々木君はこういった. 「スェーデンのような高緯度の所では, 月 が北の空に出る.」この話はそれだけで, 私は北極圏では白夜に太陽が沈 まないように, 月も地平線に沿って回るからだろうと考えたが, それ以上 はまたいつか考えようという程度であった.

このブログの最初にあった図は, 東京あたりの北緯35°くらいを想定して 描いてある. それより赤道に近い方では, 地軸がもっと寝て, 満月の軌道 はもっと起きる. 反対に極に近くを考えると, 地軸は立って, 満月の軌道 は水平になり, 春分・秋分の満月の軌道は水平線のところに来る. 冬至の 満月は地上23.5°のあたりで東→南→西→北→と周回す るようになるから, 北の空を通ることもあるわけだ.

佐々木君は大学院が終わると, 北海道大学に就職し, 理学部化学科の教授 だった. 北大を定年で退職し, 他の私大に移ったが, 間も無く病没し た.

佐々木君は学部生の頃にも高橋研に出入りし, パラメトロン回路にも興味 をもっていた. 佐々木君の考えた閾値の論理回路の驚くべき式は, Knuth がThe Art of Computer Programming, 第4A巻の7.1.15項の「対称的Boole 関数」に紹介した.

マージソート

昨今のことは存ぜぬが, かつての計算機科学科での必須の科目は, 計算機 のアーキテクチャとアルゴリズムであった. しかし最近のパソコンの中は 当時のアーキテクチャとは掛け離れて複雑になり, プログラムを書いて計 算するよりも, 情報端末として使われている. アルゴリズムを勉強しても 使う機会はない.

ところで突然マージソートのプログラムを書くことになった. マージソート に関しては大体のことは理解している心算りであったが, TAOCPの5.2.4項 を眺めたら, 意外と巧妙なことをやっているらしい. しかし, TAOCPのプ ログラムは例によって自然言語で書いてあり, 読む気にもならない. 私と しては, 再帰プログラムでさっさと書けるのではないかと思い立ち, やっ てみたのがこのブログである.

とりさえずマージソートの核心部分のおさらいから. 下の図にはA, B, Cの 3つの図がある, 似たような絵はTAOCPにあり, データの数値もそこから 借用した.

まずAを見よう. 左方のs0とs1は既にソート済みの数列である. どちらも 左の小さい値のある方が先頭である. これをマージして右方の数列t0を作る. まず先頭の503と83を比べ, 小さい方, 87をs1から外して, それまで空の列で あったt0に移す. s1は512が先頭になる. 再び先頭同士を比べ, 小さい方, 503をt0へ移す, 今度は703と512を比べ, 512がt0へ移る. 次はs1の 677がt1へ行き, s1は空になる. このように比べている数列の一方が空に なったら, もう一方の残りの数列をすべてt0へ移す. するとs0, s1は共に 空になり, すべての数はソートされてt0へ移動した. マージソートは元の 数列のソートされたと見える部分列を2つずつマージし, 最後は1つの数列に なるまで繰り返す.

その手順の最初が図Bで, 最後が図Cである. まずBから. 上の列が与えられた ソートすべきデータ列である. 下の列は部分的にマージされた列で, 元の数列と同じ 長さの領域を別に用意する.

マージされる2つの数列は, 上の数列の両端から, s0とs1で示すように内側に向って 段々大きくなるように並んでいると考える. ソートされた列の終は次の数値が今の数値より小さく なるところである. 図では縦線で示した. 上の数列の左からをs0, 右からをs1とし, これをマージして用意した領域の左からt0として詰める.

s0とs1とでt0が出来たら, 続いてs0の内側のs2と, s1の内側のs3をマージし, t0を作った領域の, 今度は右端からt1として詰めていく. こうして行くとそのうち, 比べる べき数が同じ場所になり, tの領域も中程に1個所, 空地が残った状態になり, その残った数を空地に移すとこの回のマージは終了する.

次はこうして得られたtの数列を出発側にし, また同じようなマージを繰り返す. その 内, 図Cのように, 同じ数(この図では中程の908)を比べる状態になり, tの 領域の唯一の空地が右端に残り, そこに908を移動してソートは終了する.

後はこれを正直にプログラムにすることだ. Schemeで書いたのがこれである. 関数名がalgorithm524nなのは, TAOCPのAlgorithm524Nを書き換えたから である. 引数のrはソートすべき数列.

すぐにtwowaymergesortというサブルーチンがある. これは図BやCの1回のsからtへの マージである. 引数のrsはs(source)の数列. その長さがnで, その長さのt(target)の 数列の場所. Schemeではリストの要素を取るのに, list-refのように長い関数名を 書かなくてはならないので, 数列rsの添字iを取り出す関数rや, sのi番目をtのl 番目に移動するmoveなど, よく使う関数(手続き)を用意する. (MIT Scheme にはどういう訳かlist-set!がないので, 自分で定義して使っている.)

let文は次に使うsの左の場所の添字i, 右の添字j, tの左の添字k, 右の添字lを 初期化する.

次にtoleftとtorightがあるが, 新しい領域の左にあるtへ移すか右にあるtの 移すかで使い分ける. その夫々にfromleftとfromrightがある. これは 比較の一方が空になった後, もう一方の残りをどちらのsから移動するか で使い分ける.

(define (algorithm524n r)
 (define (twowaymergesort rs)
 (let* ((n (length rs))(rt (make-list n '())))
  (define (r i) (list-ref rs i))
  (define (move l i) (list-set! rt l (r i)))
  (define (less? a b) (let ((rsa (r a)) (rsb (r b)))
    (< rsa rsb)))
  (let ((i 0) (j (- n 1)) (k 0) (l (- n 1)))
    (define (toleft)
      (define (fromright) 
        (move k j)(set! j (- j 1))(set! k (+ k 1))
        (if (less? j (+ j 1))(toright)(fromright)))
      (define (fromleft)
        (move k i)(set! i (+ i 1))(set! k (+ k 1))
        (if (less? (- i 1) i)(fromleft)(toright)))
      (if (= i j) (begin (move k i)
        (if (< k (- n 1))(twowaymergesort rt) rt))
        (if (less? i j) 
          (begin (move k i)
            (set! i (+ i 1)) (set! k (+  k 1))
	    (if (less? (- i 1) i)
	      (toleft) (fromright)))
	  (begin (move k j)
	    (set! j (- j 1)) (set! k (+ k 1))
	    (if (less? j (+ j 1))
	      (fromleft) (toleft))))))
    (define (toright) 
      (define (fromright) 
        (move l j)(set! j (- j 1))(set! l (- l 1))
        (if (less? j (+ j 1))(toleft)(fromright)))
      (define (fromleft) 
        (move l i)(set! i (+ i 1))(set! l (- l 1))
        (if (less? (- i 1) i)(fromleft)(toleft)))
          (if (= i j) (begin (move l i)
	    (twowaymergesort rt))
      (if (less? i j)
	  (begin (move l i)
	    (set! i (+ i 1)) (set! l (- l 1))
	    (if (less? (- i 1) i)
	       (toright) (fromright)))
	  (begin (move l j)
	    (set! j (- j 1)) (set! l (- l 1))
	    (if (less? j (+ j 1))
	      (fromleft) (toright))))))
    (toleft))))
(twowaymergesort r))

(algorithm524n
  '(503  87 512  61 908 170 897 275
    653 426 154 509 612 677 765 703))

TAOCPには, merge sortは1945年頃, John von Neumannが提案したと 書いてある.

干支の計算

「今年はうま年だ」のように, なに年かが気になるのは, 現代人では正月 の頃だけであろう. また「初午」, 「土用の丑の日」, 「三の酉」のよう に日付けにも十二支が対応していることも常識であろう. 十二支は曜日の ように, 正確に順番に対応して来たので, 利用できる場面が多い.

歴史上の事件ではその年の干支を付けて, 「壬申の乱」とか「辛亥革命」 とか, 干支で識別されることがある.

しかし最近は「今年の干支はひのえうま」や, 「庚申(こうしん)の日」と か十干の方まで気にする人は絶滅危惧種かも知れない.

絶滅寸前かも知れないが, 私は時々「この前の甲子(きのえね)は何年だった か」のように, 年の干支が知りたくなることがある. また青木信仰「時と 暦」によると, 日本書紀には漏刻(水を使った時計)の使用開始が671年の「夏 四月丁卯朔辛卯」と書いてあるそうで, これも日付けにしてみたい.

干支は表Aの十干(天干ともいう)と, 表Bの十二支(地支ともいう)を順に組 み合わせたものである. 10と12は互いに素ではなく, 最大公約数が2, 最 小公倍数が60なので, 組合せは表Cの様に60個になる.

十干の方は「こう, おつ, へい, てい」位までは音読みで知っている人は いるであろう. 十二支は訓で読むと「ね, うし, とら, ..., いぬ, い」で, 大 体の人は知っているに違いない.

百人一首には喜撰法師の「わが庵は都のたつみしかぞ住む 世をうぢ山と 人はいふなり」があり, それをもじって「わが庵は都のたつ, み, うま, ひつじ, さる, とり, いぬ, い, ね,うし, とら, う」と本歌取りした人 がいた.

しかし干支の組合せは音読みの方が普通なので, この機会に表の音読み を覚えのはどうか. 次の音読みは旧かな遣いで書いてあり, 実際には 「甲」, 「庚」は「こう」,「丑」は「ちゅう」, 「卯」は「ぼう」と読 む.

さらに十干や十二支の番号も知っていると役立つ. 十干は「こう, おつ」, 「へい, てい」の ように2つずつ組にして覚えると, その組が「き, ひ, つち, か, みづ」に対応する. また十二支は「し, ちう, いん, ばう」, 「しん, し, ご, び」, 「しん, ゆう, じゅつ, がい」 と4つずつにすると, それぞれの組の先頭の年がうるう年になるので都合がよい.

さてこれからが表題の「干支の計算」である. 表Cの干支の一覧が 手元にあり, 今年の干支が「丙午」であることを知っていれば, その前の「甲子」を知るには, 2026年の「丙午」 がn=42で, 「甲子」はn=0だから, 2026-42=1984と判る.

漏刻の使用開始は, 朔(4月1日)の干支, 丁卯がn=3で, 開始日の 干支, 辛卯がn=27だから, 1+(27-3)=25日と判る.

しかし, 表Cのような一覧表なしで, 干支の番号(0<=n<60)を知りたい時 はどうするか. もちろん, ここでは十干の番号(0<=s<10)と十二支の番号 (0<=b<12)は判っているとする. つまり「丙」のs=2, 「午」のb=6 から「丙午」のn=42を知る方法の話である.

下の図を見て欲しい. 縦に3つの図がある. 上の図は, 横軸のnに対して, 赤で示すs, 青で示すbの対応を示す. 赤線, 青線の左端の塗り潰しの 丸は閉区間, 白抜きの丸は開区間を示す. つまりn=10 の時のsは一番左の赤線の値ではなく, その次の赤線の値, 0であることを 示す.

中の図は, s-bを示す. 10<=n<12でs-bが負になっているが, これを 12だけ上に移動すると10になり, 12<=n<20の線と一緒になって, 10<=n<20の区間でs-bの値が2になる. 他の負の値の場所も同様になる. それが下の図である. これで判明したことは, sとbからnを作るに は, この値の2を5倍してsを足せばよい. Schemeで書くと

(define (nval s b)
  (let ((d (- s b)))
    (+ s (* 5 (if (< d 0) (+ d 12) d))))))

丙午の場合(s=2, b=6)なら

  (nval 2 6) => 42

Schemeの関数(modulo n d)は, 0<=r<dの剰余を返すから

(define (nval s b)
  (+ s (* 5 (modulo (- s b) 12))))

と書ける.

Schemeでなくても, Pythonのように負の分子に正の剰余を返す言語なら

def nval(s,b):
    return s+5*((s-b)%12)

でもこの方式だと, 丁卯(s=3, b=3)と辛卯 (s=7, b=3)の差は, それぞれのnを 計算してから差をとることになる. sとbのそれぞれの 差からnの差が得られると嬉しい. そこでそれぞれの差分dsとdbの差とnの差の 関係を作ってみると, 次の表のようであり, これは何と, sと bからnを得るのと同じであった. 考えてみるとまぁ そうなるのは当然である.

ds=7-3=4, db=3-3=0 だから

4+5*((4-0)%12)=24

である.

これで今回の話題は完了した. 以下は蛇足である.

いつか書いたかも知れないが, グレゴリオ暦の曜日を知るには, その日のユリウス日JDに 1を足して7で割った剰余をとると, (日曜を0とする)曜日の番号を得る. この方式が干支に 使えないかとやってみると, JDから1を引いて10で割った 剰余は(甲を0とする)十干の番号になる. さらにJDに1を足して12で割った剰余は(子を0 とする)十二支の番号になる. 2025年12月21日は日曜で甲子であった. その日のJDは 2461031. それを7, 10, 12で割った剰余は6, 1, 11である. JDから干支を 計算するには, s=(JD-1)mod 10. b=(JD+1)mod 12から, このブログの式 n=s+(s-b+12)mod 12 *5で 計算してもよいが, JDを60で割った剰余が11なので, JD-11を60で割るのが 簡単である.

整理すると
曜日 (JD+1)%7
十干 (JD-1)%10
十二支 (JD+1)%12
干支 (JD-11)%60
だけ覚えればよい. 偶然だが足したり引いたるがすべて1なのが覚え易い. 私のiPhone には, 個人用電卓があり, JDを計算する機能もあるから, この情報は貴重である.

干支とは関係ないが, 暦の最重要な計算は復活祭である. puzzle and paradoxにある アルゴリズムでは, 今年は4月5日. ということは, ConwayのいうDoomsdayは今年は 土曜日だ. この話題はまたいつか.

曜日の計算

我が家の書斎にのJohn Conwayへのインタビュー(Scientific American, 1999年4月号)の古いコピーがあった. それによると, Conwayは日付から曜 日を算出するのが得意であったらしい. もちろんそれには練習もした. 計 算機にログインした時, ランダムに日付を提示するようにし, その曜日を 暗算して練習したという. 2秒以内で正解が得られるようになったと書い てあったから驚きである.

立教大学の島内先生も曜日の計算は素早く, その話はいつぞやのこのブロ グにすでに記述した通りである.

私も曜日の高速計算を隠し芸にしたいが, どうすればいいか. いつもはブ ログにも何度か書いた例の 2,5,5,...の表を使う. 当日と同じ年内のことだ と, Doomsday方式も便利で, 今年(2025年)は「金曜」を覚えておくだけで よい. Conwayの曜日計算はDoomsdayの方式に違いない. その辺のことは, Winning Waysの第4巻にある.

問題はその年の定数, 若しくはDoomsdayを得ることにある. 私のは島内さ んと同様で, 西暦年数の上2桁をc, 下2桁をbとした時, (b+b//4)にcが19 なら2を引き, 20なら4を足して7の法をとっていた. 例えば2025年なら, c=20, b=25だから(25+floor(25/4)+4)%7=(25+6+4)%7=0. しかしこれはb が大きいと, 大数を扱うことになり面倒だ.

Winning WaysにあるConwayの方法は意外である. bを12で切り捨て除算し, (b=25ならb//12=2)その商(この例なら2)に剰余(同1)を足し, さらにその 剰余を4で切り捨て除算した商(同0)を足すものである((2+1+0)%7=3 それ にc=20の4を足して7で法をとるから0).

私のやり方と随分違う算法に見えるが, 取り敢えず私にやり方(Aとする) とConwayの方法(Bとする)で, bの0から35までの, 7で法をとる前の値を計 算し, 12個ずつ並べてみると

A  0  1  2  3  5  6  7  8 10 11 12 13
  15 16 17 18 20 21 22 23 25 26 27 28
  30 31 32 33 35 36 37 38 40 41 42 43 
  ...
B  0  1  2  3  5  6  7  8 10 11 12 13
   1  2  3  4  6  7  8  9 11 12 13 14
   2  3  4  5  7  8  9 10 12 13 14 15
  ...

つまりAでは2行目, 3行目は1行目に15,30を加えたもにであるのに対し, B では1,2を加えたものになっている. 15,30は7を法とすると1,2だから, 7 で法をとった結果は同じになるのであった. Aではbが99だと99+24=123まで 大きくなるのに, Bだと8+3=11にしかならない.

というわけで, また我々には12の倍数なら8倍くらいまでは常識だから, 私も最近はBを愛用している. しかし, 12で割った剰余に, さらにそれを4 で切り捨て除算をしたものを足すのは, 計算するよりそれの7での法をとっ た

0  1  2  3  5  6  0  1  3  4  5  6

を足す方が簡単のように思い, 実際の計算では表を使う.

従って, グレゴリオ暦のc*100+b年m月d日の曜日は,

西暦年数の上2桁cの定数: [4,2,0,-2][c % 4]
西暦年数の下2桁bの定数: 
  b//12+[0,1,2,3,5,6,0,1,3,4,5,6][b%12]
月mの定数: 
  [2,5,5,8,3,6,1,4,7,2,5,0][m-1] 
  ただし閏年の1,2月は1,4
日dの定数; d

これ(上からそれぞれをc項, b項, m項, d項とする)を加算して7の法をと るという計算になる. cに関する項は実用的にはcは19か20なので, -2と4 だけを覚えておけばよい.

また, bの項で使う表をbs, mの項でのそれをmsということにする.
この12要素の2つの数表bsとmsを覚えるのはそれ程面倒ではない. どちら も簡単な規則で出来ているから, すぐに構築できる. (msでは4月と9月で 7の法を取る前の値にしてある.)

ちょっと練習問題を. グレゴリオ暦になった最初の日, 1582年10月15日. ものの本には金曜とある.

c=15だからcの項は-2.
12*6=72だから, 82//12=6 82%12=10, bの項は6+5=11
10月だからmの項は2
-2+11+2+15=26 26%7=5で金曜であった. 

母親が水曜であることだけ覚えていたGaussの誕生日 1777年4月30日.

c=17だからcの項は2
77//12=6 77%12=5 bs[5]=6,
4月のmsの値は8
2+6+6+8+30=52 52%7=3 で水曜であった.

曜日の高速計算のキーポイントはこの2つの表引きが瞬時に出来ることである.

bsの表では12の値は4つずつ3組になっていて, 左の4個は見出しを二進法 にした時, 4のビットも8のビットも立っておらず, この組は0,1,2,3につ いて値も0,1,2,3だから簡単である. 中央の4個の組は4のビットが立って おり, 4,5は1ずつ多い5,6; 6,7は0,1, 右の4個は8のビットがあり, 値は5 ずつ少ない.

msの表についていえば, カレンダーの曜日は3ヶ月で大体元に戻ることを認識 するのが重要である. つまりその3ヶ月が, 30日, 31日, 30日で合計が 91日だと, 丁度7の倍数だから, 4ヶ月目は3ヶ月の最初と同じになる. 4月と 7月が同じ曜日関係になる.

90日だと曜日が1日戻り, 92日だと1日進む. 平年の1,2,3月は90日だから, 4月は1月より1日戻り, 4,5,6月は91日だから4月と7月は曜日が同じだ.

Doomsdayの規則には1月と2月への言及はないが, 私は1,2月は平年なら 10,11月と同じ, 閏年なら7,8月と同じと思っているが, 平年の1,2月の2,5 は10,11月の2,5と同じで月の大小も同じである. また閏年の1,2月の1,4は 7,8月の1,4なのと同じで,こういう事に気付いているとmsの表も覚え易い.

最近はこういう考えで曜日の高速計算に挑んでいるが, 私自身が歳を取ったことも 言い訳にして, 仲々高速にはならない. しかし結構楽しんではいる.

ユリウス日

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

ユリウス日

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

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

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

ユリウス日

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