Rubicキューブのシミュレータ

Rubicキューブのどの面をどう回せといわれても, 手でやっていれば, 手順が複雑になってくるとつい間違えてしまうのは誰にも経験があろう.

その点, シミュレータはプログラムをしっかり書けば何回でも正しく処理できる.

Winning WaysにはElena Conwayが好んだパターンを作る手順が書いてある. どんなパターンか興味があったが, シミュレータを使えばそのパターンを作るのは簡単であった.

Aは4 Windows. Bは6 Windows. CはChequers. DはHarlequin. EはStripey. FはZigzag. Gは4 Crosses. Hは6 Crossesというらしい.

手順は省略する. Winning Ways, Vol. 4, p. 876参照.

Rubicキューブのシミュレータ

情報処理学会誌にRubicキューブの話を書いたのは2005年7月号であった. その話題は島内先生の名著(島内剛一:ルービック・キューブ免許皆伝, 日本評論社(1981))の置換をHaskellで記述することであった. (この本はルービック・キューブと数学パズル, 日本評論社(2008)として再刊された.)

最近Winning Waysを見ていたら, Rubik's Hungarian Cubeなる話題があり, いろいろ書いてあるが結構ややこしい. 我が家にもRubicキューブはあったが, 最近は見当たらないので, 私のことだからPostScript(とProcessing)でシミュレータを書いた.

それを使ってWinnig Waysの方法を試みているのだが, なかなかうまくはいかない.

今回はとりあえずシミュレータの話だ.

学会誌にあった絵だが, 下の図でAは魔方体を横からみた図で, この方向を南(S)とする. 上はトップ(T)とする. 右にわずかにみえる面は東(E)である. 当然みえない右の方に西(W), 向こう側に北(N), 下にボトム(B)の面がある.

BはTの面を回転しているところで, 時計廻りをt(あるいはt+), 反時計廻りをt-という. CはSの面の回転で, s(s+)とs-になる. DはE面の回転e(e+)とe-である. (どの図もちょうど45度回ったところで, 回転方向は分からない.)

今回のシミュレータはカラーで, 上の図では見えなかった3面が右につなげて描いてある.

Bの図では左半分の上段の色が右回転した(t)後であることを示している. C, Dはsとeの結果である.

この下の図も(図の記号は違うが)学会誌の図で, 島内本にある基本操作
E) 単純3角形 9a
F) 隣辺向き替え24a
G) 巡回 28a
H) ねじり 33a
を示す. そのシミュレータ版が次だ. 54の面には次のように0から53の番号が振ってあり, リストに各面の色が入っている. 回転によりその色を入れ替えている. 例えば回転tについては

((0 2 8 6) (1 5 7 3) (15 24 33 42) (12 21 30 39) (9 18 27 36))

のような置換のリストがある. つまり面0は2へ行き, 2は8へ行き, 8は6へ行き, 6は0へ行くを示す. t-はこのリストを逆に読む. まだ前途遼遠だが, 出来たところまでをまずアップロードしよう.

MacMahonのSuperdomino

MacManonの3色四角形Superdominoにはジグソーパズル風の変形がある. つまり黒は黒と接するが, 白は灰と接するのである. それでも解はある.

白と灰の間の線を凹凸にするとジグソーパズルになるのである. そのジグソーdominoを前回のブログの3色四角形superdominoの図に対応させて描くと次のようになる.


The Winning Waysにはそれを利用したJohn Conwayの1968年のクリスマスカードの絵がある. (ウェブを探したらやっとここに見つかった.)

とりあえず制約条件をすこし修正してプログラムを走らせると, 次の解がまず得られた. この境界をジグソーに対応させたのがその下の図だ.


このジグソーパズルが通常のそれと違うのは, 絵にまったく頼らずに完成することである. それも多くの解(余詰め?)が存在する. ユニーク解にするためにはやはり絵を描くことになろう. それでCheshire Catの絵に上の解を重ねてみたのが下の図だ.


1年ほど前, IBMカードを作ったcraftroboを使って切り出せば, ジグソーパズルが作れるはずだが, こんなに優しいジグソーパズルは作るまでもあるまい.

MacMahonのSuperdomino

The Winning WaysにMacManonのSuperdominoという話題がある. 正n角形で辺と中心で出来る二等辺三角形をm色に塗り分けたものである. m色n角形superdominoという.

4色三角形superdominoは24種, 3色四角形superdominoも24種類ある.



