五角ミノ十二面体

英語でquintominal dodecahedraというクイズである. 十二面体のクイズといえば, Hamilton順路が有名だが, これは全然別だ. 今年4月に新型コロナで死去したJohn Conwayが1958年頃考えたらしい.

Conway他の書いた"Winning Ways Vol.4"の表紙に下のような興味ある図があった.


正五角形の5辺を5色(例えば赤青橙緑黒)で塗ると, 回転, 裏返しは同じと見た時, 12通りの五角形が得られる. (5!/5/2=12) 図の下にある12の五角形だ. この一組をquintomino(「五角ミノ」と訳すのは如何?)という. 一方, 正十二面体は12の正五角形か構成されるから, 辺が同じ色になるように, 各面に五角ミノを1枚づつ貼れるかというのが問題である.

Winning Waysの表紙は, 12の五角ミノにIを除いてAからMまで名をつけ, 解答を描いたものである. 各五角ミノの色の配置は以下の通り.

A=01234 B=01243 C=01324 D=01342 E=01423 F=01432
G=02134 H=02143 J=02314 K=02413 L=03124 M=03214

答をいうとこれは可能であり, 等価性を考慮すると3通りの解があるという. 表紙の3つの絵はAを底面にしたその3通りを示している.

KnuthのTAOCP, 第4巻, 第5分冊はDancing Linksが話題で, 膨大な数の演習問題があるが, その中にquintominal dodecahedron問題を解けというのがあった(ex 7221-136)ので, 挑戦してみたというのが今回のブログである.

正十二面体は下左のようなもので, こういう立体の絵では扱いにくいから, 通常は下右のような形相図(高木貞治 数学小景)を使う. 数学小景には正十二面体を床に置き, 上面の正五角形の少し上の光源で照すと, こういう影ができると書いてある.


正十二面体には, 頂点(青)と辺(赤)と面(緑)に番号を付けたが, 同じ番号を面0を底にした形相図に付けると下のようになる. 右端は面にも番号を付けたものだ. 0番の面がないのは, それは周囲の大きい五角だからだ.


ペントミノのとか8クィーンのような, ピースを置けるかどうかの問題は, 現在置いてあるピースの状態と, これから置くピースの集合を引数として受取り, あるピースがさらに置けるなら, その状態と今残るピースを新しい引数にして自分を再帰呼出しする関数を書くと解けるのが普通である.

以下のSchemeのプログラムもその方針で出来ている. colorsは12の五角ミノのリスト. facesは各々の面を作る辺の番号のリストである. allcsは色のリストを回転したり裏返したりしたもののリストを作る.

(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)))))))  

状態stateは辺の番号順の30の色のリストで, 先頭には最初の五角ミノAを固定し, 未定の25の場所は()にする. 数のリストである状態を数字列に変換したり, 逆に変換したりするのがcompressとexpandである.

(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 '())))

次のtryが再帰呼出しの関数で, 状態stateと残りの五角ミノrestを受け取る. matchはその状態でposの場所がcolの色と一致しているか()かであれば真を返す. 要するに置けるということである. fillcはposの場所をcolの色にした新しい状態を返す. 新しい状態は再帰呼出しから戻るときに捨てるから, 置いたピースを戻す仕事はしない.

ifから後がtryの本体で, 残りの五角ミノがなければ((length(rest))=0なら)解があったとして出力する. 五角ミノの1個をとってcolとし, 残りをnewrestにし, 今度置く面をposにし, colを回転, 反転してcolsにし, そのそれぞれcについて置けるようなら再帰呼び出しする.

(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)))

Aが置いてあるので0を除き, restを(1 2 ... 11)としてtryを呼ぶと60個数の解が次のように得られる.

