2022年9月20日火曜日

満月の十五夜

今年は9月10日が中秋の名月, つまりお月見であった. 幸いよく晴れて, 満月が堪能出来た.

ところで新聞が「今年の中秋の名月は満月」と報じた. 私はもちろん旧暦の15日が満月にならない方が普通で, たまには満月になることも知っていたから, この新聞記事をみて何とも思わなかったが, どの程度すれすれに満月なのかと天文年鑑を見ると満月は18時59分であった.

旧暦の15日が満月にならないのは, 新月の時点のある日を旧暦の1日とするからである. 仮に1朔望月を29.5とし, 月齢14.75を満月とすれば, ある日の朝0.25日までに新月があればその日が1日で, 15日の晩には月齢が14.75に達する. そうでないと, 月齢が14.75になるのは, 旧暦16日になる.

ところが月の公転には遅速があり, 満月は太陽と月の黄経の差が180度の時のことだから, 月齢14.75からプラスマイナス1日くらいずれ得る.

そこで新月から満月, 満月から新月の経過時間の変化を見てみたいと思った. 2013年から2022年の天文年鑑から, 新月と満月の日時を書き出し, パソコンに入力して新月と次の満月, 満月と次の新月の経過時間を計算した.

10年間にある朔望月は, 19年7閏法を考えると, 120+3.5回だから, 日時データは247, 間隔は246あった. 最大値は15.601, 最小値は13.907, 平均値は14.767.

最後の10個は
... 14.32 15.278 14.071 15.497 13.958
15.579 14.009 15.497 14.216 15.256
である.

それを絵にしたのが下だ. 赤は新月から満月, 青は満月から新月までの経過時間である. 縦軸の単位は日.

図を見ると, 時々赤と青が同じになり, つまり新月から満月までと満月から新月までが同時間になり, その中間は一方は増えて減り, 他方が減って増える.

その理由は多分こうであろう. 月の公転軌道も楕円であり, 長軸の一方が近地点, 他方が遠地点である. 近地点では公転速度が速く, 遠地点では遲いのは常識だ.

近地点と遠地点の近くに新月と満月の場所があれば, 新月から満月も満月から新月もほぼ同時間になる. しかし, 新月と満月がその中間くらいにあると, 近地点側を通る方は経過時間が短かく, 遠地点側は長い. それがこの図の謎解きであろう.

2022年の新月と満月の日の月の地心距離を天文年鑑から調べ, 新月から満月までやその逆の経過時間がどの地心距離からどの地心距離であったかの図を描いてみた. それを下に示す. 横軸は経過時間(単位は日)で, 縦軸は地心距離(単位は万km).

これを見ると, 上の推論が正しいことが分る. 左上から右下へ来る矢印は, 遠地点から近地点までの経過で14.8日くらい掛ることで, 15.7日くらいの水平の矢印は, 中間から中間への遠地点側を通る経過, 14.0日くらいの水平の矢印は, 近地点側を通る経過である.

十五夜が満月になるというのには, こういう仕掛けがあったわけだ. その後, 国立天文台暦計算室のページを見付けた.

2022年9月8日木曜日

切頭八面体

久し振りに切頭八面体について書きたい. TAOCPにおやという記述があった.

a: 切頭八面体の基本的な置き方は, (htmlでは式が書き難いから言葉で述べると) 0,1,2,3の, 同じもののない3つ組である.

b: 体心立方格子のVolonoi領域である.

bはなるほどそうらしいと思うが, そしてそういう言い方もあるのかと思うが, aはすぐには分らなかった. その組を作ってみると,
((0 1 2) (0 1 3) (0 2 1) (0 2 3) (0 3 1) (0 3 2) 
(1 0 2) (1 0 3) (1 2 0) (1 2 3) (1 3 0) (1 3 2) 
(2 0 1) (2 0 3) (2 1 0) (2 1 3) (2 3 0) (2 3 1) 
(3 0 1) (3 0 2) (3 1 0) (3 1 2) (3 2 0) (3 2 1))
になり, たしかに24個あって, 切頭八面体の頂点の数と同じである.

