polycube_packing/

polycube_packing.rs

//! The following program finds all ways to pack 25 [Y pentacubes] into
//! a $5\times5\times5$ box. (Exercise 340 in Section 7.2.2.1 of Knuth's
//! [_The Art of Computer Programming_ **4B** (2022)][taocp4b], Part 2,
//! illustrates the 29 pentacubes.) We formulate this task as an exact
//! cover problem, with one item for each of the $125$ positions to cover,
//! and one option for each [legal placement] of a piece. The program can
//! be readily adapted to fill an arbitrary shape with a given set of
//! polycubes.
//!
//! Chapter 8 of R. Honsberger's book [_Mathematical Gems II_][mgems] (1976)
//! provides a good introduction to the techniques for solving polycube packing
//! puzzles.
//!
//! [Y pentacube]: https://en.wikipedia.org/wiki/Polycube
//! [taocp4b]: https://www-cs-faculty.stanford.edu/~knuth/taocp.html#vol4
//! [legal placement]: `is_in_bounds`
//! [mgems]: https://bookstore.ams.org/dol-2

use exact_covers::{DlSolver, Solver};
use smallvec::SmallVec;
use std::collections::HashSet;
use std::iter;
use std::ops::ControlFlow;

/// A point in the cubic lattice.
type Position = (i8, i8, i8);

/// A solid object formed by joining $1\times 1\times 1$ cubies face to face.
#[derive(Debug, Eq, PartialEq, Hash, Clone)]
struct Polycube(
    /// The positions occupied by the polycube.
    SmallVec<Position, { Self::INLINE_CAP }>,
);

impl Polycube {
    /// If the number of cubies in this polycube exceeds this threshold,
    /// the buffer containing their positions is allocated on the heap.
    const INLINE_CAP: usize = 5; // In this way the array of positions takes 15 bytes,
                                 // which is one less than the space required
                                 // for a 64-bit pointer and a capacity field.

    /// Creates a polycube with one or more cubies, without copying any data.
    ///
    /// This is the `const` version of [`Polycube::from`].
    pub const fn from_const(positions: [Position; Self::INLINE_CAP]) -> Self {
        Self(SmallVec::from_buf(positions))
    }

    /// Applies a transformation to the polycube.
    pub fn transform(&self, f: impl FnMut(Position) -> Position) -> Self {
        Self(self.0.iter().copied().map(f).collect())
    }

    /// Returns the translate of the polycube whose cubie coordinates are
    /// nonnegative and as small as possible.
    ///
    /// This shape is called the _aspect_ of the polycube. Section 1.3 of
    /// the book [_Tilings and Patterns_][tap] by B. Grünbaum, G. C. Shephard
    /// (W. H. Freeman, 1987) discusses the number of distinct aspects
    /// of the tiles in particular classes of tessellations.
    ///
    /// [tap]: https://dl.acm.org/doi/10.5555/19304
    pub fn aspect(&self) -> Self {
        let x_min = self.0.iter().map(|(x, _, _)| x).min().unwrap();
        let y_min = self.0.iter().map(|(_, y, _)| y).min().unwrap();
        let z_min = self.0.iter().map(|(_, _, z)| z).min().unwrap();
        self.transform(|(x, y, z)| (x - x_min, y - y_min, z - z_min))
    }

    /// Constructs the set of aspects corresponding to all 3D-rotations of
    /// the polycube.
    ///
    /// D. E. Knuth called these the _base placements_ of a polycube in
    /// exercise 7.2.2.1–266 of [_The Art of Computer Programming_ **4B**,
    /// Part 2][taocp] (Addison-Wesley, 2022).
    ///
    /// [taocp]: https://www-cs-faculty.stanford.edu/~knuth/taocp.html#vol4
    pub fn base_placements(&self) -> HashSet<Polycube> {
        // The group of symmetries of a polycube is isomorphic to a subgroup of
        // the 24 3D-rotations of a cube, which is generated by the following
        // two elementary transformations:
        // 1. $90^\circ$ rotation about the $z$-axis, $(x,y,z)\mapsto(y,x_\text{max}-x,z)$.
        // 2. Reflection about the $x=y=z$ diagonal: $(x,y,z)\mapsto(y,z,x)$.
        // To see why, use the fact that the rotational symmetries of a cube
        // form a group that is isomorphic to the symmetric group $S_4$ on 4
        // elements. (Notice that the rotations of a cube permute its diagonals.)
        // In this way transformations 1 and 2 correspond respectively to
        // the permutations $\pi=(1234)$ and $\sigma=(142)$ under some
        // appropriate naming of the cube vertices. And $\{\pi,\sigma\}$ is
        // a generator of $S_4$.
        let mut placements = HashSet::with_capacity(24);
        let mut to_visit = Vec::with_capacity(8);
        to_visit.push(self.aspect());
        while let Some(shape) = to_visit.pop() {
            let x_max = shape.0.iter().map(|(x, _, _)| x).max().unwrap();
            let rotation = shape.transform(|(x, y, z)| (y, x_max - x, z));
            if !placements.contains(&rotation) {
                to_visit.push(rotation);
            }

            let reflection = shape.transform(|(x, y, z)| (y, z, x));
            if !placements.contains(&reflection) {
                to_visit.push(reflection);
            }

            placements.insert(shape);
        }
        placements
    }

