Xenharmonic Wiki

Code for billiard scales: Difference between revisions

← Older editNewer edit →
VisualWikitext
Inthar (talk | contribs)
autopatrolled, emailconfirmed
15,668 edits
Replaced content with "{{delete|Now part of Inthar's ternary app}}"
Tag: Replaced
Inthar (talk | contribs)
autopatrolled, emailconfirmed
15,668 edits
Undo revision 220438 by Inthar (talk)
Line 1: Line 1:
{{delete|Now part of Inthar's ternary app}}
== Ternary billiard scales ==
To run this, you need Rust installed.
# In any folder you want, create a directory <code>billiard-scales</code>.
# Using the command line, run <code>cargo init</code> in the new directory <code>billiard-scales</code>.
# Create or copy the following .rs files in <code>billiard-scales/src</code>. Copy the contents of the <code>Cargo.toml</code> provided into the new <code>Cargo.toml</code>
# Run <code>cargo run --release</code>. The script should notify you when it has finished running.
# Open the output file <code>billiard_scales_mediawiki</code> for the result.
=== plane_geometry.rs ===
<syntaxhighlight lang="rs">
// Geometry of points and lines in the plane, implemented with rational numbers. The implementation is janky, though.
// Can be used to enumerate ternary billiard scales.
use num_traits::sign::Signed;
use std::fmt;
use std::iter::zip;
use std::ops::Add;
use std::ops::Neg;
use std::ops::Sub;
 
use num_rational::Rational32 as r32;
 
use crate::helpers::gcd;
 
#[derive(Copy, Clone, Debug, Hash, PartialEq, Eq)]
/// The possible configurations of a point and an oriented line on the plane.
pub enum PointLineConfiguration {
    Left,      // The point is on the left half-plane.
    Right,    // The point is on the right half-plane.
    OnTheLine, // The point is on the line.
}
 
fn is_between_r32(a: r32, b: r32, c: r32) -> bool {
    if b == c {
        a == b
    } else if b < c {
        b < a && a < c
    } else {
        c < a && a < b
    }
}
 
/// A rational number or unsigned infinity; represents a valid slope for a line in Q^2.
#[derive(Copy, Clone, Debug)]
pub struct Slope {
    numer: i32,
    denom: i32,
}
 
impl PartialEq for Slope {
    fn eq(&self, other: &Self) -> bool {
        self.numer * other.denom == other.numer * self.denom
    }
}
 
impl fmt::Display for Slope {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "{}/{}", self.numer, self.denom)
    }
}
 
impl Slope {
    /// Converts a `Slope` to a rational; returns Err if `denom` is 0.
    pub fn as_rational(&self) -> Result<r32, String> {
        if self.denom != 0 {
            Ok(r32::new(self.numer, self.denom))
        } else {
            Err("Attempted to convert a Slope with `denom` == 0 into a Rational64".to_string())
        }
    }
 
    /// Converts the projective rational to a rational; panics if the denominator is zero.
    pub fn as_rational_raw(&self) -> r32 {
        r32::new(self.numer, self.denom)
    }
 
    /// Create a new reduced Slope; panics if both numer and denom are 0.
    pub fn new(n: i32, d: i32) -> Result<Self, String> {
        let (mut numer, mut denom) = (n, d);
        if numer == 0 && denom == 0 {
            Err("Attempted to create a Slope with both `numer` and `denom` equal to 0".to_string())
        } else {
            if denom < 0 {
                numer = -numer;
                denom = -denom;
            }
            let d = gcd(numer as u32, denom as u32) as i32;
            if d > 1 {
                numer /= d;
                denom /= d;
            }
            Ok(Self { numer, denom })
        }
    }
 
    /// Create a new Slope without checking that `numer` == `denom` == 0, or reducing the rational afterwards.
    pub fn raw(numer: i32, denom: i32) -> Self {
        Self { numer, denom }
    }
 
    /// Create a new `Slope` from a rational.
    pub fn rational(rat: r32) -> Self {
        Self {
            numer: *rat.numer(),
            denom: *rat.denom(),
        }
    }
 
    /// Create a new `Slope` from an integer.
    pub fn integer(n: i32) -> Self {
        // use 'raw' since an integer is always reduced
        Self::raw(n, 1)
    }
 
    /// Shorthand for slope 0.
    pub fn zero() -> Self {
        Self::raw(0, 1)
    }
 
    /// Shorthand for slope 1.
    pub fn one() -> Self {
        Self::raw(1, 1)
    }
 
    /// Shorthand for infinite slope.
    pub fn infinity() -> Self {
        Self::raw(1, 0)
    }
}
 
/// A finite point; a member of Q^2.
#[derive(Copy, Clone, Debug, Hash, PartialEq, Eq)]
pub struct Point {
    // The x-coordinate.
    x: r32,
    // The y-coordinate.
    y: r32,
}
 
impl Add for Point {
    type Output = Self;
 
    fn add(self, other: Self) -> Self {
        Self {
            x: self.x + other.x,
            y: self.y + other.y,
        }
    }
}
 
impl Sub for Point {
    type Output = Self;
 
    fn sub(self, other: Self) -> Self {
        Self {
            x: self.x - other.x,
            y: self.y - other.y,
        }
    }
}
 
impl Neg for Point {
    type Output = Self;
 
    fn neg(self) -> Self {
        Self {
            x: -self.x,
            y: -self.y,
        }
    }
}
 
