Gluing tetrahedra together

(This challenge exists to extend sequence A276272 in the On-Line Encyclopedia of Integer Sequences, and perhaps create a new OEIS sequence1.)

This is a , which will have you write code to compute as many terms of this sequence as possible.

Background

A polytet is a kind of polyform which is constructed by gluing regular tetrahedra together face-to-face in such a way that none of the interiors of the tetrahedra overlap. These are counted up to rotation in 3-space but not reflection, so if a particular shape is chiral it is counted twice: once for itself and once for its mirror image.

• $$\A267272(1) = 1\$$ because there is one such configuration you can make out one tetrahedron.
• $$\A267272(2) = 1\$$ because due to symmetry of the tetrahedra, gluing to a face is the same as any other face. This gives you the triangular bipyramid.
• $$\A267272(3) = 1\$$ because the triangular bipyramid itself is symmetric, so adding a tetrahedron to any face results in a biaugmented tetrahedron.
• $$\A267272(4) = 4\$$ because there are three distinct ways to add a tetrahedron to the biaugmented tetrahedron, and one is chiral (so it's counted twice). These are shown in the animated GIFs above.
• $$\A267272(5) = 10\$$ as shown in George Sicherman's catalog. Three are chiral and four are achiral.

Rules

Run your code for as long as you'd like. The winner of this challenge will be the user who posts the most terms of the sequence, along with their code. If two users post the same number of terms, then whoever posts their last term earliest wins.

(Once more terms are computed, then whoever computed the code can add the new terms to the OEIS or I can add them and list the contributor if he or she desires.)

1 A sequence which (presumably) isn't already in the OEIS is the number of polytets up to reflection (i.e. chiral compounds are counted once, not twice). The first few values of this alternative sequence can be found in the George Sicherman link.

• If I'm interpreting this correctly, I believe it's possible to keep all coordinates to rational numbers, which will lessen the burden of symbolic computation. Oct 27, 2020 at 1:03
• @Bubbler—this is correct! It might be useful to start with the tetrahedron in $\mathbb{Z}^4$ with vertices $(1,0,0,0)$, $(0,1,0,0)$, $(0,0,1,0)$, and $(0,0,0,1)$. Oct 27, 2020 at 1:15
• I'm currently working on a solution, first disregarding the collisions between tetrahedra. The bad part is that it'll give me correct answers only up to $A267272(5)$ until I implement collisions D: Oct 28, 2020 at 8:23
• @Bubbler—perhaps that's why they stopped there. If you write a program that can output the coordinates of the tetrahedra (even in 4D), then I can transform the coordinates, and create a visualization in Mathematica to visually check for collisions. Oct 28, 2020 at 16:48
• I just want to add that the oeis sequence should be 276272 not 267272. The link you have works, but the sequence name is wrong in the post. Apr 22, 2021 at 1:04

Rust, 11 (hopefully correct) terms


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
}
}

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 A267272(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.