(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"

後ではこの文字列のリストをallsolsとして使う.

ところで, Conwayの解は3通りであった. 上の60個の中には, 色の置換えについて等価のものが沢山あるということだ. 次に等価を発見する作業をやってみよう.

上の最初の解"012343421042143042321400312301"は形相図に書き入れてみると, 下左のようになる. その3と4を交換すると中の図になる. 周囲は01234ではなくなるが, よく見ると下の五角形(面の番号1)の周囲が01234である. そこで辺の接続関係を保ったまま, これが周囲(面0)に来るように図を書き直すと右のようになる. この列を圧縮したのが"012343421420130204031343214201"で, 先程の解の2番目であった. だからこの2つの解は等価であった. これを上の全ての解のすべての置換について実行し, 同じになったものをまとめると3通りになる.


以下はそのプログラムだ.

repはcolのprm0をprm1で置き換える. つまりcolの(0 1 2 3 4)を(0 1 2 4 3)で置き換える(3と4を交換する).

  (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))

(0 1 2 3 4 3 4 2 1 0 4 2 1 4 3
0 4 2 3 2 1 4 0 0 3 1 2 3 0 1))
は
(0 1 2 4 3 4 3 2 1 0 3 2 1 3 4
0 3 2 4 2 1 3 0 0 4 1 2 4 0 1)
になって左と中の絵に対応する.

remapは置き換えた色のリストを, (0 1 2 3 4)が先頭になるように(面0になるように)位置を変更した辺のリストを作る. newposとしよう. これを基本位置という. つまり右の絵にする. これは多少面倒だ.

まず後述の関数newfaceを置換後の色のリストに使うと, colorsに対応する五角ミノがいる位置が得られる. このリストをnewcolとすると, newcolは

((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))

これを見ると五角ミノ(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)にいった.

newposを作る方針は次のようだ. 基本位置に来る五角ミノの現在位置を設定する. (0 1 2 3 4)は(0 8 7 6 5)にいたので, newposは(0 8 7 6 5 ...)となる(左の絵).


先頭の(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を置く.

(0 8 7 6 5 4 ...) まで判明した(中の絵). すると後は判明した場所を隣同士にもつ五角ミノを探す作業を続ける.

まず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が担当する.)

(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)))

prmは(0 1 2 3 4)の順列. prem0はその先頭(0 1 2 3 4). equivは等価の解をまとめたリストのリスト. doループで実行が始まる.

doループではallsolsから解を1個とりcolとし, equivの表のどれかに既にあればなにもしない. どれにもないなら, そrwれを1個もつリストをequivに追加し, そのリストをeqlistとする. colを120通りの色の置換えをし, 置き換えたnewcolがqlistになければ追加する.

(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))))))))))

それぞれの等価グループの長さを出力し, 各グループを昇順に並べ替え出力する.