impl Point {
    pub fn new(x: r32, y: r32) -> Self {
        Self { x, y }
    }
 
    pub fn zero() -> Self {
        Self {
            x: r32::new(0, 1),
            y: r32::new(0, 1),
        }
    }
 
    pub fn scalar_mult(lambda: r32, v: Point) -> Self {
        Self {
            x: lambda * v.x,
            y: lambda * v.y,
        }
    }
 
    pub fn average(points: &Vec<Self>) -> Self {
        let mut sum = Self::zero();
        for p in points {
            sum = sum + *p;
        }
        Self::scalar_mult(r32::new(1, 1) / r32::new(points.len() as i32, 1), sum)
    }
 
    pub fn theta(v: Point) -> f64 {
        (*v.y.numer() as f64 / *v.y.denom() as f64).atan2(*v.x.numer() as f64 / *v.x.denom() as f64)
    }
 
    /// Return a copy of `points` sorted using the direction they make from the centroid, according to theta(v-centroid) in (-pi, pi].
    fn ordered_ccw(mut points: Vec<Self>) -> Vec<Self> {
        let centroid = Self::average(&points);
        points.sort_by(|a, b| {
            Self::theta(*a - centroid)
                .partial_cmp(&Self::theta(*b - centroid))
                .unwrap()
        });
        points
    }
 
    pub fn midpoint(p1: &Self, p2: &Self) -> Self {
        Self {
            x: (p1.x + p2.x) / r32::new(2, 1),
            y: (p1.y + p2.y) / r32::new(2, 1),
        }
    }
    // Use the slope-intercept formula including the case where the slope is infinite.
    pub fn are_collinear(p1: Self, p2: Self, p3: Self) -> bool {
        (p1.y - p3.y) * (p1.x - p2.x) == (p1.y - p2.y) * (p1.x - p3.x)
    }
}
 
impl fmt::Display for Point {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "({}, {})", self.x, self.y)
    }
}
 
/// An oriented line in Q^2 of rational or infinite slope.
#[derive(Copy, Clone, Debug)]
pub struct Line {
    point1: Point,
    point2: Point,
}
 
impl PartialEq for Line {
    fn eq(&self, other: &Self) -> bool {
        self.slope() == other.slope() && self.has_point(other.point())
    }
}
 
impl Line {
    pub fn slope(&self) -> Slope {
        Slope::new(
            (*(self.point2.y - self.point1.y).numer()) * (*(self.point2.x - self.point1.x).denom()),
            (*(self.point2.x - self.point1.x).numer()) * (*(self.point2.y - self.point1.y).denom()),
        )
        .unwrap()
    }
 
    pub fn point(&self) -> Point {
        self.point1
    }
 
    pub fn from_slope_and_point(sl: &Slope, point: &Point) -> Self {
        Self {
            point1: *point,
            point2: *point + Point::new(r32::new_raw(sl.denom, 1), r32::new_raw(sl.numer, 1)),
        }
    }
    pub fn from_points(point1: &Point, point2: &Point) -> Result<Self, String> {
        if point1 == point2 {
            Err("The points {point1} and {point2} are equal".to_string())
        } else {
            Ok(Line {
                point1: *point1,
                point2: *point2,
            })
        }
    }
 
    /// Test whether a line is vertical.
    fn is_vertical(&self) -> bool {
        self.point1.x == self.point2.x
    }
 
    /// Test whether `p` satisfies the equation of the line.
    fn has_point(&self, p: Point) -> bool {
        (p.y - self.point().y) * r32::new_raw(self.slope().denom, 1)
            == r32::new_raw(self.slope().numer, 1) * (p.x - self.point().x)
    }
 
    /// Determine whether `p1` and `p2` are on the same side of the line.
    fn on_same_side(&self, p1: Point, p2: Point) -> bool {
        if self.is_vertical() {
            (p1.x < self.point1.x && p2.x < self.point1.x)
                || (p1.x > self.point1.x && p2.x > self.point1.x)
        } else {
            // both `p1` and `p2` are above `self`
            ((p1.y - self.point1.y)*r32::new_raw(self.slope().denom,1) > r32::new_raw(self.slope().numer,1)*(p1.x - self.point1.x) &&
                (p2.y - self.point1.y)*r32::new_raw(self.slope().denom,1) > r32::new_raw(self.slope().numer,1)*(p2.x - self.point1.x))
            || // ... or both are below `self`
            ((p1.y - self.point1.y)*r32::new_raw(self.slope().denom,1) < r32::new_raw(self.slope().numer,1)*(p1.x - self.point1.x) &&
                (p2.y - self.point1.y)*r32::new_raw(self.slope().denom,1) < r32::new_raw(self.slope().numer,1)*(p2.x - self.point1.x))
        }
    }
 
