The website packomania gives packing results for a number of different situations, to high precision, but not exact . In this notebook, the results from NMinimize showing which circles touch the square and which circles touch each other are used to find an exact solution for packing three circles with radii Range[3]^(-1/2) into a square .

Approximate Solution

In[]:=
packNM[radii_]:=Block{n=Length[radii],vars,x,y,insideCons,nooverlapCons},​​​​vars=Flatten@Join[Table[{x[i],y[i]},{i,n}],{s}];​​​​insideCons=Flatten@Table[{x[i]>=radii[[i]],x[i]<=s-radii[[i]],y[i]>=radii[[i]],y[i]<=s-radii[[i]]},{i,n}];​​​​nooverlapCons=Flatten@Table[
2
(x[i]-x[j])
+
2
(y[i]-y[j])
>=
2
(radii[[i]]+radii[[j]])
,{i,n-1},{j,i+1,n}];​​​​NMinimize[{s,insideCons,nooverlapCons,y[1]>y[2],y[1]>y[3]},vars]​​​​
In[]:=
plt[sln_,radii_]:=​​Block[{n=Length[radii],cntrs},​​Show[​​Graphics@Line[{{0,0},{0,s},{s,s},{s,0},{0,0}}]/.sln[[2]],​​Graphics[Point[Table[{x[i],y[i]},{i,n}]/.sln[[2]]]],​​Graphics@Table[Circle[{x[i],y[i]},radii[[i]]],{i,n}]/.sln[[2]],​​ImageSize->Small,Frame->True​​]​​]
In[]:=
slnNM[3]=packNM
-1/2
Range[3]

Out[]=
{3.01851,{x[1]1.80001,y[1]2.01851,x[2]0.707107,y[2]0.707107,x[3]2.44116,y[3]0.57735,s3.01851}}
In[]:=
pltslnNM[3],
-1/2
Range[3]

Out[]=

Exact Solution

In[]:=
radii=
-1/2
Range[3]
;
In[]:=
y[1]=s-radii[[1]];
In[]:=
x[2]=y[2]=radii[[2]];
In[]:=
x[3]=s-radii[[3]];
In[]:=
y[3]=radii[[3]];
In[]:=
distSqr[pt1_List,pt2_List]:=(pt1-pt2).(pt1-pt2)
test model with NSolve
In[]:=
AbsoluteTiming[slnNSolve=NSolve[​​{​​distSqr[{x[1],y[1]},{x[2],y[2]}]==
2
(radii[[1]]+radii[[2]])
,​​distSqr[{x[1],y[1]},{x[3],y[3]}]==
2
(radii[[1]]+radii[[3]])
,s>3​​},{s,x[1]}]]
Out[]=
{0.0759369,{{s3.01851,x[1]1.80001}}}
use Solve
In[]:=
AbsoluteTiming[slnSolve=Solve[​​{​​distSqr[{x[1],y[1]},{x[2],y[2]}]==
2
(radii[[1]]+radii[[2]])
,​​distSqr[{x[1],y[1]},{x[3],y[3]}]==
2
(radii[[1]]+radii[[3]])
,s>3​​},{s,x[1]}]]
Out[]=
0.164861,s
3.02
…
,x[1]
1.80
…

In[]:=
N[slnSolve,20]
Out[]=
{{s3.0185131882907988217,x[1.0000000000000000000]1.8000051285999809214}}
In[]:=
plt[Join[s/.slnSolve,slnSolve],radii]
Out[]=
In[]:=
slnSolve//InputForm
Out[]//InputForm=
{{s -> Root[{-2 + #1^2 & , -3 + #2^2 & , 49 + 20*#1*#2 - 240*#3 - 216*#1*#3 -
176*#2*#3 - 96*#1*#2*#3 + 660*#3^2 + 288*#1*#3^2 + 192*#2*#3^2 +
168*#1*#2*#3^2 - 288*#3^3 - 288*#1*#3^3 - 192*#2*#3^3 + 180*#3^4 & },
{2, 2, 4}], x[1] -> Root[{-2 + #1^2 & , -3 + #2^2 & ,
49 + 20*#1*#2 - 240*#3 - 216*#1*#3 - 176*#2*#3 - 96*#1*#2*#3 + 660*#3^2 +
288*#1*#3^2 + 192*#2*#3^2 + 168*#1*#2*#3^2 - 288*#3^3 - 288*#1*#3^3 -
192*#2*#3^3 + 180*#3^4 & , -11*#1 - 6*#2 + 106*#3 + 48*#1*#3 + 32*#2*#3 +
12*#1*#2*#3 - 96*#3^2 - 66*#1*#3^2 - 44*#2*#3^2 + 60*#3^3 + 8*#4 & },
{2, 2, 4, 1}]}}
The ToRadicals result is too long to show here.

Compare to Packomania Result

Packomania scales the circles rather than the square, so a comparison is made between the inverse of the square side and the packomania result.
In[]:=
1
s
-0.331288928562270941669800278213/.slnSolve
Out[]=
{1.8834862646763788×
-14
10
}

CITE THIS NOTEBOOK

Exact packing of 3 circles with radii Range[3]^(-1/2) in a square​
by Frank Kampas​
Wolfram Community, STAFF PICKS, June 1, 2026
https://community.wolfram.com/groups/-/m/t/3724903