(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

Winningやtaocpについてはまた述べる.

プログラミングコンテスト

一週間くらい前にプログラミングコンテストがあったらしい. ツィッター のキーワードHHKBを見たら, コンテストの話題が沢山あった.

私はプログラムを急いで書くコンテストは嫌いだ. プログラムはゆっくり と考えながら書くものとDijkstraは主張する. 私もまったく同感で, ずい ぶん前に大学対抗プログラミングコンテスト(ICPC)のオフィシャルを2年 ほど引き受けたが, プログラムをあわてて書き, テストの入力データに対 して正しい数値だけ得られれば完成で, 完成までの時間の短いのを上位 とする風潮に賛成できず, オフィシャルを辞退した.

しかし世の中ではどんな問題でコンテストをやっているかとつい問題を覗いたのが 運の尽き. 解いてみたくなり, プログラムに取り憑かれた.

やったのはsquareで, 二次元空間の整数値の格子点に角を持ち, x, y軸に 平行な辺を持つ正方形の入れ子問題である. すなわち, 一辺の長さnの正方形N の内部に, 一辺の長さがそれぞれaとbの正方形AとBを重さならないように置く 置き方が何通りあるかを計算する.

取り敢えず例題にあった n=4, a=2, b=2で考える. aもbも2でややこしい が, まぁよかろう.

n×nの正方形Nの中に, 2×2の正方形Bの置き方は, 下の図のよう に3×3で9通りある. Bの左の線はNの左端から, Nの右端からBの幅bを 引いたところまで置けて, 両端を含めるから n-b+1個. 今の数値では 4-2+1=3で, これが横方向も縦方向もあるので9通りになる.


縦, 横をそれぞれl₀, l₁とする長方形内のBの置き方は,
bs(l₀,l₁)=(l₀-b+1)×(l₁-b+1)
だから, 全体に置くのはbs(n,n), 今はbs(4,4)=9である. 今後この値をtで表す.

2×2の正方形Aの置き方も, 下の図の実線の枠が示すようにで9通りである. その各の位置で, Bの一部でも隠す勢力範囲は, 灰色の部分で, この勢力範囲内の Bの数をt, つまり9から引くと, そのAの場合の数が得られる. これをすべての Aの場所で計算すれば答が得られる.


ところで, 横方向だけ考えて, 勢力範囲の左右の幅は, Aの左端の, Nの左端からの距離をxとすると, 次の下左図の太い縦線のようになる. この下限をy₀, 上限をy₁とすると, 赤線の範囲になる. 右は一般の場合の図だ.


この図から分るように

y0(x) = max (x-b+1, 0),
y1(x) = min (x+a-(b+1), n)

だから, Aの左下の座標を(x₀, x₁)とすれば, 勢力範囲は 横がy₁(x₀)-y₀(x₀), 縦がy₁(x₁)-y₀(x₁)で, この中に入るBの数をtから引くことになる. つまり t-s(y₁(x₀)-y₀(x₀),y₁(x₁)-y₀(x₁)).

私は通常, Lispの一方言のSchemeでプログラムを書くので

(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)))
))))))

実行すると

(square 4 2 2) => 32
(square 3 1 2) => 20

となる. Schemeが苦手な人のためにPythonのプログラムを示すと

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))

実行すると

20
20
32

プログラムは一旦うまく動くようになると, 改良の糸口が見えだすものだ. 仮に勢力範囲の図で, 灰色の部分の面積を足すとすると, 灰色の横の長さを足し合わせ, 縦の長さも足し合わせて, 和どうしを掛ければよいようだ. いちいちtから引く代わりに, Aの置き方の数にtを掛けたものから, 勢力範囲の 総面積を引けばよい.

勢力範囲の幅は, 図のy₀, y₁の間の縦の太線を足せばよい. 両端は 少し短いが, 破線の部分も入れてしまえば, 標準の縦の長さ a+2b-2に太線の本数 n-a+1を 掛ければよい. 破線の部分は0からb-1までの和で, 左下と右上の2箇所にあるから, (b-1)b である.

そうそう, 勢力範囲lに入るBの数がl-b+1だったように, 太線の長さからもb+1を 引く必要がある. (プログラムの(* s (- b 1)) のところ.)

下のプログラムで, tは前述の通り, sは太線の数, pが太線の総延長である.

(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

Python版は

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))

実行結果は

409369707

まずまず楽しんだ.

連分数と近似分数

前回のブログ, 連分数と近似分数でeの近似分数 878/323をStern Brocot木から求めたが, 連分数との関係については書かなかった.

容易に想像できるように, LとRの列RL0RLR2LRL4RLR6L ...は, 1年前のブログにあるeの連分数展開[2,1,2,1,1,4,1,1,6,1,1,8]と同じパターンである.

昔のブログの最後の方の中間近似分数にも同じ分数が現れる.

「コンピュータの数学」には無理数αからLとRの列を求める式もある.

if α < 1 then (output(L); α:=&alpha/(1-&alpha))
           else (output(R); α:=α-1)

途中のαと一緒にLとRを計算すると

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']