    /// Depending on the orientation of the line, determine whether the point `p` is to the left, to the right, or on the line.
    pub fn point_line_config(&self, p: Point) -> PointLineConfiguration {
        if self.point2.x != self.point1.x {
            // For non-vertical lines
            if ((p.y - self.point().y) * r32::new_raw(self.slope().denom, 1)
                - r32::new_raw(self.slope().numer, 1) * (p.x - self.point().x))
                * (self.point2.x - self.point1.x).signum()
                > r32::new_raw(0, 1)
            {
                PointLineConfiguration::Left
            } else if ((p.y - self.point().y) * r32::new_raw(self.slope().denom, 1)
                - r32::new_raw(self.slope().numer, 1) * (p.x - self.point().x))
                * (self.point2.x - self.point1.x).signum()
                == r32::new_raw(0, 1)
            {
                PointLineConfiguration::OnTheLine
            } else {
                PointLineConfiguration::Right
            }
        } else {
            // For vertical lines
            // If point2.y < point1.y, then the line is oriented downwards, so multiply by -1 to compensate.
            if (p.x - self.point().x) * (self.point2.y - self.point1.y) < r32::new_raw(0, 1) {
                // upwards vertical line and point to the left of it, or downwards vertical line and point to the right of it
                PointLineConfiguration::Left
            } else if (p.x - self.point().x) * (self.point2.y - self.point1.y) == r32::new_raw(0, 1)
            {
                PointLineConfiguration::OnTheLine
            } else {
                PointLineConfiguration::Right
            }
        }
    }
}
 
impl fmt::Display for Line {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        write!(f, "(slope: {}, point: {})", self.slope(), self.point1)
    }
}
 
/// For two non-parallel ones, find the intersection; for two parallel lines, return `None`.
pub fn intersection(l1: &Line, l2: &Line) -> Option<Point> {
    // We only call as_rational_raw() on non-infinite slopes so it's ok
    let (x1, y1, x2, y2) = (l1.point().x, l1.point().y, l2.point().x, l2.point().y);
    match (l1.slope(), l2.slope()) {
        // Equal slope
        (s1, s2) if s1 == s2 => None,
        // Unequal slopes, s1 is infinite
        (s1, s2) if s1 == Slope::infinity() && s2 != Slope::infinity() => Some(Point {
            x: x1,
            y: s2.as_rational_raw() * (x1 - x2) + y2,
        }),
        // Unequal slopes, s2 is infinite
        (s1, s2) if s1 != Slope::infinity() && s2 == Slope::infinity() => Some(Point {
            x: x2,
            y: s1.as_rational_raw() * (x2 - x1) + y1,
        }),
        // Unequal slopes, neither infinite
        (s1, s2) => Some(Point {
            x: (y1 - s1.as_rational_raw() * x1 + s2.as_rational_raw() * x2 - y2)
                / (s2.as_rational_raw() - s1.as_rational_raw()),
            y: (s1.as_rational_raw() * s2.as_rational_raw() * (x2 - x1)
                + s2.as_rational_raw() * y1
                - s1.as_rational_raw() * y2)
                / (s2.as_rational_raw() - s1.as_rational_raw()),
        }),
    }
}
 
#[derive(Clone, Debug)]
pub struct ConvexPolygon {
    vertices: Vec<Point>, // Assumes that the vertices are ordered in clockwise or counterclockwise order.
}
 
impl PartialEq for ConvexPolygon {
    fn eq(&self, other: &Self) -> bool {
        // If all the points in both polygons are equal then they are equal; we ensure this by enforcing a standard order when a `ConvexPolygon` is created.
        if other.vertices.len() != self.vertices.len() {
            return false;
        }
        for (a, b) in zip(self.vertices.iter(), other.vertices.iter()) {
            if a != b {
                return false;
            }
        }
        true
    }
}
 
impl ConvexPolygon {
    // Assumes that the vertices form a convex polygon. Does not check if the `Vec` of vertices consists of distinct elements.
    pub fn new(vertices: &[Point]) -> Option<Self> {
        if vertices.len() >= 3 {
            Some(Self {
                vertices: Point::ordered_ccw(vertices.to_vec()), // This will work for vertices of a convex polygon.
            })
        } else {
            None
        }
    }
 
    pub fn centroid(&self) -> Point {
        Point::average(&self.vertices)
    }
 
    pub fn has_point_inside(&self, point: Point) -> bool {
        // Do `point` and edge_i+1 x edge_i+2 lie on the same side of edge_i?
        for i in 0..self.vertices.len() - 2 {
            if !Line::from_points(&self.vertices[i], &self.vertices[i + 1])
                .unwrap()
                .on_same_side(point, self.vertices[i + 2])
            {
                return false;
            }
        }
        // Check last two triples.
        if !Line::from_points(
            &self.vertices[self.vertices.len() - 2],
            &self.vertices[self.vertices.len() - 1],
        )
        .unwrap()
        .on_same_side(point, self.vertices[0])
        {
            false
        } else {
            Line::from_points(&self.vertices[self.vertices.len() - 1], &self.vertices[0])
                .unwrap()
                .on_same_side(point, self.vertices[1])
        }
    }
 
    /// A vector storing the edges of the polygon in order.
    fn lines_for_edges(&self) -> Vec<Line> {
        let mut result: Vec<Line> = (0..self.vertices.len() - 1)
            .map(|i| Line::from_points(&self.vertices[i], &self.vertices[i + 1]).unwrap())
            .collect();
        result.push(
            Line::from_points(&self.vertices[self.vertices.len() - 1], &self.vertices[0]).unwrap(),
        );
        result
    }
 