5色五角形superdominoは, 各色を1回ずつ使い裏返しは同じとみると12種だ. この12種で辺の両側が同じ色になるように正十二面体に貼ることができるかがクイズである.

五角形は置いておき, 三角形では辺の長さが三角の辺の2倍の正六角形の中に, 周囲を同一色にし, 内側の相対する辺が同じ色になるように詰めるのがクイズ. 四角形では4×6の長方形で同様にするというのがクイズである.

四角形の, 周囲を黒にした解の一例は以下だ.



十年ほど前, 情報処理学会誌にプログラム・プロムナードの連載があり, 私はそこで「計算機用ジグソーパズル」という記事を書いた. それ以外にも探索問題のプログラムはなんども書いたからちょっとやってみたが, 最初の解が出てくるまでは存外大変であった.

プログラムを書くには24個の箱とその辺を下の図のように番号をつける.



制約としては箱0の辺0の色は2, 箱0の辺3の色は2, 箱1の辺0の色は2, 箱2の辺3の色は箱0の辺1の色, ...のようになる.

そして4色四角形のsuperdominoのプールから, 必要に応じて回転しながら, 制約に合うものを探すことになる.

しかし闇雲に初めても最初の解もなかなか現れない. ある週末, プログラムを走らせたまま帰宅したら, 次の週の初め, 夛くの解が出ていたから, 方針としては良かったのだが, これではあまりにも手抜きであったのでなるべく準備をしてからやりなおすことにした.

左上の箱から番号順に詰めていくとすると, まず辺0がg₀, 辺3がg₃, 辺1と辺2はどうでもよいdomino の集合を*g₀??g₃とする. ?は0,1,2だから要素が9個の集合が出来る. 箱5,11,17,23用には辺1が2だから要素数3の*g₀2?g₃. 箱18から22用には辺2が2だから要素数3の*g₀?2g₃. 箱18用には辺1と2が2だから要素数1のg₀22g₃を用意する.

*2222は(((2 2 2 2) (23 0)))
*02?2は(((0 2 0 2) (13 0)) ((0 2 1 2) (15 0)) ((0 2 2 2) (17 0)))
で, 要素は((u r d l) (i j))の形である. u, r, d, lは上右下左の色(0,1,2), iはdomino番号(0〜23), jは90度の回転した数(0〜3)である.

プログラムの中心は次のようだ.

(define (test n ps is js)
 (define (up ps) (caddr (list-ref ps 5)))
 (define (left ps) (cadar ps))
 (define (try *????)
  (for-each (lambda (x)
   (let ((p (car x)) (i (caadr x)) (j (cadadr x)))
    (if (not (member i is))
     (test (+ n 1) (cons p ps) (cons i is) (cons j js)))))
      *????))
 (case n
  ((0) (try *2??2))
  ((1 2 3 4) (case (left ps) ((0) (try *2??0))
                             ((1) (try *2??1))
                             ((2) (try *2??2))))
  ((5) (case (left ps) ((0) (try *22?0))
                       ((1) (try *22?1))
                       ((2) (try *22?2))))
  ((6 12) (case (up ps) ((0) (try *0??2))
                        ((1) (try *1??2))
            ((2) (try *2??2))))

あとは同様なのでしばらく省略.

  ((24) (begin (display ps) (newline) (display is) (newline)
   (display js)) (newline))))