    /// Returns the number of cubies in the polycube.
    pub fn cubie_count(&self) -> usize {
        self.0.len()
    }
}

impl From<&[Position]> for Polycube {
    /// Creates a polycube with one or more cubies.
    fn from(positions: &[Position]) -> Self {
        assert!(!positions.is_empty(), "a polycube has one or more cubies");
        Self(SmallVec::from_slice(positions))
    }
}

impl From<Vec<Position>> for Polycube {
    /// Creates a polycube with one or more cubies.
    fn from(positions: Vec<Position>) -> Self {
        assert!(!positions.is_empty(), "a polycube has one or more cubies");
        Self(positions.into())
    }
}

/// The Y pentacube, with the vector normal to its bottom face—the face that is
/// farthest-apart from the pentacube's center of mass—pointing downwards and
/// the cubie not in the long bar appearing in the up-west position.
const Y: Polycube = Polycube::from_const([(0, 2, 0), (1, 0, 0), (1, 1, 0), (1, 2, 0), (1, 3, 0)]);

/// Returns `true` if and only if the given polycube lies inside the
/// $l\times m\times n$ cuboid cornered at the origin.
fn is_in_bounds(polycube: &Polycube, l: i8, m: i8, n: i8) -> bool {
    let x_min = *polycube.0.iter().map(|(x, _, _)| x).min().unwrap();
    let y_min = *polycube.0.iter().map(|(_, y, _)| y).min().unwrap();
    let z_min = *polycube.0.iter().map(|(_, _, z)| z).min().unwrap();
    let x_max = *polycube.0.iter().map(|(x, _, _)| x).max().unwrap();
    let y_max = *polycube.0.iter().map(|(_, y, _)| y).max().unwrap();
    let z_max = *polycube.0.iter().map(|(_, _, z)| z).max().unwrap();
    0 <= x_min && x_max < l && 0 <= y_min && y_max < m && 0 <= z_min && z_max < n
}

fn main() {
    // Define the items of the exact cover problem, namely the $lmn$ cells of
    // the $l\times m\times n$ cuboid.
    let (l, m, n) = (5i8, 5i8, 5i8);
    let positions = (0..l)
        .flat_map(|x| iter::repeat(x).zip(0..m))
        .flat_map(|y| iter::repeat(y).zip(0..n))
        .map(|((x, y), z)| (x, y, z))
        .collect::<Vec<_>>();

    let mut solver: DlSolver<Position, ()> = DlSolver::new(&positions, &[]);
    // For each base placement $P$ of the Y pentacube, and for each offset
    // $(x_0,y_0,z_0)$ such that the piece $P'=P+(x_0,y_0,z_0)$ lies within the
    // $l\times m\times n$ cuboid cornered at the origin, define an option whose
    // primary items are the cells of $P'$. We break symmetry by constraining
    // a particular cubie of the first pentacube to be at $L_\infty$ distance
    // $\le2$ from the origin.
    let placements = Y.base_placements();
    let mut first = true;
    for placement in placements {
        for (x_0, y_0, z_0) in &positions {
            if first && (*x_0 > 2 || *y_0 > 2 || *z_0 > 2) {
                continue;
            }
            let shifted = placement.transform(|(x, y, z)| (x + x_0, y + y_0, z + z_0));
            if is_in_bounds(&shifted, l, m, n) {
                solver.add_option(&shifted.0, []);
            }
        }
        first = false;
    }

    let mut count = 0;
    let mut polycube = Vec::with_capacity(Y.cubie_count());
    solver.solve(|mut solution| {
        // Print the solution, which consists of a set of 25 pentacubes.
        print!("[");
        while solution.next(&mut polycube) {
            print!("[");
            if let Some(((&last, _), elements)) = polycube.split_last() {
                for (&cubie, _) in elements {
                    print!("[{},{},{}],", cubie.0, cubie.1, cubie.2);
                }
                print!("[{},{},{}]", last.0, last.1, last.2);
            }
            print!("],");
        }
        println!("]");
        count += 1;
        ControlFlow::Continue(())
    });
    println!("found {count} packings");
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn y_base_placements() {
        let placements = Y.base_placements();
        assert_eq!(placements.len(), 24, "the Y pentacube has 24 3D-rotations");
        for placement in &placements {
            assert_eq!(
                placement,
                &placement.aspect(),
                "a base placement is an aspect"
            );
        }
    }
}