    /// Subdivide a polygon with the given line and return the two pieces.
    pub fn subdivide(&self, line: &Line) -> (Option<ConvexPolygon>, Option<ConvexPolygon>) {
        // Collect vertices that fall to the left, on the line, and to the right, respectively.
        // The left and right collections are missing two vertices; we need to append them to make them `ConvexPolygon`s.
        let mut vertices_left: Vec<Point> = self
            .vertices
            .iter()
            .filter(|&x| line.point_line_config(*x) == PointLineConfiguration::Left)
            .cloned()
            .collect();
        let mut vertices_right: Vec<Point> = self
            .vertices
            .iter()
            .filter(|&x| line.point_line_config(*x) == PointLineConfiguration::Right)
            .cloned()
            .collect();
        let vertices_on_line: Vec<Point> = self
            .vertices
            .iter()
            .filter(|&x| line.point_line_config(*x) == PointLineConfiguration::OnTheLine)
            .cloned()
            .collect();
        if !vertices_left.is_empty() || !vertices_right.is_empty() {
            // Add any vertices on the line both to the new left polygon and the new right polygon.
            for point in &vertices_on_line {
                vertices_left.push(*point);
                vertices_right.push(*point);
            }
            // Look for intersection points on edges with the line.
            for edge in &self.lines_for_edges() {
                if let Some(intsn) = intersection(edge, line)
                    && is_between_r32(intsn.x, edge.point1.x, edge.point2.x)
                    && is_between_r32(intsn.y, edge.point1.y, edge.point2.y)
                {
                    vertices_left.push(intsn);
                    vertices_right.push(intsn);
                }
            }
            // At most two vertices are added, so if there was no vertex to one side originally, that side doesn't have a polygon.
        }
        (
            ConvexPolygon::new(&vertices_left),
            ConvexPolygon::new(&vertices_right),
        )
    }
}
</syntaxhighlight>
=== billiard.rs ===
<syntaxhighlight lang="rs">
//! Algorithm for generating ternary [billiard scales](https://en.xen.wiki/w/Hypercubic_billiard_word).
//!
//! This module implements a geometric approach to scale generation: a "billiard ball"
//! bounces inside a 3D cube, and the sequence of faces it hits determines
//! the scale's step pattern. A 2D projection of the cube is used to simplify the computation.
//!
//! # Algorithm
//!
//! 1. Project the unit cube in R³ via a map with `(a, b, c)` in its kernel
//! 2. Partition the projected hexagon by constraint planes
//! 3. For each region, trace a trajectory and record which coordinate changes
//! 4. The sequence of changes (0=L, 1=m, 2=s) forms the scale word
//!
//! # Examples
//!
//! ```
//! use ternary::billiard::billiard_scales;
//!
//! // Generate all billiard scales with signature 5L 2m 2s
//! let scales = billiard_scales(5, 2, 2);
//! assert!(!scales.is_empty());
//!
//! // Each scale has the correct number of steps
//! for scale in &scales {
//!    assert_eq!(scale.len(), 9);  // 5 + 2 + 2 = 9 notes
//! }
//! ```
//!
//! # Properties
//!
//! Billiard scales have special structural properties:
//! - They are always [deletion-MOS](https://en.xen.wiki/w/Deletion-MOS)
//!
//! Not all ternary scales are billiard scales — this is a strict subset.
 
use core::fmt;
 
use num_rational::Rational32 as r32;
 
use crate::helpers::gcd;
use crate::plane_geometry::{ConvexPolygon, Line, Point, PointLineConfiguration, Slope};
use crate::words::{Letter, least_mode};
 
// ERRORS
 
#[derive(Copy, Clone, PartialEq, Eq, Debug)]
/// Error type for illegal billiard trajectory states.
pub enum BadBilliardState {
    // Projected cube is invalid.
    ProjectedCubeInvalid((u32, u32, u32)),
    // Point is not inside the projected cube,
    PointNotInCube(Point),
    // Point is inside the projected cube but falls on a projected constraint plane
    PointOnConstraintPlane(Point),
}
 
impl fmt::Display for BadBilliardState {
    fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
        match self {
            Self::ProjectedCubeInvalid((a, b, c)) => {
                write!(f, "Projected cube for signature {a}L{b}m{c}s is invalid")
            }
            Self::PointNotInCube(pt) => write!(f, "Point {} is not inside the projected cube", pt),
            Self::PointOnConstraintPlane(pt) => write!(
                f,
                "Point {} is inside the projected cube but falls on a projected constraint plane",
                pt
            ),
        }
    }
}
 
impl std::error::Error for BadBilliardState {}
 
fn projected_constraint_planes(a: u32, b: u32, c: u32) -> (Vec<Line>, Vec<Line>, Vec<Line>) {
    let first = (0..b + c + 1)
        .map(|i| {
            Line::from_slope_and_point(
                &Slope::raw(0, 1),
                &Point::new(r32::new(0, 1), r32::new(-(b as i32) + (i as i32), c as i32)),
            )
        })
        .collect(); // horizontal lines going up (Left)
    let second = (0..a + c + 1)
        .map(|j| {
            Line::from_slope_and_point(
                &Slope::raw(1, 0),
                &Point::new(r32::new((c as i32) - (j as i32), c as i32), r32::new(0, 1)),
            )
        })
        .collect(); // vertical lines going left (Left)
    let third = (0..a + b + 1)
        .map(|k| {
            Line::from_slope_and_point(
                &Slope::raw(b as i32, a as i32),
                &Point::new(r32::new(0, 1), r32::new(-(b as i32) + (k as i32), a as i32)),
            )
        })
        .collect(); // lines of pos. slope going up (Left)
    (first, second, third)
}
 
