2021年6月7日月曜日

Piの1000桁

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

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

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

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

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

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

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

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

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

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

...

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

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

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

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

2021年5月25日火曜日

Piの1000桁

1970年頃のことだ. HITAC 10, FACOM R, MACC 7など, 1語が16ビットの4K 語のミニコンが多く使われていた. これらのミニコンは, 正面に押しボタ ンのようなスイッチが並んでいて, 電源を入れた後, スイッチを操作し, 数語の「初期入力プログラム」を手動で入力する. プログラム全体が入力 し終わると初期入力プログラムが動き出し, 紙テープのプログラムを読み 込めた.

ところがこのメーカー製プログラムの手の入力が面倒で, 誰もがもっと短 いプログラムに出来ないか工夫を凝らし, 徐々に短いプログラムが出現し た. これが一番短いと思っていると, 他から更に短いのが出来たという話 が来る. するとまた知恵を絞り, その語数のプログラムが作れるのであっ た.

つまり「何語で出来る」という情報が重要で, それを聞くとその語数のプ ログラムに向って注力できる.

最近ツイッターを眺めていたら, Automatic Computer@nt_a_c さんの

Automatic Computer@nt_a_c·3月18日
PC-1で円周率、400桁は出せた。もう数日早くできればちょうど良かったのだが。後はどのくらいコードを削れるか。

Automatic Computer@nt_a_c·3月24日
円周率、できた。1000桁、印刷と合わせて3回読み込み、合わせて40分くらい。https://automaticcomputer.github.io/PC1onUnity/pi....

というツイートがあった. 60年も前の計算機 PC-1でプログラムを書く奇 特な人がまだいるかと驚く一方, どうすればあの容量のメモリーでπが1000桁計算 できるか私も考え始めた.

πの1000桁については, このブログでも2012年11月5日に書いている. 簡単に計算するには やはりMachinの式でやることになろう.

π/4 = 4 tan-1(1/5)- tan-1(1/239)

tan-1(1/5)=1/5-1/(3·53)+1/(5·55) -1/(7·57)+...

tan-1(1/239)=1/239-1/(3·2393)+1/(5·239 5)-1/(7·2397)+...

PC-1は1語が18ビット. 計算は2語をつなげた長語36ビットを使う. 十進10 桁で1長語を使うから, 1000桁計算には100長語が必要で, eの1000桁では 300番地から500番地の100長語を作業場所に使っている. 512短語のPC-1で は100長語はどうしても2組しかとれないから,それで上の計算をしなけれ ばならない.

普通に考えると16×tan-1(1/5)を計算するには, atanを置く場所tsを0にし, 多倍長の16/5を作る. これをssとする. それを1で割り(その商をqsとする), qsをtsに足す. ssを25で割りssにし, そのssを3でわりqsとし, tsからqsを引く. ...

これをやるとtan-1(1/5)を計算するだけで100長語が3組必要になりすでに 破綻する.

tsはtan-1(1/5)とtan-1(1/239)で共用出来る. ssは 1/5が済んだ後1/239用に使えばよい.

残る手段はqsを作ってから, 下の桁からtsに足したりtsから引いたりするのを なんとかすることだ.

我々が小学校以来習った加減算は下の桁から始めるものであった. それに対して 上級になって習った算盤では上の桁から始めた. 小学生のころ, 上の桁から加減算 が出来るのが不思議であったが, そのうち繰上げをローカルに処理しながら下へ進んで いるからだと納得した. さらに一度繰上げが通過した桁に, その後から再び繰上げ が来ることはないことも証明できた. 基数が10から2^35と大きくなると, 繰上げが何桁も伝搬する確率は小さくなるが, 無視できない. い

要するにqsの結果はいらない. 各々の桁の商はqsに足すか引くかし, 剰余は次の 桁の除算に使うと作業場所は要らなくなり, なんとかなりそうだ.

そこでSchemeを使ってシミュレートしたのが次のプログラムである.

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に足すか引くかする.
(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))))
この出力の始めの方は次のようだ.

1415926535 8979323846 2643383279 5028841971
6939937510 5820974944 5923078164 628620899
8628034825 3421170679 8214808651 3282306647
938446095 5058223172 5359408128 4811174502
8410270193 8521105559 6446229489 5493038196
...
10^10進法の各桁の後に空白があるのは, 十進の一番上が0なのが分るためだ.

これで一応は計算出来るが, 考えてみるに, ssは徐々に小くなっていく. その途中にずーっと先頭から除算を始めるのはやはり無駄である. つまり 除算を開始する位置を制御すれば計算は早く終わる. 下はそのように修正 したプログラムで, 変数bが除算開始位置を覚えている.

(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 " "))
上のプログラムをPC-1のプログラムにしたのが次だ. 作業番地tsはeの1000桁と同様, 300番地以降を使う. プログラムの最初にあるインタールードでクリアしておく.