それぞれの3つ組を(x, y, z)と思い, 高さ(z) 別に絵を描くと, 下の図の上段のようになる. z=0, つまり底の面に点が6個あるから, これは六角形を床に置いた図らしいと分る. 六角形にしてはいびつだが, それは直交 座標のためらしいから, x軸とy軸の角度を60度にしてもう一度描いたのが中段の図である.

なるほど六角形が現れる.

そこでx軸とz軸, y軸とz軸の角度も60度にしたグラフ用紙を作って, 上に点を 書き込むことにした. つまり, x軸y軸z軸は正四面体の3本の辺に沿っている. それが上の赤緑青の線である. 赤はz=0の平面座標で, その三角形の中心を起点とした緑のz=1の平面座標が出来, その三角形の中心を起点とした青の z=2の平面座標が出来る. 同様にしてz=3の座標が出来るが, それはz=0のと一致する.

黒線で描いたのが六角形を底にした切頭八面体である.

これを確認するため, 別の切頭八面体の絵もある. 左上Aは, 私が通常に描く切頭八面体 である. その中心を通る直交するx軸(赤), y軸(緑), z軸(青)も描いてある.

私の描画プログラムは, 頂点の座標を座標軸に平行移動したり, 座標軸に沿って回転 したり出来る. BはAの図を, z軸について45度回転したものだ. Cはそれを右下の 5,12の辺が, xyz軸の中心に来るように移動したもの. Dはそれを右下の六角形が 底に来るように, x軸について回転したものである. 最後のEは, 座標軸の中心に 切頭八面体の中心が来るように移動した.

この時, 各頂点の直交座標での位置は次の通りである. (左は頂点の番号順, 右は高さ別)
     0 (-2.121 0.408 -0.577)      9 (-1.414 0.0 1.732)	
     1 (-1.414 0.0 -1.732)	  11 (-0.707 1.225 1.732)	
     2 (-1.414 1.633 -0.577)	  16 (-0.707 -1.225 1.732)	
     3 (-0.707 1.225 -1.732)	  18 (0.707 1.225 1.732)	
     4 (-2.121 -0.408 0.577)	  20 (0.707 -1.225 1.732)	
     5 (-0.707 -1.225 -1.732)	  22 (1.414 0.0 1.732)	
     6 (-0.707 2.041 0.577)	  1 (-1.414 0.0 -1.732)	
     7 (0.707 1.225 -1.732)	  3 (-0.707 1.225 -1.732)	
     8 (-1.414 -1.633 0.577)	  5 (-0.707 -1.225 -1.732)	
     9 (-1.414 0.0 1.732)	  7 (0.707 1.225 -1.732)	
     10 (-0.707 -2.041 -0.577)	  12 (0.707 -1.225 -1.732)	
     11 (-0.707 1.225 1.732)	  14 (1.414 0.0 -1.732)	
     12 (0.707 -1.225 -1.732)	  4 (-2.121 -0.408 0.577)	
     13 (0.707 2.041 0.577)	  6 (-0.707 2.041 0.577)	
     14 (1.414 0.0 -1.732)	  8 (-1.414 -1.633 0.577)	
     15 (1.414 1.633 -0.577)	  13 (0.707 2.041 0.577)	
     16 (-0.707 -1.225 1.732)	  21 (1.414 -1.633 0.577)	
     17 (0.707 -2.041 -0.577)	  23 (2.121 -0.408 0.577)	
     18 (0.707 1.225 1.732)	  0 (-2.121 0.408 -0.577)	
     19 (2.121 0.408 -0.577)	  2 (-1.414 1.633 -0.577)	
     20 (0.707 -1.225 1.732)	  10 (-0.707 -2.041 -0.577)	
     21 (1.414 -1.633 0.577)	  15 (1.414 1.633 -0.577)	
     22 (1.414 0.0 1.732)	  17 (0.707 -2.041 -0.577)	
     23 (2.121 -0.408 0.577))	  19 (2.121 0.408 -0.577)    
