| Bytes | Lang | Time | Link |
|---|---|---|---|
| 016 | C++ GCC | 251003T124500Z | Deadcode |
| 011 | Rust | 201108T231337Z | wastl |
C++ (GCC), 12 14 16 terms
#include <stdio.h>
#include <iostream>
#include <string.h>
#include <array>
#include <vector>
#include <unordered_set>
#include <chrono>
//#define USE_GMP
#ifdef USE_GMP
#include <gmp.h>
# if GMP_NUMB_BITS != 64
# error This is hard-coded for 64-bit limbs
# endif
#endif
auto startTime = std::chrono::steady_clock::now();
typedef uint8_t TetIndex;
typedef uint8_t TetIndexFace; // lowest 2 bits are used for a face index
#ifdef USE_GMP
class Tetrahedron
{
public:
mpz_t t[4][3];
Tetrahedron(const Tetrahedron &_t) : Tetrahedron()
{
for (int p=0; p<4; p++)
for (int d=0; d<3; d++)
mpz_set(t[p][d], _t.t[p][d]);
}
Tetrahedron()
{
for (int p=0; p<4; p++)
for (int d=0; d<3; d++)
mpz_init(t[p][d]);
}
~Tetrahedron()
{
for (int p=0; p<4; p++)
for (int d=0; d<3; d++)
mpz_clear(t[p][d]);
}
};
#else
typedef __int128 Coord;
typedef std::array<Coord, 3> Coord3;
typedef std::array<Coord3, 4> Tetrahedron;
#endif
class Tet
{
void initFaces()
{
faceAttached[0] = NULL; // t[0],t[1],t[2]
faceAttached[1] = NULL; // t[0],t[1],t[3]
faceAttached[2] = NULL; // t[0],t[2],t[3]
faceAttached[3] = NULL; // t[1],t[2],t[3]
}
public:
Tetrahedron t;
Tet *faceAttached[4];
TetIndex index; // 1-based; 0=unassigned
Tet( ) : t( ) {initFaces();}
Tet(const Tetrahedron &t) : t(t) {initFaces();}
void assignIndex(TetIndex &nextIndex)
{
if (index == 0)
index = nextIndex++;
}
};
class Polytet : public std::vector<Tet>
{
public:
TetIndex nextIndex;
void resetIndexing(size_t first) // This function should not be called if the polytet hasn't yet been populated
{
for (auto t=begin(); t!=end(); ++t)
t->index = 0;
(*this)[first].index = 1;
nextIndex = 2;
}
};
// vertex indices of faces with identical chirality
static int tetrahedronFaces[4][4] =
{
{0, 1, 2, 3},
{0, 3, 1, 2},
{0, 2, 3, 1},
{1, 3, 2, 0},
};
static int tetrahedronEdges[6][2] =
{
{0, 1},
{1, 2},
{2, 0},
{0, 3},
{1, 3},
{2, 3},
};
#ifndef USE_GMP
Coord3 operator+(const Coord3 &a, const Coord3 &b)
{
Coord3 c;
c[0] = a[0] + b[0];
c[1] = a[1] + b[1];
c[2] = a[2] + b[2];
return c;
}
Coord3 operator-(const Coord3 &a, const Coord3 &b)
{
Coord3 c;
c[0] = a[0] - b[0];
c[1] = a[1] - b[1];
c[2] = a[2] - b[2];
return c;
}
Coord3 operator*(const Coord3 &a, const Coord b)
{
Coord3 c;
c[0] = a[0] * b;
c[1] = a[1] * b;
c[2] = a[2] * b;
return c;
}
Coord dot(const Coord3 &a, const Coord3 &b)
{
return a[0]*b[0] + a[1]*b[1] + a[2]*b[2];
}
#endif
class TetrahedronOverlap
{
const Tetrahedron *a, *b;
public:
void setA(const Tetrahedron &x) {a = &x;}
void setB(const Tetrahedron &x) {b = &x;}
#ifdef USE_GMP
private:
mpz_t intersectNumerator, intersectDenominator, uNumerator, vNumerator, uvDenominator, uvNumeratorSum;
mpz_t center[3], normal[3], tmp[3], p0p1[3], intersectionPoint[3], delta[3], edge1[3], edge2[3];
mpz_t triangle[3][3];
void dot(mpz_t &result, const mpz_t _a[3], const mpz_t _b[3])
{
mpz_mul (result, _a[0], _b[0]);
mpz_addmul(result, _a[1], _b[1]);
mpz_addmul(result, _a[2], _b[2]);
}
public:
TetrahedronOverlap()
{
mpz_inits(intersectNumerator, intersectDenominator, uNumerator, vNumerator, uvDenominator, uvNumeratorSum, NULL);
for (int d=0; d<3; d++)
{
mpz_init(center[d]);
mpz_init(normal[d]);
mpz_init(tmp[d]);
mpz_init(p0p1[d]);
mpz_init(intersectionPoint[d]);
mpz_init(delta[d]);
mpz_init(edge1[d]);
mpz_init(edge2[d]);
}
for (int p=0; p<3; p++)
for (int d=0; d<3; d++)
mpz_init(triangle[p][d]);
}
~TetrahedronOverlap()
{
mpz_clears(intersectNumerator, intersectDenominator, uNumerator, vNumerator, uvDenominator, uvNumeratorSum, NULL);
for (int d=0; d<3; d++)
{
mpz_clear(center[d]);
mpz_clear(normal[d]);
mpz_clear(tmp[d]);
mpz_clear(p0p1[d]);
mpz_clear(intersectionPoint[d]);
mpz_clear(delta[d]);
mpz_clear(edge1[d]);
mpz_clear(edge2[d]);
}
for (int p=0; p<3; p++)
for (int d=0; d<3; d++)
mpz_clear(triangle[p][d]);
}
bool operator()()
{
// Take advantage of the fact that the two tetrahedrons are regular and congruent, and
// just check if any edge of tetrahedron "a" intersects with any face of tetrahedron "b".
// Don't count it if only the endpoint of an edge intersects.
for (int edgeNum=0; edgeNum<6; edgeNum++)
{
const mpz_t *p0 = a->t[tetrahedronEdges[edgeNum][0]];
const mpz_t *p1 = a->t[tetrahedronEdges[edgeNum][1]];
for (int faceNum=0; faceNum<4; faceNum++)
{
const mpz_t *normalizedTetrahedron[4][3]; // first 3 points are the face, and the 4th point is for calculating the normal
for (int i=0; i<4; i++)
for (int d=0; d<3; d++)
normalizedTetrahedron[i][d] = &b->t[tetrahedronFaces[faceNum][i]][d];
// Center coordinates will be multiplied by 3 compared to original coordinates.
// Get center of face by averaging its vertices' coordinates; the
// division by 3 is implied by omitting the multiplication by 3.
for (int d=0; d<3; d++)
mpz_set(center[d], *(normalizedTetrahedron[0][d]));
for (int p=1; p<3; p++)
for (int d=0; d<3; d++)
mpz_add(center[d], center[d], *(normalizedTetrahedron[p][d]));
for (int d=0; d<3; d++)
{
mpz_neg(normal[d], center[d]);
mpz_addmul_ui(normal[d], *(normalizedTetrahedron[3][d]), 3);
}
for (int d=0; d<3; d++)
{
mpz_sub(tmp [d], *(normalizedTetrahedron[0][d]), p0[d]);
mpz_sub(p0p1[d], p1[d], p0[d]);
}
dot(intersectNumerator , normal, tmp);
dot(intersectDenominator, normal, p0p1);
int cmp = mpz_cmp_ui(intersectDenominator, 0);
if (cmp == 0)
continue; // edge is parallel to face, which we don't count as an overlap
if (cmp < 0)
{
mpz_neg(intersectNumerator , intersectNumerator );
mpz_neg(intersectDenominator, intersectDenominator);
}
if (mpz_cmp_ui(intersectNumerator, 0) <= 0 || mpz_cmp(intersectNumerator, intersectDenominator) >= 0)
continue;
// These coordinates are all multiplied by intersectDenominator
for (int d=0; d<3; d++)
{
mpz_mul(intersectionPoint[d], p0[d], intersectDenominator);
mpz_addmul(intersectionPoint[d], p0p1[d], intersectNumerator);
}
for (int i=0; i<3; i++)
for (int d=0; d<3; d++)
mpz_mul(triangle[i][d], *(normalizedTetrahedron[i][d]), intersectDenominator);
// Check if the intersection point is inside the triangle
for (int d=0; d<3; d++)
{
mpz_sub(delta[d], intersectionPoint[d], triangle[0][d]);
mpz_sub(edge1[d], triangle[1][d] , triangle[0][d]);
mpz_sub(edge2[d], triangle[2][d] , triangle[0][d]);
}
mpz_mul (uNumerator, delta[1], edge2[0]);
mpz_submul(uNumerator, delta[0], edge2[1]);
mpz_mul (vNumerator, delta[0], edge1[1]);
mpz_submul(vNumerator, delta[1], edge1[0]);
mpz_mul (uvDenominator, edge1[1], edge2[0]);
mpz_submul(uvDenominator, edge1[0], edge2[1]);
cmp = mpz_cmp_ui(uvDenominator, 0);
if (cmp == 0)
continue;
if (cmp < 0)
{
mpz_neg(uNumerator, uNumerator);
mpz_neg(vNumerator, vNumerator);
mpz_neg(uvDenominator, uvDenominator);
}
if (mpz_cmp_ui(uNumerator, 0) <= 0 || mpz_cmp_ui(vNumerator, 0) <= 0)
continue;
mpz_add(uvNumeratorSum, uNumerator, vNumerator);
if (mpz_cmp(uvNumeratorSum, uvDenominator) < 0)
return true;
}
}
return false;
}
#else
public:
bool operator()()
{
// Take advantage of the fact that the two tetrahedrons are regular and congruent, and
// just check if any edge of tetrahedron "a" intersects with any face of tetrahedron "b".
// Don't count it if only the endpoint of an edge intersects.
for (int edgeNum=0; edgeNum<6; edgeNum++)
{
Coord3 p0 = (*a)[tetrahedronEdges[edgeNum][0]];
Coord3 p1 = (*a)[tetrahedronEdges[edgeNum][1]];
for (int faceNum=0; faceNum<4; faceNum++)
{
Tetrahedron normalizedTetrahedron; // first 3 points are the face, and the 4th point is for calculating the normal
for (int i=0; i<4; i++)
normalizedTetrahedron[i] = (*b)[tetrahedronFaces[faceNum][i]];
Coord3 center = {{0, 0, 0}}; // multiplied by 3 compared to original coordinates
// Get center of face by averaging its vertices' coordinates; the
// division by 3 is implied by omitting the multiplication by 3.
for (int p=0; p<3; p++)
for (int d=0; d<3; d++)
center[d] += normalizedTetrahedron[p][d];
Coord3 normal;
for (int d=0; d<3; d++)
normal[d] = normalizedTetrahedron[3][d] * 3 - center[d];
Coord intersectNumerator = dot(normal, normalizedTetrahedron[0] - p0);
Coord intersectDenominator = dot(normal, p1 - p0);
if (intersectDenominator == 0)
continue; // edge is parallel to face, which we don't count as an overlap
if (intersectDenominator < 0)
{
intersectNumerator = -intersectNumerator;
intersectDenominator = -intersectDenominator;
}
if (intersectNumerator <= 0 || intersectNumerator >= intersectDenominator)
continue;
// These coordinates are all multiplied by intersectDenominator
Coord3 intersectionPoint = p0 * intersectDenominator + (p1 - p0) * intersectNumerator;
Coord3 triangle[3];
for (int i=0; i<3; i++)
triangle[i] = normalizedTetrahedron[i] * intersectDenominator;
// Check if the intersection point is inside the triangle
Coord3 delta = intersectionPoint - triangle[0];
Coord3 edge1 = triangle[1] - triangle[0];
Coord3 edge2 = triangle[2] - triangle[0];
Coord uNumerator = delta[1]*edge2[0] - delta[0]*edge2[1];
Coord vNumerator = delta[0]*edge1[1] - delta[1]*edge1[0];
Coord uvDenominator = edge1[1]*edge2[0] - edge1[0]*edge2[1];
if (uvDenominator == 0)
continue;
if (uvDenominator < 0)
{
uNumerator = -uNumerator;
vNumerator = -vNumerator;
uvDenominator = -uvDenominator;
}
if (uNumerator <= 0 || vNumerator <= 0)
continue;
if (uNumerator + vNumerator < uvDenominator)
return true;
}
}
return false;
}
#endif
};
void attachNewTet(Tet &t, Tet &tetToAttachTo, const int faceNum)
{
#ifdef USE_GMP
mpz_t *newVertex = t.t.t[0];
// Get center of face by averaging its vertices' coordinates.
for (int d=0; d<3; d++)
mpz_set(newVertex[d], tetToAttachTo.t.t[tetrahedronFaces[faceNum][0]][d]);
for (int p=1; p<3; p++)
for (int d=0; d<3; d++)
mpz_add(newVertex[d], newVertex[d], tetToAttachTo.t.t[tetrahedronFaces[faceNum][p]][d]);
// Finalize the new vertex
for (int d=0; d<3; d++)
{
mpz_div_ui(newVertex[d], newVertex[d], 3);
mpz_mul_ui(newVertex[d], newVertex[d], 2);
mpz_sub (newVertex[d], newVertex[d], tetToAttachTo.t.t[3 - faceNum][d]);
}
#else
Coord3 &newVertex = t.t[0];
newVertex = {{0, 0, 0}};
// Get center of face by averaging its vertices' coordinates.
for (int p=0; p<4; p++)
{
if (p == 3 - faceNum)
continue;
for (int d=0; d<3; d++)
newVertex[d] += tetToAttachTo.t[p][d];
}
// Finalize the new vertex
for (int d=0; d<3; d++)
newVertex[d] = newVertex[d]/3 * 2 - tetToAttachTo.t[3 - faceNum][d];
#endif
// Copy the other vertices
for (int p=0; p<3; p++)
{
int p1 = tetrahedronFaces[faceNum][p];
for (int d=0; d<3; d++)
#ifdef USE_GMP
mpz_set(t.t.t[1+p][d], tetToAttachTo.t.t[p1][d]);
#else
t.t[1+p][d] = tetToAttachTo.t[p1][d];
#endif
}
t.faceAttached[0] = NULL;
t.faceAttached[1] = NULL;
t.faceAttached[2] = NULL;
t.faceAttached[3] = &tetToAttachTo;
tetToAttachTo.faceAttached[faceNum] = &t;
}
// First two tetrahedrons are implied. Each element is a subsequent tetrahedron, with the value indicating where
// it's attached. The lower 2 bits indicate which face (can only have 3 different values, because at least 1 face
// will always already be attached). The remaining bits indicate which tetrahedron (which can never be zero,
// because that one is attached implicitly).
class CompressedPolytet : public std::vector<TetIndexFace>
{
public:
void append(Polytet &polytet, Tet &tetToCompress, int vertexMap[4], int faceRotation)
// indices of vertexMap[] are compressed-output vertices; elements of vertexMap[] are the original vertices of tetToCompress
{
tetToCompress.assignIndex(polytet.nextIndex);
for (int _faceNum=0; _faceNum<3; _faceNum++)
{
int rotatedFaceNum = (_faceNum + faceRotation) % 3;
int faceNum = 3 - vertexMap[3 - rotatedFaceNum];
if (!tetToCompress.faceAttached[faceNum])
continue;
Tet *attachedTet = tetToCompress.faceAttached[faceNum];
attachedTet->assignIndex(polytet.nextIndex);
push_back((((TetIndexFace)(tetToCompress.index - 1)) << 2) + _faceNum);
int attachedFace = 0;
while (attachedTet->faceAttached[attachedFace] != &tetToCompress)
attachedFace++;
int vertexMap2[4];
int rotation = 0;
while (vertexMap[tetrahedronFaces[rotatedFaceNum][rotation]] != tetrahedronFaces[faceNum][0])
rotation++;
vertexMap2[1] = tetrahedronFaces[attachedFace][0];
vertexMap2[3] = tetrahedronFaces[attachedFace][1];
vertexMap2[2] = tetrahedronFaces[attachedFace][2];
vertexMap2[0] = tetrahedronFaces[attachedFace][3];
/*if (faceNum == 3)
{
int tmp = vertexMap2[3];
vertexMap2[3] = vertexMap2[2];
vertexMap2[2] = tmp;
}*/
append(polytet, *attachedTet, vertexMap2, rotation);
}
}
};
namespace std
{
template<>
struct hash<CompressedPolytet>
{
std::size_t operator()(const CompressedPolytet &polytet) const noexcept
{
std::size_t seed = polytet.size();
for (auto i=polytet.cbegin(); i!=polytet.cend(); ++i)
seed ^= std::hash<uint32_t>{}(*i) + 0x9e3779b9 + (seed << 6) + (seed >> 2);
return seed;
}
};
};
int main(int argc, char *argv[])
{
#ifdef USE_GMP
Tetrahedron start;
for (int d=0; d<3; d++)
{
mpz_set_si(start.t[0][d], -9);
for (int p=1; p<4; p++)
mpz_set_si(start.t[p][d], 9);
}
for (int p=1; p<4; p++)
mpz_set_si(start.t[p][p-1], -9);
#else
static Tetrahedron start =
{{
{{-9,-9,-9}},
{{-9, 9, 9}},
{{ 9,-9, 9}},
{{ 9, 9,-9}}
}};
#endif
TetrahedronOverlap overlap;
auto *polytets = new std::unordered_set<CompressedPolytet>;
polytets->insert(CompressedPolytet()); // add empty vector as the starter polytet (meaning it has two tetrahedrons)
size_t prevPolytetCount = 0;
size_t blahNum = 0;
for (int tetCount=1;;)
{
auto currentTime = std::chrono::steady_clock::now();
size_t polytetCount = polytets->size();
std::cout << tetCount << ": " << polytetCount << " [" << std::chrono::duration_cast<std::chrono::milliseconds>(currentTime - startTime).count() << " ms]" << std::endl;
if (prevPolytetCount > polytetCount)
{
std::cerr << "Quit due to apparent overflow" << std::endl;
break;
}
prevPolytetCount = polytetCount;
if (++tetCount <= 2)
continue;
/*if (tetCount > 6)
break;*/
Polytet polytet;
polytet.reserve(tetCount); // Important, to ensure pointers don't change
Tet &t0 = polytet.emplace_back(start);
attachNewTet(polytet.emplace_back(), t0, 3);
auto *newPolytets = new std::unordered_set<CompressedPolytet>;
polytet.resize(tetCount);
for (auto basePolytet=polytets->cbegin(); basePolytet!=polytets->cend(); ++basePolytet)
{
int tetNumToUncompress = 2;
for (auto elementToUncompress=basePolytet->cbegin(); elementToUncompress!=basePolytet->cend(); ++elementToUncompress)
{
int faceNum = *elementToUncompress & 3;
int tetNumToAttachTo = *elementToUncompress >> 2;
Tet &tetToAttachTo = polytet[tetNumToAttachTo];
attachNewTet(polytet[tetNumToUncompress++], tetToAttachTo, faceNum);
}
if (tetNumToUncompress != tetCount - 1)
{
std::cerr << "Error! Got " << tetNumToUncompress << ", expected " << tetCount - 1 << std::endl;
exit(-1);
}
Tet &newTet = polytet[tetCount - 1];
for (int tetNumToAttachTo = 0; tetNumToAttachTo < tetCount-1; tetNumToAttachTo++)
{
Tet &tetToAttachTo = polytet[tetNumToAttachTo];
for (int faceNum=0; faceNum<3; faceNum++) // skip last face because it's always already attached
{
if (tetToAttachTo.faceAttached[faceNum])
continue;
attachNewTet(newTet, tetToAttachTo, faceNum);
// Canonicalize the rotation of this new polytet in compressed form, so that it can be compared against others
bool haveRunningLeast = false;
CompressedPolytet runningLeastPolytet;
Tet *t = &polytet[1];
for (int i=0; i<tetCount; i++)
{
int attachedFace;
int vertexMap[4];
{
Tet &singlyAttachedTet = polytet[i];
for (int j=0; j<3; j++)
if (singlyAttachedTet.faceAttached[j])
goto skipThisTet; // not a singly attached tet
t = singlyAttachedTet.faceAttached[3];
attachedFace = 0;
while (t->faceAttached[attachedFace] != &singlyAttachedTet)
attachedFace++;
static int vertexMapTable[4][4] =
{
{3, 0, 2, 1},
{2, 0, 1, 3},
{1, 0, 3, 2},
{0, 1, 2, 3},
};
memcpy(vertexMap, vertexMapTable[attachedFace], sizeof(vertexMap));
}
for (int rotationStep=0; rotationStep<3; rotationStep++)
{
polytet.resetIndexing(i);
CompressedPolytet newRotatedPolytet;
newRotatedPolytet.reserve(tetCount - 2);
newRotatedPolytet.append(polytet, *t, vertexMap, rotationStep);
// Update the running "least" rotation
if (!haveRunningLeast ||
std::lexicographical_compare(
newRotatedPolytet .begin(), newRotatedPolytet .end(),
runningLeastPolytet.begin(), runningLeastPolytet.end()))
{
haveRunningLeast = true;
runningLeastPolytet = newRotatedPolytet;
}
}
skipThisTet:;
}
if (auto [insertedItem, wasInserted] = newPolytets->emplace(runningLeastPolytet); wasInserted)
{
// Check for overlap between this newly attached tetrahedron and the existing ones,
// and defer this until after the deduplication, to save a lot of time
overlap.setA(newTet.t);
for (auto tetCheckIntersection=polytet.cbegin(); tetCheckIntersection!=polytet.cend(); ++tetCheckIntersection)
{
if (&*tetCheckIntersection == &tetToAttachTo || &*tetCheckIntersection == &newTet)
continue; // skip this check for speed (it'll always be false anyway)
overlap.setB(tetCheckIntersection->t);
if (overlap())
{
newPolytets->erase(insertedItem);
break;
}
}
}
tetToAttachTo.faceAttached[faceNum] = NULL;
}
}
polytet[1].faceAttached[0] = NULL;
polytet[1].faceAttached[1] = NULL;
polytet[1].faceAttached[2] = NULL;
}
delete polytets;
polytets = newPolytets;
for (int p=0; p<4; p++)
for (int d=0; d<3; d++)
#ifdef USE_GMP
mpz_mul_ui(start.t[p][d], start.t[p][d], 3);
#else
start[p][d] *= 3;
#endif
}
return 0;
}
This matches the output of wastl's program up to 12 terms, and does so much faster (62 65 2680 times faster on my machine).
Certain design elements were quite natural and are shared between both our programs – namely, the use of 128-bit integer coordinates, the coordinates of the starting tetrahedron, and multiplying all coordinates by 3 for each successive term. The remaining details are implemented from scratch though, without consulting wastl's program.
The fact that both versions of my program agree with wastl's up to \$A276272(12)\$ strongly suggests that they're all correct up to that point (and beyond it by using larger integers), and that the OEIS entry can finally be updated.
My previous version used affine transformations to rotate each polytet for comparison and deduplication. I wanted to keep that version up, rather than relegating it to the edit history, but the 65536 character limit prevented that.
This new version does its rotations entirely while remaining in tree representation, using Cartesian coordinates only to check for overlaps between tetrahedrons.
Compared to the previous affine transform version, it's about 40 times faster and uses 32 times less RAM (27 times less RAM than wastl's program).
It does deduplication by serializing the tree representation with every possible choice of singly-attached tetrahedron as its root node, doing 3 rotations of each, and choosing the serialization that is lexicographically least for hash-table comparison against polytets that were previously identified as unique using the same method.
As such, it's conceptually similar to the previous version – both take advantage of the fact that all polytets have some tetrahedrons that are attached on only one face, because due to a tetrahedron's dihedral angle being irrational, a polytet can never curve around and attach to one of its own faces.
So, the computationally expensive affine transformation is replaced with a tree re-rooting. This operation now being very fast, it can do a new optimization: Delaying the tetrahedron overlap/collision checking until after deduplication, so that a provisionally new polytet won't be checked for overlaps just to find that it's a duplicate.
The timings shown in the output are cumulative.
Output in USE_GMP (GNU Multiple Precision Arithmetic Library) mode:
1: 1 [0 ms]
2: 1 [0 ms]
3: 1 [0 ms]
4: 4 [0 ms]
5: 10 [0 ms]
6: 39 [2 ms]
7: 164 [13 ms]
8: 767 [75 ms]
9: 3656 [400 ms]
10: 18186 [2347 ms]
11: 91532 [13428 ms]
12: 468203 [78222 ms]
13: 2417722 [459991 ms]
14: 12595984 [2762236 ms]
15: 66068726 [16173756 ms]
16: 348619611 [98104807 ms]
In this USE_GMP version, it is no longer possible for the calculations with Cartesian coordinates to overflow. As such, it could in theory keep going up to \$A276272(63)\$ given enough time and RAM, at which point it would no longer be able to serialize polytet trees at 1 byte per tetrahedron.
There's still plenty of room for further speed-ups, such as adding multithreading, and applying various more specific optimizations. Disk storage and merge-sorting could be used to do away with the RAM limit, allowing time to become the only limit.
Also, the STL implementation of hash tables is extremely inefficient in its RAM usage. I'm pretty sure that RAM usage could be dramatically decreased just by implementing my own hash table / unordered set, and plan to do just that to reach \$A276272(17)\$ in <64 GB without using disk storage (which, in the current implementation, would require about 184 GB RAM; the calculation of \$A276272(16)\$ peaked at 34.9 GB RAM usage).
Output in basic 128-bit integer mode, which overflows starting at \$A276272(13)\$:
1: 1 [0 ms]
2: 1 [0 ms]
3: 1 [0 ms]
4: 4 [0 ms]
5: 10 [0 ms]
6: 39 [0 ms]
7: 164 [2 ms]
8: 767 [12 ms]
9: 3656 [77 ms]
10: 18186 [447 ms]
11: 91532 [2672 ms]
12: 468203 [16470 ms]
13: 2417697 [101231 ms]
14: 349124 [549295 ms]
Output from the affine transform version
This version used about 20% more RAM than wastl's program.
In USE_GMP mode:
1: 1 [0 ms]
2: 1 [0 ms]
3: 1 [1 ms]
4: 4 [2 ms]
5: 10 [12 ms]
6: 39 [50 ms]
7: 164 [270 ms]
8: 767 [1556 ms]
9: 3656 [9505 ms]
10: 18186 [59459 ms]
11: 91532 [384340 ms]
12: 468203 [2399032 ms]
13: 2417722 [15053858 ms]
14: 12595984 [101318974 ms]
This version did not use GMP to solve for the affine transform, only to apply it after solving for it. Nor did it use GMP for storing the coordinates of vertices. As such, overflow would eventually happen. It still reaches \$A276272(14)\$ correctly though, at which point it used 42.3 GB RAM (getting close to my machine's 64 GB), so I can't run it beyond that point (nor would I want to, given how much faster the new tree version is).
In basic 128-bit integer mode, which overflows starting at \$A276272(13)\$:
1: 1 [0 ms]
2: 1 [0 ms]
3: 1 [0 ms]
4: 4 [0 ms]
5: 10 [2 ms]
6: 39 [10 ms]
7: 164 [53 ms]
8: 767 [296 ms]
9: 3656 [2350 ms]
10: 18186 [15845 ms]
11: 91532 [106165 ms]
12: 468203 [683449 ms]
13: 2417721 [4291118 ms]
14: 663520 [4734026 ms]
Rust, 11 (hopefully correct) terms
use ::std::ops::{Add, Sub, Mul};
use ::std::rc::Rc;
use ::std::hash::{Hash, Hasher};
use ::std::collections::HashSet;
use ::std::iter;
use ::std::fmt::{self, Formatter, Display};
use ::std::time::Instant;
type Coord = i128;
#[derive(Copy, Clone, PartialEq, Eq, Debug)]
struct Vec3([Coord; 3]);
impl Vec3 {
fn dot(self, other: Vec3) -> Coord {
self.0.iter().zip(other.0.iter()).map(|(a, b)| a * b).sum()
}
}
impl Display for Vec3 {
fn fmt(&self, f: &mut Formatter) -> fmt::Result {
write!(f, "({}, {}, {})", self.0[0], self.0[1], self.0[2])
}
}
impl Mul<Vec3> for Coord {
type Output = Vec3;
fn mul(self, mut vec: Vec3) -> Vec3 {
for i in &mut vec.0 {
*i *= self;
}
vec
}
}
impl Add for Vec3 {
type Output = Vec3;
fn add(mut self, other: Vec3) -> Vec3 {
for (i, n) in self.0.iter_mut().enumerate() {
*n += other.0[i];
}
self
}
}
impl Sub for Vec3 {
type Output = Vec3;
fn sub(self, other: Vec3) -> Vec3 {
self + (-1 * other)
}
}
#[derive(Clone, Debug)]
struct Tetrahedron([Vec3; 4]);
impl Default for Tetrahedron {
fn default() -> Tetrahedron {
Tetrahedron([
Vec3([-1, -1, -1]),
Vec3([-1, 1, 1]),
Vec3([ 1, -1, 1]),
Vec3([ 1, 1, -1]),
])
}
}
impl Display for Tetrahedron {
fn fmt(&self, f: &mut Formatter) -> fmt::Result {
write!(f, "Tetrahedron({}, {}, {}, {})", self.0[0], self.0[1], self.0[2], self.0[3])
}
}
impl Tetrahedron {
fn collides(&self, other: &Tetrahedron) -> bool {
let mut othervecs = [3 * self.0[1], 3 * self.0[2], 3 * self.0[3]];
let sum = self.0[0] + self.0[1] + self.0[2] + self.0[3];
let mut same = 0;
for (i, &vec) in self.0.iter().enumerate() {
let sum = sum - vec;
let vec = 3 * vec;
let through = vec - sum;
for (j, &u) in other.0.iter().enumerate() {
let u = 3 * u;
if u == vec { same += 1 }
let up = (sum - u).dot(through);
for &v in &other.0[j+1..] {
let v = 3 * v;
let edge = v - u;
let ep = edge.dot(through);
if up.signum() != ep.signum() || up.abs() >= ep.abs() {
continue
}
let intersection = ep * u + up * edge;
if othervecs.iter().enumerate().all(|(i, &ov)| {
let mid = othervecs[(i+1)%3] + othervecs[(i+2)%3];
let ov = 2 * ov - mid;
(2 * intersection - ep * mid).dot(ep * ov) > 0
}) {
return true
}
}
}
if i != 3 { othervecs[i] = vec; }
}
if same == 4 { panic!("EQUAL TETRAHEDRA IN .collides()"); }
false
}
// Mirroring also scales by 3
fn mirror(&self, i: usize) -> Tetrahedron {
let mut sum = Vec3([0; 3]);
for (j, &vec) in self.0.iter().enumerate() {
if j != i { sum = sum + vec; }
}
let mut copy = self.clone();
copy.scale();
copy.0[i] = 2 * sum - copy.0[i];
copy.swap(i, 3);
copy.rotate_left(i);
copy
}
fn scale(&mut self) {
for vec in self.0.iter_mut() {
*vec = 3 * *vec;
}
}
fn rotate_left(&mut self, n: usize) {
self.0[..3].rotate_left(n)
}
fn swap(&mut self, a: usize, b: usize) {
self.0.swap(a, b)
}
fn reverse(&mut self) {
self.0.reverse()
}
}
type EndpointIter<'a> = Box<dyn Iterator<Item = Endpoint> + 'a>;
#[derive(Debug, Clone)]
struct Endpoint {
tree: Rc<TetraTree>,
hedron: Rc<Tetrahedron>,
}
impl Endpoint {
fn iter_endpoints(&self) -> EndpointIter {
let mut hedron = Rc::clone(&self.hedron);
//Rc::make_mut(&mut hedron).swap(0, 3);
Rc::make_mut(&mut hedron).reverse();
Box::new(self.clone().into_iter_directions()
.chain(self.tree.iter_endpoints(
Rc::new(TetraTree {
subtrees: [None, None, None],
hedron,
})
))
)
}
fn into_iter_directions(self) -> EndpointIter<'static> {
Box::new(iter::successors(Some(self), |this| {
let mut this = this.clone();
Rc::make_mut(&mut this.tree).rotate_left(1);
Rc::make_mut(&mut this.hedron).0[1..].rotate_left(1);
Some(this)
}).take(3))
}
fn iter_extensions(&self) -> EndpointIter {
let mut tree = self.tree.as_ref().clone();
tree.scale();
let mut hedron = self.hedron.as_ref().clone();
//hedron.swap(0, 3);
hedron.reverse();
hedron.scale();
Box::new(Endpoint {
tree: Rc::new(tree),
hedron: Rc::clone(&self.hedron),
}.into_iter_directions().filter_map(|endpoint| {
let mut new = endpoint.hedron.mirror(3);
new.swap(0, 3);
if endpoint.tree.collides(&new) {
None
} else {
let mut hedron = endpoint.hedron;
Rc::make_mut(&mut hedron).scale();
Some(Endpoint {
tree: Rc::new(TetraTree {
subtrees: [Some(endpoint.tree), None, None],
hedron,
}),
hedron: Rc::new(new),
})
}
}).chain(self.tree.iter_extensions(Rc::new(TetraTree {
subtrees: [None, None, None],
hedron: Rc::new(hedron),
}))))
}
fn iter_tetrahedra<'a>(&'a self)
-> Box<dyn Iterator<Item = &Tetrahedron> + 'a> {
Box::new(
iter::once(self.hedron.as_ref()).chain(self.tree.iter_tetrahedra())
)
}
}
impl Default for Endpoint {
fn default() -> Endpoint {
Endpoint::from(Tetrahedron::default())
}
}
impl From<Tetrahedron> for Endpoint {
fn from(mut hedron: Tetrahedron) -> Endpoint {
let mirrored = hedron.mirror(0);
hedron.scale();
Endpoint {
tree: Rc::new(TetraTree {
subtrees: [None, None, None],
hedron: Rc::new(mirrored),
}),
hedron: Rc::new(hedron),
}
}
}
impl Hash for Endpoint {
fn hash<H: Hasher>(&self, hasher: &mut H) {
let mut stuff: Vec<_> = self.tree.hash_helper(1).collect();
stuff.push(self.tree.len());
stuff.sort();
stuff.hash(hasher);
}
}
impl PartialEq for Endpoint {
fn eq(&self, other: &Endpoint) -> bool {
self.iter_endpoints().any(|ep| ep.tree == other.tree)
}
}
impl Eq for Endpoint {}
#[derive(Debug, Clone)]
struct TetraTree {
subtrees: [Option<Rc<TetraTree>>; 3],
hedron: Rc<Tetrahedron>,
}
impl TetraTree {
fn iter_endpoints<'x>(&'x self, behind: Rc<TetraTree>) -> EndpointIter<'x> {
let mut iterator = self.subtrees.iter().enumerate()
.filter_map(|(i, opt)| opt.as_ref().map(|some| (i, some)));
if let Some(first) = iterator.next() {
let closure = move |(i, subtree): (usize, &'x Rc<TetraTree>)| -> EndpointIter<'x> {
let mut behind = behind.as_ref().clone();
behind.rotate_left(3-i);
let mut this = self.clone();
this.rotate_left(i);
Rc::make_mut(&mut this.hedron).reverse();
this.subtrees.swap(1, 2);
this.subtrees[0] = Some(Rc::new(behind));
for (j, subtree) in this.subtrees.iter_mut().enumerate().skip(1)
.filter_map(|(j, s)| s.as_mut().map(|s| (j, s))) {
Rc::make_mut(subtree).rotate_left(j);
}
subtree.iter_endpoints(Rc::new(this))
};
Box::new(closure(first).chain(iterator.flat_map(closure)))
} else {
let mut hedron = self.hedron.as_ref().clone();
//hedron.swap(0, 3);
hedron.reverse();
Endpoint {
tree: Rc::clone(&behind),
hedron: Rc::new(hedron),
}.into_iter_directions()
}
}
fn rotate_left(&mut self, i: usize) {
self.subtrees.rotate_left(i);
Rc::make_mut(&mut self.hedron).rotate_left(i);
}
fn collides(&self, hedron: &Tetrahedron) -> bool {
self.hedron.collides(hedron) ||
self.subtrees.iter()
.filter_map(Option::as_ref)
.any(|subtree| subtree.collides(hedron))
}
fn scale(&mut self) {
for subtree in self.subtrees.iter_mut().filter_map(Option::as_mut) {
Rc::make_mut(subtree).scale();
}
Rc::make_mut(&mut self.hedron).scale();
}
fn iter_extensions(&self, behind: Rc<TetraTree>) -> EndpointIter {
Box::new(self.subtrees.iter().enumerate().flat_map(move |(i, next)| {
let mut behind = behind.as_ref().clone();
behind.rotate_left(3-i);
let mut this = self.clone();
this.rotate_left(i);
let hedron = Rc::make_mut(&mut this.hedron);
hedron.reverse();
this.subtrees.swap(1, 2);
this.subtrees[0] = Some(Rc::new(behind));
for (j, subtree) in this.subtrees.iter_mut().enumerate().skip(1)
.filter_map(|(j, s)| s.as_mut().map(|s| (j, s))) {
let subtree = Rc::make_mut(subtree);
subtree.scale();
subtree.rotate_left(j);
}
if let Some(next) = next {
hedron.scale();
next.iter_extensions(Rc::new(this))
} else {
let mut mirrored = hedron.mirror(3);
mirrored.swap(0, 3);
hedron.scale();
if this.collides(&mirrored) {
Box::new(iter::empty()) as EndpointIter
} else {
Box::new(iter::once(Endpoint {
tree: Rc::new(this),
hedron: Rc::new(mirrored),
}))
}
}
}))
}
fn iter_tetrahedra<'a>(&'a self)
-> Box<dyn Iterator<Item = &Tetrahedron> + 'a> {
Box::new(iter::once(self.hedron.as_ref()).chain(
self.subtrees.iter()
.filter_map(Option::as_deref)
.flat_map(TetraTree::iter_tetrahedra)
))
}
fn len(&self) -> usize {
1usize + self.subtrees.iter()
.filter_map(Option::as_deref).map(TetraTree::len).sum::<usize>()
}
fn hash_helper<'a>(&'a self, behind: usize)
-> Box<dyn Iterator<Item = usize> + 'a> {
let sub: Vec<_> = self.subtrees.iter().filter_map(Option::as_deref)
.map(|a| (a, a.len())).collect();
let sum = sub.iter()
.map(|&(_, len)| len).sum::<usize>() + behind + 1;
Box::new(iter::once(behind).chain(
sub.into_iter().flat_map(move |(sub, len)|
iter::once(len).chain(sub.hash_helper(sum - len)))
))
}
}
impl PartialEq for TetraTree {
fn eq(&self, rhs: &TetraTree) -> bool {
self.subtrees == rhs.subtrees
}
}
impl Eq for TetraTree {}
fn main() {
let verbose = std::env::args().skip(1).any(|arg| arg == "-v");
let begin = Instant::now();
println!("1: 1 [{}ms]", begin.elapsed().as_millis());
if verbose {
println!("{}\n--", Tetrahedron::default());
}
let mut polytets = HashSet::new();
polytets.insert(Endpoint::default());
for i in 2.. {
println!("{}: {} [{}ms]", i, polytets.len(), begin.elapsed().as_millis());
if verbose {
for polytet in &polytets {
for hedron in polytet.iter_tetrahedra() {
println!("{}", hedron);
}
println!("--");
}
}
polytets = polytets.iter()
.flat_map(Endpoint::iter_extensions)
.collect();
}
}
Try it online! The footer contains some compatability implementations because TIO's Rust is a little old. (This has already been reported and added to the list.)
Output:
1: 1 [0ms]
2: 1 [0ms]
3: 1 [2ms]
4: 4 [8ms]
5: 10 [35ms]
6: 39 [139ms]
7: 164 [738ms]
8: 767 [4328ms]
9: 3656 [31298ms]
10: 18186 [287871ms]
11: 91532 [3154716ms]
As you can see, it also reports how long it took. (These are cumulated times.) You can use the -v option to list all tetrahedra for each polytet.
I don't really have a good way of verifying the results beyond A276272(5). I hope it works, but I'm not sure.
The idea is that we store the polytet as a tree of tetrahedra that also encodes orientation. But for collision detection we need actual tetrahedra. We start with the tetrahedron with vertices (-1, -1, -1), (-1, 1, 1), (1, -1, 1), (1, 1, -1) and scale all tetrahedra by 3 for every term. This avoids the need for fractions.
Can probably be made faster, but I don't know how.