NO_WAIT

主にプログラミング

デュードニーのThe Haberdasher's PuzzleをJavaScript(D3)とSageを利用して作る

The Haberdasher's Puzzle という名のパズルの存在を知り、JavaScriptで再現してみたくなりました。 このパズルはHenry Ernest Dudeneyというイギリスのパズル作家による有名なパズルだそうで、 正三角形を4分割し、分割されたパーツを配置し直すことで正方形ができてしまうというなかなか巧妙なパズルです。 "The Haberdasher's Puzzle"や「デュードニー分割」などで検索すれば図はいくらでも出てくるのですが、 やはり実際に動かしてみたいところです。 考案したDudeney氏はマホガニーの板で実際に動くものを作って実演したそうですが、 ここではプログラミングによって動くものを作ってみました。

正三角形の断片をドラッグすることで黒丸で表された蝶番を中心に回転します。 蝶番をドラッグすると図形を変形させずに移動できます。 このパズルの面白いところは、断片をバラバラにせず蝶番で1本につないだまま、正方形と正三角形の間の変形が可能な点です。

1.f:id:shinaisan:20140906212845p:plain 2.f:id:shinaisan:20140906212850p:plain 3.f:id:shinaisan:20140906212859p:plain 4.f:id:shinaisan:20140906212904p:plain

三角形の分割位置の計算は下記のようになります。かなり面倒です。 定規とコンパスでがんばって作図してみるのも面白いかもしれません。

分割位置の座標の計算(Sage使用)

正三角形\(ABC\)を分割していきます。ここでは点\(A\)を原点とし、下図中各点の座標を点\(A\)からの相対位置として求めます。 少しわかりにくいですが、x軸(水平方向)は右向きを正、y軸(垂直方向)は下向きを正としています。

f:id:shinaisan:20140115204228p:plain

まず正三角形の各辺の長さを\(l\)とします。

var("l")
AB = BC = CA = l

辺\(AB\)、辺\(BC\)の中点をそれぞれ\(D\)、\(E\)とします。

AD = AB/2
BE = BC/2

以下コード中では、点\(A\)から点\(X\)の相対位置をvAXと表すことにします。 点\(B\)、\(C\)、\(D\)、\(E\)の座標をそれぞれ ベクトルvAB, vAC, vAD, vAEで表現します:

AE = AB*sin(pi/3)
vAB = vector([-BE, AE])
vAC = vector([BE, AE])
vAD = vAB/2
vAE = vector([0, AE])

線分\(AE\)を延長して\(AF\)を作ります。このとき\(EF = BE\)となるようにします。

EF = BE
AF = AE + EF
vAF = vector([0, AF])

線分\(AF\)の中点を\(G\)とします。

AG = AF/2
GE = AE - AG
vAG = vector([0, AG])

次に点\(G\)を中心に半径\(AG\)の円を描きます。 この円と点\(E\)から点\(B\)へ延ばした半直線との交点を\(H\)とします。 この\(H\)の座標vAHのx成分を変数xとおいて円の方程式eqを立てます。

var("x y") # yは後で使います。
eq = (x**2 + GE**2 == AG**2)

これをsolvexに関して解き、負の解を選択します。点\(B\)がx座標で負の位置にあるためです。

EH = - next(s.rhs() for s in solve(eq, x) if s.rhs().substitute(l = 1) < 0)
vAH = vector([-EH, AE])

点\(H\)の座標は

\[ \left(-\frac{\sqrt{\sqrt{3}}\,l}{2},\,\frac{\sqrt{3} l}{2} \right) \]

となります。

さらに点\(E\)を中心に半径\(EH\)の円を描きます。 この円と辺\(AC\)との交点を\(J\)とします。 \(J\)の座標vAJを変数x yとおいて今度は円の方程式を含む連立方程式を解きます。 次のeq1は線分\(AC\)と\(AJ\)が平行であるという条件を表し、 eq2は円の方程式を表します。

eq1 = (x*vAC[1] == y*vAC[0])
eq2 = (x**2 + (y - AE)**2 == EH**2)

solvexに関して解きます。

S = solve([eq1, eq2], x, y)
vAJ = next(vector([sx.rhs(), sy.rhs()]) for [sx, sy] in S if sx.rhs().substitute(l = 1) < 1/2)

点\(J\)の座標は

\[ \left(-\frac{1}{8} \, \sqrt{4 \, \sqrt{3} - 3} l + \frac{3}{8} \, l,\,-\frac{1}{8} \, {\left(\sqrt{4 \, \sqrt{3} - 3} \sqrt{3} - 3 \, \sqrt{3}\right)} l\right) \]