ssは0番地から200番地までを使う. つまりイニシアルオーダーR0は壊す. 計算だけで記憶場所を使うので, 結果tsの出力は別のプログラムにする.

まだコメントが少ないが, 徐々に追加するつもりだ.

      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,      z2の除算ループ. 最初の剰余を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.
このプログラムは私のPC-1シミュレータに追加してある. シミュレータを開き, R0 を読み込み, Place Tapeをクリック. メニューの下方にあるpi1000cmpを クリックすると計算のプログラムが設定される. そこでFree Runをクリック; さらにInit Startをクリックするとプログラムが走り出す.

計算は私のMacBook Pro環境では約3分半で終り, シミュレータは停止する.

結果の出力には作業場所 を同じにしてあるので, eの1000桁の出力プログラムを多少修正して使う.

計算プログラムがR0領域を破壊したから, 出力プログラムの読込み前にもう一度 R0をロードする. 後はPlace TapeからPi1000prtを選び, そのプログラムを 読込み起動すればπの1000桁が出力される.

久し振りにPC-1のプログラムを書いたが, インタプリタで必要な情報を取り出し たり, Schemeのプログラムで予備テストが出来たり, 便利な世の中である.

2021年3月29日月曜日

PC-1シミュレータ

このブログにパラメトロン計算機 PC-1のシミュレータの話を書いたのはそろそろ10年前の ことだ. そのシミュレータはProcessingで書いたものだが, その後シミュレータの 置いてある計算機が引退したりで, 長いあいだ動かなかった. ブログを見た方には申し訳 ないことをした.

ずーっとシミュレータの書き直しを考えていたが果せず, 最近やっとJavaScriptで書き, どうやら無事に動くようなので, 公開することにした.

シミュレータは https://www.iijlab.net/~ew/pc1sim2/pc1sim.html にあるので試してみて欲しい.

シミュレータの見掛けは以前のままである. もともとPC-1が出来たころ, 近隣の 計算機と競争していた, 自然対数の底eの値を1000桁計算するプログラムを シミュレートするのが目的だったので, そのプログラムに特化してはあるが, 他の プログラムも少数ながら動く.

シミュレータを開くと次のような画面になる. 右手に3x3の押しボタンがある.



最初にやることは初期入力ルーチンR0の読込みで, まず左下のPlace R0をクリックする. R0のテープが現れたら, 次にその上のInit Loadをクリック. これでR0が読み込まれる.

今度は中下のPlace e1000をクリック. e1000桁のテープが現れる. そこで 中央のFree Runを押し, 続けて中央上のInit Startをクリック. これでe1000桁の テープがR0で読み込まれ, 計算が始まる. 程なくして結果の印字が始まり, 計算が 終了する.

他に用意してあるテープは,

factorize32 231より少し小さい整数の素因数分解
eratosthenes エラトステネスの篩
v2test フーリエ変換
va1test 素数を法とする行列式の計算

で, これらを走らせるには右下のPlace Tapeをクリックする. そうすると右上に プログラム名が表示されるから, それをクリックする. テープが現れたら中央の Free Run, その上のInit Startをクリック. なかにはテープ読込み中に停止する ものもあるが, その時は中央のFree Runと, 右上のRestartをクリックする.

2020年12月10日木曜日

無限クイーン

TAOCPの演習問題7221-38に, 無限クイーン列<qn>の効率的計算法を考案せよというのがある.

面白い解答なので, その説明をしたい.

二進の配列 a,b,cを使う. cは正と負の添字がある.
G1. [初期化.] r←0, s←1, t←0, n←0とする. ( 1≤k≤nについて, qkの計算が済んだ.)
G2. [qn≤nを調べる.] (この時点で1≤k<sについて ak=1でa2=0; また-r<k≤tについてck=1 でc-r=ct+1=0; 各ベクターにはn個の1がある.) n←n+1, k←sとする.
G3. [見つけた?] ak=bk+n=ck-n=0なら G5へ進む. そうでないなら, k←k+1とし, k≤n-rならこのステップを繰り返す.
G4. [qn>nを作る.] t←t+1; qn←n+t; ak←b2n+t←ct←1とし, G2へ戻る.
G5. [qn≤nを作る.] qn←k; ak←bk+n←ck-n←1とする. k=sなら, as=0になるまでs←s+1を繰り返す. k=n-rなら c-r=0になるまでr←r+1を繰り返す. G2へ戻る.

このプログラムを忠実に辿ると, 大体こういうことをやっていると判明する.

TAOCPのプログラムでは, チェス盤の行も列も1から始めるが, 私はやはり0から 始めたい. 前回の無限クイーンの最後にあるように, 解の列は

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

