tag:blogger.com,1999:blog-88513448908512823732024-03-13T11:07:48.272+09:00パラメトロン計算機和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.comBlogger374125tag:blogger.com,1999:blog-8851344890851282373.post-45408177174124304952022-09-20T16:10:00.001+09:002022-09-25T14:27:47.884+09:00満月の十五夜<html>
<body>
今年は9月10日が中秋の名月, つまりお月見であった. 幸いよく晴れて, 満月が堪能出来た.<br><br>
ところで新聞が「今年の中秋の名月は満月」と報じた. 私はもちろん旧暦の15日が満月にならない方が普通で, たまには満月になることも知っていたから, この新聞記事をみて何とも思わなかったが, どの程度すれすれに満月なのかと天文年鑑を見ると満月は18時59分であった.
<br><br>
旧暦の15日が満月にならないのは, 新月の時点のある日を旧暦の1日とするからである. 仮に1朔望月を29.5とし, 月齢14.75を満月とすれば, ある日の朝0.25日までに新月があればその日が1日で, 15日の晩には月齢が14.75に達する. そうでないと, 月齢が14.75になるのは, 旧暦16日になる.<br><br>
ところが月の公転には遅速があり, 満月は太陽と月の黄経の差が180度の時のことだから, 月齢14.75からプラスマイナス1日くらいずれ得る.<br><br>
そこで新月から満月, 満月から新月の経過時間の変化を見てみたいと思った. 2013年から2022年の天文年鑑から, 新月と満月の日時を書き出し, パソコンに入力して新月と次の満月, 満月と次の新月の経過時間を計算した.<br><br>
10年間にある朔望月は, 19年7閏法を考えると, 120+3.5回だから, 日時データは247, 間隔は246あった. 最大値は15.601, 最小値は13.907, 平均値は14.767. <br><br>
最後の10個は<br>
... 14.32 15.278 14.071 15.497 13.958 <br>15.579 14.009 15.497 14.216
15.256<br> である. <br><br>
それを絵にしたのが下だ. 赤は新月から満月, 青は満月から新月までの経過時間である. 縦軸の単位は日.<br><br>
<div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj0D91Nq3cDju_5OOTd75IO9JuiEr5GrdqQcTWj_D33ucsmTkzbMiMEsFgzZKJSW7_MkoBJmF3Cx0jG0OaGINEmHWVrCOqAMKE_kZK7LdUZVetemNYhr1P-dAziRtoiT5HDvY01MCPFML9i3G2KwPwdrQqm-9zdSiJpp3Ug04VLWn33OU03ZaAa8C-log/s1198/new-moon-fig0.png" style="display: block; padding: 1em 0; text-align: center; "><img alt="" border="0" width="400" data-original-height="230" data-original-width="1198" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj0D91Nq3cDju_5OOTd75IO9JuiEr5GrdqQcTWj_D33ucsmTkzbMiMEsFgzZKJSW7_MkoBJmF3Cx0jG0OaGINEmHWVrCOqAMKE_kZK7LdUZVetemNYhr1P-dAziRtoiT5HDvY01MCPFML9i3G2KwPwdrQqm-9zdSiJpp3Ug04VLWn33OU03ZaAa8C-log/s320/new-moon-fig0.png"/></a></div>
図を見ると, 時々赤と青が同じになり, つまり新月から満月までと満月から新月までが同時間になり, その中間は一方は増えて減り, 他方が減って増える.<br><br>
その理由は多分こうであろう. 月の公転軌道も楕円であり, 長軸の一方が近地点, 他方が遠地点である. 近地点では公転速度が速く, 遠地点では遲いのは常識だ. <br><br>
近地点と遠地点の近くに新月と満月の場所があれば, 新月から満月も満月から新月もほぼ同時間になる. しかし, 新月と満月がその中間くらいにあると, 近地点側を通る方は経過時間が短かく, 遠地点側は長い. それがこの図の謎解きであろう.<br><br>
2022年の新月と満月の日の月の地心距離を天文年鑑から調べ, 新月から満月までやその逆の経過時間がどの地心距離からどの地心距離であったかの図を描いてみた. それを下に示す. 横軸は経過時間(単位は日)で, 縦軸は地心距離(単位は万km). <br><br>
<div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiO7LagwhipA7nWcQJ31ewdc4mdZNm4i5DdMffFd7DXO6PsEkYVOIjcgUuo0-2QBdCd7RIFEVMOBt2dMK6hzCmPUj3QPUxOe5ID57dncQ_VRGrsnYRKG_bGR3gVh-CJeMkcp_wLuZtY6MAMQ3WZBrn_aCAfxmF-H4v0wUZ5YqdFQGDO4Q6vk6Cfe5VmcQ/s846/new-moon-fig1.png" style="display: block; padding: 1em 0; text-align: center; "><img alt="" border="0" width="400" data-original-height="344" data-original-width="846" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiO7LagwhipA7nWcQJ31ewdc4mdZNm4i5DdMffFd7DXO6PsEkYVOIjcgUuo0-2QBdCd7RIFEVMOBt2dMK6hzCmPUj3QPUxOe5ID57dncQ_VRGrsnYRKG_bGR3gVh-CJeMkcp_wLuZtY6MAMQ3WZBrn_aCAfxmF-H4v0wUZ5YqdFQGDO4Q6vk6Cfe5VmcQ/s320/new-moon-fig1.png"/></a></div>
これを見ると, 上の推論が正しいことが分る. 左上から右下へ来る矢印は, 遠地点から近地点までの経過で14.8日くらい掛ることで, 15.7日くらいの水平の矢印は, 中間から中間への遠地点側を通る経過, 14.0日くらいの水平の矢印は, 近地点側を通る経過である.<br><br>
十五夜が満月になるというのには, こういう仕掛けがあったわけだ.
その後, 国立天文台暦計算室の<a href="https://eco.mtk.nao.ac.jp/koyomi/wiki/B7EEA4CECBFEA4C1B7E7A4B12FBCFEB4FCA4CECAD1C6B0.html">ページ</a>を見付けた.
</body</html>
和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-181588227720588212022-09-08T15:43:00.001+09:002022-09-08T15:46:11.590+09:00切頭八面体久し振りに切頭八面体について書きたい. TAOCPにおやという記述があった.<br><br>
a: 切頭八面体の基本的な置き方は, (htmlでは式が書き難いから言葉で述べると)
0,1,2,3の, 同じもののない3つ組である.<br><br>
b: 体心立方格子のVolonoi領域である.<br><br>
bはなるほどそうらしいと思うが, そしてそういう言い方もあるのかと思うが,
aはすぐには分らなかった. その組を作ってみると,
<pre>
((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))
</pre>
になり, たしかに24個あって, 切頭八面体の頂点の数と同じである. <br><br>
それぞれの3つ組を(<it>x</it>, <it>y</it>, <it>z</it>)と思い, 高さ(<it>z</it>)
別に絵を描くと, 下の図の上段のようになる. z=0, つまり底の面に点が6個あるから,
これは六角形を床に置いた図らしいと分る. 六角形にしてはいびつだが, それは直交
座標のためらしいから, x軸とy軸の角度を60度にしてもう一度描いたのが中段の図である.<br><br>
なるほど六角形が現れる.<br><br>
<div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiMRX1vLZaHDZuCHgepJzG0wlqldTFRb_T2WBVki_0F2J-wvZx3t_uAbfAlCLGTR-eBVJTjYulzBuT8b6Tf-ueOpSr8qn5RHFlAEkkfGCFfi2naZAffr6-ZBdXbwm2H9Uy3VcriCehE96CugNOwHVLaGG7uL73mbL1SVekda1O12VAQHvtFzQKUPb2Kjw/s1340/truncoct0.png" style="display: block; padding: 1em 0; text-align: center; "><img alt="" border="0" height="450" data-original-height="1340" data-original-width="1232" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEiMRX1vLZaHDZuCHgepJzG0wlqldTFRb_T2WBVki_0F2J-wvZx3t_uAbfAlCLGTR-eBVJTjYulzBuT8b6Tf-ueOpSr8qn5RHFlAEkkfGCFfi2naZAffr6-ZBdXbwm2H9Uy3VcriCehE96CugNOwHVLaGG7uL73mbL1SVekda1O12VAQHvtFzQKUPb2Kjw/s320/truncoct0.png"/></a></div>
そこでx軸とz軸, y軸とz軸の角度も60度にしたグラフ用紙を作って, 上に点を
書き込むことにした. つまり, x軸y軸z軸は正四面体の3本の辺に沿っている.
それが上の赤緑青の線である. 赤はz=0の平面座標で,
その三角形の中心を起点とした緑のz=1の平面座標が出来, その三角形の中心を起点とした青の
z=2の平面座標が出来る. 同様にしてz=3の座標が出来るが, それはz=0のと一致する.<br><br>
黒線で描いたのが六角形を底にした切頭八面体である. <br><br>
<!--ps/polyhedron/truncoct.ps -->
これを確認するため, 別の切頭八面体の絵もある. 左上Aは, 私が通常に描く切頭八面体
である. その中心を通る直交するx軸(赤), y軸(緑), z軸(青)も描いてある.<br><br>
<div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi6vHrdP8y5xXXVz40pYIfBYClvOvCTIibZK9toT9EnpDFsimOyFFRo7BJs-8GvJE9BCUWrawfnbxUe75qyaVgmveRr82A_gtFQZ83Jnl9fugGNaS5aZGZL7GNHuk4CBMLgdLWSAGD2ICF6gyZJslwzCj_JEN_Y_85LL3zaHP5vbc4pLPLjs8sw7fjEVg/s1266/truncoct1.png" style="display: block; padding: 1em 0; text-align: center; "><img alt="" border="0" width="400" data-original-height="920" data-original-width="1266" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEi6vHrdP8y5xXXVz40pYIfBYClvOvCTIibZK9toT9EnpDFsimOyFFRo7BJs-8GvJE9BCUWrawfnbxUe75qyaVgmveRr82A_gtFQZ83Jnl9fugGNaS5aZGZL7GNHuk4CBMLgdLWSAGD2ICF6gyZJslwzCj_JEN_Y_85LL3zaHP5vbc4pLPLjs8sw7fjEVg/s320/truncoct1.png"/></a></div>
私の描画プログラムは, 頂点の座標を座標軸に平行移動したり, 座標軸に沿って回転
したり出来る. BはAの図を, z軸について45度回転したものだ. Cはそれを右下の
5,12の辺が, xyz軸の中心に来るように移動したもの. Dはそれを右下の六角形が
底に来るように, x軸について回転したものである. 最後のEは, 座標軸の中心に
切頭八面体の中心が来るように移動した.<br><br>
この時, 各頂点の直交座標での位置は次の通りである. (左は頂点の番号順, 右は高さ別)
<pre>
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)
</pre>
各頂点の直交座標(x,y,z)を, 斜方座標(a,b,c)に変換するには, 例えば
<pre>
(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)))
</pre>
しかしこれで変換すると0.707や2.121が出るので, 正規化すると, 先程の24個の
直交座標は
<pre>
((-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))
</pre>
なるほど, 最初に書いた座標は0,1,2,3であったが, 今回のは. 切頭八面体の中心に座標の
中心を置いたので, -3,-1,1,3で構成されていた.
</body>
</html>
和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-72410993122340486912022-09-06T17:21:00.002+09:002022-09-06T20:02:24.899+09:00Wolframの対数表学生の頃から気になっていたことがある. 高木貞治先生の「解析概論」のp.215の脚注に
「Wolframノ表ニハ1000以下ノ素数ノ自然対数ノ50桁ノ表ガ掲ゲラレテヰル.
コノ表ハ既ニ少年がうすガ愛用シタモノデアル.」と書いてあったことだ.
また, 岩波文庫の「近世数学史談」のp.65に
「Wolframの計算した10000以下の素数及び若干の特殊の数の自然対数の表が
附いていたのである.」と書いてあるのも見付けた. 1000までの素数か10000までの
素数かに違いはあるが, とにかくそういう対数表が存在していたらしい.<br><br>
学生の頃にそんな数表を探すのは夢のまた夢であったが, 長生きはするもので, この歳に
なってインターネットを検索していたら, 関係のあるウェブページがあった.<br><br>
そのひとつは<a href="https://hal.inria.fr/hal-01615064/document">
A reconstruction of the Mathematical Tables Project’s table of natural logarithms</a>という文献で, その図1は, Wolframが計算したものを, Schulzeが1778年(Gaussの生まれた翌年)に
出版したものの最後のページの写真. これを見て, 1000までの方が違っているのが分かった.
しかも最後は10009で精度は48桁であった. また6個の素数には対数が記入
されていず, Wolframが病気で計算できなかったらしい.
その図2はVegaが欠けている値を計算して
追加して1794年に出版したものの最後のページの写真がある. いずれにしても,
Gaussが所持していたWolframの数表はこういうものであったらしい.<br><br>
もうひとつは, <a href="https://www.ams.org/journals/mcom/1950-04-032/S0025-5718-1950-0037814-2/S0025-5718-1950-0037814-2.pdf">
New Information Concerning Isaac Wolfram's Life and Calculations</a>である.
これによると, Isaac Wolframはオランダの砲兵の士官であった.<br><br>
Wolframは何年もかけてこの数表を作ったらしいが, 現代は便利な世の中で, 自然対数も
簡単に得ることが出来る<br><br>
そのひとつは, 別のWolframだが, <a href="https://www.wolframalpha.com/">
WolframAlpha</a>である. このページを開くと入力用窓があり,
そこに「ln(2)」といれると, 2の自然対数が64桁くらい表示される.
0.6931471805599453094172321214581765680755001343602552541206800094
私もこの数表を作ってみたくなった. しかし. WolframAlphaは手動でしか使えないから,
別の方法が必要だ. あれこれ探していたら, Python3でも値が得られることが分った.
<pre>
from decimal import *
getcontext().prec = 55
Decimal(2).ln() =>
Decimal('0.6931471805599453094172321214581765680755001343602552541')
</pre>
これを使うと, Wolframの対数表みたいなものはすぐ作れる. 10000までの素数は
1229個あるが, 最初の20個の素数と, 10000未満の最後の19個の素数に対して
作た対数表はこんなようだ.
<div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgaNpcjdm7TbRAq3zGMHHF_HvJWqaGNTZpZ6LkkQP0PsBfyAAsgCc2BUDuwJp-3SBVkwGN-aQTon_cSWPrYuUsCAzQPGNyPyC7Otk-3ce_vi0qod09lKYTlnq9PyyKaIbag7KTCRLkTlkkR_o4lbXfiHvgA7AcI20md4Ymhu2GD5nYquaG3yAAkmfJc-g/s1146/natlogtab0.png" style="display: block; padding: 1em 0; text-align: center; "><img alt="" border="0" height="500" data-original-height="1146" data-original-width="884" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEgaNpcjdm7TbRAq3zGMHHF_HvJWqaGNTZpZ6LkkQP0PsBfyAAsgCc2BUDuwJp-3SBVkwGN-aQTon_cSWPrYuUsCAzQPGNyPyC7Otk-3ce_vi0qod09lKYTlnq9PyyKaIbag7KTCRLkTlkkR_o4lbXfiHvgA7AcI20md4Ymhu2GD5nYquaG3yAAkmfJc-g/s320/natlogtab0.png"/></a></div>
2から9973までの素数に対する自然対数の表は,
<a href="https://www.iijlab.net/~ew/natlogtab.pdf">私の計算機</a>にある.
和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-52844480293855959122022-04-12T15:56:00.000+09:002022-04-12T15:56:06.216+09:00mex関数<html>
<body>
mex関数についてブログに全回書いたのは, もう2年も前だった. その頃, Schemeでの実装も
出来ていたが, なんとなく気に入らなかった.<br>
最近, もう少し気持ちよくしたいと思い, やった書き直したので, ブログに記録して
おく.
まずinsertion. <tt>(ins n t)</tt>は, nを, 木tに挿入する関数である.
<pre>
(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)))))))
</pre>
一方, これを使う<tt>(mex ms)</tt>は次のようだ.
<pre>
(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))))))
</pre>
全回と似たような図だが,
<tt> 8 4 12 2 6 10 14 7 3 5 1 11 13 15 9 </tt>
の順に挿入した時の木構造の移り代りの図を示す.
<div class="separator" style="clear: both;"><a href="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj-rCsXqvj49cr9373wYYvV4dNlsWCz2j9n-hexxrZvfgpkQ4HtPZBX27VDt0v-wG7gsaKdQgHawzyh8ZvGnZ5Jh9o5suam8YS_FmpHZaJ7BKZ1UU-lDkcq0nMXBOK4nHpwXtt4Dp4FihV8mhnNtPtHBHjeVVGs3QQQBN954WG0-uCIgVSl-jpnmnW60w/s1560/mexfig2.png" style="display: block; padding: 1em 0; text-align: center; "><img alt="" border="0" width="400" data-original-height="846" data-original-width="1560" src="https://blogger.googleusercontent.com/img/b/R29vZ2xl/AVvXsEj-rCsXqvj49cr9373wYYvV4dNlsWCz2j9n-hexxrZvfgpkQ4HtPZBX27VDt0v-wG7gsaKdQgHawzyh8ZvGnZ5Jh9o5suam8YS_FmpHZaJ7BKZ1UU-lDkcq0nMXBOK4nHpwXtt4Dp4FihV8mhnNtPtHBHjeVVGs3QQQBN954WG0-uCIgVSl-jpnmnW60w/s320/mexfig2.png"/></a></div>
左上の A [8,8] 8 は, 左のAが図の記号, 右の8が今回挿入したn,
中央が(最初は空だった)木にnが挿入された, 新しい木 のように見る.
</body>
</html>
和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-73870967714326872072021-10-31T17:19:00.003+09:002022-05-29T20:57:56.818+09:00クイーン支配問題<html>
<body>
前回のクイーン支配問題のブログは, 8×8の大きさのチェス盤に, 5個の
クイーンを置き, 盤面のすべての場所が支配できる配置の数を計算するものであった.
4860という解が得られたが, どう考えても, 90度の回転や, 左右の反転で同じになる
仲間が沢山含まれている筈である.<br><br>
そこで次の課題は, 本質的に異なる配置の数を知ることである. 前回のすべての配置の
それぞれに, 回転や反転を施し, それで一緒になる配置を仲間とするプログラムを書き,
実行すると, 638という答が出た.<br><br>
クイーン支配の配置は, かなりランダムらしいから, 回転で仲間が4倍になり, 反転で
また2倍になって, 仲間は8人がほとんどになりそうである. しかし, 反転するとどうなるか.
左右反転では, 偶数の列に奇数のクイーンを置くから, クイーンの奇数個を中央の列に置き,
残りの偶数個を, 左右対称の位置に置けばいいわけだが, 中央の列が存在しないから,
左右や上下の反転で一緒になるものはない. 一方, 対角線で反転することを考えると,
対角線には中央があるから, そこに奇数個を置くと, 対称になり得る. この対角線の
中央を通る, もう一方の対角線で, もう一回反転できるかというと, 対角線の長さ
が8という偶数なので, そういう反転はあり得ない.<br><br>
従って, 4人の仲間のグループと, 8人の仲間のグループとで, グループは638あり,
メンバーの総数は, 4860人ということになる.<br><br>
4人と8人のグループはそれぞれ何個あるかは, いわゆる鶴亀算で得られる.
a+b=638, 4a+8b=4860 だから a=61, b=577である.<br><br>
そこで, その638の配置をすべて図示したのが次である. 1行に16個, 縦に20行あって,
最初のページに320個, 次のページに318並んでいる. 638個の配置には, 0番から637番
の番号があり, 最初のページの左上が0, そこから右へ1,2,3,...と続く.
次のページの最下行の右端が, 637番である.<br><br>
<div class="separator" style="clear: both;"><a href="https://1.bp.blogspot.com/-5sdfARSlbVA/YX4BZMuP0zI/AAAAAAAAAgQ/BMluXIUXPQUrjjYnCCHEiXAXyteDnR-JwCLcBGAsYHQ/s792/sols0.png" style="display: block; padding: 1em 0; text-align: center; "><img alt="" border="0" height="530" data-original-height="792" data-original-width="500" src="https://1.bp.blogspot.com/-5sdfARSlbVA/YX4BZMuP0zI/AAAAAAAAAgQ/BMluXIUXPQUrjjYnCCHEiXAXyteDnR-JwCLcBGAsYHQ/s600/sols0.png"/></a></div>
<div class="separator" style="clear: both;"><a href="https://1.bp.blogspot.com/-FabiVykHUW8/YX4C9KV3psI/AAAAAAAAAgY/bqT3y3xFKYkMU3bUzgCIqo-kqqQIiSbIQCLcBGAsYHQ/s792/sols1.png" style="display: block; padding: 1em 0; text-align: center; "><img alt="" border="0" height="530" data-original-height="792" data-original-width="612" src="https://1.bp.blogspot.com/-FabiVykHUW8/YX4C9KV3psI/AAAAAAAAAgY/bqT3y3xFKYkMU3bUzgCIqo-kqqQIiSbIQCLcBGAsYHQ/s600/sols1.png"/></a></div>
これらの図は, 小さ過ぎるので, それぞれの図で, 各クイーンがどのように全盤を支配する
かを示す図を描いた. 下の20の図は, 上の図の左端の配置の一つ置きに描いたもので,
塗りつぶしの丸がクイーン, それから四方八方に出る破線が支配領域である. 8×8
のすべての盤面が破線で蔽われているのが分る. <br><br>
<div class="separator" style="clear: both;"><a href="https://1.bp.blogspot.com/-ux7mde12l9M/YX4LG4yMVZI/AAAAAAAAAgo/OdJ8lO3oS4wg7bp9WjSU2h0AILBiPGOswCLcBGAsYHQ/s792/queendom.png" style="display: block; padding: 1em 0; text-align: center; "><img alt="" border="0" height="580" data-original-height="792" data-original-width="612" src="https://1.bp.blogspot.com/-ux7mde12l9M/YX4LG4yMVZI/AAAAAAAAAgo/OdJ8lO3oS4wg7bp9WjSU2h0AILBiPGOswCLcBGAsYHQ/s600/queendom.png"/></a></div>
最後に4人仲間の配置はどれかという話題である. 合格発表のような下の番号が,
2回対称軸を持つ61の4人組である. <br><br>
(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)<br><br>
この表の0, 15, 30, 45番目を選び(20, 183, 289, 524番になる), 上と同じように
示したのが, 下の図である. 対角線で反転するのが分る. 中のふたつは面白い.
<div class="separator" style="clear: both;"><a href="https://1.bp.blogspot.com/-ds5FpMPcOnw/YX5PNkHPcPI/AAAAAAAAAgw/nfT3CJt666UYj5Vdkif9Hr-ly0_eq-jZwCLcBGAsYHQ/s1982/group4.png" style="display: block; padding: 1em 0; text-align: center; "><img alt="" border="0" width="420" data-original-height="494" data-original-width="1982" src="https://1.bp.blogspot.com/-ds5FpMPcOnw/YX5PNkHPcPI/AAAAAAAAAgw/nfT3CJt666UYj5Vdkif9Hr-ly0_eq-jZwCLcBGAsYHQ/s400/group4.png"/></a></div>
</body>
</html>
和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-44804592582668706712021-10-27T11:16:00.001+09:002022-05-29T20:55:59.035+09:00クイーン支配問題<html>
<body>
クイーン支配問題
The Art of Computer Programmingに, クイーン支配問題(queens domination
probmel)という話題があった.<br><br>
8×8の盤に, 5個のクイーンを, 盤のすべての場所を支配するように置く
という問題である. 1863年の本にあるそうだから, 古典的問題である.<br><br>
この解は, 4860通りあるらしいので, 自分でも計算することにした.<br><br>
私のことだから, Schemeを使う. またSchemeを使うなら, bit-stringが役立つ.
<br><br>
プログラムを書いたら, 思ったより簡単で, 1時間はかからなかった. こんな具合だ.
<br><br>
<pre>
(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))))))))))))
</pre>
盤のますには, 左上から右上にかけ, そして上から下へ, 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とする.
<br><br>
先頭の関数, r,c,d,uは, 行n, 列n, 下り対角線n, 上り対角線nにいるクイーンが支配する
ますの番号を返す. <br><br>
bssは, ますnのクイーンが支配するますをビット0, それ以外をビット1にした, 64ビット
のビット列の配列である。<br><br>
これだけ用意し, 後は5個のクイーンを, 全面に置きながら, 支配場所を集め, 5個置いた
時に, すべてのます(に対応するビット)が0なら, countを1増やす.<br><br>
やってみたら, 18秒くらいで, 4860が得られた.
</body>
</html>
和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-54171012772614462512021-06-07T11:44:00.002+09:002021-06-07T20:32:59.250+09:00Piの1000桁前回のブログのPiの1000桁をもう少し詳しく説明したい. Piの値の計算に Machinの式<br><br>
π/4 = 4 tan<sup>-1</sup>(1/5)- tan<sup>-1</sup>(1/239)<br>
tan<sup>-1</sup>(1/5)=1/5-1/(3·5<sup>3</sup>)+1/(5·5<sup>5</sup>)
-1/(7·5<sup>7</sup>)+...<br>
tan<sup>-1</sup>(1/239)=1/239-1/(3·239<sup>3</sup>)+1/(5·239<sup>
5</sup
>)-1/(7·239<sup>7</sup>)+...<br /><br />
を使った. これは絵を描いてみるとこういうことだ.
<div class="separator" style="clear: both;"><a href="https://1.bp.blogspot.com/-SWn0fqlyEP8/YLl5GzagFJI/AAAAAAAAAeE/-B6LOMZhxMgbqevaqMxIdk2bW2eywETmACLcBGAsYHQ/s842/pifig.png" style="display: block; padding: 1em 0; text-align: center; "><img alt="" border="0" width="320" data-original-height="840" data-original-width="842" src="https://1.bp.blogspot.com/-SWn0fqlyEP8/YLl5GzagFJI/AAAAAAAAAeE/-B6LOMZhxMgbqevaqMxIdk2bW2eywETmACLcBGAsYHQ/s320/pifig.png"/></a></div>
円弧の中に底辺:高さが5:1の直角三角形を描く. その中心角をαとする. つまり
α=tan<sup>-1</sup>(1/5). αは11.3°くらいだ.<br><br>
αの斜辺の上にまた今の直角三角形を描く. それを3回繰り返すと4αの
角度が得られ, 上の値を使えば4α=45.2°だ. 一方π/4は45°だから,
今度は底辺:高さ=239:1の直角三角形(赤線)を角度を引く方法に描くとその斜辺が
丁度45°になるというのがMachinの式であった. 計算すると次のようになる.
<div class="separator" style="clear: both;"><a href="https://1.bp.blogspot.com/-Pxpsj-HDlcI/YLl55Ad3FaI/AAAAAAAAAeM/1kFh5LmDT-s7aWnqtZG1lAbBg1m2Pe5IwCLcBGAsYHQ/s720/pieq.png" style="display: block; padding: 1em 0; text-align: center; "><img alt="" border="0" height="320" data-original-height="720" data-original-width="640" src="https://1.bp.blogspot.com/-Pxpsj-HDlcI/YLl55Ad3FaI/AAAAAAAAAeM/1kFh5LmDT-s7aWnqtZG1lAbBg1m2Pe5IwCLcBGAsYHQ/s320/pieq.png"/></a></div>
一番上はtanの加法公式である. これを使いαから2α, 4α を
計算すると4α=120/119が得られる. この120と119の差が1なのが 重要だ.
4αとπとの差の角度βのtanはまた加法公式を使うと1/239 になるわけだ.
一番下のようにMachinの式になる.<br><br>
Machinが分ったから計算に進む. 最初のPiの1000桁のブログでやったように,
40桁計算する様子を見るには, 小数点を40桁右に移動して長い整数(bignum)で計算する.
つまり1/5とせず, ss=10<sup>40</sup>/5 から始める.<br><br>
それを1で割ってtsに足し, ss=ss/(5×5)にしてから ss/3をtsから引き,
ss=ss/(5×5)にしてから ss/5をtsに足し, ... を計算する. しかし最初だけ
5で割り, 後は5×5で割るなら, 最初の分子を5倍しておき, いつも5×5
で割るのがプログラムは簡単だ.<br><br>
さらにMachinの式はπ/4を計算するが, 我々はπが欲しいから, 最初のss
は 4×4×5×10<sup>40</sup>とする.<br><br>
従って
<pre>
(set! ts 0)
(set! ss (* 4 4 5 (expt 10 40)))
</pre>
続けて
<pre>
(set! ss (quotient ss 25))
(set! ts (+ ts (quotient ss 1)))
(set! ss (quotient ss 25))
(set! ts (- ts (quotient ss 3)))
...
</pre>
これをssが0になるまで実行する. 初期化の後 doループを使い
<pre>
(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))))
</pre>
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になったかを知るためだ.<br><br>
((if ff + -) ts (quotient ss n))は, ffが#tなら
(+ ts (quotient...)), #fなら(- ts (quotient ...))のことである.
つまり1回置きに加算と減算を繰り返す.<br><br>
初期化したssとループでのssとtsを出力すると
<pre>
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
</pre>
最後の61はdoループが返したnの値である. 下の方の23や0は元々は右端に
揃っていたのだが, ブログはそういう風には出してくれない. <br><br>
同じようにtan<sup>-1</sup>1/239も計算すると
<pre>
(n 17)
ss 1
ts 31415926535897932384626433832795028841973
(n 19)
ss 0
ts 31415926535897932384626433832795028841973
21
</pre>
と円周率が見えてくる. <br><br>
余談だが, 私の恩師 故高橋秀俊先生は中学生の時円周率を筆算で50桁計算されたそうだ.
40桁でもn=57まで計算するのだから, 50桁はまだまだ先までで大変だっただろう.<br><br>
次はこの計算を有限語長ので分割しながら実行するわけだが, それは次のブログにしよう.
和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-317905278627000312021-05-25T15:57:00.003+09:002021-06-14T13:26:17.575+09:00Piの1000桁1970年頃のことだ. HITAC 10, FACOM R, MACC 7など, 1語が16ビットの4K
語のミニコンが多く使われていた. これらのミニコンは, 正面に押しボタ
ンのようなスイッチが並んでいて, 電源を入れた後, スイッチを操作し,
数語の「初期入力プログラム」を手動で入力する. プログラム全体が入力
し終わると初期入力プログラムが動き出し, 紙テープのプログラムを読み
込めた.<br><br>
ところがこのメーカー製プログラムの手の入力が面倒で, 誰もがもっと短
いプログラムに出来ないか工夫を凝らし, 徐々に短いプログラムが出現し
た. これが一番短いと思っていると, 他から更に短いのが出来たという話
が来る. するとまた知恵を絞り, その語数のプログラムが作れるのであっ
た.<br><br>
つまり「何語で出来る」という情報が重要で, それを聞くとその語数のプ
ログラムに向って注力できる.<br><br>
最近ツイッターを眺めていたら, Automatic Computer@nt_a_c さんの<br><br>
Automatic Computer@nt_a_c·3月18日<br>
PC-1で円周率、400桁は出せた。もう数日早くできればちょうど良かったのだが。後はどのくらいコードを削れるか。<br><br>
Automatic Computer@nt_a_c·3月24日<br>
円周率、できた。1000桁、印刷と合わせて3回読み込み、合わせて40分くらい。https://automaticcomputer.github.io/PC1onUnity/pi....<br><br>
というツイートがあった. 60年も前の計算機 PC-1でプログラムを書く奇
特な人がまだいるかと驚く一方, どうすればあの容量のメモリーでπが1000桁計算
できるか私も考え始めた.<br><br>
πの1000桁については, このブログでも2012年11月5日に書いている. 簡単に計算するには
やはりMachinの式でやることになろう.<br><br>
π/4 = 4 tan<sup>-1</sup>(1/5)- tan<sup>-1</sup>(1/239)<br><br>
tan<sup>-1</sup>(1/5)=1/5-1/(3·5<sup>3</sup>)+1/(5·5<sup>5</sup>)
-1/(7·5<sup>7</sup>)+...<br><br>
tan<sup>-1</sup>(1/239)=1/239-1/(3·239<sup>3</sup>)+1/(5·239<sup>
5</sup>)-1/(7·239<sup>7</sup>)+...<br><br>
PC-1は1語が18ビット. 計算は2語をつなげた長語36ビットを使う. 十進10
桁で1長語を使うから, 1000桁計算には100長語が必要で, eの1000桁では
300番地から500番地の100長語を作業場所に使っている. 512短語のPC-1で
は100長語はどうしても2組しかとれないから,それで上の計算をしなけれ
ばならない.<br><br>
普通に考えると16×tan<sup>-1</sup>(1/5)を計算するには, atanを置く場所tsを0にし,
多倍長の16/5を作る. これをssとする. それを1で割り(その商をqsとする), qsをtsに足す.
ssを25で割りssにし, そのssを3でわりqsとし, tsからqsを引く. ...<br><br>
これをやるとtan<sup>-1</sup>(1/5)を計算するだけで100長語が3組必要になりすでに
破綻する.<br><br>
tsはtan<sup>-1</sup>(1/5)とtan<sup>-1</sup>(1/239)で共用出来る. ssは
1/5が済んだ後1/239用に使えばよい.<br><br>
残る手段はqsを作ってから, 下の桁からtsに足したりtsから引いたりするのを
なんとかすることだ. <br><br>
我々が小学校以来習った加減算は下の桁から始めるものであった. それに対して
上級になって習った算盤では上の桁から始めた. 小学生のころ, 上の桁から加減算
が出来るのが不思議であったが, そのうち繰上げをローカルに処理しながら下へ進んで
いるからだと納得した. さらに一度繰上げが通過した桁に, その後から再び繰上げ
が来ることはないことも証明できた. 基数が10から2^35と大きくなると,
繰上げが何桁も伝搬する確率は小さくなるが, 無視できない. い<br><br>
要するにqsの結果はいらない. 各々の桁の商はqsに足すか引くかし, 剰余は次の
桁の除算に使うと作業場所は要らなくなり, なんとかなりそうだ.<br><br>
そこでSchemeを使ってシミュレートしたのが次のプログラムである.<br><br>
lnは長語数. b35, d10は二進35桁, 十進10桁の定数. tsはatanの値を蓄積する
場所, ssはatan(1/z)の計算で1/z, 1/z^3, 1/z^5,...を保持する場所である.
(zは5と239)
ldivはssをd = z^2で次々割る. ldiv+- はssをn(1,3,5,..)で割りながら,
部分商をtsに足すか引くかする.
<!-- scheme/pitest4.scm -->
<pre>
(define ln 100)
(define b35 (expt 2 35))
(define d10 (expt 10 10))
(define (ldiv ss d)
(let ((r 0) (n 0))
(do ((i 0 (+ i 1))) ((= i ln))
(set! n (+ (fix:lsh r 35) (list-ref ss i)))
(list-set! ss i (quotient n d))
(set! r (modulo n d)))))
(define (ldiv+- ts ss ff d)
(let ((r 0) (n 0) (s 0) (c 0))
(do ((i 0 (+ i 1))) ((= i ln))
(set! n (+ (fix:lsh r 35) (list-ref ss i)))
(set! r (modulo n d))
(set! s ((if ff + -) (list-ref ts i)
(quotient n d)))
(list-set! ts i (fix:and s (- b35 1)))
(set! c (fix:and s b35))
(do ((j (- i 1) (- j 1))) ((= c 0))
(set! s ((if ff + -) (list-ref ts j) 1))
(list-set! ts j (fix:and s (- b35 1)))
(set! c (fix:and s b35))))))
(define (lmul ns m)
(let ((e 0) (p 0))
(do ((i (- ln 1) (- i 1))) ((= i 0) e)
(set! p (+ (* (list-ref ns i) m) e))
(list-set! ns i (modulo p b35))
(set! e (quotient p b35)))))
(define ts (make-list ln 0))
(define ss (make-list ln 0)) (list-set! ss 0 80)
(do ((n 1 (+ n 2)) (ff #t (not ff))) ((> n 1433))
(ldiv ss (* 5 5))
(ldiv+- ts ss ff n))
(define ss (make-list ln 0))(list-set! ss 0 (* 4 239))
(do ((n 1 (+ n 2)) (ff #f (not ff))) ((> n 423))
(ldiv ss (* 239 239))
(ldiv+- ts ss ff n))
(do ((n 0 (+ n 1))) ((= n 100))
(display (list (lmul ts d10))))
</pre>
この出力の始めの方は次のようだ. <br><br>
<pre>
1415926535 8979323846 2643383279 5028841971
6939937510 5820974944 5923078164 628620899
8628034825 3421170679 8214808651 3282306647
938446095 5058223172 5359408128 4811174502
8410270193 8521105559 6446229489 5493038196
...
</pre>
10^10進法の各桁の後に空白があるのは, 十進の一番上が0なのが分るためだ. <br><br>
これで一応は計算出来るが, 考えてみるに, ssは徐々に小くなっていく.
その途中にずーっと先頭から除算を始めるのはやはり無駄である. つまり
除算を開始する位置を制御すれば計算は早く終わる. 下はそのように修正
したプログラムで, 変数bが除算開始位置を覚えている.<br><br>
<!-- scheme/pitest5.scm -->
<pre>
(define b 0)
(define (ldiv ss d)
(let ((r 0) (n 0) (q 0))
(do ((i b (+ i 1))) ((= i ln))
(set! n (+ (fix:lsh r 35) (list-ref ss i)))
(set! q (quotient n d))
(list-set! ss i q)
(if (= q 0) (set! b (+ i 1)))
(set! r (modulo n d)))))
(define (ldiv+- ts ss ff d)
(let ((r 0) (n 0) (s 0) (c 0))
(do ((i b (+ i 1))) ((= i ln))
(set! n (+ (fix:lsh r 35) (list-ref ss i)))
(set! r (modulo n d))
(set! s ((if ff + -)
(list-ref ts i) (quotient n d)))
(list-set! ts i (fix:and s (- b35 1)))
(set! c (fix:and s b35))
(do ((j (- i 1) (- j 1))) ((= c 0))
(set! s ((if ff + -) (list-ref ts j) 1))
(list-set! ts j (fix:and s (- b35 1)))
(set! c (fix:and s b35))))))
(define (lmul ns m)
(let ((e 0) (p 0))
(do ((i (- ln 1) (- i 1))) ((= i 0) e)
(set! p (+ (* (list-ref ns i) m) e))
(list-set! ns i (modulo p b35))
(set! e (quotient p b35)))))
(define start (real-time-clock))
(define ts (make-list ln 0))
(define ss (make-list ln 0)) (list-set! ss 0 80)
(set! b 0)
(do ((n 1 (+ n 2)) (ff #t (not ff))) ((> n 1433))
(ldiv ss (* 5 5)) (ldiv+- ts ss ff n))
(define ss (make-list ln 0)) (list-set! ss 0 (* 4 239))
(set! b 0)
(do ((n 1 (+ n 2)) (ff #f (not ff))) ((> n 423))
(ldiv ss (* 239 239)) (ldiv+- ts ss ff n))
(do ((n 0 (+ n 1))) ((= n 100))
(display (lmul ts d10)) (display " "))
</pre>
上のプログラムをPC-1のプログラムにしたのが次だ. 作業番地tsはeの1000桁と同様,
300番地以降を使う. プログラムの最初にあるインタールードでクリアしておく.<br><br>
ssは0番地から200番地までを使う. つまりイニシアルオーダーR0は壊す.
計算だけで記憶場所を使うので, 結果tsの出力は別のプログラムにする. <br><br>
まだコメントが少ないが, 徐々に追加するつもりだ.<br><br>
<!-- scheme/pc1sim/pi1000-6 -->
<pre>
0p=200,
0 0p=196,0p:p2, インタールード 0をaccへ
1 tl(300), 作業領域ts(300から498)をクリア
2 p1r,
3 a40, 2を足す
4 x1r,
5 s8r, 終了か
6 z40, R0へ戻る
7 jlr, ループする
8 500,jlp.
0 0p=196,0n=276,0m=496,0p:a7n,
1 x16r, atan サブルーチン 帰り番地を入れる
2 p4n, 作業場所ssをクリア
3 x7r,
4 s19n, 終了?
5 z11r,
6 p4n, 0をaccに
7 tl( ), 作業場所へ
8 p7r, アドレスを増やす
9 a7n,
10 jl3r,
11 pl10n, ssの最初の語に数を入れる
12 tl0,
13 pl4n, atan 計算ループの回数
14 tl2n,
15 sl8n, その回数が0になる迄
16 zl( ), 帰り番地
17 q4n, z<sup>2</sub>の除算ループ. 最初の剰余をRregに
18 p22n,
19 x23r, ssの語のアドレスを順に作る
20 x25r,
21 s19n, 終り?
22 z34r,
23 ql( ), 被除数をaccとRregに
24 dl12n, z^2による除算
25 tl( ), 商を格納
26 sl4n, 商が零かテスト 1を引き零なら負にする
27 kl31r,
28 p22n, 除算開始の位置b
29 a7n, bを2増やす
30 x22n,
31 p23r,
32 a7n, ssのアドレスを増やす
33 jl19r,
34 q4n,
35 p22n, ssをnで除し, tsに加減する
36 x43r,
37 a20n,
38 x46r,
39 x49r,
40 x52r,
41 s21n, ループの終か
42 z70r,
43 ql( ), 被除数を作り
44 dl2n, nで割る.
45 tln, 商を一旦しまう
46 pl( ), ts
47 aln, aln<->sln 商を足すか引く
48 k51r, sign bitが立っていたら繰上げ処理へ
49 tl( ), 商をしまってこの語は終り
50 jl67r,
51 cl14n, sign bitを消して
52 tl( ), しまう
53 p46r, 繰上げ処理
54 s7n,
55 x58r, 繰上げ先の語のアドレス
56 x61r,
57 x64r,
58 pl( ),
59 al4n, aln<->sln 繰上げを足すか引く
60 k63r, sign bitがあれば繰上げ
61 tl( ), 繰上げがなければ処理は終り
62 jl67r,
63 cl14n, sign bitを消す
64 tl( ), しまう
65 p58r, 繰上げ先のアドレスを進める
66 jl54r,
67 p43r, ssの語のアドレスを進める
68 a7n,
69 jl36r,
70 p79r, バイナリスイッチ切替え
71 b47r,
72 t47r,
73 p79r,
74 b59r,
75 t59r,
76 pl2n, nを増やす
77 al6n,
78 jl14r,
79 i, a xor s = i
0 0n:,, データ領域
2 ,,
4 ,1,
6 ,2,
8 ,1435,
10 ,80, =5*16 atan(1/5)で使う
12 ,25,
14 131071,262143, マスクパターン
16 425,956, 956=239*4
18 57121,194, 57121=239*239
20 300,494,
22 0,
0 0m:it, 主ルーチン
1 jlp, atan(1/5)
2 p16n, atan(1/239)用に
3 t9n, データを入替え
4 p17n,
5 t11n,
6 p18n,
7 t13n,
8 p4n,
9 t22n,
10 p4n,
11 t22n,
12 it,
13 jlp, atan(1/239)
14 0, 計算終り 停止する
jlm.
</pre>
このプログラムは私のPC-1シミュレータに追加してある. シミュレータを開き, R0
を読み込み, Place Tapeをクリック. メニューの下方にあるpi1000cmpを
クリックすると計算のプログラムが設定される. そこでFree Runをクリック;
さらにInit Startをクリックするとプログラムが走り出す.<br><br>
計算は私のMacBook Pro環境では約3分半で終り, シミュレータは停止する.<br><br>
結果の出力には作業場所
を同じにしてあるので, eの1000桁の出力プログラムを多少修正して使う. <br><br>
計算プログラムがR0領域を破壊したから, 出力プログラムの読込み前にもう一度
R0をロードする. 後はPlace TapeからPi1000prtを選び, そのプログラムを
読込み起動すればπの1000桁が出力される.<br><br>
久し振りにPC-1のプログラムを書いたが, インタプリタで必要な情報を取り出し
たり, Schemeのプログラムで予備テストが出来たり, 便利な世の中である.
和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-5801227873097392502021-03-29T19:02:00.000+09:002021-03-29T19:02:21.073+09:00PC-1シミュレータこのブログにパラメトロン計算機 PC-1のシミュレータの話を書いたのはそろそろ10年前の
ことだ. そのシミュレータはProcessingで書いたものだが, その後シミュレータの
置いてある計算機が引退したりで, 長いあいだ動かなかった. ブログを見た方には申し訳
ないことをした.<br><br>
ずーっとシミュレータの書き直しを考えていたが果せず, 最近やっとJavaScriptで書き,
どうやら無事に動くようなので, 公開することにした. <br><br>
シミュレータは<a href="https://www.iijlab.net/~ew/pc1sim2/pc1sim.html">
https://www.iijlab.net/~ew/pc1sim2/pc1sim.html</a>
にあるので試してみて欲しい.<br><br>
シミュレータの見掛けは以前のままである. もともとPC-1が出来たころ, 近隣の
計算機と競争していた, 自然対数の底eの値を1000桁計算するプログラムを
シミュレートするのが目的だったので, そのプログラムに特化してはあるが, 他の
プログラムも少数ながら動く.<br><br>
シミュレータを開くと次のような画面になる. 右手に3x3の押しボタンがある.<br><br>
<img src="http://www.iijlab.net/~ew/pc1sim2/pc1simface.png" width="400">
<br><br>
最初にやることは初期入力ルーチンR0の読込みで, まず左下のPlace R0をクリックする.
R0のテープが現れたら, 次にその上のInit Loadをクリック. これでR0が読み込まれる.<br><br>
今度は中下のPlace e1000をクリック. e1000桁のテープが現れる. そこで
中央のFree Runを押し, 続けて中央上のInit Startをクリック. これでe1000桁の
テープがR0で読み込まれ, 計算が始まる. 程なくして結果の印字が始まり, 計算が
終了する.<br><br>
他に用意してあるテープは,<br><br>
factorize32 2<sup>31</sup>より少し小さい整数の素因数分解<br>
eratosthenes エラトステネスの篩<br>
v2test フーリエ変換<br>
va1test 素数を法とする行列式の計算<br><br>
で, これらを走らせるには右下のPlace Tapeをクリックする. そうすると右上に
プログラム名が表示されるから, それをクリックする. テープが現れたら中央の
Free Run, その上のInit Startをクリック. なかにはテープ読込み中に停止する
ものもあるが, その時は中央のFree Runと, 右上のRestartをクリックする.<br><br>
</body>
</html>
和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-81125950521956713652020-12-10T14:37:00.003+09:002020-12-10T20:08:18.682+09:00無限クイーン<html>
<body>
TAOCPの演習問題7221-38に, 無限クイーン列<q<sub>n</sub>>の効率的計算法を考案せよというのがある.<br><br>
面白い解答なので, その説明をしたい.<br><br>
二進の配列 a,b,cを使う. cは正と負の添字がある.<br>
<b>G1.</b> [初期化.] r←0, s←1, t←0, n←0とする. (
1≤k≤nについて, q<sub>k</sub>の計算が済んだ.)<br>
<b>G2.</b> [q<sub>n</sub>≤nを調べる.] (この時点で1≤k<sについて
a<sub>k</sub>=1でa<sub>2</sub>=0; また-r<k≤tについてc<sub>k</sub>=1
でc<sub>-r</sub>=c<sub>t+1</sub>=0; 各ベクターにはn個の1がある.)
n←n+1, k←sとする.<br>
<b>G3.</b> [見つけた?] a<sub>k</sub>=b<sub>k+n</sub>=c<sub>k-n</sub>=0なら
G5へ進む. そうでないなら, k←k+1とし, k≤n-rならこのステップを繰り返す.<br>
<b>G4.</b> [q<sub>n</sub>>nを作る.] t←t+1; q<sub>n</sub>←n+t;
a<sub>k</sub>←b<sub>2n+t</sub>←c<sub>t</sub>←1とし,
G2へ戻る. <br>
<b>G5.</b> [q<sub>n</sub>≤nを作る.] q<sub>n</sub>←k;
a<sub>k</sub>←b<sub>k+n</sub>←c<sub>k-n</sub>←1とする.
k=sなら, a<sub>s</sub>=0になるまでs←s+1を繰り返す. k=n-rなら
c<sub>-r</sub>=0になるまでr←r+1を繰り返す. G2へ戻る. <br><br>
このプログラムを忠実に辿ると, 大体こういうことをやっていると判明する.<br><br>
TAOCPのプログラムでは, チェス盤の行も列も1から始めるが, 私はやはり0から
始めたい. 前回の無限クイーンの最後にあるように, 解の列は<br><br>
0 2 4 1 3 8 10 12 14 5 7 18 6 21 9 24 26 28 30 11 13 34
36 38 40 15 17 44 16 47
<br><br>
のように始まる. その始めの方の置きかたを見ると(下の図 前回の図では, 行は下に
延びていたのに, 今回は逆に上に延びるように描いて申し訳けない.)
<br>
<img src="http://www.iijlab.net/~ew/images/infqueens/infqueens3.png" width="400">
<br>
まず列n=0は行0にクイーンが置ける. それが左下の黒丸だ. そうするとどの行にはもう置けない
ことを示す配列aの0を1に, 行-列の対角線に置けないことを示す配列cの0-0=0を1にする.
図ではそれを黒丸のすぐ右の青小丸とすぐ右上の赤小丸で示す. 配列bも0+0=0が
1になるがそれは図には示さない. <br><br>
続いてn=1にクイーンを置くのだが, 青線が示すようにこの行には置けず, 赤の対角線が
示すようにこの斜線にも置けない. よって赤斜線の上のq=2に置くことになる. それが
(1,2)の黒丸で, その右の青小丸と右上の赤小丸で配列の設定も分る.<br><br>
n=2の列は, 行0は青線でだめ, 行2は赤線でだめ, その間の行1は, すぐ左上にクイーン
がいてだめ(これは配列bで分る). 行4も(1,2)からの斜線のためだめで, 結局
行4に置く. <br><br>
このようにして置いていくのだが, 青線はこれより下はすべて詰っている; 2本の赤線は
この間はすべて詰っていることを示す. したがって行0の青線は, 列3の行1に置くと
行2まで詰ったので, ここからは行2から横に進む. <br><br>
赤の斜線も同様に出来ている. 従ってある列でクイーンが置けるか調べるには, 青線路
の上から始め, 下の赤線の下までである. ここまでで配列a, b, cを調べ, 置けない
時は, 上の赤線の上に置くことになる.<br><br>
この方針で私流に書いたプログラムが次だ. <br>
<pre>
(define qs (make-list 30))
(define as (make-list 50 0))
(define bs (make-list 80 0))
(define cs (make-list 40 0))
(define coff 20)
(define (a k) (list-ref as k))
(define (b k) (list-ref bs k))
(define (c k) (list-ref cs (+ k coff)))
(define (a! k) (list-set! as k 1))
(define (b! k) (list-set! bs k 1))
(define (c! k) (list-set! cs (+ k coff) 1))
(define (place k) (list-set! qs n k)
(a! k) (b! (+ k n)) (c! (- k n)))
(define n 0) (define s -1) (define t 0) (define r 0)
(define k)
(while (< n 30) (set! k (+ s 1))
(while
(and (< k (+ n r))
(or (= (a k) 1) (= (b (+ k n)) 1) (= (c (- k n)) 1)))
(set! k (+ k 1)))
(if (and (= (a k) 0) (= (b (+ k n)) 0)
(= (c (- k n)) 0))
(begin (place k)
(while (= (a (+ s 1)) 1) (set! s (+ s 1)))
(while (= (c (- r 1)) 1) (set! r (- r 1))))
(begin (place (+ n t 1)) (set! t (+ t 1))))
(set! n (+ n 1)))
qs
=>(0 2 4 1 3 8 10 12 14 5 7 18 6 21 9 24 26 28 30 11 13
34 36 38 40 15 17 44 16 47)
</pre>
変数sは青線を, rは下の赤線を, tは右の赤線を示す. これは最初のTAOCPのプログラムと
同じである. 変数kの使い方も同様. <br><br>
下の図は, 上のと同じだが, nの範囲を広げた. また配列を調べた位置を白丸で示した.
白丸の数が少ないことに注目して欲しい. TAOCPのプログラムは判定の方法が違うので,
白丸の数はこれよりかなり多い. <br><br>
<img src="http://www.iijlab.net/~ew/images/infqueens/infqueens4.png" width="400">
<br>
</body></html>
和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-14182628215202892742020-11-22T11:42:00.003+09:002021-04-23T08:12:19.385+09:00mex関数<html>
<body>
アルゴリズムを読んでいると, mexに出会うことがある. mexとは
minimum excludantのことで, 最小の排除されたもの, つまり
引数の集合の最小不在者である. 例えば<br>
<pre>
(mex '()) => 0
(mex '(1 2)) => 0
(mex '(0 2)) => 1
(mex '(0 1)) => 2
</pre>
前回のブログ, 無限クイーンで説明したbit stringのある処理系では,<br>
<pre>
(define (mex ns)
(if (null? ns) 0
(let* ((l (+ (apply max ns) 2))
(abs (make-bit-string l #t)))
(for-each (lambda (n) (bit-string-clear! bs n)) ns)
(bit-substring-find-next-set-bit bs 0 l))))
</pre>
nsが(0 1 1 2 3 5 8 13)の時, <t>(mex ns)</t>では. nsのmaxが13だからlは15.
従って1が15個のbit stringが出来, 最初の(右端の)ビットのビット番号は0, 最後の(左端の)
ビット番号は14だ.<br>
<pre>
#*111111111111111
</pre>
次に位置0, 1, 1, 2,..., 13のビットをクリアして0にすると
<pre>
#*101111011010000
</pre>
になる. そして0から15の範囲で最初の1の位置を探すと4であり,
<pre>
(mex ns) => 4
</pre>
上のmex関数はこのように動作するが, bit stringの最初の1を探すような飛道具を
使っていて, 快しとしない向きには, 以下の二進木(二分木)を使うアルゴリズムはどうだろうか.<br><br>
二進木の節点は, 左部分木, その節点の情報, 右部分木 で構成される. 今の
アルゴリズムの場合, その節点の情報は範囲[l, r]で, mexをとる集合nsに
l<=k<=r要素があることを示す. 二進木に含まれる節点の範囲には重なりがないだけ
ではなく, 隣同士の範囲の間には隙間があるようになっている.<br><br>
隣同士の範囲の間に新しい要素が現れてその隙間を埋めると, 二つの範囲は合体して
一つの範囲になる.
<pre>
ns=(12 1 3 3 4 4 5 11 6 6 10 14 15 13 4 7)
</pre>
として, 二進木の出来かたを下の図に示す.<br><br>
<img src="http://www.iijlab.net/~ew/images/binarymex.png" width="400">
最初の12が来た時に, 範囲[12,12]が出来る(図の右上.) その右の12は, この時の
入力が12であったことを示す.<br><br>
次に1が来ると, 1は既にある範囲[12,12]の中でもないし, 両端の続きでもなく離れ
ているので, 別の範囲[1,1]を作り, [12,12]の節点の左部分木とする(先程の図の下.) <br><br>
その次の3は, 新しい範囲[3,3]になって, [12,12]の左の[1,1]の右に繋る.<br><br>
次の3は[3,3]の範囲にあるから, 何もしない.<br><br>
4は[3,3]の右に接しているから, 範囲[3,3]が[3,4]になる. しかし, その右部分木
の最も左に5がないので, 合体は生じず, 仕事はここで終わる. <br><br>
図の1列目が終り, 中央の列の上に戻り, 5が来ると, 範囲[3,4]が[3,5]になる.<br><br>
11が来ると, [12,12]が[11,12]になり, ここでも合体しない.<br><br>
更に進んで右の列になり, 14,15,13と来ると, 13は[10,12]の
右を13にし, 右の節点に[14,15]があるのでこの二つが合体して,
[10,13]の範囲が出来る.<br><br>
こうしてnsを全て処理した二進木が右下の図になる.<br><br>
この二進木を中順序で辿り, 最初に来る節点の範囲[1,1]を見る.<br><br>
この左が0でなければ, mexは0, 0ならば右プラス1がmexである.<br><br>
二進木の操作になれていれば, プログラムを書くのは造作ない.<br><br>
なお, nsが空なら, 木も空で, mexは0である.<br><br>
P.S. Knuth先生へのメイルでこの実装にも触れたら返事に<br>
It is attractive. I don't recall seeing it before.<br>
と書いてあった.
</body></html>
和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-50337709654210252362020-11-10T15:34:00.002+09:002020-12-10T14:48:37.830+09:00無限クイーンTAOCPのV4F5に,
<pre>
1,3,5,2,4,9,11,13,15,6,8,19,7,22,10,25,27,29,31,12,14,35,37,39,
41,16,18,45,...
</pre>
という数列があり, これは「無限クイーン」の辞書式順で最小の解とあった(演習問題 7221-42). これを計算してみたというのが今回の話題だ. <br><br>
情報科学の標準問題のひとつに8クイーンがある. チェス盤の上に, 飛車と角行を併せて動くクイーンを8個, 互いに当らないように置く方法を求める問題だ. その8を無限大にしたのが「無限クイーン」である. つまり縦横とも半無限に広いチェス盤に, 孫悟空の分身の術のように無限に現れるクイーンを置くのである.<br><br>
釈迦に説法かも知れないが, 順序としてまず8クイーンの解法を復習しよう. 次の図を見て欲しい.
<br><br>
8×8のチェス盤の列にも行にも0から7の番号をつける. そして左端の列0からクイーンの
置けるところを探す. 既に置いてあるクイーンと当らないようにということは, 置く場合,
その行には他のクイーンはない; その場所の右上がりの筋には他のクイーンはいない.
その場所の右下がりの筋には他のクイーンはいない. という条件が揃えばそこに
新しいクイーンが置けて, さらに右の列の検討が始まる.<br><br>
その列のどの行にも上の条件の場所がなければ, 1列戻って, その列のクイーンを更に下に
置けるかどうか調べる. こうして8個置くことが出来れば解である. 一つの解が見付かっても,
まだ他の解を探すのが全解探索だ.<br><br>
上のような戦術で8クイーンを探すプログラムを示すと次のようになる. <br><br>
この行, この右上がりの筋, この右下がりの筋にまだクイーンが置けることを示す配列
as, bs, csを用意する. 真理値が入る. 最初はどこにでも置けるから, 初期値は真だ.
<br><br>
列の添字をm, 行の添字をnとすると, 右上がりの筋はm+n=一定, 右下がりはm-n=一定で,
mもnも0から7の値をとるから, m+nの値は0から14, m-nの値は-7から7である.
配列の添字は0以上なので, 右下がりはm-n+7の場所に置く.<br><br>
各列の置き方を調べる関数が col m で, その下請けの各行の置き方を調べるのが row n
である.<br><br>
<img src="http://www.iijlab.net/~ew/images/infqueens/eightqueens.png" width="400">
<pre>
(define qs (make-list 8 '()))
(define as (make-list 8 #t))
(define bs (make-list 15 #t))
(define cs (make-list 15 #t))
(define (col m)
(define (row n)
(if (< n 8) (begin
(if (and (list-ref as n) (list-ref bs (+ m n))
(list-ref cs (- m n -7)))
(begin (list-set! qs m n) (list-set! as n #f)
(list-set! bs (+ m n) #f)
(list-set! cs (- m n -7) #f) (col (+ m 1))
(list-set! cs (- m n -7) #t)
(list-set! bs (+ m n) #t) (list-set! as n #t)
(list-set! qs m '())))
(row (+ n 1)))))
(if (< m 8) (row 0) (display qs)))
(col 0)
</pre>
このプログラムはバックトラックを繰り返す. 最初の解が得られるまでのバックトラック
の様子を見ると下のようだ. <br><br>
<img src="http://www.iijlab.net/~ew/images/infqueens/backtrack.png" width="400">
<br><br>
そうして得られた92通りの解答は次の通り.
<pre>
(0 4 7 5 2 6 1 3)(0 5 7 2 6 3 1 4)(0 6 3 5 7 1 4 2)
(0 6 4 7 1 3 5 2)(1 3 5 7 2 0 6 4)(1 4 6 0 2 7 5 3)
(1 4 6 3 0 7 5 2)(1 5 0 6 3 7 2 4)(1 5 7 2 0 3 6 4)
(1 6 2 5 7 4 0 3)(1 6 4 7 0 3 5 2)(1 7 5 0 2 4 6 3)
(2 0 6 4 7 1 3 5)(2 4 1 7 0 6 3 5)(2 4 1 7 5 3 6 0)
(2 4 6 0 3 1 7 5)(2 4 7 3 0 6 1 5)(2 5 1 4 7 0 6 3)
(2 5 1 6 0 3 7 4)(2 5 1 6 4 0 7 3)(2 5 3 0 7 4 6 1)
(2 5 3 1 7 4 6 0)(2 5 7 0 3 6 4 1)(2 5 7 0 4 6 1 3)
(2 5 7 1 3 0 6 4)(2 6 1 7 4 0 3 5)(2 6 1 7 5 3 0 4)
(2 7 3 6 0 5 1 4)(3 0 4 7 1 6 2 5)(3 0 4 7 5 2 6 1)
(3 1 4 7 5 0 2 6)(3 1 6 2 5 7 0 4)(3 1 6 2 5 7 4 0)
(3 1 6 4 0 7 5 2)(3 1 7 4 6 0 2 5)(3 1 7 5 0 2 4 6)
(3 5 0 4 1 7 2 6)(3 5 7 1 6 0 2 4)(3 5 7 2 0 6 4 1)
(3 6 0 7 4 1 5 2)(3 6 2 7 1 4 0 5)(3 6 4 1 5 0 2 7)
(3 6 4 2 0 5 7 1)(3 7 0 2 5 1 6 4)(3 7 0 4 6 1 5 2)
(3 7 4 2 0 6 1 5)(4 0 3 5 7 1 6 2)(4 0 7 3 1 6 2 5)
(4 0 7 5 2 6 1 3)(4 1 3 5 7 2 0 6)(4 1 3 6 2 7 5 0)
(4 1 5 0 6 3 7 2)(4 1 7 0 3 6 2 5)(4 2 0 5 7 1 3 6)
(4 2 0 6 1 7 5 3)(4 2 7 3 6 0 5 1)(4 6 0 2 7 5 3 1)
(4 6 0 3 1 7 5 2)(4 6 1 3 7 0 2 5)(4 6 1 5 2 0 3 7)
(4 6 1 5 2 0 7 3)(4 6 3 0 2 7 5 1)(4 7 3 0 2 5 1 6)
(4 7 3 0 6 1 5 2)(5 0 4 1 7 2 6 3)(5 1 6 0 2 4 7 3)
(5 1 6 0 3 7 4 2)(5 2 0 6 4 7 1 3)(5 2 0 7 3 1 6 4)
(5 2 0 7 4 1 3 6)(5 2 4 6 0 3 1 7)(5 2 4 7 0 3 1 6)
(5 2 6 1 3 7 0 4)(5 2 6 1 7 4 0 3)(5 2 6 3 0 7 1 4)
(5 3 0 4 7 1 6 2)(5 3 1 7 4 6 0 2)(5 3 6 0 2 4 1 7)
(5 3 6 0 7 1 4 2)(5 7 1 3 0 6 4 2)(6 0 2 7 5 3 1 4)
(6 1 3 0 7 4 2 5)(6 1 5 2 0 3 7 4)(6 2 0 5 7 4 1 3)
(6 2 7 1 4 0 5 3)(6 3 1 4 7 0 2 5)(6 3 1 7 5 0 2 4)
(6 4 2 0 5 7 1 3)(7 1 3 0 6 4 2 5)(7 1 4 2 0 6 3 5)
(7 2 0 5 1 4 6 3)(7 3 0 2 5 1 6 4)
</pre>
ところで, MIT Schemeにはbit stringというデータの型があり, Eratosthenesの篩の
ような真理値の配列に便利である. 今回使う関数を説明すると<br>
(make-bit-string 8 #t)<br>
各ビットが1の8ビットのbit stringを生成する. #fなら0になる.<br>
(bit-string-ref bs k)<br>
bit string bsのk番目のビットを調べる. 0なら#f, 1なら#tが返る.<br>
(bit-string-set! as k)<br>
bit string asのk番目のビットを1にする.<br>
(bit-string-clear! as k)<br>
bit string asのk番目のビットを0にする. bit stringにはある範囲で
最初の1の位置を探す飛道具がある. <br>
(bit-string-find-next-set-bit as k 8)<br>
bit string asの[k .. 8)の間で, 最初の1であるビットの位置を返す. 1がなければ#fを返す.<br>
(bit-string-append as bs)<br>
bit string asとbsを接続した新しいbit stringを返す. as内のビット番号は不変. 新しい
bsの部分の番号は前の番号にasの長さを足したものになる.
<br><br>
bit stringを使った8クイーンのプログラムは次のようだ. as, bs, csはどの位置にも
クイーンが置けるから#t, つまり1にしておく. bit stringからビットを取り出すのは
(b k)はbsのkビットを返す. (as k)はasのkからの最初の#tの位置を返す.
ビットの設定は(a~ k)はasのkビットを0に, (a! k)はasのkビットを1にする.
<pre>
(define qs (make-list 8 '()))
(define as (make-bit-string 8 #t))
(define bs (make-bit-string 15 #t))
(define cs (make-bit-string 15 #t))
(define (a k) (bit-substring-find-next-set-bit as k 8))
(define (b k) (bit-string-ref bs k))
(define (c k) (bit-string-ref cs (+ k 7)))
(define (a! k) (bit-string-set! as k))
(define (b! k) (bit-string-set! bs k))
(define (c! k) (bit-string-set! cs (+ k 7)))
(define (a~ k) (bit-string-clear! as k))
(define (b~ k) (bit-string-clear! bs k))
(define (c~ k) (bit-string-clear! cs (+ k 7)))
(define (col m)
(define (row n)
(let ((n1 (a n)))
(if n1 (begin
(if (and (b (+ m n1)) (c (- m n1))) (begin
(list-set! qs m n1) (a~ n1) (b~ (+ m n1))
(c~ (- m n1)) (col (+ m 1)) (c! (- m n1))
(b! (+ m n1)) (a! n1) (list-set! qs m '())))
(row (+ n1 1))))))
(if (< m 8) (row 0) (display qs)))
(col 0)
</pre>
これで様子が分ったから, 無限クイーンに取り掛かろう. 問題はcのbit stringで,
0がbit stringの途中にあるので, オフセットを足すことだ.
無限クイーンでいいのは, 列の中で次に置ける位置はどんどん先まで探せばよく, バックトラック
しないことだ. 従って#tに戻す(a! k)のような関数はない. 無限といっても無限まで
計算することはないので, 列の最大番号はきめることにした. mmaxの値である.
するとas, bs, csなどの長さも余裕をもって決められる.
<pre>
(define mmax 256)
(define nmax (* mmax 2))
(define qmax mmax)
(define amax nmax)
(define bmax (+ mmax nmax))
(define cmax (+ mmax nmax))
(define qs (make-list qmax '()))
(define as (make-bit-string amax #t))
(define bs (make-bit-string bmax #t))
(define cs (make-bit-string cmax #t))
(define (a k)
(bit-substring-find-next-set-bit as k amax))
(define (b k) (bit-string-ref bs k))
(define (c k) (bit-string-ref cs (+ k nmax)))
(define (a~ k) (bit-string-clear! as k))
(define (b~ k) (bit-string-clear! bs k))
(define (c~ k) (bit-string-clear! cs (+ k nmax)))
(define (q! m n) (list-set! qs m n))
(define (col m)
(define (row n)
(let ((n1 (a n)))
(if (and (b (+ m n1)) (c (- m n1)))
(begin (q! m n1) (a~ n1) (b~ (+ m n1))
(c~ (- m n1)) (col (+ m 1)))
(row (+ n1 1)))))
(if (< m mmax) (row 0) (display (take m qs))))
(col 0)
</pre>
実行結果を下に示す.
<pre>
(0 2 4 1 3 8 10 12 14 5 7 18 6 21 9 24 26 28 30 11 13 34
36 38 40 15 17 44 16 47 19 50 52 20 55 57 59 22 62 23 65
27 25 69 71 73 75 77 29 31 81 83 85 32 88 33 91 37 35 95
97 99 101 39 104 106 41 109 42 112 43 115 117 119 45 122
46 49 126 48 129 131 133 135 51 53 139 141 143 54 56 147
149 151 58 154 156 158 60 161 61 64 165 63 168 170 172
66 175 177 67 180 68 183 185 70 74 189 191 72 194 196 76
199 201 203 205 78 80 209 211 79 82 215 217 84 220 222
224 226 86 229 87 90 233 89 236 238 240 242 92 94 246 93
249 96 252 254 256 98 259 261 100 264 266 268 102 271
103 274 107 105 278 280 282 284 108 110 288 290 292 111
113 296 298 300 114 116 304 306 308 310 118 120 314 316
318 121 123 322 324 326 124 329 125 128 333 127 336 338
340 342 130 132 346 348 350 352 134 136 356 358 360 137
363 138 366 142 140 370 372 374 376 378 144 146 382 384
145 148 388 390 150 393 395 397 399 152 402 153 405 407
155 159 411 413)
</pre>
<br><br>
TAOCPには<i>q<sub>n</sub></i>は<i>nφ</i>+<i>O</i>(1)か
<i>n/φ</i>+<i>O</i>(1)と書いてあったので, 図にしてみた.
赤線が黄金比の傾きの斜線だ.
<br><br>
<img src="http://www.iijlab.net/~ew/images/infqueens/infqueens.png" width="400">
<br><br>
上のプログラムは, as, bs, などが最初から固定長であった. これを最初はある程度に
とり, 必要に応じて長くするように書いたのが次のプログラムだ. 上のと違う部分だけ
書いておく.<br><br>
<pre>
(define mmax 8) (define nmax (* mmax 2))
(define qmax nmax) (define amax nmax)
(define bmax (+ mmax nmax)) (define cmax (+ mmax nmax))
(define coff nmax)
(define qs (make-list qmax '()))
(define as (make-bit-string amax #t))
(define bs (make-bit-string bmax #t))
(define cs (make-bit-string cmax #t))
(define (a k) (let
((k1 (bit-substring-find-next-set-bit as k amax)))
(if k1 k1 (begin
(set! as (bit-string-append as
(make-bit-string amax #t)))
(set! amax (* amax 2)) (a k)))))
(define (b k) (if (>= k bmax) (begin
(set! bs (bit-string-append bs
(make-bit-string bmax #t)))
(set! bmax (* bmax 2))))
(bit-string-ref bs k))
(define (c k)
(if (or (>= (+ k coff) cmax) (< (+ k coff) 0)) (begin
(set! cs (bit-string-append
(bit-string-append (make-bit-string (/ cmax 2) #t)
cs) (make-bit-string (/ cmax 2) #t)))
(set! coff (+ coff (/ cmax 2)))
(set! cmax (* cmax 2))))
(bit-string-ref cs (+ k coff)))
(define (q! m n) (if (>= m qmax) (begin
(set! qs (append qs (make-list qmax '())))
(set! qmax (* qmax 2))))
(list-set! qs m n))
この下 (a~ k), (b~ k)からは同じ
</pre>
</body>
</html>
和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-57655992356381472792020-10-19T17:01:00.003+09:002020-12-10T14:53:02.152+09:00五角ミノ十二面体<html>
<body>
英語でquintominal dodecahedraというクイズである. 十二面体のクイズといえば, Hamilton順路が有名だが, これは全然別だ. 今年4月に新型コロナで死去したJohn Conwayが1958年頃考えたらしい.<br><br>
Conway他の書いた"Winning Ways Vol.4"の表紙に下のような興味ある図があった.<br><br>
<img src="http://www.iijlab.net/~ew/images/quintomino/quintominofig0.png" width=400><br>
正五角形の5辺を5色(例えば赤青橙緑黒)で塗ると, 回転, 裏返しは同じと見た時, 12通りの五角形が得られる. (5!/5/2=12) 図の下にある12の五角形だ. この一組をquintomino(「五角ミノ」と訳すのは如何?)という. 一方, 正十二面体は12の正五角形か構成されるから, 辺が同じ色になるように, 各面に五角ミノを1枚づつ貼れるかというのが問題である.<br><br>
Winning Waysの表紙は, 12の五角ミノにIを除いてAからMまで名をつけ, 解答を描いたものである. 各五角ミノの色の配置は以下の通り.<br><br>
A=01234 B=01243 C=01324 D=01342 E=01423 F=01432<br>
G=02134 H=02143 J=02314 K=02413 L=03124 M=03214<br><br>
答をいうとこれは可能であり, 等価性を考慮すると3通りの解があるという. 表紙の3つの絵はAを底面にしたその3通りを示している.<br><br>
KnuthのTAOCP, 第4巻, 第5分冊はDancing Linksが話題で, 膨大な数の演習問題があるが, その中にquintominal dodecahedron問題を解けというのがあった(ex 7221-136)ので, 挑戦してみたというのが今回のブログである.<br><br>
正十二面体は下左のようなもので, こういう立体の絵では扱いにくいから, 通常は下右のような形相図(高木貞治 数学小景)を使う. 数学小景には正十二面体を床に置き, 上面の正五角形の少し上の光源で照すと, こういう影ができると書いてある.<br><br>
<img src="http://www.iijlab.net/~ew/images/quintomino/quintominofig1.png" width=401><br>
正十二面体には, 頂点(青)と辺(赤)と面(緑)に番号を付けたが, 同じ番号を面0を底にした形相図に付けると下のようになる. 右端は面にも番号を付けたものだ. 0番の面がないのは, それは周囲の大きい五角だからだ.<br><br>
<img src="http://www.iijlab.net/~ew/images/quintomino/quintominofig2.png" width=402><br>
ペントミノのとか8クィーンのような, ピースを置けるかどうかの問題は, 現在置いてあるピースの状態と, これから置くピースの集合を引数として受取り, あるピースがさらに置けるなら, その状態と今残るピースを新しい引数にして自分を再帰呼出しする関数を書くと解けるのが普通である.<br><br>
以下のSchemeのプログラムもその方針で出来ている. colorsは12の五角ミノのリスト. facesは各々の面を作る辺の番号のリストである. allcsは色のリストを回転したり裏返したりしたもののリストを作る.<br><br>
<pre>
(define colors '
((0 1 2 3 4)(0 1 2 4 3)(0 1 3 2 4)(0 1 3 4 2)(0 1 4 2 3)
(0 1 4 3 2)(0 2 1 3 4)(0 2 1 4 3)(0 2 3 1 4)(0 2 4 1 3)
(0 3 1 2 4)(0 3 2 1 4)))
(define faces '
((0 1 2 3 4)(0 5 6 7 8)(1 9 10 11 5)(2 12 13 14 9)(3 15 16 17 12)
(4 8 18 19 15)(7 20 21 22 18)(6 11 23 24 20)(10 14 25 26 23)
(13 17 27 28 25)(16 19 22 29 27)(21 24 26 28 29)))
(define (allcs cs)
(let ((ccs '()))
(do ((i 0 (+ i 1))) ((= i 5))
(set! ccs (cons cs ccs))
(set! cs (append (cdr cs) (list (car cs)))))
(set! cs (reverse cs))
(do ((i 0 (+ i 1))) ((= i 5) (reverse ccs))
(set! ccs (cons cs ccs))
(set! cs (append (cdr cs) (list (car cs)))))))
</pre>
状態stateは辺の番号順の30の色のリストで, 先頭には最初の五角ミノAを固定し, 未定の25の場所は()にする. 数のリストである状態を数字列に変換したり, 逆に変換したりするのがcompressとexpandである.<br><br>
<pre>
(define (compress col)
(list->string (map (lambda (n) (integer->char (+ n 48))) col)))
(define (expand col)
(map (lambda (n) (- (char->integer n) 48)) (string->list col)))
(define state (append '(0 1 2 3 4) (make-list 25 '())))
</pre>
次のtryが再帰呼出しの関数で, 状態stateと残りの五角ミノrestを受け取る. matchはその状態でposの場所がcolの色と一致しているか()かであれば真を返す. 要するに置けるということである. fillcはposの場所をcolの色にした新しい状態を返す. 新しい状態は再帰呼出しから戻るときに捨てるから, 置いたピースを戻す仕事はしない.<br><br>
ifから後がtryの本体で, 残りの五角ミノがなければ((length(rest))=0なら)解があったとして出力する. 五角ミノの1個をとってcolとし, 残りをnewrestにし, 今度置く面をposにし, colを回転, 反転してcolsにし, そのそれぞれcについて置けるようなら再帰呼び出しする.<br><br>
<pre>
(define count 0)
(define (try state rest)
(define (match state col pos)
(every (lambda (c p) (let ((d (list-ref state p)))
(or (null? d) (= d c)))) col pos))
(define (fillc state col pos)
(let ((newstate (list-copy state)))
(for-each (lambda (c p)
(list-set! newstate p c)) col pos) newstate))
(if (= (length rest) 0)
(begin (display (compress state)) (newline)
(set! count (+ count 1)))
(for-each (lambda (n)
(let* ((col (list-ref colors n))
(newrest (delete n rest))
(pos (list-ref faces (- 12 (length rest))))
(cols (allcs col)))
(for-each (lambda (c)
(if (match state c pos)
(try (fillc state c pos) newrest))) cols))) rest)))
</pre>
Aが置いてあるので0を除き, restを(1 2 ... 11)としてtryを呼ぶと60個数の解が次のように得られる.<br><br>
<pre>
(newline)
(try state (range 1 12))
(display count)
"012343421042143042321400312301"
"012343421402013214300141324302"
"012343421420130204031343214201"
"012344231023134204031240413201"
"012344231320410012234101334402"
"012342431340401201031242330241"
"012342431430031214031242304210"
"012342431430130240031422341120"
"012342431430130204032141324102"
"012342431430130024231042314021"
"012342431403103204031242014231"
"012343241024134042321203001143"
"012343241024143240301020314321"
"012343241024143042320211304310"
"012343241024143042233010114032"
"012343241024413012321200341304"
"012343241420103042321204310341"
"012342341034431021231300224401"
"012342341340140042321202413301"
"012342341430130204032314112204"
"012342341304041241300212134203"
"012344312032431021314200124403"
"012344312320410102034321234401"
"012344312302014124300241423301"
"012342143034134240100232310241"
"012342143340041241102033402231"
"012342143304401012120232334410"
"012342143304014214100232331204"
"012342143304041241013122034023"
"012342143304041241100323224301"
"012342143304041142200132334102"
"012342143403031241100232443201"
"012343412042413120300241320341"
"012343412402103024310241324310"
"012343412402013142300421342301"
"012343412402013124033241024031"
"012343412402031124300243124103"
"012343412420013124302041304321"
"012344132032431012134020342014"
"012344132032431120034210324410"
"012344132032431021130244320041"
"012344132032134024134200321104"
"012344132023431021132404003421"
"012344132302401021134203024431"
"012342413340410102022431332401"
"012342413430130204013242114203"
"012342413403013124200341224301"
"012343142042431021133200423401"
"012343142024143042312100314302"
"012343142402013124300214321304"
"012344213023413102024231004431"
"012344213320401102024233140134"
"012344213320410102023241443301"
"012344213320410120204031334421"
"012344213320014142024231330041"
"012344213320140240014231334210"
"012344123032431201014130324402"
"012344123023143042124300214301"
"012344123320410102023414231403"
"012344123302041241100434332201"
</pre>
後ではこの文字列のリストをallsolsとして使う.<br><br>
ところで, Conwayの解は3通りであった. 上の60個の中には, 色の置換えについて等価のものが沢山あるということだ. 次に等価を発見する作業をやってみよう.<br><br>
上の最初の解"012343421042143042321400312301"は形相図に書き入れてみると, 下左のようになる. その3と4を交換すると中の図になる. 周囲は01234ではなくなるが, よく見ると下の五角形(面の番号1)の周囲が01234である. そこで辺の接続関係を保ったまま, これが周囲(面0)に来るように図を書き直すと右のようになる. この列を圧縮したのが"012343421420130204031343214201"で, 先程の解の2番目であった. だからこの2つの解は等価であった. これを上の全ての解のすべての置換について実行し, 同じになったものをまとめると3通りになる.<br><br>
<img src="http://www.iijlab.net/~ew/images/quintomino/quintominofig3.png" width=404><br>
以下はそのプログラムだ.<br><br>
repはcolのprm0をprm1で置き換える. つまりcolの(0 1 2 3 4)を(0 1 2 4 3)で置き換える(3と4を交換する).<br><br>
<pre>
(define (rep prm0 prm1 col) ;col expanded form
(define (pair as bs)
(if (null? as) as
(cons (list (car as) (car bs))
(pair (cdr as) (cdr bs)))))
(define (subst ps ls)
(if (null? ls) ls
(cons (cadr (assoc (car ls) ps))
(subst ps (cdr ls)))))
(subst (pair prm0 prm1) col))
</pre><br>
(0 1 2 3 4 3 4 2 1 0 4 2 1 4 3<br>
0 4 2 3 2 1 4 0 0 3 1 2 3 0 1))<br>
は<br>
(0 1 2 4 3 4 3 2 1 0 3 2 1 3 4<br>
0 3 2 4 2 1 3 0 0 4 1 2 4 0 1)<br>
になって左と中の絵に対応する.<br><br>
remapは置き換えた色のリストを, (0 1 2 3 4)が先頭になるように(面0になるように)位置を変更した辺のリストを作る. newposとしよう. これを基本位置という. つまり右の絵にする. これは多少面倒だ.<br><br>
まず後述の関数newfaceを置換後の色のリストに使うと, colorsに対応する五角ミノがいる位置が得られる. このリストをnewcolとすると, newcolは<br><br>
<pre>
((0 8 7 6 5) (0 1 2 3 4) (28 25 13 17 27)
(28 29 21 24 26) (9 1 5 11 10) (22 29 27 16 19)
(9 2 12 13 14) (23 26 25 14 10) (23 11 6 20 24)
(15 19 18 8 4) (22 21 20 7 18) (15 16 17 12 3))
</pre>
これを見ると五角ミノ(0 1 2 3 4)は(0 1 2 3 4)にいたのに(0 8 7 6 5)に, (0 1 2 4 3)は(0 8 7 6 5)にいたのに(0 1 2 3 4)に, (0 1 3 2 4)は(1 9 10 11 5)にいたのが(28 25 13 17 27)に, (0 1 3 4 2)は(22 29 27 16 19)にいたのに(28 29 21 24 26)にいった.<br><br>
newposを作る方針は次のようだ. 基本位置に来る五角ミノの現在位置を設定する. (0 1 2 3 4)は(0 8 7 6 5)にいたので, newposは(0 8 7 6 5 ...)となる(左の絵).<br><br>
<img src="http://www.iijlab.net/~ew/images/quintomino/quintominofig4.png" width=404><br>
先頭の(0 8 7 6 5)の0と1の位置にある0の辺と8の辺の隣の面を探す. つまり(cdr newcol)中の0と8を持つリストを探すと(0 1 2 3 4)と(15 19 18 8 4)が見付かるが, これらに共通の4が辺0と8から別れるもう1本の辺と分る. これはnewposの5の位置なので, そこに4を置く.<br><br>
(0 8 7 6 5 4 ...) まで判明した(中の絵). すると後は判明した場所を隣同士にもつ五角ミノを探す作業を続ける.<br><br>
まずnewposの(0 5)の位置の(0 4)を持つ五角ミノ, (0 1)の位置の(0 8)を持つ五角ミノを探すと
(0 1 2 3 4)と(15 19 18 8 4)が見付かるので, それらの入る(6 7 8)の位置に(2 3 4), (11 10 9)の位置に(15 19 18)が入り, 右の絵になる. 以後同様にして各辺が決る. (consfの下請け関数cfsが担当する.)<br><br>
<pre>
(define (remap newcol)
(compress (map (lambda (n) (list-ref newcol n))
(consf (newface newcol)))))
(define (newface col)
(define (sel ps ls)
(map (lambda (p) (list-ref ls p)) ps))
(append-map (lambda (c)
(filter (lambda (f) (equal? c (sel f col)))
(apply append (map (lambda (a) (allcs a)) faces))))
colors))
(define (consf col)
(let ((newpos (append (car col) (make-list 25 '()))))
(define (csf a b c d e)
(define (ls n f m)
(list-set! newpos n (list-ref f m)))
(define (search a b)
(car (filter (lambda (c) (subset? (list a b) c))
col)))
(define (reshape a b l)
(define (left a l)
(cadr (member a (cons (car l) (reverse l)))))
(if (= b (left a l)) (set! l (reverse l)))
(let ((p (elemindex a l)))
(append (drop p l) (take p l))))
(let* ((fa (list-ref newpos a))
(fb (list-ref newpos b))
(ff (reshape fa fb (search fa fb))))
(ls c ff 2) (ls d ff 3) (ls e ff 4)))
(let* ((f (car col)) (f0 (car f)) (f1 (cadr f))
(g (car (filter (lambda (a) (member f0 a))
(cdr col))))
(h (car (filter (lambda (a) (member f1 a))
(cdr col))))
(n (car (intersection g h))))
(list-set! newpos 5 n)
(csf 0 5 6 7 8) (csf 1 5 11 10 9) (csf 2 9 14 13 12)
(csf 3 12 17 16 15) (csf 4 15 19 18 8)
(csf 7 18 22 21 20) (csf 6 11 23 24 20)
(csf 10 14 25 26 23) (csf 13 17 27 28 25)
(csf 16 19 22 29 27)
newpos)))
</pre>
prmは(0 1 2 3 4)の順列. prem0はその先頭(0 1 2 3 4). equivは等価の解をまとめたリストのリスト. doループで実行が始まる.<br><br>
doループではallsolsから解を1個とりcolとし, equivの表のどれかに既にあればなにもしない. どれにもないなら, そrwれを1個もつリストをequivに追加し, そのリストをeqlistとする. colを120通りの色の置換えをし, 置き換えたnewcolがqlistになければ追加する.<br><br>
<pre>
(define prm (permutation (range 0 5)))
(define prm0 (list-ref prm 0))
(define equiv '())
(do ((j 0 (+ j 1))) ((= j 60))
(let* ((col (list-ref allsols j)))
(if (every (lambda (a) (not (member col a))) equiv)
(begin (set! equiv (cons (list col) equiv))
(let ((eqlist (filter (lambda (a) (member col a))
equiv)))
(do ((i 1 (+ i 1))) ((= i 120))
(let ((newcol (remap (rep prm0 (list-ref prm i)
(expand col)))))
(if (not (member newcol (car eqlist)))
(set-cdr! (last-pair (car eqlist))
(list newcol))))))))))
</pre>
それぞれの等価グループの長さを出力し, 各グループを昇順に並べ替え出力する.<br><br>
<pre>
(display (map length equiv))
(define (formprint ls)
(for-each (lambda (a) (display a) (newline)) ls))
(for-each (lambda (a) (newline)
(formprint (sort a string<?))) equiv)
(20 10 30)
012342143034134240100232310241 ;20 taocp B
012342143304041241013122034023
012342143304401012120232334410
012342143340041241102033402231
012342341034431021231300224401
012342341340140042321202413301
012342341430130204032314112204
012342431340401201031242330241
012342431430130240031422341120
012343241024134042321203001143
012343241024143042233010114032
012344123023143042124300214301
012344123032431201014130324402
012344123320410102023414231403 ;winning C
012344132023431021132404003421
012344132032431012134020342014
012344213320140240014231334210
012344213320401102024233140134
012344231023134204031240413201
012344231320410012234101334402
012342413403013124200341224301 ;10
012343142402013124300214321304 ;taocp C
012343412042413120300241320341
012343412402013124033241024031
012343412402013142300421342301
012343412402031124300243124103
012343412402103024310241324310
012343412420013124302041304321
012343421402013214300141324302
012344312302014124300241423301 ;winning A
012342143304014214100232331204 ;30
012342143304041142200132334102 ;winning B
012342143304041241100323224301
012342143403031241100232443201
012342341304041241300212134203
012342413340410102022431332401
012342413430130204013242114203
012342431403103204031242014231
012342431430031214031242304210
012342431430130024231042314021
012342431430130204032141324102
012343142024143042312100314302
012343142042431021133200423401
012343241024143042320211304310
012343241024143240301020314321
012343241024413012321200341304
012343241420103042321204310341
012343421042143042321400312301
012343421420130204031343214201
012344123302041241100434332201
012344132032134024134200321104
012344132032431021130244320041
012344132032431120034210324410 ;taocp A
012344132302401021134203024431
012344213023413102024231004431
012344213320014142024231330041
012344213320410102023241443301
012344213320410120204031334421
012344312032431021314200124403
012344312320410102034321234401
</pre>
Winningやtaocpについてはまた述べる.
</body>
</html>
和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-14822902350048159612020-10-18T11:17:00.003+09:002020-12-10T14:54:25.594+09:00プログラミングコンテスト一週間くらい前にプログラミングコンテストがあったらしい. ツィッター
のキーワードHHKBを見たら, コンテストの話題が沢山あった.<br><br>
私はプログラムを急いで書くコンテストは嫌いだ. プログラムはゆっくり
と考えながら書くものとDijkstraは主張する. 私もまったく同感で, ずい
ぶん前に大学対抗プログラミングコンテスト(ICPC)のオフィシャルを2年
ほど引き受けたが, プログラムをあわてて書き, テストの入力データに対
して正しい数値だけ得られれば完成で, 完成までの時間の短いのを上位
とする風潮に賛成できず, オフィシャルを辞退した.<br><br>
しかし世の中ではどんな問題でコンテストをやっているかとつい問題を覗いたのが
運の尽き. 解いてみたくなり, プログラムに取り憑かれた.<br><br>
やったのは<a href="https://atcoder.jp/contests/hhkb2020/tasks/hhkb2020_d">square</a>で, 二次元空間の整数値の格子点に角を持ち, x, y軸に
平行な辺を持つ正方形の入れ子問題である. すなわち, 一辺の長さnの正方形N
の内部に, 一辺の長さがそれぞれaとbの正方形AとBを重さならないように置く
置き方が何通りあるかを計算する.
<br><br>
取り敢えず例題にあった n=4, a=2, b=2で考える. aもbも2でややこしい
が, まぁよかろう. <br><br>
n×nの正方形Nの中に, 2×2の正方形Bの置き方は, 下の図のよう
に3×3で9通りある. Bの左の線はNの左端から, Nの右端からBの幅bを
引いたところまで置けて, 両端を含めるから n-b+1個. 今の数値では
4-2+1=3で, これが横方向も縦方向もあるので9通りになる.<br><br>
<img src="http://www.iijlab.net/~ew/images/square/square0.png" width=400><br>
縦, 横をそれぞれl<sub>0</sub>, l<sub>1</sub>とする長方形内のBの置き方は,<br>
bs(l<sub>0</sub>,l<sub>1</sub>)=(l<sub>0</sub>-b+1)×(l<sub>1</sub>-b+1)<br>
だから, 全体に置くのはbs(n,n), 今はbs(4,4)=9である. 今後この値をtで表す. <br><br>
2×2の正方形Aの置き方も, 下の図の実線の枠が示すようにで9通りである.
その各の位置で, Bの一部でも隠す勢力範囲は, 灰色の部分で, この勢力範囲内の
Bの数をt, つまり9から引くと, そのAの場合の数が得られる. これをすべての
Aの場所で計算すれば答が得られる. <br><br>
<img src="http://www.iijlab.net/~ew/images/square/square1.png" width="400"><br>
ところで, 横方向だけ考えて, 勢力範囲の左右の幅は, Aの左端の, Nの左端からの距離をxとすると, 次の下左図の太い縦線のようになる. この下限をy<sub>0</sub>, 上限をy<sub>1</sub>とすると,
赤線の範囲になる. 右は一般の場合の図だ.<br><br>
<img src="http://www.iijlab.net/~ew/images/square/square2.png" width="400"><br>
この図から分るように<br><br>
<pre>
y<sub>0</sub>(x) = max (x-b+1, 0),
y<sub>1</sub>(x) = min (x+a-(b+1), n)
</pre>
だから, Aの左下の座標を(x<sub>0</sub>, x<sub>1</sub>)とすれば, 勢力範囲は
横がy<sub>1</sub>(x<sub>0</sub>)-y<sub>0</sub>(x<sub>0</sub>),
縦がy<sub>1</sub>(x<sub>1</sub>)-y<sub>0</sub>(x<sub>1</sub>)で, この中に入るBの数をtから引くことになる.
つまり t-s(y<sub>1</sub>(x<sub>0</sub>)-y<sub>0</sub>(x<sub>0</sub>),y<sub>1</sub>(x<sub>1</sub>)-y<sub>0</sub>(x<sub>1</sub>)).<br><br>
私は通常, Lispの一方言のSchemeでプログラムを書くので
<pre>
(define (square n a b)
(define (y0 x) (max (- x (- b 1)) 0))
(define (y1 x) (min (+ x a (- b 1)) n))
(define (bs l0 l1) (* (- l0 b -1) (- l1 b -1)))
(let ((t (bs n n)))
(do ((x0 0 (+ x0 1)) (sum 0)) ((> x0 (- n a)) sum)
(do ((x1 0 (+ x1 1))) ((> x1 (- n a)))
(set! sum (+ sum
(- t (bs (- (y1 x0) (y0 x0)) (- (y1 x1) (y0 x1)))
))))))
</pre>
実行すると<br>
<pre>
(square 4 2 2) => 32
(square 3 1 2) => 20
</pre>
となる. Schemeが苦手な人のためにPythonのプログラムを示すと
<pre>
def square (n, a, b):
def y0 (x):
return max(x-(b-1),0)
def y1 (x):
return min(x+a+(b-1),n)
def bs (l0,l1):
return (l0-b+1)*(l1-b+1)
t=bs(n,n)
sum=0
for x0 in range(0,n-a+1):
for x1 in range(0,n-a+1):
sum = sum + \
t - bs (y1(x0)-y0(x0), y1(x1)-y0(x1))
return sum
print (square(3,1,2))
print (square(3,2,1))
print (square(4,2,2))
</pre>
実行すると
<pre>
20
20
32
</pre>
</pre><br>
プログラムは一旦うまく動くようになると, 改良の糸口が見えだすものだ.
仮に勢力範囲の図で, 灰色の部分の面積を足すとすると, 灰色の横の長さを足し合わせ,
縦の長さも足し合わせて, 和どうしを掛ければよいようだ.
いちいちtから引く代わりに, Aの置き方の数にtを掛けたものから, 勢力範囲の
総面積を引けばよい.<br><br>
勢力範囲の幅は, 図のy<sub>0</sub>, y<sub>1</sub>の間の縦の太線を足せばよい. 両端は
少し短いが, 破線の部分も入れてしまえば, 標準の縦の長さ a+2b-2に太線の本数 n-a+1を
掛ければよい. 破線の部分は0からb-1までの和で, 左下と右上の2箇所にあるから,
(b-1)b である.<br><br>
そうそう, 勢力範囲lに入るBの数がl-b+1だったように, 太線の長さからもb+1を
引く必要がある. (プログラムの(* s (- b 1)) のところ.)<br><br>
下のプログラムで, tは前述の通り, sは太線の数, pが太線の総延長である. <br><br>
<pre>
(define (square n a b)
(let* ((t (* (- n b -1) (- n b -1)))
(s (- n a -1))
(p (- (* (+ a b b -2) s) (* b (- b 1)) (* s (- b 1)))))
(modulo (- (* t s s) (* p p)) 1000000007)))
(square 331895368 154715807 13941326) => 409369707
</pre>
Python版は
<pre>
def square(n, a, b):
t = (n-b+1) ** 2
s = n - a + 1
p = (a+2*b-2)*s - b*(b-1) - s*(b-1)
return (t*s*s - p*p) % 1000000007
print (square(331895368,154715807,13941326))
</pre>
実行結果は
<pre>
409369707
</pre>
まずまず楽しんだ.
和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-43361044790426502032020-10-12T10:48:00.003+09:002020-12-10T14:56:10.285+09:00連分数と近似分数前回のブログ, 連分数と近似分数でeの近似分数 878/323をStern Brocot木から求めたが, 連分数との関係については書かなかった.<br><br>
容易に想像できるように, LとRの列RL0RLR2LRL4RLR6L ...は, 1年前のブログにあるeの連分数展開[2,1,2,1,1,4,1,1,6,1,1,8]と同じパターンである.<br><br>
昔のブログの最後の方の中間近似分数にも同じ分数が現れる.<br><br>
「コンピュータの数学」には無理数αからLとRの列を求める式もある.<br><br>
<pre>
if α < 1 then (output(L); α:=&alpha/(1-&alpha))
else (output(R); α:=α-1)
</pre>
途中のαと一緒にLとRを計算すると<br><br>
<pre>
from math import e
print(e)
alpha=e
lri=[]
def sb():
global alpha
if alpha<1:
lri.append("L")
alpha=alpha/(1-alpha)
else:
lri.append("R")
alpha=alpha-1
print(alpha)
for i in range(16):
sb()
print(lri)
2.718281828459045
1.718281828459045
0.7182818284590451
2.5496467783038432
1.5496467783038432
0.5496467783038432
1.2204792856454296
0.22047928564542962
0.2828395468977148
0.39438809777395056
0.651222501282252
1.8671574389873726
0.8671574389873726
6.5277079301786785
5.5277079301786785
4.5277079301786785
3.5277079301786785
['R','R','L','R','R','L','R','L','L','L','L','R','L',
'R','R','R']
</pre>
つまり連分数で1なら文字は1回, 2なら2回 繰り返すのであり, Stern Brocot木による計算は連分数展開の近似を, 中間近似分数まで含めて得ていたのであった.<br><br>
和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-91125674963306125442020-10-11T16:32:00.003+09:002020-10-11T17:05:01.935+09:00連分数と近似分数<p> 昨年6月のブログで, eの近似分数 878/323の求め方を話題にした. 先月出版された「コンピュータの数学 第2版」を拾い読みしていたら, Stern Brocot木の関連でeの近似分数を求める話があった.<br /><br />Stern Brocot木は, 互いに素なmとnの分数を大きさの順に並べたもので, 次のような形である. ある場所の値は, その左上の一番近い祖先m, nと右上の一番近い祖先m', n'から m+m'/n+n' で得られる.<br /></p><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-OomvFUBDi6k/X4KzYRyjH2I/AAAAAAAAAaw/-JWqCMiqX-0ChvT2Rs8b5D38tSvRe8KGACLcBGAsYHQ/s1018/sternbrocot.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="872" data-original-width="1018" height="274" src="https://1.bp.blogspot.com/-OomvFUBDi6k/X4KzYRyjH2I/AAAAAAAAAaw/-JWqCMiqX-0ChvT2Rs8b5D38tSvRe8KGACLcBGAsYHQ/w364-h274/sternbrocot.png" width="364" /></a></div><p><br />最上部の1/1からLでは左下へ, Rで右下へ辿ると, 色々な分数が一意的に表せる. そしてeは<br /><br />e=RL<sup>0</sup>RLR<sup>2</sup>LRL<sup>4</sup>RLR<sup>6</sup>L ...<br /><br />つまりRL<sup>0</sup>R, LR<sup>2</sup>L, RL<sup>4</sup>R, LR<sup>6</sup>L, ...<br /><br />であって, Eulerが24歳の時に発見したと書いてあった.<br /><br />早速計算するプログラムをSchemeで書いた.<br /></p><p>(l)は左下へ, (r)は右下へ進む関数. lmnは左祖先, rmnは右祖先, cmnは現在の分数, nmnは次の世代の分数である. <br /></p><p>(define (show) (display cmn))<br />(define lmn '(0 1)) (define cmn '(1 1)) (define rmn '(1 0))<br />(define (reset!) (set! lmn '(0 1)) (set! cmn '(1 1))<br /> (set! rmn '(1 0)) (show))<br />(define (l) (let ((nmn (map + lmn cmn)))<br /> (set! rmn cmn) (set! cmn nmn) (show)))<br />(define (r) (let ((nmn (map + cmn rmn)))<br /> (set! lmn cmn) (set! cmn nmn) (show)))<br />(for-each (lambda (c)<br /> (cond ((eq? c 'l) (l))<br /> ((eq? c 'r) (r))<br /> (else (reset!))))<br /> '(c r r l r r l r l l l l r l r r r))<br /><br />実行すると, </p><p></p><p>(1 1)(2 1)(3 1)(5 2)(8 3)(11 4)(19 7)(30 11)(49 18)(68 25)</p><p>(87 32)(106 39)(193 71)(299 110)(492 181)(685 252)</p><p>(878 323)</p><p><br />すげー.<br /><br /></p><br />和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-47711803820975166272020-10-01T17:21:00.002+09:002020-12-10T14:59:31.607+09:00XYプロッター<p> IEEE Specturm 09.20を眺めていたら, AxidrawというXYプロッターの記事があった. それに下のような図があったので, その機構を考えてみた.<br /><a href="https://1.bp.blogspot.com/-Xw751X4p0do/X3WRJFO13RI/AAAAAAAAAac/HIX3cjTpi8A3CHKNCKsRTn8ENsWmmaR_QCLcBGAsYHQ/s1078/xyplotter.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="772" data-original-width="1078" src="https://1.bp.blogspot.com/-Xw751X4p0do/X3WRJFO13RI/AAAAAAAAAac/HIX3cjTpi8A3CHKNCKsRTn8ENsWmmaR_QCLcBGAsYHQ/s320/xyplotter.png" width="320" /></a><br />中央に左右に広がる基盤があり, これは固定されている. その左端と右端にそれぞれモーターと同軸で回転するプーリーがある(L, R). 縦に長い板は, 縦方向を保ったまま, 上下左右に移動し, 下に自由に回転するプーリーMがあって, その辺りに上げ下げできる筆記具がある. それでXYプロッターになるわけだ.<br /><br />3個のプーリーは, 中央の4個の補助プーリーを介して1本のベルトで結ばれており, ベルトの最後は, 縦の板の上端である. <br /><br />まずLを固定し, Rの周囲が長さaだけ右回転したとする. その際, 縦の板が上下だけに動いたとすると, 上端に左側のベルトもaだけ弛むが, Lが回転しない約束なので, その弛みは上下の板を右へ移動して吸収することになる.<br /><br />というわけで, Rがaだけ回転すると, 上下の板は右へa/2, 下へa/2だけ移動する. そうすると, 右上部分のベルトは, x1もy1もa/2だけ短くなり, Rはaだけ回転し, 右下部分ではx1はa/2短く, y0はa/2長くなるので, Rから出たaはMを回って左下へ行く. <br /><br />左上はy1がa/2短く, x0がa/2長くなるから, 予定通りLは回転しない. 左下はx0もy0も伸びるが, その分だけMから貰うので万事うまくいく.<br /><br />つまりRの右回転は上下板を右下へ移動させ, 左回転は左上へ移動させる. 同様にLの左回転は上下板を左下へ移動させ, 右回転は右上は移動させる.<br /><br />左右だけの移動, 上下だけの移動は, LとRの回転を合成すればよい.<br /><br />なるほどうまい仕掛けであった. <br /></p><br /><br />和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-68680428219986121892020-08-29T16:41:00.005+09:002021-10-01T14:05:40.675+09:00恒星時計<p> 恒星時計(こうせいじ-けい)を実装した. <a href="http://www.iijlab.net/~ew/siderealclock12.html" target="_blank">ここ</a>からダウンロードできる.</p><p> ダウンロードした時から走行し, 下の上の図のような形である. 右上のディジタル表示は現在の太陽時の年月日と時刻(時, 分, 秒). 左上の表示は現在の(地方)恒星時の時, 分, 秒である.<br /><br />中央の時計がその恒星時をアナログ的に表示する. 恒星時は0時から24時まであるが, 時計の文字盤には時の表示が12個しかなく, この時計の時針は恒星時の12時間で1回転する. その代わり, 恒星時の正時に文字盤の時針と反対の位置の文字が12増えたり12減ったりする. <br /><br />この図では, 時針は9時と10時の間にあり, そこから左周りに6時間は過去の時刻を, 右周りに6時間は未来の6時間を示す.<br /><br />時計の周囲の四角形の内部をクリックすると, 時計は太陽時の秒が繰り上がった時点で停止し, その時刻での恒星時がミリ秒まで表示される. 下の時計の図は停止中の様子を示す. もう一度クリックすると走行を再開する. <br /></p><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-78tU44xvGVU/X0n7tjbuJPI/AAAAAAAAAYE/86j6fXAfnkcyssrNKPbSKXyIZ4CKZbREgCLcBGAsYHQ/s1024/siderealtimeclock0.png" style="margin-left: 1em; margin-right: 1em;"></a><br /></div><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-IWMRp8UYnrw/X0n8wOAQX9I/AAAAAAAAAYc/tZCykwfnuqsSMepkZSEyZmc3cnYMIPG-QCLcBGAsYHQ/s1024/siderealtimeclock0.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="1024" data-original-width="1020" height="328" src="https://1.bp.blogspot.com/-IWMRp8UYnrw/X0n8wOAQX9I/AAAAAAAAAYc/tZCykwfnuqsSMepkZSEyZmc3cnYMIPG-QCLcBGAsYHQ/w326-h328/siderealtimeclock0.png" width="326" /></a></div><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-FpSgGAct5OM/X0n8eJn1rRI/AAAAAAAAAYU/XXHHswkSJSwh_6_HDoVCuU7fnq1B6zZ2ACLcBGAsYHQ/s1030/siderealtimeclock1.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="1030" data-original-width="1024" height="328" src="https://1.bp.blogspot.com/-FpSgGAct5OM/X0n8eJn1rRI/AAAAAAAAAYU/XXHHswkSJSwh_6_HDoVCuU7fnq1B6zZ2ACLcBGAsYHQ/w326-h328/siderealtimeclock1.png" width="326" /></a></div><p></p><p>時計の盤面に, 北極星を中心とする北天の星図が表示してある. この星図は星座早見盤の自作<a href="http://www.jodrellbank.net/wp-content/uploads/2014/10/Make-a-planisphere.pdf" target="_blank">キット</a>から使わせて頂いた. この星図は恒星時の24時間で左周りに1回転する. 通常の時計の6時相当する真下を時角0h, 3時に相当する右を時角6h, 真上を12h, 左を18hとして北天が日周運動する. <br /><br />星図に薄く天の川が表示されている. 天の川が天の北極に近いのが赤経0hの方向で, 中央からその方向に辿ると天の川の中にカシオペアのWがあり, さらに外周へ探すとペガススの四角がある. その少し外の緑の点が春分点(Υ)で, 春分点の時角と恒星時が一致しているのが分る.<br /><br />この時計の星図の下2/3くらいの横長の楕円が, その時刻に見える星空になる. <br /><br />次の大きな円は星図の元の図である. 右の方の3本の左向き矢印の先の星が, 上からうお(Pisces)のλとω, くじら座(Cetus)のιである. その下の図がうお座とくじら座で, これから見ると春分点はこれら3つの星の中央にあり, 恒星時計の文字盤の春分点をその位置にマークした. <br /><br /><br /></p><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-oD7HB7OviFM/X0n-y0ea6HI/AAAAAAAAAYo/NFdqisUz188fXF_Rrfi-dG2ykwgFlx3gwCLcBGAsYHQ/s1498/northstar.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="1498" data-original-width="1492" height="328" src="https://1.bp.blogspot.com/-oD7HB7OviFM/X0n-y0ea6HI/AAAAAAAAAYo/NFdqisUz188fXF_Rrfi-dG2ykwgFlx3gwCLcBGAsYHQ/w326-h328/northstar.png" width="326" /></a></div><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-Nqz3n7qnHww/X0n_xzP50cI/AAAAAAAAAYw/6HyGDNg1DakXn_yeA91JqC9_fYgNo326QCLcBGAsYHQ/s1990/pisces.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="1584" data-original-width="1990" height="326" src="https://1.bp.blogspot.com/-Nqz3n7qnHww/X0n_xzP50cI/AAAAAAAAAYw/6HyGDNg1DakXn_yeA91JqC9_fYgNo326QCLcBGAsYHQ/w410-h326/pisces.png" width="410" /></a></div><br /><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-RaW28Jeuy7c/X0n_9dbUU7I/AAAAAAAAAY0/fJ9zw_jwxq08FdAV-2J4mVMRDY_VhZP1QCLcBGAsYHQ/s1990/cetus.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="1554" data-original-width="1990" height="320" src="https://1.bp.blogspot.com/-RaW28Jeuy7c/X0n_9dbUU7I/AAAAAAAAAY0/fJ9zw_jwxq08FdAV-2J4mVMRDY_VhZP1QCLcBGAsYHQ/w410-h320/cetus.png" width="410" /></a></div><p>この時計では, ある太陽時の日時, 時刻の恒星時を知ることが出来る. 時計の枠内でキーボードから数字を入力すると時計が停止し, その数字を空白で区切りながら, 年, 月, 日, 時, 分, 秒と入力すると, 対応する恒星時が示される.<br /><br />最初の時計の図の時刻では, 2020⨆8⨆8⨆11⨆53⨆9と入力すると, 左上に恒星時が現れ, 文字盤にはその時刻に現れる星空が表示される. ミリ秒の値は少し違うが.<br /><br />恒星時の計算は暗算としてはこのようにする. 3月21日頃の春分には太陽が春分点にいるから, 昼頃, 太陽の南中時が恒星時の0h. それに太陽時の正午からの差を加える. 太陽の赤経は, 1ヶ月に2h増えるから, 例えば8月8日, 立秋の頃なら, 正午に太陽が南中するとき, 春分点は4ヶ月半近く, 9hだけ時角が離れる. つまり正午が9hである. 従って8月8日, 21時だと9h+9h=18hとおよその見当がつく.<br /><br />しかしこの方法は時計には使えない. 理科年表や天文年鑑には, その年の毎日の世界時0時のグリニジ視恒星時の表があり, インターネットで探すと多くの関連ページが見付かる.<br /><br />国立天文台の グリニジ恒星時計算<a href="https://eco.mtk.nao.ac.jp/cgi-bin/koyomi/cande/gst.cgi" target="_blank">HP</a>で2020年1月1日から世界時0時の恒星時を1ヶ月表示したのの一部が次の図である. (出力時刻もGMT)</p><div class="separator" style="clear: both; text-align: center;"><a href="https://1.bp.blogspot.com/-OM0rHkvf6P8/X0oBC3tDkUI/AAAAAAAAAZE/es9qvQCpsiI_e6FvXgsrEdjn5Zgsjd10gCLcBGAsYHQ/s1596/gsttable.png" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="1596" data-original-width="1394" height="512" src="https://1.bp.blogspot.com/-OM0rHkvf6P8/X0oBC3tDkUI/AAAAAAAAAZE/es9qvQCpsiI_e6FvXgsrEdjn5Zgsjd10gCLcBGAsYHQ/w447-h512/gsttable.png" width="447" /></a></div><br /><p>理科年表では0.1秒, 天文年鑑では0.1分まで表示されているのが, このHPでは0.001秒の精度で表示される.<br /><br />インターネットにあるいくつかの恒星時計算アルゴリズムで, 2020年1月1日世界時0時の恒星時を計算して見ると, どういうわけかこの表の値にならないので, この値が得られるような数値を自分で計算した.<br /><br />あるユリウス日が<i>jd</i><span style="font-size: xx-small;"> 0</span>の日時での恒星時<i>st </i><span style="font-size: xx-small;">0</span>が既知の時, ユリウス日が<i>jd</i>での日時での恒星時<i>st</i>は<br /> </p><p><i>st</i>=<i>st</i><span style="font-size: xx-small;"> 0</span>+1.0027378<span style="font-size: small;">×</span>(<i>jd</i>-<i>jd</i><span style="font-size: xx-small;"> 0</span>)<br /><br /><br />だから, <i>jd</i>と<i>st</i>のペアがあれば, ある<i>jd</i><span style="font-size: xx-small;"> 0</span>での<i>st</i><span style="font-size: xx-small;"> 0</span>が計算できる.<br /> </p><p><i>jd</i><span style="font-size: xx-small;"> 0</span>としては, 何かの恒星時計算アルゴリズムにあった切断(?truncated)ユリウス日の起算時点(<i>jd</i>=2440000.5)を使うことにした. つまり1月中の0時の<i>jd</i>と対応する<i>st</i>から31個の<br /> </p><p><i>st</i><span style="font-size: xx-small;"> 0</span>=<i>st</i>-1.0027378×(<i>jd</i>-2440000.5)<br /><br />を求めて平均し, <i>st</i><span style="font-size: xx-small;"> 0</span>を作る. 恒星時の表はミリ秒のリスト<br />(226 778 330 882 435 990 548 109 672 237 802 365 925 481 34 <br />585 136 688 243 801 362 925 487 49 609 167 721 274 824 374 924))<br /><br />を用意し, ミリ秒が減った時は前日の秒に237秒を, 増えた時は236秒を足し, ミリ秒その日の値としながら, 1月の恒星時を入力した. そして定数<i>st</i><span style="font-size: xx-small;">0</span>を得たが, それで12月の恒星時を計算すると, 天文台の表と多少ずれるので, 12月についても同様に計算し, 両者の定数の平均を使うことにした. つまり<br /> </p><p><i>st</i>=0.6712386796084648+(1.00273791×(<i>jd</i>-2440000.5))<br /><br />を使う. <br /><br />この値を使う, Pythonで書いた恒星時計算プログラムは次の通り. <br /><br />def jd(y,m,d):<br /> def leap(y):<br /> return y%400==0 if y%100==0 else y%4==0<br /> return y*365+y//4-y//100+y//400-(489*(13-m)+26)//16\<br /> +((1 if leap(y)else 2)if m<3 else 0)+d+1721426<br />def sidereal(y,mo,d,h,mi,s):<br /> def frac(a):<br /> return a-int(a)<br /> ut=jd(y,mo,d)+(h-9+(mi+s/60)/60)/24-0.5 #-9: to GMT <br /> st=0.6712386796084648+1.00273791*(ut-2440000.5)+9/24</p><p> #+9/24: to JST <br /> h=frac(st)*24<br /> m=frac(h)*60<br /> s=frac(m)*60<br /> ms=frac(s)<br /> return [int(h),int(m),int(s),ms]<br />print(sidereal(2020,1,1,9,0,0))<br /><br /><br />出力は[15, 40, 28, 0.20889443391934037]で, 前掲の1月1日の値 6h40m28.226sに9時間足したのにかなり合っている. <br /><br />このような恒星時計が, 宮沢賢治の「銀河鉄道の夜」に出て来る 白鳥, わし, さそり, ケンタウルスなどの駅のホームや操車場やアルビレオ観測所で動いていると想像すると楽しい.<br /><br />「ケンタウルス, 露をふらせ.」<br /><br /></p>和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-90304290649104627442020-08-26T14:01:00.003+09:002020-12-10T15:02:27.827+09:00老子<p> 夏目漱石が護国寺で運慶を見たとか(<a href="https://www.aozora.gr.jp/cards/000148/files/799_14972.html" target="_blank">「夢十夜」の「第六夜」</a>), 芥川龍之介が市電から寒山, 拾得を見とか(<a href="https://www.aozora.gr.jp/cards/000879/files/3809_26610.html" target="_blank">「寒山拾得」</a>)いう話があるが, 寺田寅彦が<a href="https://www.aozora.gr.jp/cards/000042/files/4359_9286.html" target="_blank">「電車で老子に会った話」</a>は趣きが異なる.<br /><br />つまり寺田先生は「インゼル・ビュフェライ叢書」にアレクサンダー・ウラールが訳した<a href="https://www.amazon.de/Bahn-rechte-Lao-Tse-Insel-B%C3%BCcherei/dp/3458089918" target="_blank">『老子』</a>を見付けて購入し, 電車の中で読み始めたら面白くなり, 数回の車中で読了したという. 日本語の解説より分りやすかったらしい.<br /><br />私は老子は近寄り難く避けていたが, 昨年末, 中国深圳で少し話をしたら, 謝礼として為書きのある書を貰った. それに「大道至簡」と書いてあり, 調べたところ老子の言葉らしいことが判明. それではというので老子を読むことにした.<br /><a href="https://1.bp.blogspot.com/-VhpkfXCSbxk/X0Xq0wrWcyI/AAAAAAAAAXw/S07QQ7SP9x0t547Qu92wq9eu_055k3YWQCLcBGAsYHQ/s2016/C10.jpg" style="margin-left: 1em; margin-right: 1em;"><img border="0" data-original-height="1344" data-original-width="2016" height="342" src="https://1.bp.blogspot.com/-VhpkfXCSbxk/X0Xq0wrWcyI/AAAAAAAAAXw/S07QQ7SP9x0t547Qu92wq9eu_055k3YWQCLcBGAsYHQ/w512-h342/C10.jpg" width="512" /></a><br />早速読んだのは<br /><br />池田知久 「老子 全訳注」 講談社学術文庫<br /><br />だが,<br /><br />保立道久 「現代語訳」 老子 ちくま新書<br /><br />もインターネットに<br /><br />保立道久の研究雑記 <br /><a href="https://note.com/michihisahotate/n/n4f5c6abee089" target="_blank">『老子』、現代語訳・原文・読み下し</a><br /><br />があったので, それも参考にする. <br /><br /><br />老子は第1章から第37章までの道経と, 第38章から第81章までの徳経の, 計81章で構成される. 今回は夫々の最初の章の最初の文を例にする.<br /><br />池田本では, 第1章の最初は<br /><br />道可道也、非恆道也。名可名也、非恆名也。<br /><br />道の道とす可きは、恒の道に非ざるなり。名の名とす可きは、恒の名に非ざるなり。<br /><br />道というものがあれこれ唱えられている中で、道として定立することのできるわたしの道も、恒常不変の真の道であるとは言えない。道についての名(理論)があれこれ説かれている中で、道の名(理論)として表現することのできるわたしの名(理論)も、恒常不変の真の名(理論)であるとは言えない。<br /><br />一方, 保立本は<br /><br />道可道也、非恒道也。名可名也、非恒名也。<br /><br />道の道(ゆ)くべきは、恒なる道に非(あら)ざるなり。名の名づくべきは、恒なる名に非ざるなり。<br /><br />普通に行く道と、ここでいう「恒なる道」はまったく違うものだ。普通に名づけることができる名と、ここでいう「恒なる名」もまったく違う。<br /><br />どちらも, 特に後者は短いだけあって, 禅問答のようだ. <br /><br />次に第38章は, 池田本で<br /><br />上徳不徳、是以有徳。下徳不失徳、是以徳无。<br /><br />上徳は徳ならず、是を以て徳有り。下徳は徳を失わず、是を以て徳无し。<br /><br />そもそも最上の徳は、世間的な徳とは正反対である。だからこそ、真の徳がある。下等な徳は、世間的な徳を捨てることが出来ない。だからこそ、真の徳がないのだ。<br /><br />保立本は<br /><br />上徳、不徳是以有徳。下徳、不失徳是以無徳。<br /><br />上徳は、徳ならずして是(ここ)を以て徳あり。下徳は、徳を失わずして是を以て徳なし。<br /><br />「道」から発した最上の「徳(いきおい)」(上徳)は徳(はたらき)がないようにみえて大きな徳(いきおい)があり、そうでない「下徳」は徳(はたらき)を失っていないようだが実は徳(いきおい)がなくなっている。<br /><br />同じような感じだ. もともと道とか徳とかが抽象的だから, 説明が難しくて当然であろう.<br /><br />寺田先生がドイツ語訳で分ったような気分になったというので, 英訳を探すと<br /><br /><a href="https://terebess.hu/english/tao/gia.html#Kap01" target="_blank">The Tao Te Ching by Lao Tzu</a><br />馮家福 Jane English<br /><br />があった. 表題のTao Te ChingはTao=道, Te=徳, Ching=経のことだ.<br /><br />道の説明は次のようだ. <br /><br />One<br /><br />The Tao that can be told is not the eternal Tao.<br />The name that can be named is not the eternal name.<br /><br />また徳の方は<br /><br />Thirty-eight<br /><br />A truly good man is not aware of his goodness,<br />And is therefore good.<br />A foolish man tries to be good,<br />And is therefore not good.<br /><br />上徳と下徳がgood manとfoolish manになっているのが凄い. なるほど日本語の説明より具体的である.<br /><br />インターネットでさらに探すと, 寺田先生が読んだAlexander Ularの訳の一部があった. <br /><br />Die Bahn und der rechte Weg des Lao-Tse<br />Alexander Ular<br /><br />DER ERSTE SPRUCH<br /><br />Die Bahn der Bahnen ist nicht die Alltagsbahn;<br />Der Name der Namen ist nicht der Alltagsname.<br /><br />38章は掲載されていなかった.<br /><br />ドイツ語を和訳するのは簡単だ. 複数をたちと訳せば<br /><br />道たちの道は毎日の道ではない<br />名たちの名は毎日の名ではない<br /><br />でもこれが何を言おうとしているかは, やはり分らない.<br /><br />結果的には池田本の訳が一番しっくりする. つまりドイツ語が分り易いというのも章によるのではないかということである. 沢山の訳を比較すれば, 段々と理解しやすくなるという当然の結論に達する. <br /><br /><br />ところで論語だが, 大学の先生方がよくいうのが<br /><br />学而不思則罔 思而不学則殆<br /><br />学びて思わざれば則ち罔(くら)し、思いて学ばざれば則ち殆(あやう)し。<br /><br />これは罔や殆が理解を妨げている. 私の持っている論語の英訳を見ると<br /><br />The Analects of Confucius<br /><br />The Master said: `Learing without thinking is useless. Thinking without learning is dangerous.'<br /><br />のように, uselessとdangerousになっている.<br /><br />つまり, 数学の公式を学んでも, それを十分に理解しなければ, 使えない. また思い付いた名案も, その基礎が分っていないと, 迷案になるというようなことだ.<br /><br />この場合は英訳が分り易い. <br /><br /><br /></p><br /><br />和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-66029819357747648132019-06-17T21:19:00.000+09:002019-06-24T19:27:59.088+09:00SDGsロゴの錯視貧困なし, 飢餓ゼロ, 健康と幸福, 高度の教育,...など, 2030年までの17達成目標を掲げる国連のSDGs(持続可能な開発目標 Sustainable Development Goals)のロゴはなかなか華麗だ. 早速PostScriptで描いてみた.<br />
<br />
<br />
<a href="https://www.iijlab.net/~ew/images/sdgslogo1.png" imageanchor="1" style="clear: left; float: left; margin-bottom: 1em; margin-right: 1em;"><img border="0" src="https://www.iijlab.net/~ew/images/sdgslogo1.png" width="400" /></a><br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
<br />
カラフルな17個の扇形で囲まれていて, その隣り同志の扇形の間には白い間隔がある. この間隔は等幅のようでもあるが, なんとなく内側の方が広く見えるから不思議だ. 扇形が内側で狭くなるのにつられて, 間隔は逆に内側が広く感じるのであろうか. これを取り敢えずSDGsロゴの錯視といおう. <br />
<br />
そこだけ取り出して描いたのがこの図だ.<br />
<br />
<br />
<img src="https://www.iijlab.net/~ew/images/sdgslogo2.png" width="200" /><br />
<br />
<br />
これだけ眺めると, 外縁での間隔<i>g<sub>o</sub></i>と, 内縁での間隔<i>g<sub>i</sub></i>は別に違っては見えないが, 上の図ではたしかに錯視が起きる.<br />
<br />
北岡さんの<a href="http://www.psy.ritsumei.ac.jp/~akitaoka/catalog.html">錯視のカタログ<t>http://www.psy.ritsumei.ac.jp/~akitaoka/catalog.html</t></a>を見たが, こんなに簡単なのは掲載されていないようだ. <br />
<br />
間隔の幅Gを変えたり, 内縁半径Rを変えたりして, 錯視の起きる範囲をしらべたいと思い作った<a href="https://www.iijlab.net/~ew/sdgslogo.html">ウェブアプリ</a>が下だ. 左下は+, -でRを加減する. Rをクリックすると標準に戻る. 右下はGについての同様なメニューである.<br />
<br />
<img src="https://www.iijlab.net/~ew/images/sdgslogo.png" width="400" /><br />
<br />
錯視の起きる起きないの境界は微妙なので, はっきりしたことはいえないが, 錯視はほどほどの間隔の時に起き易いらしい. このウェブアプリは出来たばかりで枯れていないので, うまく動作しない時はご容赦を乞う. <br />
<br />
<br />和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-82123033512643795672019-06-15T13:26:00.000+09:002019-06-15T13:29:09.034+09:00Tパズル昨年夏のプログラミング・シンポジウムは諏訪のかんぽの宿で開催された. その客室にはお茶などと一緒にTパズルが置いてあり, それをいじった泊り客が気に入って宿の売店で購入するのを期待している. <br />
<br />
<img src="https://www.iijlab.net/~ew/images/Tpuzzle/houseT.png" width="400" /><br />
<br />
<br />
Tパズルは上の図の左の家型の五角形を4分割したなんの変哲もないピースを組み合わせ, 指示されたいろいろな形を構成するパズルである. 最初の指示が図の右のT型を作るので, Tパズルというらしい. 元の家型の時のピースの縦横の辺が, T型では斜めになっているので, この問題は最初からいささか意地が悪い. <br />
<br />
インターネットで探すと, 下の図のように様々な問題があるので, ちょっとやってみたくなり, 百円ショップで小学生の使う工作用紙を購入し, 自作した. <br />
<br />
<img src="https://www.iijlab.net/~ew/images/Tpuzzle/Tpuzprob.png" width="400" /><br />
<br />
自作するには下のような図を描き, 同じ形の表と裏を貼り合せた. この升目が工作用紙の1cmである. 各辺の内側の数値は, T型の足の幅を基本単位とするその辺の長さで, 出来上がったピースにも同様に記入してある. この図では各ピースが繋って描いてあるが, 作業に便利なように, 本当はピースの間に隙間があった. <br />
<br />
<img src="https://www.iijlab.net/~ew/images/Tpuzzle/pattern.png" width="400" /><br />
<br />
実際に出来たものは下の写真の通りで, 旅館で販売しているものより遥かに大きい. 市販のものは辺に寸法が書いてなく, いわばアナログ仕様だが, それでは√2と1.5の区別もつき難そうだ. 私のは寸法が記入してあるので間違うことはない. <br />
<br />
<img src="https://www.iijlab.net/~ew/images/Tpuzzle/Tpuzpiece.png" width="300" /><br />
<br />
<br />
<br />
さて, 私のやり方をちょっと種明かししよう. 上にあった問題の1番のT型と12番の家型を例にする. <br />
<br />
型の図をPCに読み込み, 適当に拡大する. 角に順に番号をつける. マウスを動かし, それぞれの角の座標を読み取る. 座標はそれほどの精度を要求しない. それが下の2枚の図だ. <br />
<img src="https://www.iijlab.net/~ew/images/Tpuzzle/caption1.png" width="400" /><br />
<br />
<img src="https://www.iijlab.net/~ew/images/Tpuzzle/caption12.png" width="400" /><br />
<br />
次にこの形を角を頂点とする三角形に分割する. T型では((0 6 7) (0 1 6) (4 5 6) (4 6 1) (4 1 2) (4 2 3)), 家型では((1 2 3) (1 3 4) (1 4 0))とした. そして夫々の型で, 座標と三角形の頂点の組をSchemeのプログラムで読み込む. その値を確認するため, プログラムはSchemeのgraphics機能を利用して出力する. それが次の図である. T型の足など多少怪しいがまぁよしとしよう. <br />
<br />
<img src="https://www.iijlab.net/~ew/images/Tpuzzle/T1skelton.png" width="400" /><br />
<br />
<img src="https://www.iijlab.net/~ew/images/Tpuzzle/T12skelton.png" width="400" /><br />
<br />
続いてこの座標系において相隣る2点間の距離を計算する. それぞれの点のx, y座標が既知だから, Pythagorasの定理で距離は分る. 更にこの座標系での各々の三角形の面積を求め, 形全体の面積を算出する. <br />
<br />
<i>x</i><sub>0</sub>,<i>y</i><sub>0</sub>; <i>x</i><sub>1</sub>,<i>y</i><sub>1</sub>; <i>x</i><sub>2</sub>,<i>y</i><sub>2</sub>を頂点とする三角形の面積は, 行列式<br />
<br />
|<i>x</i><sub>0</sub> <i>y</i><sub>0</sub> 1|<br />
|<i>x</i><sub>1</sub> <i>y</i><sub>1</sub> 1|<br />
|<i>x</i><sub>2</sub> <i>y</i><sub>2</sub> 1|<br />
<br />
の半分であることは高校生のころ, 解析幾何の授業で学んだのでそれを使う. ただ三角形の頂点を時計廻りか反時計廻りかで辿ると面積が正負になり, いつも怪しくなる. そういう時は(0,0) (1,0) (0,1)の三角形で試みると反時計廻りで正になると分るが, プログラムでは絶対値をとってから合計している.<br />
<br />
一方, 4個のピースの面積の和は6平方基本単位であるから, 今得られた面積を6で除して平方根をとると, それが座標系の長さと基本単位の比で, それぞれの辺が基本単位に変換できる.<br />
<br />
これでT型を計算すると<br />
(2.8999986101772106 .8652158634699881 .9960165420683874 2.935646424994157 1.0222409451321806 .9698622511287074 2.8779826845214163 1.0169988220903106)<br />
が出てきて, 番号0の角からの順の辺の長さが, 3,1,1,3,1,1,3,1であることが分る. (0.865を1にするには勇気がいるが.)<br />
<br />
また家型の方は<br />
(1.418992855950106 2.022013340831029 1.993599932770006 1.4095980475184093 2.7873862791620128)<br />
だから, √2, 2,2,√2, 2√2らしいとなる. <br />
<br />
ここまで来たら大体の問題はわけもなく解けるが, 実際にやってみると案外微妙な値になる型もあり(1.5か√2かなど), ここに述べた方法が特効薬であるわけでもない. 元々の問題の図もさほど精密ではなく, また中にはこれはどうかと思う問題もあったが, 目くじらを立てることもあるまい.<br />
<br />和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com2tag:blogger.com,1999:blog-8851344890851282373.post-53549145459202221642019-06-13T13:29:00.002+09:002020-10-12T10:19:39.478+09:00連分数と近似分数<div class="separator" style="clear: both; text-align: center;">
</div>
<html><br /> <body><br /> </body></html><br />Martin GardnerがそのThe Unexpected Hanging(ネットで見つかる)の40ページで, 自然対数の底<i>e</i>を近似する, 分母分子とも3桁以内の分数は何かと問うている. その答は覚えやすい878/323 (=2.718266253869969...)だが, その計算法が今回の話題である.<br /><br />Gardnerは自分では説明せず, <br />George Chrystal, Algebra An Elementary Text-Book for the Higher Classes of Secondary Schools and for Colleges, Part II. 1906<br />を見よとしか書いていない. この古典は有名らしく, ネットで探すと読むことができる. (648ページもある!)<br /><t>http://www.archive.org/details/algebraelementar02chryuoft.pdf</t><br /><br />しかしこれは, 高木貞治先生の「初等整数論講義」(手元のは第2版だが)の「中間近似分数」の章にそのものずばりの解説がある. <br /><br />円周率πの同様な近似値355/113は遥かに有名であった. この分数を得るには連分数を使う.<br /><br />普通の連分数は<br /><br />
<img br="" src="https://www.iijlab.net/~ew/images/contfrac/cf0.png" width="400" /><br /><br />の左上のような形である. この元のTexプログラムは<br /><pre>
\[a_0+{b_1\over\displaystyle a_1+
{\strut b_2\over\displaystyle a_2+
{\strut b_3\over\displaystyle\ddots+
{\strut b_n\over a_n+\ddots
}}}}\]
\[a_0+\frac{b_1}{a_1}\:\raisebox{-1.3ex}{+}\:\frac{b_2}{a_2}
\:\raisebox{-1.3ex}{+}\:\raisebox{-1.ex}{\ldots}\,
\raisebox{-1.3ex}{+}\:\frac{b_n}{a_n}
\raisebox{-1.3ex}{+}\:\raisebox{-1.ex}{\ldots}\]</pre>
<br />この式にある<i>b<sub>i</sub></i>を1, <i>a<sub>i</sub></i>を正の整数<i>k<sub>i</sub></i>にしたものを単純(simple)連分数とか正則(regular)連分数といい, 右上のように書き, 物の本にはこれだけを扱ううものが少くない.<br /><br />こういう表示法は空間をとるので, 通常は下の左や右のように書く. <i>k<sub>i</sub></i>を部分商という.<br /><br />こういう連分数から, 878/323や355/113のような近似分数が得られる.<br /><br />円周率πの連分数展開で得られる部分商は7, 15, 1, 292, 1,...<br />で, この値を順に求める方法が「初等整数論講義」に載っている. これを見ると連分数は, Euclidの互除法に似ていることが分る.<br /><br /><img src="https://www.iijlab.net/~ew/images/contfrac/contfracdivision.png" width="400" /><br /><br />部分商<i>k</i>の値から, <br /><br /><img src="https://www.iijlab.net/~ew/images/contfrac/cf1.png" width="400" /><br />の関係で<i>p, q</i>が計算できる. πの場合の<i>k, p, q</i>は下の表の通り. <i>p/q</i>は連分数の展開を途中で止めたときの分数の値であり, 下の計算のようになる. 各行の右は, 左の分数の値とπとの差で, 見て分るように, 正負が交互に現れる.<br /><br />この表の下の方に, 355/113がある. この近似分数はこうして発見された.<br /><br /><img br="" src="https://www.iijlab.net/~ew/images/contfrac/cf2.png" width="400" /><br /><i>e</i>については, 今回はPythonで書いて<br /><br /><pre>
import math
x=math.e
ai=[]
while len(ai)<12:
a=math.floor(x)
ai.append(a)
x=1/(x-a)
print(ai)</12:></pre>
<br />すると得られた結果は<br />
<br />[2, 1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8]<br /><br /><br />πと同じような表を作ると<br /><br /><img src="https://www.iijlab.net/~ew/images/contfrac/cf3.png" width="400" /><br /><br />しかしこの計算には, 問題の878/323がない. そこで「中間近似分数」が登場する.<br /><br /> <i>p</i><sub>5</sub>=106, <i>p</i><sub>6</sub>=193,<i>p</i><sub>7</sub>=1264<br /><i>q</i><sub>5</sub>=39, <i>q</i><sub>6</sub>=71,<i>q</i><sub>7</sub>=465<br />の間には<br /><i>p</i><sub>7</sub>-<i>p</i><sub>5</sub>=<i>p</i><sub>6</sub>×<i>k</i><sub>6</sub><br /><i>q</i><sub>7</sub>-<i>q</i><sub>5</sub>=<i>q</i><sub>6</sub>×<i>k</i><sub>6</sub><br /><br />の関係があるから, <i>p</i>の方は106に193を6回足すと1264になり, <i>q</i>は39に71を6回足すと465になる. この足し算の途中の値<i>r<sub>λ</sub></i>, <i>s<sub>λ</sub></i>を分子, 分母にして計算すると, <i>p</i><sub>5</sub>/<i>q</i><sub>5</sub>から<i>p</i><sub>7</sub>/<i>q</i><sub>7</sub>の間の近似値が次々と得られる. これを中間近似分数(intermediate convergents)という. <br /><br />この計算が次の表である. <br /> <br /> <img src="https://www.iijlab.net/~ew/images/contfrac/cf4.png" width="400" /><br /> <br />そしてこの中に, 問題の878/323があったということだ. <br /><br /><br /><br /><br /><br /><br />和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-83986649251776683252019-04-22T19:59:00.000+09:002019-04-22T19:59:41.857+09:0015パズル最近, 高木貞治先生の「数学小景」を眺めていたら, その中に「十五の駒遊び」という題で15パズルの話があった. <br /><br />ところで, 最近知った15パズルの面白い解き方は, 福岡教育大学の藤本さんの提案する「回転型操作」によるものだ.<br /><br />2件の論文がある.<br /><br />情報処理学会研究報告にある「スライドパズルにおける回転型操作とアクセスビリティ」と, 近着の数式処理にある「Loop generatorによる15パズルの最適アルゴリズムとGod's numberについて」.<br /><br />私は明解な手順には関心があるが, 最短手順には興味がないので, 後の論文は眺めた程度であった. 前の論文も解法は数式処理しシステムGAPを使って記述してあるので, GAPに馴染まぬ凡人にはとんと理解し得ぬ. 以下では私がSchemeで書いたプログラムの考え方を述べる. <br /><br />この回転操作(英語ではloop generatorというらしい)は, 前回の15パズルのブログと似たやり方だが, 経路の単純なのが取り柄だ. 回転操作の前と後では, 空白は右下隅(15)に置くことにする. (イタリックの数字は場所の番地である.) 回転操作にはa, b, cがあり, それぞれ下の図の赤, 緑, 青の線のようにこまを回す. もちろん反対方向にも回す(回転数にマイナスを付けて示す). 例えばaを1回というと, <i>0</i>にあったこまは<i>4</i>へ移動, 2回というと<i>8</i>へ移動する. bを1回では, <i>5</i>のこまは<i>9</i>へ, 2回では<i>13</i>へ進む. 逆回転はは回数を負にする. aを-2回というと, <i>0</i>が<i>2</i>へ行く. <br /><br /><image src="http://www.iijlab.net/~ew/images/15puzzle/base4.png" width="120"><br /><br />回転経路が重なるのは, 図の<i>7, 11, 13, 14</i>の場所で, ここを経由して回転だけを使い, こまを任意の場所から他の任意の場所へ移すのである. <br /> <br />例えば<i>5</i>にあるこまを<i>0</i>へ移すには, aを4回して移動先<i>0</i>を<i>13</i>へ移動し, bを2回で移動元を<i>13</i>へ移し, aを-4回して最終目的地<i>0</i>へ動かす. <br /><br />移動元と移動先が同じループに属しているといささか面倒だ. 例えば<i>4</i>から<i>1</i>へ移動するには,どちらもaのループにあるから, あらかじめ<i>4</i>を内側のループ, bかcのaと重ならない場所へ移動しておく. つまりaを3回で, <i>4</i>のものを<i>13</i>へ; 次にbを-1回で<i>9</i>へ; 移動先の<i>1</i>は, a3回で<i>8</i>へ来ているから, 後2回で<i>13</i>へ進める. そこでbを1回行い, <i>4</i>にあったものを<i>13</i>へ移す. それからaを-5回行うと, <i>I1</i>の位置へ行くわけである. <br /><br />こういう手順が判明したから, aのループにある, 上端と左端にあるべき7個(4,3,2,1,,5,9,13)のこまを目的位置へ移す. この時はb, cのループの場所は自由に使ってよい. その後, bのループの5個(8,7,6,10,14)を揃える. この作業場所はcの10である. bが揃うと残りはcのループの<i>10, 11, 14</i>の番地に11, 12, 15のこまが残るだけになり, もとの配置のパリティが合っていれば, c1回かc-1回で完成する.<br /><br />プログラムは後で示すことにし, 実際に走ったところを見よう. 下の図の左上(A)が最初の状態で, 空白はすでに右下<i>15</i>にある. </image><br />
<br />こま4はbのループにあるから, (a 7)(b 2)(a -7)を行うと, (a 7)で4の目的地<i>3</i>が<i>13</i>へ行き, (b 2)で4が<i>13</i>へ行き, (a -7)で4が<i>3</i>へ移る. (B)になる.<br />
<br /><image src="http://www.iijlab.net/~ew/images/15puzzle/base5.png" width="400"><br /><br />こま3は元々は目的地<i>2</i>にあったのだが, 4の移動のとばっちりで<i>0</i>にいる. そこで(a 4)(b -1)(a -4)(a 6)(b 1)(a -6)を実行する. (a 4)で3を<i>13</i>へ, (b -1)で<i>9</i>へ, (a -4)でaループを一旦元へ戻す. 改めて<i>2</i>の位置を(a 6)で<i>13</i>へ運び, (b 1)で3をそこに置き, (a -6)で(C)のように3が収まる. (a -4) (a 6)と続けるのは無駄だが, プログラムが明瞭になるからこのままにしてある. <br /><br />こま2は3に連られて目的地に来ていたからなにもしない. (D).<br /><br />次は1の番だ. またaのループにあるから, (a 3)(b -1)(a -3)(a 4)(b 1)(a -4)の前半で<i>9</i>の位置に置き, 後半で<i>0</i>へ移す. (E).<br /><br />5はaループだが同時にbループでもある. この場合は, bの中へ入れればよいが, 私のプログラムでは(b 2)で<i>5</i>へ入れることにしている. それから(a 3)(b 2)(a -3)で(F).<br /><br />9は(a 1)(b -1)(a -1)(a 2)(b 1)(a -2)とする. (G).<br /> <br />次の13もa, cループだから, 一旦<i>10</i>へ移す. (c 1). その後(a 2)(c 1)(a -2) で(H). ここまででaループの7個が終わる. 以後はaループの7個には影響しないように注意. <br /><br />(H)を見ると8はb, cループにあるから, (c -1)で<i>10</i>へ避難. それから(b 5)(c 1)(b -5). (I)になる. <br /> <br />7はcループにあるから簡単だ. 6の場所を(b 4)で<i>14</i>へ. それから(c 1)(b -4)で(J). 幸運にもついでに6, 10も揃うから(K),(L)も同じ図だ. <br /><br />14はb, cループ上にいるから(c -1)で<i>10</i>へ引上げ, (b 1)(c 1)(b -1).<br /><br />最後は11, 12, 15のcループ. 11が<i>10</i>にあれば何もしない. <i>11</i>にあれば(c 1), <i>14</i>にあれば(c -1)だ. (N).<br /><br />全体の手は<br /><pre>
(a 7)(b 2)(a -7);4
(a 4)(b -1)(a -4)(a 6)(b 1)(a -6);3 2
(a 3)(b -1)(a -3)(a 4)(b 1)(a -4);1
(b 2)(a 3)(b 2)(a -3);5
(a 1)(b -1)(a -1)(a 2)(b 1)(a -2);9
(c 1)(a 2)(c 1)(a -2);13
(c -1)(b 5)(c 1)(b -5);8
(b 4)(c 1)(b -4);7 6 10
(c -1)(b 1)(c 1)(b -1);14
(c -1);11 12 15</pre>
<br />であった.<br /> <br />一方, Schemeによるプログラムはこんな具合いだ. <br /><pre>
(define bs (range 0 16))</pre>
<br />は盤面のこまのリスト. <br /><pre>
(define atab '(
(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14)
(4 0 1 2 8 5 6 3 12 9 10 7 13 14 11)
(8 4 0 1 12 5 6 2 13 9 10 3 14 11 7)
(12 8 4 0 13 5 6 1 14 9 10 2 11 7 3)
(13 12 8 4 14 5 6 0 11 9 10 1 7 3 2)
(14 13 12 8 11 5 6 4 7 9 10 0 3 2 1)
(11 14 13 12 7 5 6 8 3 9 10 4 2 1 0)
(7 11 14 13 3 5 6 12 2 9 10 8 1 0 4)
(3 7 11 14 2 5 6 13 1 9 10 12 0 4 8)
(2 3 7 11 1 5 6 14 0 9 10 13 4 8 12)
(1 2 3 7 0 5 6 11 4 9 10 14 8 12 13)))
(define btab '(
(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14)
(0 1 2 3 4 9 5 6 8 13 10 7 12 14 11)
(0 1 2 3 4 13 9 5 8 14 10 6 12 11 7)
(0 1 2 3 4 14 13 9 8 11 10 5 12 7 6)
(0 1 2 3 4 11 14 13 8 7 10 9 12 6 5)
(0 1 2 3 4 7 11 14 8 6 10 13 12 5 9)
(0 1 2 3 4 6 7 11 8 5 10 14 12 9 13)))
(define ctab'(
(0 1 2 3 4 5 6 7 8 9 10 11 12 13 14)
(0 1 2 3 4 5 6 7 8 9 14 10 12 13 11)
(0 1 2 3 4 5 6 7 8 9 11 14 12 13 10)))</pre>
<br />はa, b, cそれぞれの回転で, <i>0</i>,...のこまがどこへ行くかを示す. atab, つまりaの表は, 最初の行が(a 0)に対応. 同じ場所に留まる. 次の行(4 0 1...)はaの1回転で, <i>0</i>は<i>4</i>へ, <i>1</i>は<i>0</i>へ移ることを示す. 最後の行(atab[10])は-1回のものだ. b, cについても同様. <br /><pre>
(define (rotate tab)
(let ((b (make-list 16 0)))
(do ((i 0 (+ i 1))) ((= i 15))
(list-set! b (list-ref tab i) (list-ref bs i)))
(set! bs (list-copy b)) 'ok))
(define (a n) (rotate (list-ref atab (modulo n 11))))
(define (b n) (rotate (list-ref btab (modulo n 7))))
(define (c n) (rotate (list-ref ctab (modulo n 3))))</pre>
<br />(a n)がaのループをn回実行する. aの表から使うべき行をとりだし, rotateへ. rotateは16個のリストbを用意し, bs内のこまの番号をbへ移して最後にbをbsへ戻す.<br /><pre>
(define ps0 '(4 3 2 1 5 9 13))
(define qs0 '(3 2 1 0 4 8 12))
(define ps1 '(8 7 6 10 14))
(define qs1 '(7 6 5 9 13))</pre>
<br />ps0はaループのこまの番号, qs0はその目的地. ps1とqs1はbループのものである. 次がいよいよ解くプログラムである. <br /><pre>
(define (solve2)
(do ((i 0 (+ i 1))) ((= i 7))
(let* ((p (list-ref ps0 i)) (q (list-ref qs0 i))
(r (elemindex p bs)) )
(if (not (= r q)) (begin (set! r
(case r
((7) (b 2) 5)
((11) (c 1) 10)
((14) (c -1) 10)
((13) (b -2) 5)
((0 1 2 3 4 8 12) (let ((l (length (member r qs0))))
(a l) (b -1) (a (- l)) 9))
(else r)))
(if (= r 10)
(let ((l (- 8 i))) (a l) (c 1) (a (- l)))
(let ((l (- 7 i))) (a l)
(b (length (member r '(6 5 9)))) (a (- l))))))))</pre>
<br />ここまでがaループのプログラムで, 7個についてdoループを回す. letへ来て, pはこまの番号, qは目的地, rは現在地である. r=qなら何もしない. そうでないならcase式へ来て, r=7なら(b 2)を実行, こまを5へ, という風に読む. rが0,1,2,3,4,8,12なら, rがqs0にある位置から後方の長さlを知り, (a l) (b -1) (a (- l))を行ない, こまを9へ移す. 上のいずれでもなければelseへ来て, 回転はせず, rはそのまま. 上のbeginの次の(set! r ...)でrを更新する.<br /><br />結局移動するこまはcの<i>10</i>かbの<i>5,6,9</i>のどこかにいるから, それらの情報を使い, 目的地へ移動する. <br /><br />bループのプログラムが以下だが, 説明は同様だ. <br /><pre>
(do ((i 0 (+ i 1))) ((= i 5))
(let* ((p (list-ref ps1 i)) (q (list-ref qs1 i))
(r (elemindex p bs)))
(if (not (= r q)) (begin
(case r
((11) (c 1))
((14) (c -1))
((7 6 5 9 13) (let ((l (length (member r qs1))))
(b l) (c -1) (b (- l)))))
(let ((l (- 5 i))) (b l) (c 1) (b (- l)))))))
(case (elemindex 11 bs)
((11) (c 1))
((14) (c -1)))
(drawbs))</pre>
<br />という次第で無事に並べ替られる. なかなか楽しい.<br /></image>和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-34363045505699992842019-03-11T13:00:00.000+09:002019-03-12T07:51:50.130+09:00Lyonsの7テストMartin Gardnerは面白い話を沢山書いた. 最近, Lyonsの7テストというものを読んだ. ニューヨークの精神神経科医 L. Vosburgh Lyonsが考案したという. <br />
<br />
7で割ったときの剰余を計算するのは面倒だが, Gardnerによるとこれが一番効率がいいそうだ.<br />
<br />
まず例から. 14桁くらいの整数を考える. n=12340123456789とでもするか. 右から2桁ずつ区切り, 12 34 01 23 45 67 89とし, それぞれの区画の2桁の数の7での剰余をとる. つまり上の例では<br />
<br />
5 6 1 2 3 4 5<br />
<br />
が得られる. これを3個ずつ揃え<br />
<pre></pre>
3 4 5<br />
6 1 2<br />
5<br />
ーーーーー<br />
<br />
<br />
それぞれをmod 7で足す. 2 5 5になるわけだ. 右の2個を55のように読み, 左の2個を25のように読み, 右の7の剰余6から左の7の剰余4を引くと元の数の7の剰余2が求まる.<br />
<br />
なぜうまくいくか. とりあえずn=123456のような6桁だけで考えてみる. <br />
<br />
2桁ずつをa, b, cとすると n=10000a+100b+c. n mod 7が欲しいから<br />
<br />
10000 mod 7 = 4, 100 mod 7 = 2を入れて<br />
<br />
n mod 7 = 4a+2b+c. こうして右の2桁 10b+c から 左の2桁 10a+b を引いたから,<br />
<br />
10b+c - (10a+b)=-10a+9b+c. -10=4 mod 7, 9=2 mod 7だから<br />
<br />
4a+2b+cになり目出度し目出度し.<br />
<br />
6桁より大きいときは n=10<sup>6</sup>a+bだが 10000=4 mod 7, 100=2 mod 7 だったから1000000=8 mod 7=1 mod 7で, n=a+bと同じになり, 6桁ずつ束にして扱えばよいことがわかる.<br />
<br />
なるほど.<br />
<br />
この方法が優れている理由を冷静に考えてみると, 2桁の数を7で割った剰余は簡単に得られることを利用しているからだ. 我々は7の段 7,2 14; 7,3 21; 7,4 28を諳じているから, それに70を足した84, 91, 98が7の倍数であることも承知している.<br />
<br />
3桁を7で割るのはもっと面倒だが, それが出来るなら, Gardnerがいうように7の剰余は次のように得る. <br />
<br />
n=12340123456789 を右から3桁ずつに分ける.<br />
<br />
12 340 123 456 789 それぞれの7の剰余をとる.<br />
<br />
5 4 4 1 5<br />
<br />
右から交互に 5 - 1 + 4 - 4 + 5 = 9 = 2 mod 7<br />
<br />
この方法は11での剰余なら<br />
<br />
1 10 2 5 8<br />
<br />
8 - 5 + 2 - 10 + 1 = -4 = 7 mod 11. n mod 11 = 7.<br />
<br />
13の剰余なら<br />
<br />
12 2 6 1 9<br />
<br />
9 - 1 + 6 - 2 + 12 = 24 = 11 mod 13. n mod 13 = 11.<br />
<br />
というふうに使えるが, 11や13の除数を得るのに計算機を使うなら最初から計算機を使えばよいから, そんなにメリットはない.<br />
<br />
ところで今の手品の秘密は 7×11×13=1001 にある.<br />
<br />
<br />和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0tag:blogger.com,1999:blog-8851344890851282373.post-56857264592046587952019-02-10T11:57:00.000+09:002019-02-10T11:57:03.366+09:00中学入試問題数日前の新聞に, 開成中学の数学の入試問題が掲載されていた. そのうち1問をやってみたが, 最近の小学生はこれが解けるのかと感心した.<br /><br />下の図の左に示す直方体ABCD-EFGHがある. DC上にPを, DP=8, PC=12になるようにとる. EF上にQを, EQ=4にとる. CG上にRを, CR=9にとる.<br /><br />直方体をPQRを通る面で切り, 断面をXとする. <br /><br />そのXをABFE側から(以下南から)見た面積は228; ABCD側から(以下上から)見た面積は266.<br /><br />この時高さAE, 奥行ADを計算する.<br /><br /><br /><br />すぐ分るのはPRはXの辺であること. すると<br /><br />平行な2面と, これと交わる他の面との交線は平行であり, 3次元空間での平行線はどの方向から見ても平行に見える.<br /><br />ことから, Qを通るRPと平行な線はXの辺である. その線とAEとの交点をSとする.<br /><br /> PR || SQ.<br /><br />Xの辺には, PからADへ向って引く線のADとの交点をT, QからFGへ向った引く線のFGとの交点をUとする. PT || UQ, TS || RU.<br /><br /><image src="http://www.iijlab.net/~ew/images/exam/kaisei0.png" width="400"> <br /> <br />Xの大体の形が分ったので, 南から見た図を作ると下のようになる. 太線の内側の面積が228, 左下の三角の面積が6, 右上の三角の面積が54, 長方形の面積は228+6+54=288で幅が20だから高さは14.4となる.<br /><br /><image src="http://www.iijlab.net/~ew/images/exam/kaisei1.png" width="400"> <br /><br />次は東から見た図である. 高さが分ったから, 図のAS=11.4とRG=5.4が分る. TSとRUが平行だから, ある比例定数aに対してAT=11.4a, UG=5.4aとする.<br /><br />さらに上から見た図では, DP=8とQF=16が分っており, TPとQUが平行だから, DT=bとすると, FU=2b.<br /><br />直方体の奥行=ST+TD=FU+URだから<br /><br /> 11.4a+b=2b+5.4a<br /><br /> これから b=6a<br /><br /> 従って 奥行=17.4a<br /> 左上の三角=24a, 右下の三角=96a,<br /> 266+24a+96a=17.4a×-120<br /> 226=348a-120a=228a a=266/228,<br /><br /> 故に奥行=17.4×266/228=20.3<br /><br /> <br />とまぁすらすら書いたが, 結構図を何度も書きなおし, 沈思したことも何度かあった. 私が子供の頃はこんな問題はなくて助かったなぁ.<br /><br /></image></image>和田英一http://www.blogger.com/profile/03766161406401142226noreply@blogger.com0