testの引数はpsがそれまでのdominoの列, isはi, jsはjの列である. (test 0 '() '() '())で起動する.

プログラムをしばらく走らせると解が次々と現れる. とても全解探索する気分にならないが, Martin GardnerのNew Revised Edition Mathematical Diversionsの193ページに解の総数は12261と書いてあった.

1964年の初め, Stanford大学の計算センターで, Gary Feldmanが全解を求めるプログラムを書いた. Algolで書いたプログラムはB5000計算機で40時間かかったそうだ. 結果は8ページの"Documentation of the MacMahon Squares Problem," a Stanford Artificial Intelligence Project Memo No. 12で, 1964年1月16日に計算センターから刊行された.

なおこの話題によると2009年11月21日と29日のブログ 彩色立方体もMacMahonの考えたものらしい.

生後n日目

最近 還暦, 古希, 喜寿, 傘寿などの祝事で兄弟があつまる事が少くない. そういう折に私が「生れて何万何千日たった」と発言したので, まわりが似たような計算をしたがった.

私自身はJulian Date(DJ)が計算できる例の個人用電卓を携帯しているので, 「今日の日付を入力, JDに変換, 生れた日付を入力, JDに変換, 引き算」で簡単に得られる. もっともこれでは生れた日が第0日になり, 生れた日を第1日にする常識的な勘定と違うが, 満n日といえば良いかもしれない.

しかし一方, 「ちょうどn日になるのはいつかなぁ」にはJDだけでは片付かない. そこで誕生日とnを入力すると, 生後n日目になる日付を計算するウェブアプリを作ってみようと考えた.

アルゴリズムは私のブログ, 2011年5月23日のunix timeにあるのを使えばよい. 現在生きている人が対象なら, Julian Calendarは考えなくてよい. Gregorian Calendarが昔まで続いていたとするFixed Date(Rata Die)というのがあるから, 年月日からFixed Dateへの変換とその逆変換があれば, 出来たようなものである.

私はhtmlのウェブページは何回も書いたことがあるが, 計算するにはJavaScriptでも使うことになるであろう. O'reillyのJavaScriptの本は持っている. その例をみながらやって見ることにした.

<html>
<head>
<title>nth day after birth</title>
<script language="JavaScript">

この辺まではお作法のうちだ. JavaScriptでは整数の除算が見当たらなかったので, 整数除算の関数を作る.

function quotient(a,b){
var r=a%b;
var q=(a-r)/b;
return q;}

つぎは私のカレンダーのプログラムによく登場する前月までの日数の和のリストで, 平年用と閏年用を用意する.

var mon0=[0,31,59,90,120,151,181,212,243,273,304,334];
var mon1=[0,31,60,91,121,152,182,213,244,274,305,335];

次は年月日からFixed Dateを計算する関数rdだ.

function rd(y,m,d){
var a=(y-1)*365;
var b=quotient(y+3,4);
var c=quotient(y,400)-quotient(y,100);
var f=(((y%100)==0)?((y%400)==0):((y%4)==0))?mon1[m-1]:
  mon0[m-1];
return a+b+c+f+d-1;}

ご覧のように簡単だ.

反対にFixed Dateから年月日を求める関数grだ. ブログにあったアルゴリズムはfloor関数を使うが, JavaScriptにはないらしいので, 1で整数除算をすることにした.

function gr(rd){
var d0=rd+366-1;
var y400=quotient(d0,146097);
var d1=d0%146097;
var y100=quotient(d1/36524.25,1);
var d2=d1-146097+quotient(36524.25*(4-y100),1);
var y4=(y100==0)?quotient(d2,1461):
  (24-quotient((36523-d2)/1460.96,1));
var d3=(y100==0)?(d2%1461):(d2-quotient(y4*1460.96,1));
var y1=((y100>0)&&(y4==0))?quotient(d3,365):
  quotient(d3/365.25,1);
var d4=((y100>0)&&(y4==0))?d3%365:
  (d3-1461+quotient(365.25*(4-y1),1));
var y=y400*400+y100*100+y4*4+y1;
var leap=((y%100)==0)?((y%400)==0):((y%4)==0);
var e=leap?mon1:mon0;
var mon=0;
for(i=0;e[i]<=d4;i++){mon=mon+1;}
var d=d4-e[mon-1]+1;
document.write(y.toString(10)+" "+mon.toString(10)+" "
  +d.toString(10));
document.write("<br>");}

これでテストしてみるとうまく行く. 結果は関数の中でdocument.writeする.

ところでy,m,d,nの入力はどうするか. 読み込んだものは文字列らしい. 解説書で使えそうなものを探すと, データ間はカンマで区切り, カンマの位置をindexOfで探し, substringを取れば, それぞれのデータが得られることが分った. またそれだけだと文字列であるが, 0を引くと整数になるという奥の手も知った. こうして書いたのが下の<body>部分である.

</script>
</head>
<body>
<script language="JavaScript">
document.write("<br>");
var s,i,y,m,d,n;
s = prompt("Type in your birthyear,birthmonth,birthday and n
  as `2000,1,1,1000' without space", "");
i = s.indexOf(',');
y = s.substring(0,i);
s = s.substring(i+1,s.length);
i = s.indexOf(',');
m = s.substring(0,i);
s = s.substring(i+1,s.length);
i = s.indexOf(',');
d = s.substring(0,i);
n = s.substring(i+1,s.length);
document.write("y = " + y);document.write("<br>");
document.write("m = " + m);document.write("<br>");
document.write("d = " + d);document.write("<br>");
document.write("n = " + n);document.write("<br>");
y=y-0;
m=m-0;
d=d-0;
n=n-0;
var r=rd(y,m,d);
document.write("Your "+ n.toString(10) +"th day is ");
gr(r+n);document.write("<br>");
</script>
</body>
</html>

さて上のプログラムを呼出すと下のような画面が現れる.





例えばAlan Turingは1912年6月23日生まれ. その滿1万5千日目をみるには

1912,6,23,15000

と入力すると

y = 1912
m = 6
d = 23
n = 15000
Your 15000th day is 1953 7 18

と得られる. Turingは1954年6月7日に他界したから, 1万5千日とすこししか生きていなかったわけである.

Piの1000桁

IEEE Annals of the History of Computing July-September 2012に, 1949年のLabor Dayの休暇にENIACで円周率を2035計算したという記事があった.

円周率(Ludolph's number)といえばWilliam Shanksが1873年に707桁計算し, 自分の墓に刻んたという話が有名であったが, その値が528桁目から下は違っていて, それを見付けたのはENIACの計算によるのかと思っていたが, ウェブで調べると見つけたのは同時代だが, 電動計算機だったらしい.

今から半世紀前, 日本でもあちこちの大学や研究所で計算機が作られはじめ, デモ用に自然対数の底e(Napier's number)を1000桁計算するプログラムが作られた. 私たちもパラメトロン計算機PC-1でプログラムを作ったのは 2012年5月9日のこのブログに書いた通りだ.

ENIACの記事によれば, Machinの式で計算しそうだ. つまり

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

tan^-1(1/5)=1/5-1/(3·5³)+1/(5·5⁵)-1/(7·5⁷)+...

tan^-1(1/239)=1/239-1/(3·239³)+1/(5·239⁵)-1/(7·239⁷)+...

見て分るようにこの式はeの式

e=1+1/1!+1/2!+1/3!+...

より面倒なので, πに手を出す人はあまりいなかった. しかしtan^-1は交代級数なので, 計算は楽な筈である. 1000桁の計算ならtan^-1(1/5)は1/(2n+1·5²ⁿ⁺¹)が1/10¹⁰⁰⁰より小さくなる辺りまで計算すればよい.

手始めに100桁でやってみる.

(define c (expt 10 100))
(do ((n 1 (+ n 2))) ((< (/ c n (expt 5 n)) 1) n)) => 141
(do ((n 1 (+ n 2))) ((< (/ c n (expt 239 n)) 1) n)) => 43

これで必要な項数は判ったので, 105桁のbignumで計算する. aは tan^-1(1/5),　bはtan^-1(1/239)で, 最後にaの16倍からbの4倍を引き, 最後の5桁を捨てる. opはループで毎回0と1に切り替わりそれに従って次の項を足したり引いたりする.

(define c (expt 10 105))
(define op 0) (define a 0)
 (do ((n 1 (+ n 2))) ((= n 143))
  (set! a ((if (= op 0) + -) a
   (quotient (quotient c n) (expt 5 n))))
  (set! op (- 1 op)))
(define op 0) (define b 0)
(do ((n 1 (+ n 2))) ((= n 45))
 (set! b ((if (= op 0) + -) b
  (quotient (quotient c n) (expt 239 n))))
 (set! op (- 1 op)))
(quotient (- (* a 16) (* b 4)) 100000)

で結果は

314159265358979323846264338327950288419716939937510
 58209749445923078164062862089986280348253421170679

πの100桁, 1000桁の値はすぐに見付かる. 私は城, 牧之内「計算機械」にあったのを知っているので書棚から取り出し同じであることを確認した.

1000桁も同様に計算したが, 書くまでもあるまい.

多面体描画道楽

このタイトルのブログを最後にアップロードしたのは今年の9月9日で, 大二十面体の6枚の面に色づけした.

その頃から試みたかったのは, ある方向からの光線で陰影をつけることであった.

その前に見える面のすべてに色づけした図を描いてみるとこのようになった. 色の選び方には自信がないが.



これが出来たので, いよいよ陰影の図に取りかかる. 光の来る方向と各面の法線とのスカラー積を計算し, +1は光線に向いているから白に, -1は黒にすればよいのではと考えて次のような図が出来た.



しかし灰色が似ていてやはり面の識別は難しいというのが感想である. 枠となる正二十面体の頂点0,2,3のところだけ取り出したのがこの図である. 緑の正五角形の中が凹み, そこに5/2角錐が立っているのが見てとれればよいが, これもなかなか難しい.



何回か前のブログに書いたように, やはり難解な星形多面体である.