つまり連分数で1なら文字は1回, 2なら2回 繰り返すのであり, Stern Brocot木による計算は連分数展開の近似を, 中間近似分数まで含めて得ていたのであった.

連分数と近似分数

 昨年6月のブログで, eの近似分数 878/323の求め方を話題にした. 先月出版された「コンピュータの数学 第2版」を拾い読みしていたら, Stern Brocot木の関連でeの近似分数を求める話があった.

Stern Brocot木は, 互いに素なmとnの分数を大きさの順に並べたもので, 次のような形である. ある場所の値は, その左上の一番近い祖先m, nと右上の一番近い祖先m', n'から m+m'/n+n' で得られる.

最上部の1/1からLでは左下へ, Rで右下へ辿ると, 色々な分数が一意的に表せる. そしてeは

e=RL⁰RLR²LRL⁴RLR⁶L ...

つまりRL⁰R, LR²L, RL⁴R, LR⁶L, ...

であって, Eulerが24歳の時に発見したと書いてあった.

早速計算するプログラムをSchemeで書いた.

(l)は左下へ, (r)は右下へ進む関数. lmnは左祖先, rmnは右祖先, cmnは現在の分数, nmnは次の世代の分数である.

(define (show) (display cmn))
(define lmn '(0 1)) (define cmn '(1 1)) (define rmn '(1 0))
(define (reset!) (set! lmn '(0 1)) (set! cmn '(1 1))
  (set! rmn '(1 0)) (show))
(define (l) (let ((nmn (map + lmn cmn)))
 (set! rmn cmn) (set! cmn nmn) (show)))
(define (r) (let ((nmn (map + cmn rmn)))
 (set! lmn cmn) (set! cmn nmn) (show)))
(for-each (lambda (c)
 (cond ((eq? c 'l) (l))
       ((eq? c 'r) (r))
       (else (reset!))))
  '(c r r l r r l r l l l l r l r r r))

実行すると,

(1 1)(2 1)(3 1)(5 2)(8 3)(11 4)(19 7)(30 11)(49 18)(68 25)

(87 32)(106 39)(193 71)(299 110)(492 181)(685 252)

(878 323)

すげー.

XYプロッター

 IEEE Specturm 09.20を眺めていたら, AxidrawというXYプロッターの記事があった.  それに下のような図があったので, その機構を考えてみた.

中央に左右に広がる基盤があり, これは固定されている. その左端と右端にそれぞれモーターと同軸で回転するプーリーがある(L, R). 縦に長い板は, 縦方向を保ったまま, 上下左右に移動し, 下に自由に回転するプーリーMがあって, その辺りに上げ下げできる筆記具がある. それでXYプロッターになるわけだ.

3個のプーリーは, 中央の4個の補助プーリーを介して1本のベルトで結ばれており, ベルトの最後は, 縦の板の上端である.

まずLを固定し, Rの周囲が長さaだけ右回転したとする. その際, 縦の板が上下だけに動いたとすると, 上端に左側のベルトもaだけ弛むが, Lが回転しない約束なので, その弛みは上下の板を右へ移動して吸収することになる.

というわけで, Rがaだけ回転すると, 上下の板は右へa/2, 下へa/2だけ移動する. そうすると, 右上部分のベルトは, x1もy1もa/2だけ短くなり, Rはaだけ回転し, 右下部分ではx1はa/2短く, y0はa/2長くなるので, Rから出たaはMを回って左下へ行く.

左上はy1がa/2短く, x0がa/2長くなるから, 予定通りLは回転しない. 左下はx0もy0も伸びるが, その分だけMから貰うので万事うまくいく.

つまりRの右回転は上下板を右下へ移動させ, 左回転は左上へ移動させる. 同様にLの左回転は上下板を左下へ移動させ, 右回転は右上は移動させる.

左右だけの移動, 上下だけの移動は, LとRの回転を合成すればよい.

なるほどうまい仕掛けであった.

恒星時計

 恒星時計(こうせいじ-けい)を実装した. ここからダウンロードできる.

 ダウンロードした時から走行し, 下の上の図のような形である. 右上のディジタル表示は現在の太陽時の年月日と時刻(時, 分, 秒). 左上の表示は現在の(地方)恒星時の時, 分, 秒である.

中央の時計がその恒星時をアナログ的に表示する. 恒星時は0時から24時まであるが, 時計の文字盤には時の表示が12個しかなく, この時計の時針は恒星時の12時間で1回転する. その代わり, 恒星時の正時に文字盤の時針と反対の位置の文字が12増えたり12減ったりする.

この図では, 時針は9時と10時の間にあり, そこから左周りに6時間は過去の時刻を, 右周りに6時間は未来の6時間を示す.

時計の周囲の四角形の内部をクリックすると, 時計は太陽時の秒が繰り上がった時点で停止し, その時刻での恒星時がミリ秒まで表示される. 下の時計の図は停止中の様子を示す. もう一度クリックすると走行を再開する.

時計の盤面に, 北極星を中心とする北天の星図が表示してある. この星図は星座早見盤の自作キットから使わせて頂いた. この星図は恒星時の24時間で左周りに1回転する. 通常の時計の6時相当する真下を時角0h, 3時に相当する右を時角6h, 真上を12h, 左を18hとして北天が日周運動する.

星図に薄く天の川が表示されている. 天の川が天の北極に近いのが赤経0hの方向で, 中央からその方向に辿ると天の川の中にカシオペアのWがあり, さらに外周へ探すとペガススの四角がある. その少し外の緑の点が春分点(Υ)で, 春分点の時角と恒星時が一致しているのが分る.

この時計の星図の下2/3くらいの横長の楕円が, その時刻に見える星空になる.

次の大きな円は星図の元の図である. 右の方の3本の左向き矢印の先の星が, 上からうお(Pisces)のλとω, くじら座(Cetus)のιである. その下の図がうお座とくじら座で, これから見ると春分点はこれら3つの星の中央にあり, 恒星時計の文字盤の春分点をその位置にマークした.

この時計では, ある太陽時の日時, 時刻の恒星時を知ることが出来る. 時計の枠内でキーボードから数字を入力すると時計が停止し, その数字を空白で区切りながら, 年, 月, 日, 時, 分, 秒と入力すると, 対応する恒星時が示される.

最初の時計の図の時刻では, 2020⨆8⨆8⨆11⨆53⨆9と入力すると, 左上に恒星時が現れ, 文字盤にはその時刻に現れる星空が表示される. ミリ秒の値は少し違うが.

恒星時の計算は暗算としてはこのようにする. 3月21日頃の春分には太陽が春分点にいるから, 昼頃, 太陽の南中時が恒星時の0h. それに太陽時の正午からの差を加える. 太陽の赤経は, 1ヶ月に2h増えるから, 例えば8月8日, 立秋の頃なら, 正午に太陽が南中するとき, 春分点は4ヶ月半近く, 9hだけ時角が離れる. つまり正午が9hである. 従って8月8日, 21時だと9h+9h=18hとおよその見当がつく.

しかしこの方法は時計には使えない. 理科年表や天文年鑑には, その年の毎日の世界時0時のグリニジ視恒星時の表があり, インターネットで探すと多くの関連ページが見付かる.

国立天文台の グリニジ恒星時計算HPで2020年1月1日から世界時0時の恒星時を1ヶ月表示したのの一部が次の図である. (出力時刻もGMT)

理科年表では0.1秒, 天文年鑑では0.1分まで表示されているのが, このHPでは0.001秒の精度で表示される.

インターネットにあるいくつかの恒星時計算アルゴリズムで, 2020年1月1日世界時0時の恒星時を計算して見ると, どういうわけかこの表の値にならないので, この値が得られるような数値を自分で計算した.

あるユリウス日がjd 0の日時での恒星時st 0が既知の時, ユリウス日がjdでの日時での恒星時stは
 

st=st 0+1.0027378×(jd-jd 0)


だから, jdとstのペアがあれば, あるjd 0でのst 0が計算できる.
 

jd 0としては, 何かの恒星時計算アルゴリズムにあった切断(?truncated)ユリウス日の起算時点(jd=2440000.5)を使うことにした. つまり1月中の0時のjdと対応するstから31個の
 

st 0=st-1.0027378×(jd-2440000.5)

を求めて平均し, st 0を作る. 恒星時の表はミリ秒のリスト
(226 778 330 882 435 990 548 109 672 237 802 365 925 481 34
585 136 688 243 801 362 925 487 49 609 167 721 274 824 374 924))

を用意し, ミリ秒が減った時は前日の秒に237秒を, 増えた時は236秒を足し, ミリ秒その日の値としながら, 1月の恒星時を入力した. そして定数st0を得たが, それで12月の恒星時を計算すると, 天文台の表と多少ずれるので, 12月についても同様に計算し, 両者の定数の平均を使うことにした. つまり
 

st=0.6712386796084648+(1.00273791×(jd-2440000.5))

を使う.

この値を使う, Pythonで書いた恒星時計算プログラムは次の通り.

def jd(y,m,d):
    def leap(y):
      return y%400==0 if y%100==0 else y%4==0
    return y*365+y//4-y//100+y//400-(489*(13-m)+26)//16\
      +((1 if leap(y)else 2)if m<3 else 0)+d+1721426
def sidereal(y,mo,d,h,mi,s):
    def frac(a):
      return a-int(a)
    ut=jd(y,mo,d)+(h-9+(mi+s/60)/60)/24-0.5    #-9: to GMT   
    st=0.6712386796084648+1.00273791*(ut-2440000.5)+9/24

    #+9/24: to JST 
    h=frac(st)*24
    m=frac(h)*60
    s=frac(m)*60
    ms=frac(s)
    return [int(h),int(m),int(s),ms]
print(sidereal(2020,1,1,9,0,0))


出力は[15, 40, 28, 0.20889443391934037]で, 前掲の1月1日の値 6h40m28.226sに9時間足したのにかなり合っている.

このような恒星時計が, 宮沢賢治の「銀河鉄道の夜」に出て来る 白鳥, わし, さそり, ケンタウルスなどの駅のホームや操車場やアルビレオ観測所で動いていると想像すると楽しい.

「ケンタウルス, 露をふらせ.」

老子

 夏目漱石が護国寺で運慶を見たとか(「夢十夜」の「第六夜」), 芥川龍之介が市電から寒山, 拾得を見とか(「寒山拾得」)いう話があるが, 寺田寅彦が「電車で老子に会った話」は趣きが異なる.

つまり寺田先生は「インゼル・ビュフェライ叢書」にアレクサンダー・ウラールが訳した『老子』を見付けて購入し, 電車の中で読み始めたら面白くなり, 数回の車中で読了したという. 日本語の解説より分りやすかったらしい.

私は老子は近寄り難く避けていたが, 昨年末, 中国深圳で少し話をしたら, 謝礼として為書きのある書を貰った. それに「大道至簡」と書いてあり, 調べたところ老子の言葉らしいことが判明. それではというので老子を読むことにした.

早速読んだのは

池田知久 「老子 全訳注」 講談社学術文庫

だが,

保立道久 「現代語訳」 老子 ちくま新書

もインターネットに

保立道久の研究雑記
『老子』、現代語訳・原文・読み下し

があったので, それも参考にする.


老子は第1章から第37章までの道経と, 第38章から第81章までの徳経の, 計81章で構成される. 今回は夫々の最初の章の最初の文を例にする.

池田本では, 第1章の最初は

道可道也、非恆道也。名可名也、非恆名也。

道の道とす可きは、恒の道に非ざるなり。名の名とす可きは、恒の名に非ざるなり。

道というものがあれこれ唱えられている中で、道として定立することのできるわたしの道も、恒常不変の真の道であるとは言えない。道についての名(理論)があれこれ説かれている中で、道の名(理論)として表現することのできるわたしの名(理論)も、恒常不変の真の名(理論)であるとは言えない。

一方, 保立本は

道可道也、非恒道也。名可名也、非恒名也。

道の道(ゆ)くべきは、恒なる道に非(あら)ざるなり。名の名づくべきは、恒なる名に非ざるなり。

普通に行く道と、ここでいう「恒なる道」はまったく違うものだ。普通に名づけることができる名と、ここでいう「恒なる名」もまったく違う。

どちらも, 特に後者は短いだけあって, 禅問答のようだ.

次に第38章は, 池田本で

上徳不徳、是以有徳。下徳不失徳、是以徳无。

上徳は徳ならず、是を以て徳有り。下徳は徳を失わず、是を以て徳无し。

そもそも最上の徳は、世間的な徳とは正反対である。だからこそ、真の徳がある。下等な徳は、世間的な徳を捨てることが出来ない。だからこそ、真の徳がないのだ。

保立本は

上徳、不徳是以有徳。下徳、不失徳是以無徳。

上徳は、徳ならずして是(ここ)を以て徳あり。下徳は、徳を失わずして是を以て徳なし。

「道」から発した最上の「徳(いきおい)」（上徳）は徳(はたらき)がないようにみえて大きな徳(いきおい)があり、そうでない「下徳」は徳(はたらき)を失っていないようだが実は徳(いきおい)がなくなっている。

同じような感じだ. もともと道とか徳とかが抽象的だから, 説明が難しくて当然であろう.

寺田先生がドイツ語訳で分ったような気分になったというので, 英訳を探すと

The Tao Te Ching by Lao Tzu
馮家福 Jane English

があった. 表題のTao Te ChingはTao=道, Te=徳, Ching=経のことだ.

道の説明は次のようだ.

One

The Tao that can be told is not the eternal Tao.
The name that can be named is not the eternal name.

また徳の方は

Thirty-eight

A truly good man is not aware of his goodness,
And is therefore good.
A foolish man tries to be good,
And is therefore not good.

上徳と下徳がgood manとfoolish manになっているのが凄い. なるほど日本語の説明より具体的である.

インターネットでさらに探すと, 寺田先生が読んだAlexander Ularの訳の一部があった.

Die Bahn und der rechte Weg des Lao-Tse
Alexander Ular

DER ERSTE SPRUCH

Die Bahn der Bahnen ist nicht die Alltagsbahn;
Der Name der Namen ist nicht der Alltagsname.

38章は掲載されていなかった.

ドイツ語を和訳するのは簡単だ. 複数をたちと訳せば

道たちの道は毎日の道ではない
名たちの名は毎日の名ではない

でもこれが何を言おうとしているかは, やはり分らない.

結果的には池田本の訳が一番しっくりする. つまりドイツ語が分り易いというのも章によるのではないかということである. 沢山の訳を比較すれば, 段々と理解しやすくなるという当然の結論に達する.


ところで論語だが, 大学の先生方がよくいうのが

学而不思則罔 思而不学則殆

学びて思わざれば則ち罔(くら)し、思いて学ばざれば則ち殆(あやう)し。

これは罔や殆が理解を妨げている. 私の持っている論語の英訳を見ると

The Analects of Confucius

The Master said: `Learing without thinking is useless. Thinking without learning is dangerous.'

のように, uselessとdangerousになっている.

つまり, 数学の公式を学んでも, それを十分に理解しなければ, 使えない. また思い付いた名案も, その基礎が分っていないと, 迷案になるというようなことだ.

この場合は英訳が分り易い.