10
$\begingroup$

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?

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – Ghoster
    Aug 22, 2023 at 4:28
  • $\begingroup$ Thank you very much. You are right. $\endgroup$
    – Laurenso
    Aug 22, 2023 at 10:04

3 Answers 3

7
$\begingroup$
  • 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]]

enter image description here

{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];
adjmatrix = 
  SparseArray[{i_, 
      j_} /; (pts[[i]] - pts[[j]]) . (pts[[i]] - pts[[j]]) == 
      length^2 -> 1, {Length@pts, Length@pts}];
adjgraph = AdjacencyGraph[adjmatrix];
subgraphs = FindIsomorphicSubgraph[adjgraph, CompleteGraph[4], All];
tetrahedrons = pts[[VertexList[#]]] & /@ subgraphs;
Graphics3D[ConvexHullRegion /@ tetrahedrons, Boxed -> False]
adjgraph

enter image description here

enter image description here

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
$\endgroup$
3
  • $\begingroup$ It is interresting if we draw one regular tetrahedron by using this method $\endgroup$ Aug 23, 2023 at 0:11
  • $\begingroup$ Wouldn't it be better instead to show the AbsoluteTiming of the whole function? not just the linetetrahedrons = pts[[VertexList[#]]] & /@ subgraphs; // AbsoluteTiming? $\endgroup$
    – ydd
    Aug 23, 2023 at 7:52
  • $\begingroup$ 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. $\endgroup$
    – kirma
    Aug 25, 2023 at 5:55
7
$\begingroup$

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.

enter image description here


Old answer:

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.

$\endgroup$
4
  • 1
    $\begingroup$ (+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} $\endgroup$
    – cvgmt
    Aug 22, 2023 at 11:49
  • 1
    $\begingroup$ @kirma Your code need DeleteCases when I try Sphere[{0, 0, 0}, 9 Sqrt[3]]. The result must be smaller. There are 34 tetrahedrons. $\endgroup$
    – Laurenso
    Aug 22, 2023 at 12:27
  • $\begingroup$ @cvgmt Added a duplicate-removal step in the end. $\endgroup$
    – kirma
    Aug 22, 2023 at 12:30
  • $\begingroup$ Very clean and very fast, impressive! $\endgroup$
    – ydd
    Aug 23, 2023 at 15:13
6
$\begingroup$

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;
 
 (*side length from cvgmt's answer*)
 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

An improvement on Version 1, added some comments

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;
 
 (*side length from cvgmt's answer*)
 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}*)

enter image description here

r = 5 Sqrt[3]

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

(*{0.028064, Null}*)

enter image description here


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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.