各頂点の直交座標(x,y,z)を, 斜方座標(a,b,c)に変換するには, 例えば
(define (tcood x y z)
(let* ((c (/ (* z 3) (sqrt 6)))
       (b (/ (* (- y (/ (* c (sqrt 3)) 6)) 2) (sqrt 3)))
       (a (- x (/ b 2) (/ c 2))))
(list a b c)))
しかしこれで変換すると0.707や2.121が出るので, 正規化すると, 先程の24個の 直交座標は
((-3 1 -1) (-1 1 -3) (-3 3 -1) (-1 3 -3) (-3 -1 1) (1 -1 -3)
 (-3 3 1) (1 3 -3) (-1 -3 1) (-3 -1 3) (1 -3 -1) (-3 1 3)
 (3 -1 -3) (-1 3 1) (3 1 -3) (1 3 -1) (-1 -3 3) (3 -3 -1)
 (-1 1 3) (3 1 -1) (1 -3 3) (3 -3 1) (1 -1 3) (3 -1 1))
なるほど, 最初に書いた座標は0,1,2,3であったが, 今回のは. 切頭八面体の中心に座標の 中心を置いたので, -3,-1,1,3で構成されていた.

2022年9月6日火曜日

Wolframの対数表

学生の頃から気になっていたことがある. 高木貞治先生の「解析概論」のp.215の脚注に 「Wolframノ表ニハ1000以下ノ素数ノ自然対数ノ50桁ノ表ガ掲ゲラレテヰル. コノ表ハ既ニ少年がうすガ愛用シタモノデアル.」と書いてあったことだ. また, 岩波文庫の「近世数学史談」のp.65に 「Wolframの計算した10000以下の素数及び若干の特殊の数の自然対数の表が 附いていたのである.」と書いてあるのも見付けた. 1000までの素数か10000までの 素数かに違いはあるが, とにかくそういう対数表が存在していたらしい.

学生の頃にそんな数表を探すのは夢のまた夢であったが, 長生きはするもので, この歳に なってインターネットを検索していたら, 関係のあるウェブページがあった.

そのひとつは A reconstruction of the Mathematical Tables Project’s table of natural logarithmsという文献で, その図1は, Wolframが計算したものを, Schulzeが1778年(Gaussの生まれた翌年)に 出版したものの最後のページの写真. これを見て, 1000までの方が違っているのが分かった. しかも最後は10009で精度は48桁であった. また6個の素数には対数が記入 されていず, Wolframが病気で計算できなかったらしい. その図2はVegaが欠けている値を計算して 追加して1794年に出版したものの最後のページの写真がある. いずれにしても, Gaussが所持していたWolframの数表はこういうものであったらしい.

もうひとつは, New Information Concerning Isaac Wolfram's Life and Calculationsである. これによると, Isaac Wolframはオランダの砲兵の士官であった.

Wolframは何年もかけてこの数表を作ったらしいが, 現代は便利な世の中で, 自然対数も 簡単に得ることが出来る

そのひとつは, 別のWolframだが, WolframAlphaである. このページを開くと入力用窓があり, そこに「ln(2)」といれると, 2の自然対数が64桁くらい表示される. 0.6931471805599453094172321214581765680755001343602552541206800094 私もこの数表を作ってみたくなった. しかし. WolframAlphaは手動でしか使えないから, 別の方法が必要だ. あれこれ探していたら, Python3でも値が得られることが分った.
from decimal import *
getcontext().prec = 55
Decimal(2).ln() => 
Decimal('0.6931471805599453094172321214581765680755001343602552541')
これを使うと, Wolframの対数表みたいなものはすぐ作れる. 10000までの素数は 1229個あるが, 最初の20個の素数と, 10000未満の最後の19個の素数に対して 作た対数表はこんなようだ.
2から9973までの素数に対する自然対数の表は, 私の計算機にある.

2022年4月12日火曜日

mex関数