/// Project the unit cube in R^3 via a linear map that has (a,b,c) in its kernel.
fn projected_cube(a: u32, b: u32, c: u32) -> Option<ConvexPolygon> {
    if a > 0 && b > 0 && c > 0 {
        let projected_cube_vertices = vec![
            Point::new(r32::new(0, 1), r32::new(1, 1)),
            Point::new(r32::new(1, 1), r32::new(1, 1)),
            Point::new(r32::new(1, 1), r32::new(0, 1)),
            Point::new(
                r32::new(-(a as i32), c as i32),
                r32::new(-(b as i32), c as i32) + r32::new(1, 1),
            ),
            Point::new(
                r32::new(-(a as i32), c as i32),
                r32::new(-(b as i32), c as i32),
            ),
            Point::new(
                r32::new(-(a as i32), c as i32) + r32::new(1, 1),
                r32::new(-(b as i32), c as i32),
            ),
        ];
        ConvexPolygon::new(&projected_cube_vertices)
    } else {
        None
    }
}
 
/// Partition the projected unit cube with the projected constraint planes.
fn project_and_partition(a: u32, b: u32, c: u32) -> Option<Vec<ConvexPolygon>> {
    projected_cube(a, b, c).map(|projected_cube| {
        let lineses = projected_constraint_planes(a, b, c); // it's a collection of collections of `Line`s
        let mut first_shapes = Vec::<ConvexPolygon>::new();
        let mut second_shapes = Vec::<ConvexPolygon>::new();
        let mut third_shapes = Vec::<ConvexPolygon>::new();
        let mut i: usize = 0;
        let mut remaining: ConvexPolygon = projected_cube;
        loop {
            let (maybe_remaining, maybe_next) = remaining.subdivide(&(lineses.0)[i]);
            // The remaining portion is to the Left.
            match (maybe_remaining, maybe_next) {
                (Some(r), None) => {
                    remaining = r;
                }
                (Some(r), Some(next)) => {
                    remaining = r;
                    first_shapes.push(next);
                }
                (None, Some(next)) => {
                    // This is the last `next` you need to look at.
                    first_shapes.push(next);
                    break;
                }
                (_, _) => {
                    break;
                }
            }
            i += 1;
        }
        for polygon in first_shapes {
            let mut j: usize = 0;
            let mut remaining2: ConvexPolygon = polygon;
            loop {
                let (maybe_remaining, maybe_next) = remaining2.subdivide(&(lineses.1)[j]); // The remaining portion is to the Left.
                match (maybe_remaining, maybe_next) {
                    (Some(r), None) => {
                        remaining2 = r;
                    }
                    (Some(r), Some(next)) => {
                        remaining2 = r;
                        second_shapes.push(next);
                    }
                    (None, Some(next)) => {
                        second_shapes.push(next);
                        break;
                    }
                    (_, _) => {
                        break;
                    }
                }
                j += 1;
            }
        }
        for polygon in second_shapes {
            let mut j: usize = 0;
            let mut remaining3: ConvexPolygon = polygon;
            loop {
                let (maybe_remaining, maybe_next) = remaining3.subdivide(&(lineses.2)[j]);
                // The remaining portion is to the Left.
                match (maybe_remaining, maybe_next) {
                    (Some(r), None) => {
                        remaining3 = r;
                    }
                    (Some(r), Some(next)) => {
                        remaining3 = r;
                        third_shapes.push(next);
                    }
                    (None, Some(next)) => {
                        third_shapes.push(next);
                        break;
                    }
                    (_, _) => {
                        break;
                    }
                }
                j += 1;
            }
        }
        third_shapes
    })
}
 
/// Advance a projected point, `pt`, by one letter according to what integer coordinate it next hits as it traverses in the given velocity,
/// and return the next point and the corresponding letter, one of "0", "1", and "2".
fn advance(a: u32, b: u32, c: u32, pt: Point) -> Result<(Point, Letter), BadBilliardState> {
    if let Some(projected_cube) = projected_cube(a, b, c) {
        // Consider: _ /
        //            |
        if !projected_cube.has_point_inside(pt) {
            Err(BadBilliardState::PointNotInCube(pt))
        } else {
            let projected_x_y_equals_1_1 = Line::from_slope_and_point(
                &Slope::raw(b as i32, a as i32),
                &Point::new(r32::new(0, 1), r32::new(-(b as i32) + (a as i32), a as i32)),
            ); // the /
            let projected_x_z_equals_1_1 = Line::from_slope_and_point(
                &Slope::raw(0, 1),
                &Point::new(r32::new(0, 1), r32::new(-(b as i32) + (c as i32), c as i32)),
            ); // the _, oriented right
            let projected_y_z_equals_1_1 = Line::from_slope_and_point(
                &Slope::raw(1, 0),
                &Point::new(r32::new((c as i32) - (a as i32), c as i32), r32::new(0, 1)),
            ); // the |, oriented up
            if projected_x_y_equals_1_1.point_line_config(pt) == PointLineConfiguration::Right
                && projected_y_z_equals_1_1.point_line_config(pt) == PointLineConfiguration::Right
            {
                // Below the / and to the right of the | => x has changed, change x coordinate to 0 and project; move x coordinate to the left
                Ok((pt - Point::new(r32::new(1, 1), r32::new(0, 1)), 0)) // new point is one unit to the left
            } else if projected_x_z_equals_1_1.point_line_config(pt) == PointLineConfiguration::Left
                && projected_x_y_equals_1_1.point_line_config(pt) == PointLineConfiguration::Left
            {
                // Above the _ and to the left of the / => y has changed, change y coordinate to 0 and project; move y coordinate down
                Ok((pt - Point::new(r32::new(0, 1), r32::new(1, 1)), 1)) // new point is one unit below
            } else if projected_x_z_equals_1_1.point_line_config(pt)
                == PointLineConfiguration::Right
                && projected_y_z_equals_1_1.point_line_config(pt) == PointLineConfiguration::Left
            {
                // Below the _ and to the left of the | => z has changed, change y coordinate to 0 and project
                Ok((
                    pt + Point::new(r32::new(a as i32, c as i32), r32::new(b as i32, c as i32)),
                    2,
                )) // new point is one unit below in the z-direction, and e3 is projected to (-a/c, -b/c)
            } else {
                Err(BadBilliardState::PointOnConstraintPlane(pt))
            }
        }
    } else {
        Err(BadBilliardState::ProjectedCubeInvalid((a, b, c)))
    }
}
 
