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})
r = 5 Sqrt[3]
(*same old code as Version 1, just different r*)
({0.028064, Null})