mex関数についてブログに全回書いたのは, もう2年も前だった. その頃, Schemeでの実装も 出来ていたが, なんとなく気に入らなかった.
最近, もう少し気持ちよくしたいと思い, やった書き直したので, ブログに記録して おく. まずinsertion. (ins n t)は, nを, 木tに挿入する関数である.
(define (ins n t)  ;inserts n in a tree t, returns a new tree
  (if (null? t) (list '() (list n n) '())
      (let ((lt (car t)) (l (caadr t)) (r (cadadr t))
            (rt (caddr t)))
  (define (insr t)
    (let ((lt (car t)) (m (caadr t)) (s (cadadr t))
          (rt (caddr t)))
      (if (null? rt)
        (if (= s (- n 1)) (begin (set! l m) lt)
            (begin (set! l n) t))
        (list lt (list m s) (insr rt)))))
  (define (insl t)
    (let ((lt (car t)) (m (caadr t)) (s (cadadr t))
          (rt (caddr t)))
      (if (null? lt)
        (if (= m (+ n 1)) (begin (set! r s) rt)
            (begin (set! r n) t))
        (list (insl lt) (list m s) rt))))
  (cond ((<= l n r) t)
        ((= n (- l 1))
          (if (null? lt) (list lt (list n r) rt)
              (let ((lt (insr lt)))
                (list lt (list l r) rt))))
        ((< n (- l 1)) (list (ins n lt) (list l r) rt))
        ((= n (+ r 1))
          (if (null? rt) (list lt (list l n) rt)
              (let ((rt (insl rt)))
                (list lt (list l r) rt))))
        ((> n (+ r 1)) (list lt (list l r) (ins n rt)))))))