/// Generate all billiard scales with the given step signature.
///
/// Returns scales in canonical form (lexicographically first rotation),
/// sorted and deduplicated.
///
/// # Arguments
///
/// * `a` - Number of L (large) steps
/// * `b` - Number of m (medium) steps
/// * `c` - Number of s (small) steps
///
/// # Examples
///
/// ```
/// use ternary::billiard::billiard_scales;
///
/// // Generate billiard scales for 5L 2m 2s
/// let scales = billiard_scales(5, 2, 2);
/// assert!(!scales.is_empty());
///
/// // All scales are in canonical form (sorted)
/// for scale in &scales {
///    assert_eq!(scale.iter().filter(|&&x| x == 0).count(), 5);  // 5 L steps
/// }
/// ```
///
/// # Returns
///
/// Empty vector if the signature is invalid (any component is 0).
pub fn billiard_scales(a: usize, b: usize, c: usize) -> Vec<Vec<Letter>> {
    let (a, b, c) = (a as u32, b as u32, c as u32);
    let d = gcd(gcd(a, b), c);
    if d != 0 {
        let mut result = Vec::<Vec<Letter>>::new();
        let n = a / d + b / d + c / d;
        if let Some(regions) = project_and_partition(a / d, b / d, c / d) {
            for region in regions {
                let mut word = Vec::<Letter>::new();
                let mut next_point = region.centroid();
                let first_point = region.centroid();
                for i in 0..n {
                    let res = advance(a, b, c, next_point)
                    .expect("This should never fail as long as you begin from a point inside a polygonal region");
                    next_point = res.0;
                    debug_assert!(i == n - 1 || first_point != next_point); // There should not be a subperiod in the orbit of the "billiard ball".
                    word.push(res.1);
                }
                result.push({
                    // Repeat each word d times
                    let mut r = Vec::with_capacity(word.len() * d as usize); // Pre-allocate memory
                    for _ in 0..d {
                        r.extend_from_slice(&word); // Append a slice of the original vec
                    }
                    r
                });
                debug_assert_eq!(first_point, next_point); // The last point should come back to the first point. When i = n - 1, next_point gets set to point number n.
            }
            result = result
                .into_iter()
                .map(|scale| least_mode(&scale))
                .collect::<Vec<_>>();
            result.sort();
            result.dedup();
            result
        } else {
            vec![]
        }
    } else {
        vec![]
    }
}
</syntaxhighlight>
=== lib.rs ===
<syntaxhighlight lang="rs">use std::collections::BTreeSet;
use num_rational::Rational64 as r64;
mod plane_geometry;
use plane_geometry::*;
mod billiard;
use billiard::*;
use std::cmp::{min, max};
 
type EquivalenceRelation<T> = fn(&T, &T) -> bool;
 
/// Treating `scale` as a cyclic string (that is, "scale[i] == scale[i % scale.len()]"),
/// take a slice of length `slice_length` from `basepoint`; assumes `slice_length` <= `scale`.len()
fn slice_cyclic_string(scale: &str, basepoint: usize, slice_length: usize) -> String {
let arr = scale;
if basepoint + slice_length < scale.len() {
return String::from(&arr[basepoint..(basepoint+slice_length)]);
} else {
let mut ret = String::new();
ret.push_str(&arr[basepoint..]);
ret.push_str(&arr[..(basepoint + slice_length - scale.len())]);
return ret;
}
}
 
pub fn are_rotationally_equivalent(scale1: &String, scale2: &String)-> bool{
    if scale1.len() != scale2.len() {
        return false;
    } else {
        let length = scale1.len();
        for i in 0..length {
            if slice_cyclic_string(scale1, i, length) == *scale2 {
                return true;
            }
        }
        false
    }
}
 
/// Assuming that "`equiv` : T x T -> bool" is an equivalence relation,
/// return the equivalence class representatives `set` under `equiv`.
pub fn equivalence_class_representatives<T>(set: Vec<T>, equiv: EquivalenceRelation<T>) -> Vec<T>
where T : Clone + PartialEq {
    let mut result : Vec<T> = Vec::new();
    for i in 0..set.len() {
        let mut found_class_equivalent = false;
        'a: for j in 0..result.len(){ // All vecs in result are nonempty.
            if equiv(&result[j], &set[i]){
                found_class_equivalent = true;
                continue 'a;
            }
        }
        if found_class_equivalent == false {
            result.push(set[i].clone());
        }
    }
    result
}
 