のように始まる. その始めの方の置きかたを見ると(下の図 前回の図では, 行は下に 延びていたのに, 今回は逆に上に延びるように描いて申し訳けない.)

まず列n=0は行0にクイーンが置ける. それが左下の黒丸だ. そうするとどの行にはもう置けない ことを示す配列aの0を1に, 行-列の対角線に置けないことを示す配列cの0-0=0を1にする. 図ではそれを黒丸のすぐ右の青小丸とすぐ右上の赤小丸で示す. 配列bも0+0=0が 1になるがそれは図には示さない.

続いてn=1にクイーンを置くのだが, 青線が示すようにこの行には置けず, 赤の対角線が 示すようにこの斜線にも置けない. よって赤斜線の上のq=2に置くことになる. それが (1,2)の黒丸で, その右の青小丸と右上の赤小丸で配列の設定も分る.

n=2の列は, 行0は青線でだめ, 行2は赤線でだめ, その間の行1は, すぐ左上にクイーン がいてだめ(これは配列bで分る). 行4も(1,2)からの斜線のためだめで, 結局 行4に置く.

このようにして置いていくのだが, 青線はこれより下はすべて詰っている; 2本の赤線は この間はすべて詰っていることを示す. したがって行0の青線は, 列3の行1に置くと 行2まで詰ったので, ここからは行2から横に進む.

赤の斜線も同様に出来ている. 従ってある列でクイーンが置けるか調べるには, 青線路 の上から始め, 下の赤線の下までである. ここまでで配列a, b, cを調べ, 置けない 時は, 上の赤線の上に置くことになる.

この方針で私流に書いたプログラムが次だ.
(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)	  
変数sは青線を, rは下の赤線を, tは右の赤線を示す. これは最初のTAOCPのプログラムと 同じである. 変数kの使い方も同様.

下の図は, 上のと同じだが, nの範囲を広げた. また配列を調べた位置を白丸で示した. 白丸の数が少ないことに注目して欲しい. TAOCPのプログラムは判定の方法が違うので, 白丸の数はこれよりかなり多い.


2020年11月22日日曜日

mex関数