一方, これを使う(mex ms)は次のようだ.
(define (mex ns) (let ((t '()))
 (for-each (lambda (n) (set! t (ins n t))) ns)
  (if (null? t) 0
      (begin (while (not (null? (car t))) (set! t (car t)))
             (if (< 0 (caadr t)) 0 (+ (cadadr t) 1))))))
全回と似たような図だが, 8 4 12 2 6 10 14 7 3 5 1 11 13 15 9 の順に挿入した時の木構造の移り代りの図を示す.
左上の A [8,8] 8 は, 左のAが図の記号, 右の8が今回挿入したn, 中央が(最初は空だった)木にnが挿入された, 新しい木 のように見る.

2021年10月31日日曜日

クイーン支配問題

前回のクイーン支配問題のブログは, 8×8の大きさのチェス盤に, 5個の クイーンを置き, 盤面のすべての場所が支配できる配置の数を計算するものであった. 4860という解が得られたが, どう考えても, 90度の回転や, 左右の反転で同じになる 仲間が沢山含まれている筈である.

そこで次の課題は, 本質的に異なる配置の数を知ることである. 前回のすべての配置の それぞれに, 回転や反転を施し, それで一緒になる配置を仲間とするプログラムを書き, 実行すると, 638という答が出た.

クイーン支配の配置は, かなりランダムらしいから, 回転で仲間が4倍になり, 反転で また2倍になって, 仲間は8人がほとんどになりそうである. しかし, 反転するとどうなるか. 左右反転では, 偶数の列に奇数のクイーンを置くから, クイーンの奇数個を中央の列に置き, 残りの偶数個を, 左右対称の位置に置けばいいわけだが, 中央の列が存在しないから, 左右や上下の反転で一緒になるものはない. 一方, 対角線で反転することを考えると, 対角線には中央があるから, そこに奇数個を置くと, 対称になり得る. この対角線の 中央を通る, もう一方の対角線で, もう一回反転できるかというと, 対角線の長さ が8という偶数なので, そういう反転はあり得ない.

従って, 4人の仲間のグループと, 8人の仲間のグループとで, グループは638あり, メンバーの総数は, 4860人ということになる.

4人と8人のグループはそれぞれ何個あるかは, いわゆる鶴亀算で得られる. a+b=638, 4a+8b=4860 だから a=61, b=577である.

そこで, その638の配置をすべて図示したのが次である. 1行に16個, 縦に20行あって, 最初のページに320個, 次のページに318並んでいる. 638個の配置には, 0番から637番 の番号があり, 最初のページの左上が0, そこから右へ1,2,3,...と続く. 次のページの最下行の右端が, 637番である.

これらの図は, 小さ過ぎるので, それぞれの図で, 各クイーンがどのように全盤を支配する かを示す図を描いた. 下の20の図は, 上の図の左端の配置の一つ置きに描いたもので, 塗りつぶしの丸がクイーン, それから四方八方に出る破線が支配領域である. 8×8 のすべての盤面が破線で蔽われているのが分る.

最後に4人仲間の配置はどれかという話題である. 合格発表のような下の番号が, 2回対称軸を持つ61の4人組である.

(20 25 29 30 59 128 129 137 141 149 163 165 171 180 182 183 204 215 219 229 232 238 242 258 259 261 267 276 283 288 289 365 368 372 391 393 441 468 470 473 474 511 516 517 518 524 526 530 551 568 569 574 582 584 592 594 605 617 620 623 633)

この表の0, 15, 30, 45番目を選び(20, 183, 289, 524番になる), 上と同じように 示したのが, 下の図である. 対角線で反転するのが分る. 中のふたつは面白い.

2021年10月27日水曜日

クイーン支配問題

クイーン支配問題 The Art of Computer Programmingに, クイーン支配問題(queens domination probmel)という話題があった.

8×8の盤に, 5個のクイーンを, 盤のすべての場所を支配するように置く という問題である. 1863年の本にあるそうだから, 古典的問題である.

この解は, 4860通りあるらしいので, 自分でも計算することにした.

私のことだから, Schemeを使う. またSchemeを使うなら, bit-stringが役立つ.

プログラムを書いたら, 思ったより簡単で, 1時間はかからなかった. こんな具合だ.

(define (r n) (map (lambda (i) (+ (* n 8) i)) (range 0 8)))
(define (c n) (map (lambda (i) (+ (* i 8) n)) (range 0 8)))
(define (d n) (if (< n 0)  ;n=i-j -7<=n<8
  (map (lambda (i) (+ (* i 9) (- 0 n))) (range 0 (+ n 8)))
  (map (lambda (i) (+ (* i 9) (* n 8))) (range 0 (- 8 n)))))
(define (u n) (if (< n 7)  ;n=i+j 0<=n<15
  (map (lambda (i) (+ (* i -7) (* n 8))) (range 0 (+ n 1)))
  (map (lambda (i) (+ (* i -7) (+ n 49))) (range 0 (- 15 n)))))
(define (bs n)
  (let ((i (quotient n 8)) (j (modulo n 8))
	(b (make-bit-string 64 #t)))
    (for-each (lambda (m) (bit-string-clear! b m))
      (append (r i) (c j) (d (- i j)) (u (+ i j))))
    b))
(define bss (map bs (range 0 64)))
(do ((i 0 (+ i 1)) (count 0) (bs (make-bit-string 64 #t)))
    ((= i 60) count) (display (list 'i i))
(let ((bs (bit-string-copy
       (bit-string-and bs (list-ref bss i)))))
  (do ((j (+ i 1) (+ j 1))) ((= j 61))
   (let ((bs (bit-string-copy
          (bit-string-and bs (list-ref bss j)))))
    (do ((k (+ j 1) (+ k 1))) ((= k 62))
     (let ((bs (bit-string-copy
            (bit-string-and bs (list-ref bss k)))))
      (do ((l (+ k 1) (+ l 1))) ((= l 63))
       (let ((bs (bit-string-copy
              (bit-string-and bs (list-ref bss l)))))
	(do ((m (+ l 1) (+ m 1))) ((= m 64))
         (if (bit-string-zero?
              (bit-string-and bs (list-ref bss m)))
          (set! count (+ count 1))))))))))))
盤のますには, 左上から右上にかけ, そして上から下へ, 0から63まで番号を付ける. 行rは上が0, 下が7(その番号をiとする), 列cは左が0, 右が7(その番号をjとする), 左上から右下への対角線dは, ます7を通る(i=0,j=7)のを-7, ます0を通るのを0, ます56を通る(i=7, j=0)のを7, 左下から右上への対角線uは, ます0を通る(i=0,j=0) のを0, ます56を通る(i=7,j=0)のを7, ます63を通る(i=7,j=7)のを14とする.

先頭の関数, r,c,d,uは, 行n, 列n, 下り対角線n, 上り対角線nにいるクイーンが支配する ますの番号を返す.

bssは, ますnのクイーンが支配するますをビット0, それ以外をビット1にした, 64ビット のビット列の配列である。<br>
これだけ用意し, 後は5個のクイーンを, 全面に置きながら, 支配場所を集め, 5個置いた 時に, すべてのます(に対応するビット)が0なら, countを1増やす.

やってみたら, 18秒くらいで, 4860が得られた.

2021年6月7日月曜日

Piの1000桁

前回のブログのPiの1000桁をもう少し詳しく説明したい. Piの値の計算に Machinの式

π/4 = 4 tan-1(1/5)- tan-1(1/239)
tan-1(1/5)=1/5-1/(3·53)+1/(5·55) -1/(7·57)+...
tan-1(1/239)=1/239-1/(3·2393)+1/(5·239 5)-1/(7·2397)+...

を使った. これは絵を描いてみるとこういうことだ.
円弧の中に底辺:高さが5:1の直角三角形を描く. その中心角をαとする. つまり α=tan-1(1/5). αは11.3°くらいだ.

αの斜辺の上にまた今の直角三角形を描く. それを3回繰り返すと4αの 角度が得られ, 上の値を使えば4α=45.2°だ. 一方π/4は45°だから, 今度は底辺:高さ=239:1の直角三角形(赤線)を角度を引く方法に描くとその斜辺が 丁度45°になるというのがMachinの式であった. 計算すると次のようになる.
一番上はtanの加法公式である. これを使いαから2α, 4α を 計算すると4α=120/119が得られる. この120と119の差が1なのが 重要だ. 4αとπとの差の角度βのtanはまた加法公式を使うと1/239 になるわけだ. 一番下のようにMachinの式になる.

Machinが分ったから計算に進む. 最初のPiの1000桁のブログでやったように, 40桁計算する様子を見るには, 小数点を40桁右に移動して長い整数(bignum)で計算する. つまり1/5とせず, ss=1040/5 から始める.

それを1で割ってtsに足し, ss=ss/(5×5)にしてから ss/3をtsから引き, ss=ss/(5×5)にしてから ss/5をtsに足し, ... を計算する. しかし最初だけ 5で割り, 後は5×5で割るなら, 最初の分子を5倍しておき, いつも5×5 で割るのがプログラムは簡単だ.

さらにMachinの式はπ/4を計算するが, 我々はπが欲しいから, 最初のss は 4×4×5×1040とする.

従って
(set! ts 0)
(set! ss (* 4 4 5 (expt 10 40)))
続けて
(set! ss (quotient ss 25))
(set! ts (+ ts (quotient ss 1)))
(set! ss (quotient ss 25))
(set! ts (- ts (quotient ss 3)))
...
これをssが0になるまで実行する. 初期化の後 doループを使い
(set! ts 0)
(set! ss (* 4 4 5 (expt 10 40)))
(do ((n 1 (+ n 2)) (ff #t (not ff))) ((= ss 0) n)
  (set! ss (quotient ss 25))
  (set! ts ((if ff + -) ts (quotient ss n))))
doループは, nを1, ffを#tで始め, ループを回る度にnはn+2にし, ffは ¬ffにしてssが0になるまで, 本体を実行すると読む. つまりnは1,3,5,7,... ffは#t,#f,#t,#f,...となる. 終了条件s=0の後のnはdoループの返す値である. nをどこまで進めるとsが0になったかを知るためだ.

((if ff + -) ts (quotient ss n))は, ffが#tなら (+ ts (quotient...)), #fなら(- ts (quotient ...))のことである. つまり1回置きに加算と減算を繰り返す.

初期化したssとループでのssとtsを出力すると
ss 800000000000000000000000000000000000000000(
(n 1)
ss 32000000000000000000000000000000000000000
ts 32000000000000000000000000000000000000000
(n 3)
ss  1280000000000000000000000000000000000000
ts 31573333333333333333333333333333333333334
(n 5)
ss    51200000000000000000000000000000000000
ts 31583573333333333333333333333333333333334
(n 7)
ss     2048000000000000000000000000000000000
ts 31583280761904761904761904761904761904763

...

(n 57)
ss                                        23
ts 31583289575980921339207962431166446951613
(n 59)
ss                                         0
ts 31583289575980921339207962431166446951613
61
最後の61はdoループが返したnの値である. 下の方の23や0は元々は右端に 揃っていたのだが, ブログはそういう風には出してくれない.

同じようにtan-11/239も計算すると
(n 17)
ss                                         1
ts 31415926535897932384626433832795028841973
(n 19)
ss                                         0
ts 31415926535897932384626433832795028841973
21
と円周率が見えてくる.

余談だが, 私の恩師 故高橋秀俊先生は中学生の時円周率を筆算で50桁計算されたそうだ. 40桁でもn=57まで計算するのだから, 50桁はまだまだ先までで大変だっただろう.

次はこの計算を有限語長ので分割しながら実行するわけだが, それは次のブログにしよう.