/// Gives the greatest common denominator of the two inputs, unless that's 2^63.
/// 2^63 doesn't fit in an `i64`, so it returns -2^63, which does.
pub fn gcd(u: i64, v: i64) -> i64 {
    // `wrapping_abs` gives a number's absolute value, unless that's 2^63. 2^63
    // won't fit in `i64`, so it gives -2^63 instead.
    let mut v = v.wrapping_abs() as u64;
    if u == 0 {
        return v as i64;
    }
    let mut u = u.wrapping_abs() as u64;
    if v == 0 {
        return u as i64;
    }
 
    // `|` is bitwise OR. `trailing_zeros` quickly counts a binary number's
    // trailing zeros, giving its prime factorization's exponent on two.
    let gcd_exponent_on_two = (u | v).trailing_zeros();
 
    // `>>=` divides the left by two to the power of the right, storing that in
    // the left variable. `u` divided by its prime factorization's power of two
    // turns it odd.
    u >>= u.trailing_zeros();
    v >>= v.trailing_zeros();
 
    while u != v {
        if u < v {
            // Swap the variables' values with each other.
            core::mem::swap(&mut u, &mut v);
        }
        u -= v;
        u >>= u.trailing_zeros();
    }
 
    // `<<` multiplies the left by two to the power of the right.
    (u << gcd_exponent_on_two) as i64
}
 
fn lcm(m: i64, n: i64) -> i64 {
m*n/gcd(m,n)
}
 
/// Computes the maximum variety of a ternary scale word.
pub fn maximum_variety(scale: &str) -> usize {
    let mut result = 0;
let floor_half: usize = scale.len()/2;
for slice_length in 1..(floor_half+1) {
let mut sizes: BTreeSet<(usize, usize, usize)> = BTreeSet::new();
for basepoint in 0..(scale.len()) {
let mut size: (usize, usize, usize) = (0, 0, 0);
let sl = &slice_cyclic_string(scale, basepoint, slice_length);
let chars_in_sl = sl.chars();
for ch in chars_in_sl {
match ch {
'x' => {size.0 += 1;},
'y' => {size.1 += 1;},
'z' => {size.2 += 1;},
_ => {panic!()},
}
}
sizes.insert(size);
}
        result = max(result, sizes.len()); // update result
}
result
}
 
/// Return the balance of `s` = `s`(x, y, z), defined as max { | |w|_{x_i} - |w'|_{x_i} | : x_i is a letter of `s` and k = len(w) = len(w') }.
pub fn balance(scale: &str) -> usize {
    if scale.len() <= 1 {
        scale.len()
    } else {
        let mut result = 0;
    let floor_half: usize = scale.len()/2;
    for slice_length in 1..(floor_half+1) {
            let (mut min_x, mut max_x, mut min_y, mut max_y, mut min_z, mut max_z) = (usize::MAX, 0, usize::MAX, 0, usize::MAX, 0); // min and max for each length
    let mut sizes: BTreeSet<(usize, usize, usize)> = BTreeSet::new();
    for basepoint in 0..(scale.len()) {
    let mut size: (usize, usize, usize) = (0, 0, 0);
    let sl = &slice_cyclic_string(scale, basepoint, slice_length);
    let chars_in_sl = sl.chars();
    for ch in chars_in_sl {
    match ch {
    'x' => {size.0 += 1;},
    'y' => {size.1 += 1;},
    'z' => {size.2 += 1;},
    _ => {panic!()},
    }
    }
    sizes.insert(size);
    }
            for vector in &sizes {
                min_x = min(min_x, vector.0);
                max_x = max(max_x, vector.0);
                min_y = min(min_y, vector.1);
                max_y = max(max_y, vector.1);
                min_z = min(min_z, vector.2);
                max_z = max(max_z, vector.2);
            }
            result = max(result, max(max(max_x.abs_diff(min_x), max_y.abs_diff(min_y)), max_z.abs_diff(min_z)));
    }
        result
    }
}
 
/// Function for computing two properties, maximum variety and balance, while executing the subroutine to extract words only once at a time for both properties.
pub fn max_variety_and_balance(scale: &str) -> (usize, usize) {
    if scale.len() <= 1 {
        (scale.len(), 0)
    } else {
        let mut mv_result = 0;
        let mut balance_result = 0;
        let (mut min_x, mut max_x, mut min_y, mut max_y, mut min_z, mut max_z) = (usize::MAX, 0, usize::MAX, 0, usize::MAX, 0);
    let floor_half: usize = scale.len()/2;
    for slice_length in 1..(floor_half+1) {
    let mut sizes: BTreeSet<(usize, usize, usize)> = BTreeSet::new();
    for basepoint in 0..scale.len() {
    let mut size: (usize, usize, usize) = (0, 0, 0);
    let sl = &slice_cyclic_string(scale, basepoint, slice_length);
    let chars_in_sl = sl.chars();
    for ch in chars_in_sl {
    match ch {
    'x' => {size.0 += 1;},
    'y' => {size.1 += 1;},
    'z' => {size.2 += 1;},
    _ => {panic!()},
    }
    }
    sizes.insert(size);
    }
            for vector in &sizes {
                min_x = min(min_x, vector.0);
                max_x = max(max_x, vector.0);
                min_y = min(min_y, vector.1);
                max_y = max(max_y, vector.1);
                min_z = min(min_z, vector.2);
                max_z = max(max_z, vector.2);
            }
            mv_result = max(mv_result, sizes.len()); // update result
            balance_result = max(balance_result, max(max(max_x.abs_diff(min_x), max_y.abs_diff(min_y)), max_z.abs_diff(min_z)));
    }
    (mv_result, balance_result)
    }
}
 
