What is the exact value of the radius in the Six Disks Problem?

Solution 1:

The minimal polynomial for $r(6)$ is $$\color{red}{7841367r^{18}-33449976r^{16}+62607492r^{14}-63156942r^{12}+41451480r^{10}-19376280r^8+5156603r^6-746832r^4+54016r^2+3072}$$ I obtained this result in a marathon coding session across Chinese New Year, going well into the night of 15 February and continuing into the afternoon of 16 February.

By symmetry, we can restrict ourselves to one half of the configuration, which I will take to be the upper half-plane. Then consider the smaller discs' centres and their points of intersection (both among themselves and with the outer unit disc, centred on the origin), and the following diagram is obtained:

The black dots are the key points in the problem; we will derive expressions for their coordinates later. Each red line segment is a radius of some small disc, so all of them have the same length, which I will call $r$. The red dot is the centre of the disc in the top-left corner of the diagram, and the two segments incident to it form a diameter.

Besides $r$, we define two other variables $a$ and $b$, then set the coordinates of the two key points on the $x$-axis as $A=(-a,0)$ and $B=(b,0)$ as shown in the diagram ($a,b,r$ are positive). The coordinates of the other key points can then be derived as follows (I will use $1$ to refer to the point with coordinates $(x_1,y_1)$ and so on, while $O$ assumes its usual meaning of the origin):

• The law of cosines gives $\cos\angle1OA=\frac{1+a^2-r^2}{2a}$. Because this is also the supplement of $1$'s argument, we have $$x_1=-\frac{1+a^2-r^2}{2a}$$ $y_1$ is then straightforward: $$y_1=\sqrt{1-x_1^2}$$ We can do the same thing with $\angle2OB$ to get $$x_2=\frac{1+b^2-r^2}{2b}\qquad y_2=\sqrt{1-x_2^2}$$
• $3$ lies on the perpendicular bisector of $A$ and $B$: $$x_3=\frac{b-a}2$$ If the altitude is dropped from $3$ to the $x$-axis, with $F$ as the foot, then $\angle3FB$ is right, so the Pythagorean theorem gives $$y_3=\sqrt{r^2-\left(\frac{a+b}2\right)^2}$$
• The coordinates of $4$ can now be derived as simply as adding two vectors (in this case, $\vec{B3}$ and $\vec{B2}$): $$x_4=x_2+x_3-b\qquad y_4=y_2+y_3$$
• $\angle BO5=\angle BO1-\angle5O1$, where it is trivial to compute $\sin\angle5O1=2r\sqrt{1-r^2}$ and $\cos\angle5O1=1-2r^2$. We already worked out $\sin\angle BO1$ and $\cos\angle BO1$ as $y_1$ and $x_1$ respectively, so applying angle addition formulae gives $$x_5=\cos\angle BO5=x_1(1-2r^2)+y_1(2r\sqrt{1-r^2})\\ y_5=\sin\angle BO5=y_1(1-2r^2)-x_1(2r\sqrt{1-r^2})$$

For the given configuration to be valid, the distance between $4$ and $5$ must be $r$. This can be expressed as a constraint: $$g(a,b,r)=(x_5-x_4)^2+(y_5-y_4)^2-r^2=0$$

Express the problem of finding $r(6)$ as an optimisation problem: $$\text{minimise }f(a,b,r)=r\text{ subject to }g(a,b,r)=0$$ Now try to solve this using Lagrange multipliers. The Lagrangian is $$f-\lambda g=r-\lambda g(a,b,r)$$ and all its first partial derivatives must be zero, giving a system of four equations: $$-\lambda\frac{\partial g}{\partial a}=-\lambda\frac{\partial g}{\partial b}=1-\lambda\frac{\partial g}{\partial r}=-g=0$$ However, if $\lambda=0$ then the third equation reduces into the absurdity $1=0$. Thus $\frac{\partial g}{\partial a}=\frac{\partial g}{\partial b}=0$ and combining these with the fourth equation yields a system of three equations in three unknowns: $$\frac{\partial g}{\partial a}=\frac{\partial g}{\partial b}=g=0$$ It is this last system that we will solve.

