# How can I select four points on a sphere to make a regular tetrahedron so that its coordinates are integer numbers?

I want to select four points lie on the sphere (x-1)^2 + (y-3)^2 + (z-5)^2 = (5* Sqrt[3])^2 so that its coordinates are integer numbers to make a regular tetrahedron. I tried

ClearAll[a, b, r, c];
a = 1;
b = 3;
c = 5;
r = 5* Sqrt[3]; ss =
Subsets[{x, y, z} /.
Solve[{(x - a)^2 + (y - b)^2 + (z - c)^2 == r^2}, {x, y, z},
Integers], {4}];
list = Select[
ss, ( EuclideanDistance[#[[1]], #[[2]]] ==
EuclideanDistance[#[[1]], #[[3]]] ==
EuclideanDistance[#[[1]], #[[4]]] ==
EuclideanDistance[#[[2]], #[[4]]] ==
EuclideanDistance[#[[2]], #[[3]]] ==
EuclideanDistance[#[[3]], #[[4]]] &&
Det[{#[[1]] - #[[2]], #[[1]] - #[[3]], #[[1]] - #[[4]]}] !=
0 &) ]


About ten minutes, I can not get the result? How can I get the result?

• It appears that you are trying to find a regular tetrahedron. In general, a tetrahedron doesn’t have edges of equal length, so you have probably found many tetrahedra with integral vertices on that sphere. Commented Aug 22, 2023 at 4:28
• Thank you very much. You are right. Commented Aug 22, 2023 at 10:04

• No such regular tetrahedron,even no triangle of such regular tetrahedron.
• If there are exist such regular tetrahedron, the length of its side should be r/(Sqrt[3/2]/2). We found many lines satisfies this condition, but no triangle,so there are no regular tetrahedron.
Clear["Global*"];
{a,b,c}={1,3,5};
r = 15;
pairs = Subsets[{x, y, z} /.
Solve[{(x - a)^2 + (y - b)^2 + (z - c)^2 == r^2}, {x, y, z},
Integers], {2}];
test2[{p1_, p2_}] := (p1 - p2) . (p1 - p2) == (r/(Sqrt[3/2]/2))^2;
lines = Pick[pairs, test2 /@ pairs];
Graphics3D[Line[lines]]


{a,b,c}={1,3,5};
r = 15;
triples =
Subsets[{x, y, z} /.
Solve[{(x - a)^2 + (y - b)^2 + (z - c)^2 == r^2}, {x, y, z},
Integers], {3}];
test3[{p1_, p2_,
p3_}] := (p1 - p2) . (p1 - p2) == (p2 - p3) . (p2 - p3) == (p3 -
p1) . (p3 - p1) == (r/(Sqrt[3/2]/2))^2;
triangles = Pick[triples, test3 /@ triples]


{}

## Edit : Graph theory method

For r = 33*Sqrt[3];, there are 26*5=130 subgraphs means that there are 26*5=130 tetrahedron in the original question.

Clear["Global*"];
{a, b, c} = {1, 3, 5};
r = 33*Sqrt[3];
length = r/(Sqrt[3/2]/2);
pts = SolveValues[{(x - a)^2 + (y - b)^2 + (z - c)^2 == r^2}, {x, y,
z}, Integers];
SparseArray[{i_,
j_} /; (pts[[i]] - pts[[j]]) . (pts[[i]] - pts[[j]]) ==
length^2 -> 1, {Length@pts, Length@pts}];
tetrahedrons = pts[[VertexList[#]]] & /@ subgraphs;
Graphics3D[ConvexHullRegion /@ tetrahedrons, Boxed -> False]

adjgraph


## Edit

• Now we do not assume that we know in advance that the tetrahedron has sides of length r/(Sqrt[3/2]/2).
{a, b, c} = {1, 3, 5};
r = 33 Sqrt[3];
pts = {x, y, z} /.
Solve[{(x - a)^2 + (y - b)^2 + (z - c)^2 == r^2}, {x, y, z},
Integers];
pairs = Table[
If[i > j, {i, j} -> (pts[[i]] - pts[[j]]) . (pts[[i]] - pts[[j]]),
Nothing], {i, Length@pts}, {j, Length@pts}];
groups = GatherBy[Flatten[pairs, 1], Last];
graphs = Graph[pts, #, VertexCoordinates -> pts] & /@ Keys@groups;
subgraphs =
FindIsomorphicSubgraph[#, CompleteGraph[4], All] & /@ graphs /. {} ->
Nothing // First

Graph[VertexList@#, EdgeList@#, VertexCoordinates -> VertexList@#,
EdgeStyle -> RandomColor[]] & /@ subgraphs

• It is interresting if we draw one regular tetrahedron by using this method Commented Aug 23, 2023 at 0:11
• Wouldn't it be better instead to show the AbsoluteTiming of the whole function? not just the linetetrahedrons = pts[[VertexList[#]]] & /@ subgraphs; // AbsoluteTiming?
– ydd
Commented Aug 23, 2023 at 7:52
• FindIsomorphicSubgraph is a correct, but not a necessary method here. The graph contains at most 4-cliques, and thus more efficient FindClique can be used in this case. Commented Aug 25, 2023 at 5:55

EDIT: Improved distance computation efficiency a lot by using DistanceMatrix:

With[{r = 33 Sqrt[3]},
With[
(* Find integer-valued points on the sphere. *)
{pts = {x, y, z} /.
Solve[Element[{x, y, z}, Sphere[{1, 3, 5}, r]], Integers]},
(* Compute distance matrix for the points. *)
DistanceMatrix[pts, DistanceFunction -> SquaredEuclideanDistance] //
(* Find indices of edges which may be a part of a regular
tetrahedron in this circumsphere by length. *)
Position[(8/3) r^2] //
(* Construct graph edges with coordinates as vertices. *)
Map[UndirectedEdge @@ pts[[#]] &]]] //
(* Find all 4-cliques. *)
FindClique[Graph[#], {4}, All] &


$$r = 33\sqrt{3}$$ evaluates now in 24 milliseconds on average on my laptop.

It should be noted that FindClique returns maximal cliques. For exact 4-clique subgraphs one should normally use FindIsomorphicSubgraph but that is not necessary here since these graphs don't have larger than 4-cliques as one can have only four points at equal distances from each other in the three-dimensional Euclidean space.

EDIT: Converted the search into a graph clique problem:

With[{r = 33 Sqrt[3]},
Select[
Subsets[
(* Find integer-valued points on the sphere. *)
{x, y, z} /.
Solve[Element[{x, y, z}, Sphere[{1, 3, 5}, r]], Integers], {2}],
(* Select edges which may be a part of a regular tetrahedron in
this circumsphere by length. *)
SquaredEuclideanDistance @@ # == (8/3) r^2 &] //
(* Construct a graph and find all 4-cliques. *)
MapApply@UndirectedEdge // FindClique[#, {4}, All] &]


This is not the fastest solution but it's quite concise, and for $$r = 33\sqrt{3}$$ evaluates in 0.56 seconds on my laptop.

EDIT: A new, much more scalable approach here, the original solution under it.

With[{r = 33 Sqrt[3]},
Select[
Subsets[
(* Find integer-valued points on the sphere. *)
{x, y, z} /.
Solve[Element[{x, y, z}, Sphere[{1, 3, 5}, r]],
Integers], {2}],
(* Select edges which may be a part of a regular tetrahedron in
this circumsphere by length. *)
SquaredEuclideanDistance @@ # == Evaluate[(8/3) r^2] &]] //
Select[
(* Select from pairs of such edges so that *)
Flatten[#, 1] & /@ Subsets[#, {2}],
(* resulting vertices form a regular tetrahedron. *)
Equal @@ (SquaredEuclideanDistance @@@ Subsets[#, {2}]) &] & //
(* Delete redundant reorderings. *)
DeleteDuplicatesBy[Sort]


... this computes 130 solutions for $$r = 33\sqrt{3}$$ in less than 10 seconds on my laptop while my old code below requires almost four minutes to do it.

Solution for $$r = 99\sqrt{3}$$ below.

Brute-force search with a small tweak (find valid equilateral triangles first) is probably the most efficient solution when the sphere radius is 15 (original question):

With[
{(* Integer-valued points on the sphere. *)
onsphere =
{x, y, z} /.
Solve[Element[{x, y, z}, Sphere[{1, 3, 5}, 15]], Integers],
(* Function which returns True if all points in the list
are pairwise equidistant. *)
eqdist = Equal @@ (SquaredEuclideanDistance @@@ Subsets[#, {2}]) &},
Table[
(* Return valid tetrahedra. *)
If[eqdist[Append[i, j]],
Append[i, j],
Nothing],
(* Equilateral triangles on the sphere. *)
{i, Select[Subsets[onsphere, {3}], eqdist]},
(* All individual points to test along with the triangles. *)
{j, onsphere}]] //
(* Delete redundant reorderings. *)
Flatten[#, 1] & // DeleteDuplicatesBy[Sort]

(* {} *)


Result is empty which means there are no such solutions.

With sphere radius of $$5\sqrt{3}$$ 14 solutions are found in 0.6 seconds, matching @cvgmt's result.

• (+1) Faster. And there are essential 14 different tetrahedrons in the list. {GatherBy[list, Map@*Sort] // Length, Gather[list, RegionEqual[ConvexHullRegion[#1], ConvexHullRegion[#2]] &] // Length} Commented Aug 22, 2023 at 11:49
• @kirma Your code need DeleteCases when I try Sphere[{0, 0, 0}, 9 Sqrt[3]]. The result must be smaller. There are 34 tetrahedrons. Commented Aug 22, 2023 at 12:27
• @cvgmt Added a duplicate-removal step in the end. Commented Aug 22, 2023 at 12:30
• Very clean and very fast, impressive!
– ydd
Commented Aug 23, 2023 at 15:13

I used PowersRepresentations as Daniel Lichtblau did here as it felt natural to me when I first saw this type of problem.

With Version 3, I get the 130 solutions for $$r=33~\sqrt{3}$$ in less than 0.5 seconds on my laptop:

# Version 3

Noticing that the distances matrix in Version 2 was symmetrical, I roughly cut my time in half. I used Stelio's answer here to create the distances matrix entries below the diagonal when I needed to index on it.

r =33 Sqrt[3], using symmetry of distance matrix

Clear["Global*"];
AbsoluteTiming[center = {1, 3, 5};
r = 33 Sqrt[3];
(*get positive integer coordinates on r=15 sphere*)
nnvals = PowersRepresentations[r^2, 3, 2];
permVals = Flatten[Permutations /@ nnvals, 1];

(*multiply the coords by all possible signs*)
signs = Tuples[{-1, 1}, {3}];
alltriples = Union[Flatten[Outer[Times, signs, permVals, 1], 1]];

(*shift to center of sphere*)
alltriples = # + center & /@ alltriples;

sideLength = r/(Sqrt[3/2]/2);

(*note because of the +/- and permutations, lTrips is always even*)
lTrips = Length@alltriples;
(*take only the first half of alltriples, since the distance matrix \
is symmetrical*)
upper = Take[alltriples, lTrips/2];
(*calculate distances for first half of list to all of list*)
distances = Outer[EuclideanDistance, upper, alltriples, 1];

(*get list of coords that are sideLength away from each coord*)
verticesWithoutSelf = (Flatten@Position[#, sideLength]) & /@
distances;
(*add the coord itself to its own list*)
vertices = MapIndexed[Sort@Join[#2, #1] &, verticesWithoutSelf];

(*get tetrahedron candidates,sort,and delete duplicates*)
verts4 = Flatten[Sort@Subsets[#, {4}] & /@ vertices, 1];
verts4 = DeleteDuplicates[verts4];

(*create part of distance matrix below diagonal so we can index on \
it*)
lowerDistances =
Reverse /@ (Transpose[Reverse /@ distances]) // Transpose;
fullDistances = Join[distances, lowerDistances];

(*grab indices from full distance matrix fullDistances*)
distFun[v1_, v2_] := fullDistances[[v1, v2]];
(*get distances between each vertex of the tetrahedron candidates*)
edges = Outer[distFun[#1, #2] &, #, #, 1] & /@ verts4;
(*get unique edge lengths*)
edges = Sort[DeleteDuplicates[Flatten[#]]] & /@ edges;
(*a regular tetrahedron will only have lengths sideLength and 0 \
(length to self)*)
allSameSideLength = Flatten@Position[edges, {0, sideLength}];
(*get tetrahedron vertex coordinates*)
tetraVerts = verts4[[allSameSideLength]];
tetrahedrons = Part[alltriples, #] & /@ tetraVerts;
]

(*{0.427283, Null}*)


r = 5 Sqrt[3]

(*same Version 3 code, just different r*)

(*{0.007682, Null}*)


As a small note, you don't actually have shift alltriples to the center of the sphere, you could apply this shift later to tetrahedrons and get the same thing since integer-integer = integer. This doesn't seem to make any difference in computation time however.

# Version 2

r= 33 Sqrt[3]

Clear["Global*"];
AbsoluteTiming[
center = {1, 3, 5};
r = 33 Sqrt[3];
(*get positive integer coordinates on r=15 sphere*)
nnvals = PowersRepresentations[r^2, 3, 2];
permVals = Flatten[Permutations /@ nnvals, 1];

(*multiply the coords by all possible signs*)
signs = Tuples[{-1, 1}, {3}];
alltriples = Union[Flatten[Outer[Times, signs, permVals, 1], 1]];

(*shift to center of sphere*)
alltriples = # + center & /@ alltriples;

sideLength = r/(Sqrt[3/2]/2);

(*calculate distance between all coords on the sphere*)
distances = Outer[EuclideanDistance, alltriples, alltriples, 1];

(*get list of coords that are sideLength away from each coord*)
verticesWithoutSelf = (Flatten@Position[#, sideLength]) & /@
distances;
(*add the coord itself to its own list*)
vertices = MapIndexed[Sort@Join[#2, #1] &, verticesWithoutSelf];

(*get tetrahedron candidates, sort, and delete duplicates*)
verts4 = Flatten[Sort@Subsets[#, {4}] & /@ vertices, 1];
verts4 = DeleteDuplicates[verts4];

(*I just made this function to get distances[[v1,v2]]*)
distFun[v1_, v2_] := distances[[v1, v2]];
(*get distances between each vertex of the tetrahedron candidates*)
edges = Outer[distFun[#1, #2] &, #, #, 1] & /@ verts4;
(*get unique edge lengths*)
edges = Sort[DeleteDuplicates[Flatten[#]]] & /@ edges;

(*a regular tetrahedron will only have lengths sideLength and 0 \
(length to self)*)
allSameSideLength = Flatten@Position[edges, {0, sideLength}];
(*get tetrahedron vertex coordinates*)
tetraVerts = verts4[[allSameSideLength]];
tetrahedrons = Part[alltriples, #] & /@ tetraVerts;]
Graphics3D[ConvexHullRegion /@ tetrahedrons]

(*{0.805232, Null}*)


For anyone wanting to further optimize for speed, the bottleneck is calculating the distances between all the integer coordinates:

AbsoluteTiming[
distances = Outer[EuclideanDistance, alltriples, alltriples, 1];]
(*{0.762782, Null}*)


Also, the AbsoluteTiming of distances looks like it should be proportional to Length[alltriples]^2. Part of alltriples definition is a Permutation...this is probably not possible but, if we calculated the distances for one permutation of coordinates, could we draw conclusions about other permutations?

r = 5 Sqrt[3]

(*Same Version 2 code, just different r*)

(*{0.013104, Null}*)


# Version 1

r = 33 Sqrt[3]

 Clear["Global*"];
AbsoluteTiming[

center = {1, 3, 5};
r = 33 Sqrt[3];
nnvals = PowersRepresentations[r^2, 3, 2];
permVals = Flatten[Permutations /@ nnvals, 1];
signs = Tuples[{-1, 1}, {3}];
alltriples = Union[Flatten[Outer[Times, signs, permVals, 1], 1]];
alltriples = # + center & /@ alltriples;

sideLength = r/(Sqrt[3/2]/2);

distances = Outer[EuclideanDistance, alltriples, alltriples, 1];

tetras = Flatten[Table[
candidateVertices =
Join[{i},
Flatten[Position[distances[[i]], n_ /; n == sideLength]]];
indices = Subsets[candidateVertices, {4}];
Part[alltriples, #] & /@ indices
, {i, distances // Length}], 1];
sortedTetras = DeleteDuplicates[Sort /@ tetras];
allSameSidePositions =
Flatten[Position[
Table[Tally[
DeleteCases[Flatten[Outer[EuclideanDistance, i, i, 1]],
0]][[All, 1]], {i, sortedTetras}], n_ /; n == {sideLength}]];
tetrahedrons = sortedTetras[[allSameSidePositions]];
]
Graphics3D[ConvexHullRegion /@ tetrahedrons]
(*{2.0544, Null}*)


r = 5 Sqrt[3]

(*same old code as Version 1, just different r*)

(*{0.028064, Null}*)


• (+1) fastest. But the reason for the fast is not because of the use of PowersRepresentations. Commented Aug 22, 2023 at 23:03
• @cvgmt Can you repair your code to become faster? Commented Aug 22, 2023 at 23:35
• @cvgmt I added comments and further optimized my answer if you would like to play with it.
– ydd
Commented Aug 23, 2023 at 4:05
• @ydd Thank you for your edits. You can add an answer to this question Commented Aug 23, 2023 at 4:10
• @cvgmt I think I beat it by using DistanceMatrix and FindClique`... Commented Aug 23, 2023 at 10:21