/// Billiard scales with the given signature up to rotation.
pub fn billiard_scales(a: u32, b: u32, c: u32) -> Vec<String> {
    let d = gcd(gcd(a as i64, b as i64), c as i64) as u32; // take gcd(a, b, c)
    let mut result = Vec::<String>::new();
    let regions: Vec<ConvexPolygon> = plane_geometry::project_and_partition(a/d, b/d, c/d);
    let n = a/d + b/d + c/d;
    for region in regions {
        let mut word = String::new();
        let mut next_point = region.centroid();
        let first_point = region.centroid().clone();
        for i in 0..n {
            let res = plane_geometry::advance(a, b, c, next_point).unwrap();
            next_point = res.0;
            debug_assert!(i == n-1 || first_point != next_point); // There should not be a subperiod in the orbit of the "billiard ball".
            word.push_str(&(res.1));
        }
        result.push(std::iter::repeat(word).take(d.try_into().unwrap()).collect::<String>()); // Each word should be repeated d times.
        debug_assert_eq!(first_point, next_point); // The last point should come back to the first point. When i = n - 1, next_point gets set to point number n.
    }
    result = equivalence_class_representatives::<String>(result, are_rotationally_equivalent);
    result
}
</syntaxhighlight>
 
=== main.rs ===
<syntaxhighlight lang="rust">
use billiard_scales;
use std::env;
use std::io;
use std::fs::File;
use std::io::Write;
 
fn main() -> std::io::Result<()> {
    env::set_var("RUST_BACKTRACE", "1");
    let mut input = String::new();
    let mut f = File::create("billiard_scales_mediawiki")?;
    println!("This script will generate the list of all ternary billiard scales of the form ax by cz.");
    println!("Enter a > 0");
    std::io::stdin().read_line(&mut input).unwrap();
    let mut result_a: Result<i64, String> = match input.trim().parse::<i64>(){
        Ok(n) => if n > 0 {Ok(n)} else {Err("input must be a positive integer".to_string())}
        Err(e) => Err(e.to_string())
    };
    while let Err(ref e) = result_a {
        input.clear();
        println!("Input must be a positive integer. Please try again.");
        std::io::stdin().read_line(&mut input).unwrap();
        result_a = match input.trim().parse::<i64>(){
            Ok(n) => if n > 0 {Ok(n)} else {Err("input must be a positive integer".to_string())}
            Err(e) => Err(e.to_string())
        };
    }
    let a: u32 = result_a.expect("Unknown error") as u32;
    input.clear();
    println!("Enter b > 0");
    std::io::stdin().read_line(&mut input).unwrap();
    let mut result_b: Result<i64, String> = match input.trim().parse::<i64>(){
        Ok(n) => if n > 0 {Ok(n)} else {Err("input must be a positive integer".to_string())}
        Err(e) => Err(e.to_string())
    };
    while let Err(ref e) = result_b {
        input.clear();
        println!("Input must be a positive integer. Please try again.");
        std::io::stdin().read_line(&mut input).unwrap();
        result_b = match input.trim().parse::<i64>(){
            Ok(n) => if n > 0 {Ok(n)} else {Err("input must be a positive integer".to_string())}
            Err(e) => Err(e.to_string())
        };
    }
    let b: u32 = result_b.expect("Unknown error") as u32;
    input.clear();
    println!("Enter c > 0");
    std::io::stdin().read_line(&mut input).unwrap();
    let mut result_c: Result<i64, String> = match input.trim().parse::<i64>(){
        Ok(n) => if n > 0 {Ok(n)} else {Err("input must be a positive integer".to_string())}
        Err(e) => Err(e.to_string())
    };
    while let Err(ref e) = result_c {
        input.clear();
        println!("Input must be a positive integer. Please try again.");
        std::io::stdin().read_line(&mut input).unwrap();
        result_c = match input.trim().parse::<i64>(){
            Ok(n) => if n > 0 {Ok(n)} else {Err("input must be a positive integer".to_string())}
            Err(e) => Err(e.to_string())
        };
    }
    let c: u32 = result_c.expect("Unknown error") as u32;
    input.clear();
    let d = billiard_scales::gcd(billiard_scales::gcd(a as i64, b as i64), c as i64);
    if d > 1 {
        println!("Note: a, b, and c are not coprime. You will get scale words with {} periods.", d);
    }
    writeln!(&f, "Billiard scales of signature {}x{}y{}z", a, b, c);
    let billiard_scales = billiard_scales::billiard_scales((a as usize).try_into().unwrap(), (b as usize).try_into().unwrap(), (c as usize).try_into().unwrap());
    for scale in &billiard_scales{
        writeln!(&f, "# {}\n#* maximum variety {}, balance {}", scale, billiard_scales::maximum_variety(scale), billiard_scales::balance(scale));
    }
    println!("Done. Output has been saved in `billiard_scales_mediawiki` in the current directory.");
 
    Ok(())
}
</syntaxhighlight>
 
=== Cargo.toml ===
<syntaxhighlight lang="toml">
[package]
name = "billiard_scales"
version = "0.0.0"
edition = "2021"
 
# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
 
[dependencies]
num-rational = "0.4"
num-traits = "0.2.17"
</syntaxhighlight>
 
[[Category:Code]]
Retrieved from "https://en.xen.wiki/w/Code_for_billiard_scales"