Before we proceed with solving for the minimal polynomial of $r(6)$, we need to convert $g$ into a polynomial form, even though this introduces extra solutions. This is done by removing square roots. If $$g=p+m\sqrt q=0$$ where $p,m,q$ are arbitrary expressions, the root over $q$ can be removed by bringing $m\sqrt q$ to the RHS, squaring and moving back: $$g=p^2-m^2q=0$$ All the original solutions of $g$ are retained in this lifting process, but care must be taken to keep the degrees small. Lifting with $q=-a^4+2a^2r^2+2a^2-r^4+2r^2-1$, then $q=-b^4+2b^2r^2+2b^2-r^4+2r^2-1$, already allows several small polynomials to factor out, which I removed by comparing with numerical approximations for the optimal $r,a,b$. I then performed one last lift with $q=\sqrt{-r^2+1}\sqrt{-a^2-2ab-b^2+4r^2}$, whereupon $g$ was an irreducible polynomial of total degree $22$ with $286$ monomials. The degree of $r$ in $g$ is always even, so I substituted $s=r^2$, lowering the total degree to $18$.

The Lagrangian system was finally solved in Singular via resultant computations: $a$ was eliminated from each of the pairs $\left(g,\frac{\partial g}{\partial a}\right)$ and $\left(g,\frac{\partial g}{\partial b}\right)$ and I factored out irrelevant factors as before. The two relevant polynomials in $s$ and $b$ were then combined into a third pair and $b$ was eliminated, giving a final univariate polynomial in $s$ of degree $616$. The gold at the end of the rainbow – $s$'s true minimal polynomial – is of degree $9$, and replacing back $s=r^2$ gives the degree-$18$ minimal polynomial of $r(6)$ that I wrote at the very top of this answer: $$7841367r^{18}-33449976r^{16}+62607492r^{14}-63156942r^{12}+41451480r^{10}-19376280r^8+5156603r^6-746832r^4+54016r^2+3072$$

By eliminating $s$ instead of $b$ in the resultant of $c_1$ and $c_2$, $b$'s minimal polynomial can be obtained: $$7841367b^{18}-7610085b^{16}-71679924b^{14}-45430974b^{12}+58085154b^{10}+26617188b^8-12289505b^6-5546269b^4+2814464b^2-327680$$ By eliminating $b$ instead of $a$ when computing the first two resultants, $a$'s minimal polynomial can be obtained: $$7841367a^{18}-90811125a^{16}+344734812a^{14}-376618374a^{12}-189418782a^{10}+373328980a^8+71968935a^6-130579181a^4+27056736a^2-800000$$ $r$ is indexed as A299695 in the OEIS.

3 May 2021: Three years after this monumental computation I had learned enough of Gröbner bases and Singular's libraries to get a cleaner derivation of the minimal polynomial for $r(6)$. Instead of iteratively constructing points from a radius variable and two base points, I added more variables and eliminated some of them at once to get the same reduced $g$ as described above faster. A few aspects of the last stage were also sped up, such that the new combined code takes only $2.5$ minutes on the Python portion and $3$ minutes on the Singular portion.

This is the new model:

The constraints are $1A,B2,45,15$ with the vertex labels referring to the 2018 picture above; eliminating $t,u,v$ via the Gröbner-based procedure elim gives $g$ after one division by $r$ (this $r$ is equivalent to the $s$ of the old derivation). A resultant computation is still used for eliminating a variable (this time $b$) from $\left(g,\frac{\partial g}{\partial a}\right)$ and $\left(g,\frac{\partial g}{\partial b}\right)$, and getting $a$ and $b$'s minimal polynomials, but the final elimination for $r$ uses elim.