となります。

続いて、線分\(JC\)上に\(JK = BE\)となるように点\(K\)を取ります。

vJK = vector([BE*cos(pi/3), BE*sin(pi/3)])
vAK = vAJ + vJK

点\(K\)の座標は

\[ \left(-\frac{1}{8} \, \sqrt{4 \, \sqrt{3} - 3} l + \frac{5}{8} \, l,\,-\frac{1}{8} \, {\left(\sqrt{4 \, \sqrt{3} - 3} \sqrt{3} - 3 \, \sqrt{3}\right)} l + \frac{1}{4} \, \sqrt{3} l\right) \]

となります。

点\(D\)から線分\(JE\)に下ろした垂線の足を\(L\)とします。 \(L\)の座標vALを変数x yとおいて連立方程式を立てます。 以下、eq1は線分\(JE\)と\(DL\)の直交性を内積(dot_product)0として表した方程式で、 eq2は線分\(JE\)と\(JL\)が平行であるという条件を表した方程式です。

vJE = vAE - vAJ
eq1 = ((vector([x, y]) - vAD).dot_product(vJE) == 0)
eq2 = ((x - vAJ[0])*vJE[1] == vJE[0]*(y - vAJ[1]))

連立させてx, yに関して解きます。

S = solve([eq1, eq2], x, y)
vAL = vector([S[0][0].rhs(), S[0][1].rhs()])

点\(L\)の座標は

\[ \left(\frac{1}{16} \, {\left(\sqrt{4 \, \sqrt{3} - 3} \sqrt{3} + \sqrt{3} - 4\right)} l,\,-\frac{1}{16} \, {\left(\sqrt{4 \, \sqrt{3} - 3} - 4 \, \sqrt{3} - 3\right)} l\right) \]

となります。

最後に、点\(K\)から線分\(JE\)に下ろした垂線の足を\(M\)とします。 \(M\)の座標vAMを変数x yとおいて\(L\)と同様に求めます。

eq1 = ((vector([x, y]) - vAK).dot_product(vJE) == 0)
eq2 = ((x - vAJ[0])*vJE[1] == vJE[0]*(y - vAJ[1]))
S = solve([eq1, eq2], x, y)
vAM = vector([S[0][0].rhs(), S[0][1].rhs()])

点\(M\)の座標は

\[ \left(-\frac{1}{16} \, {\left({\left(\sqrt{3} + 2\right)} \sqrt{4 \, \sqrt{3} - 3} + \sqrt{3} - 10\right)} l,\,-\frac{1}{16} \, {\left({\left(2 \, \sqrt{3} - 1\right)} \sqrt{4 \, \sqrt{3} - 3} - 10 \, \sqrt{3} + 3\right)} l\right) \]

となります。

プログラムで使用するための数値解はlに具体的な値を代入(substitute)して n()呼び出しで得ます。n()numerical_approx()でも可です。

for v in [["AB", vAB], ["AC", vAC], ["AD", vAD], ["AE", vAE], ["AF", vAF], ["AG", vAG],
          ["AH", vAH], ["AJ", vAJ], ["AK", vAK], ["AL", vAL], ["AM", vAM]]:
    n = v[1].substitute(l = 300).n()
    print "%s:\t{x: %s, y: %s}" % (v[0], n[0], n[1])

上のJavaScriptプログラムではl = 300とした以下の結果をそのまま用いています:

AB: {x: -150.000000000000, y: 259.807621135332}
AC: {x: 150.000000000000, y: 259.807621135332}
AD: {x: -75.0000000000000, y: 129.903810567666}
AE: {x: 0.000000000000000, y: 259.807621135332}
AF: {x: 0.000000000000000, y: 409.807621135332}
AG: {x: 0.000000000000000, y: 204.903810567666}
AH: {x: -197.411101942874, y: 259.807621135332}
AJ: {x: 38.1761425074362, y: 66.1230184598694}
AK: {x: 113.176142507436, y: 196.026829027535}
AL: {x: 21.8423013377311, y: 148.991881821384}
AM: {x: 16.3338411697051, y: 176.938757773817}

TODO

プログラムでは蝶番の部分をかなり無理に実装していますが、 物理シミュレーションエンジンで、Jointをサポートしているものがあればより簡単に書けそうです。 Box2Dならb2RevoluteJointDefを用いることでシミュレートできるそうなので試してみたいです。

参考