アルゴリズムを読んでいると, mexに出会うことがある. mexとは minimum excludantのことで, 最小の排除されたもの, つまり 引数の集合の最小不在者である. 例えば
(mex '()) => 0
(mex '(1 2)) => 0
(mex '(0 2)) => 1
(mex '(0 1)) => 2
前回のブログ, 無限クイーンで説明したbit stringのある処理系では,
(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))))
nsが(0 1 1 2 3 5 8 13)の時, (mex ns)では. nsのmaxが13だからlは15. 従って1が15個のbit stringが出来, 最初の(右端の)ビットのビット番号は0, 最後の(左端の) ビット番号は14だ.
#*111111111111111
次に位置0, 1, 1, 2,..., 13のビットをクリアして0にすると
#*101111011010000
になる. そして0から15の範囲で最初の1の位置を探すと4であり,
(mex ns) => 4
上のmex関数はこのように動作するが, bit stringの最初の1を探すような飛道具を 使っていて, 快しとしない向きには, 以下の二進木(二分木)を使うアルゴリズムはどうだろうか.

二進木の節点は, 左部分木, その節点の情報, 右部分木 で構成される. 今の アルゴリズムの場合, その節点の情報は範囲[l, r]で, mexをとる集合nsに l<=k<=r要素があることを示す. 二進木に含まれる節点の範囲には重なりがないだけ ではなく, 隣同士の範囲の間には隙間があるようになっている.

隣同士の範囲の間に新しい要素が現れてその隙間を埋めると, 二つの範囲は合体して 一つの範囲になる.
			     
ns=(12 1 3 3 4 4 5 11 6 6 10 14 15 13 4 7)
として, 二進木の出来かたを下の図に示す.

最初の12が来た時に, 範囲[12,12]が出来る(図の右上.) その右の12は, この時の 入力が12であったことを示す.

次に1が来ると, 1は既にある範囲[12,12]の中でもないし, 両端の続きでもなく離れ ているので, 別の範囲[1,1]を作り, [12,12]の節点の左部分木とする(先程の図の下.)

その次の3は, 新しい範囲[3,3]になって, [12,12]の左の[1,1]の右に繋る.

次の3は[3,3]の範囲にあるから, 何もしない.

4は[3,3]の右に接しているから, 範囲[3,3]が[3,4]になる. しかし, その右部分木 の最も左に5がないので, 合体は生じず, 仕事はここで終わる.

図の1列目が終り, 中央の列の上に戻り, 5が来ると, 範囲[3,4]が[3,5]になる.

11が来ると, [12,12]が[11,12]になり, ここでも合体しない.

更に進んで右の列になり, 14,15,13と来ると, 13は[10,12]の 右を13にし, 右の節点に[14,15]があるのでこの二つが合体して, [10,13]の範囲が出来る.

こうしてnsを全て処理した二進木が右下の図になる.

この二進木を中順序で辿り, 最初に来る節点の範囲[1,1]を見る.

この左が0でなければ, mexは0, 0ならば右プラス1がmexである.

二進木の操作になれていれば, プログラムを書くのは造作ない.

なお, nsが空なら, 木も空で, mexは0である.

P.S. Knuth先生へのメイルでこの実装にも触れたら返事に
It is attractive. I don't recall seeing it before.
と書いてあった.

2020年11月10日火曜日

無限クイーン

TAOCPのV4F5に,
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,...
という数列があり, これは「無限クイーン」の辞書式順で最小の解とあった(演習問題 7221-42). これを計算してみたというのが今回の話題だ.

情報科学の標準問題のひとつに8クイーンがある. チェス盤の上に, 飛車と角行を併せて動くクイーンを8個, 互いに当らないように置く方法を求める問題だ. その8を無限大にしたのが「無限クイーン」である. つまり縦横とも半無限に広いチェス盤に, 孫悟空の分身の術のように無限に現れるクイーンを置くのである.

釈迦に説法かも知れないが, 順序としてまず8クイーンの解法を復習しよう. 次の図を見て欲しい.

8×8のチェス盤の列にも行にも0から7の番号をつける. そして左端の列0からクイーンの 置けるところを探す. 既に置いてあるクイーンと当らないようにということは, 置く場合, その行には他のクイーンはない; その場所の右上がりの筋には他のクイーンはいない. その場所の右下がりの筋には他のクイーンはいない. という条件が揃えばそこに 新しいクイーンが置けて, さらに右の列の検討が始まる.

その列のどの行にも上の条件の場所がなければ, 1列戻って, その列のクイーンを更に下に 置けるかどうか調べる. こうして8個置くことが出来れば解である. 一つの解が見付かっても, まだ他の解を探すのが全解探索だ.

上のような戦術で8クイーンを探すプログラムを示すと次のようになる.

この行, この右上がりの筋, この右下がりの筋にまだクイーンが置けることを示す配列 as, bs, csを用意する. 真理値が入る. 最初はどこにでも置けるから, 初期値は真だ.

列の添字をm, 行の添字をnとすると, 右上がりの筋はm+n=一定, 右下がりはm-n=一定で, mもnも0から7の値をとるから, m+nの値は0から14, m-nの値は-7から7である. 配列の添字は0以上なので, 右下がりはm-n+7の場所に置く.

各列の置き方を調べる関数が col m で, その下請けの各行の置き方を調べるのが row n である.

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

このプログラムはバックトラックを繰り返す. 最初の解が得られるまでのバックトラック の様子を見ると下のようだ.



そうして得られた92通りの解答は次の通り.

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

ところで, MIT Schemeにはbit stringというデータの型があり, Eratosthenesの篩の ような真理値の配列に便利である. 今回使う関数を説明すると
(make-bit-string 8 #t)
各ビットが1の8ビットのbit stringを生成する. #fなら0になる.
(bit-string-ref bs k)
bit string bsのk番目のビットを調べる. 0なら#f, 1なら#tが返る.
(bit-string-set! as k)
bit string asのk番目のビットを1にする.
(bit-string-clear! as k)
bit string asのk番目のビットを0にする. bit stringにはある範囲で 最初の1の位置を探す飛道具がある.
(bit-string-find-next-set-bit as k 8)
bit string asの[k .. 8)の間で, 最初の1であるビットの位置を返す. 1がなければ#fを返す.
(bit-string-append as bs)
bit string asとbsを接続した新しいbit stringを返す. as内のビット番号は不変. 新しい bsの部分の番号は前の番号にasの長さを足したものになる.

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にする.

(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)
これで様子が分ったから, 無限クイーンに取り掛かろう. 問題はcのbit stringで, 0がbit stringの途中にあるので, オフセットを足すことだ. 無限クイーンでいいのは, 列の中で次に置ける位置はどんどん先まで探せばよく, バックトラック しないことだ. 従って#tに戻す(a! k)のような関数はない. 無限といっても無限まで 計算することはないので, 列の最大番号はきめることにした. mmaxの値である. するとas, bs, csなどの長さも余裕をもって決められる.
(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)
実行結果を下に示す.


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


TAOCPにはqn+O(1)か n/φ+O(1)と書いてあったので, 図にしてみた. 赤線が黄金比の傾きの斜線だ.



上のプログラムは, as, bs, などが最初から固定長であった. これを最初はある程度に とり, 必要に応じて長くするように書いたのが次のプログラムだ. 上のと違う部分だけ 書いておく.

 (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)からは同じ

2020年10月19日月曜日

五角ミノ十二面体

英語で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についてはまた述べる.