Gluing tetrahedra together

Question

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

This is a code-challenge, 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.)

---

¹ 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.

Answer 1

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