User:Inthar/Code

Billiard scales

To run this, you need Rust installed.

1. In any folder you want, create a directory billiard-scales.
2. Using the command line, run cargo init in the new directory billiard-scales.
3. Create or copy the following .rs files in billiard-scales/src. Copy the contents of the Cargo.toml provided into the new Cargo.toml
4. Run cargo run --release.

planegeometry.rs

// plane_geometry.rs
// v0.0.0 by inthar
// 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 std::ops::Add;
use std::ops::Sub;
use std::ops::Neg;
use std::iter::zip;
use num_traits::sign::Signed;
use std::fmt;

use num_rational::Rational64 as r64;

#[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(a:r64, b:r64, c:r64) -> bool {
    if b == c {
        a == b
    } else if b < c {
        b < a && a < c
    } else {
        c < a && a < b
    }
}

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

/// A rational number or unsigned infinity; represents a valid slope for a line in Q^2.
#[derive(Copy, Clone, Debug, Hash)]
pub struct Slope {
    numer: i64,
    denom: i64,
}

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<r64, String> {
        if self.denom != 0 {
            Ok(r64::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) -> r64 {
        r64::new(self.numer, self.denom)
    }
    
    /// Create a new reduced Slope; panics if both numer and denom are 0.
    pub fn new(n: i64, d: i64) -> 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, denom);
            if d > 1 {
                numer = numer/d;
                denom = 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: i64, denom: i64) -> Self {
        Self {
            numer,
            denom,
        }
    }
    
    /// Create a new `Slope` from a rational.
    pub fn rational(rat: r64) -> Self {
        Self {
            numer: *rat.numer(),
            denom: *rat.denom(),
        }
    }
    
    /// Create a new `Slope` from an integer.
    pub fn integer(n: i64) -> 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)]
pub struct Point {
    // The x-coordinate.
    x: r64,
    // The y-coordinate.
    y: r64,
}

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: r64, y: r64) -> Self {
        Self {
            x,
            y,
        }
    }
    
    pub fn zero() -> Self {
        Self {
            x: r64::new(0,1),
            y: r64::new(0,1),
        }
    }
    
    pub fn scalar_mult(lambda: r64, 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(r64::new(1,1)/r64::new(points.len() as i64,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)/r64::new(2,1),
            y: (p1.y+p2.y)/r64::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, Hash)]
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: Line) -> 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: Line) -> Point {
        self.point1
    }
    
    pub fn from_slope_and_point(sl: Slope, point: Point) -> Self {
        Self {
            point1: point,
            point2: point + Point::new(r64::new_raw(sl.denom, 1), r64::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, point2})
        }
    }
    
    /// Test whether a line is vertical.
    fn is_vertical(&self) -> bool {
        self.point1.x == self.point2.x
    }

    /// Test whether a line is horizontal.
    fn is_horizontal(&self) -> bool {
        self.point1.y == self.point2.y
    }
    
    /// Test whether two lines are parallel.
    fn is_parallel(&self, other: Self) -> bool {
        self.slope() == other.slope()
    }
    
    /// Test whether `p` satisfies the equation of the line.
    fn has_point(&self, p: Point) -> bool {
        (p.y - self.point().y)*r64::new_raw(self.slope().denom,1) == r64::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)*r64::new_raw(self.slope().denom,1) > r64::new_raw(self.slope().numer,1)*(p1.x - self.point1.x) &&
                (p2.y - self.point1.y)*r64::new_raw(self.slope().denom,1) > r64::new_raw(self.slope().numer,1)*(p2.x - self.point1.x))
            || // ... or both are below `self`
            ((p1.y - self.point1.y)*r64::new_raw(self.slope().denom,1) < r64::new_raw(self.slope().numer,1)*(p1.x - self.point1.x) &&
                (p2.y - self.point1.y)*r64::new_raw(self.slope().denom,1) < r64::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)*r64::new_raw(self.slope().denom, 1) - r64::new_raw(self.slope().numer, 1)*(p.x - self.point().x)) * (self.point2.x - self.point1.x).signum() > r64::new_raw(0,1) {
                    PointLineConfiguration::Left
                } else if ((p.y - self.point().y)*r64::new_raw(self.slope().denom, 1) - r64::new_raw(self.slope().numer, 1)*(p.x - self.point().x)) * (self.point2.x - self.point1.x).signum() == r64::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) < r64::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) == r64::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> {
    let (x1, y1, x2, y2) = (l1.point().x, l1.point().y, l2.point().x, l2.point().y);
    match (l1.slope(), l2.slope()) {
        (s1, s2) if s1 == s2 => None, // Same slope
        (s1, s2) if s1 == Slope::infinity() && s2 != Slope::infinity() => Some(Point {
            x: x1,
            y: s2.as_rational_raw() * (x1 - x2) + y2,
        }),
        (s1, s2) if s1 != Slope::infinity() && s2 == Slope::infinity() => Some(Point {
            x: x2,
            y: s1.as_rational_raw() * (x2 - x1) + y1,
        }),
        (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, Hash)]
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: Vec<Point>) -> Option<Self> {
        if vertices.len() >= 3 {
            Some(Self {
                vertices: Point::ordered_ccw(vertices), // This will work for vertices of a convex polygon.
            })
        } else {
            None
        }
    }
    
    fn raw(vertices: Vec<Point>) -> Self { // Does not check for ordering.
        Self {
            vertices: Point::ordered_ccw(vertices), // This will work for vertices of a convex polygon.
        }
            
    }
    
    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]) {
            return false;
        } else if !Line::from_points(self.vertices[self.vertices.len()-1], self.vertices[0]).unwrap().on_same_side(point, self.vertices[1]) {
            return false;
        } else {
            return true;
        }
    }
    
    /// 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).into_iter().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 thhe line, and to the right, respectively. The left and right collections are missing two vertices.
        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) {
                    if is_between(intsn.x, edge.point1.x, edge.point2.x) && is_between(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))
    }
    

    fn triangle_from_lines(l1: Line, l2: Line, l3: Line) -> Option<Self> {
        if !l1.is_parallel(l2) && !l2.is_parallel(l3) && !l3.is_parallel(l1) {
            let p1 = intersection(l1, l2).unwrap();
            let p2 = intersection(l2, l3).unwrap();
            let p3 = intersection(l3, l1).unwrap();
            Self::new(vec![p1, p2, p3])
        } else {
            None
        }
    }
    
    fn parallelogram_from_lines(l1: Line, l2: Line, m1: Line, m2: Line) -> Option<Self> {
        if l1.is_parallel(l2) && l1 != l2 && m1.is_parallel(m2) && m1 != m2 {
            let p1 = intersection(l1, m1).unwrap();
            let p2 = intersection(l1, m2).unwrap();
            let p3 = intersection(l2, m1).unwrap();
            let p4 = intersection(l2, m2).unwrap();
            Self::new(vec![p1, p2, p3, p4])
        } else {
            None
        }
    }
}

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(r64::new(0, 1), r64::new(-(b as i64) + (i as i64), c as i64 ) ) )).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(r64::new((c as i64)-(j as i64), c as i64), r64::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 i64, a as i64), Point::new( r64::new(0, 1), r64::new(-(b as i64)+(k as i64), a as i64))) ).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.
pub 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(r64::new(0,1), r64::new(1,1)),  
                    Point::new(r64::new(1,1), r64::new(1,1)), 
                    Point::new(r64::new(1,1), r64::new(0,1)),
                    Point::new(r64::new(-(a as i64), c as i64), r64::new(-(b as i64), c as i64) + r64::new(1,1)), 
                    Point::new(r64::new(-(a as i64), c as i64), r64::new(-(b as i64), c as i64)), 
                    Point::new(r64::new(-(a as i64), c as i64) + r64::new(1,1), r64::new(-(b as i64), c as i64))];
        ConvexPolygon::new(projected_cube_vertices)
    } else {
        None
    }
}

/// Partition the projected unit cube with the projected constraint planes.
pub fn project_and_partition(a: u32, b: u32, c: u32) -> Vec<ConvexPolygon> {
    let projected_cube = projected_cube(a, b, c).unwrap();
    let lineses = projected_constraint_planes(a, b, c);  // it's a vector of vector 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;
    while 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 = 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 = j+1;
        }
    }
    for polygon in second_shapes {
        let mut j: usize = 0;
        let mut remaining3: ConvexPolygon = polygon;
        while 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 = 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 "x", "y", and "z".
pub fn advance(a: u32, b: u32, c: u32, pt: Point) -> Result<(Point, String), String> {
    let projected_cube = projected_cube(a, b, c).unwrap();
    // Consider: _ / 
    //            |
    if !projected_cube.has_point_inside(pt) {
        Err(format!("Point passed to `advance`, {}, is not inside the projected cube", pt))
    } else { 
        let projected_x_y_equals_1_1 = Line::from_slope_and_point(Slope::raw(b as i64, a as i64), Point::new( r64::new(0, 1), r64::new(-(b as i64)+(a as i64), a as i64))); // the /
        let projected_x_z_equals_1_1 = Line::from_slope_and_point(Slope::raw(0, 1), Point::new(r64::new(0, 1), r64::new(-(b as i64)+(c as i64), c as i64))); // the _, oriented right
        let projected_y_z_equals_1_1 = Line::from_slope_and_point(Slope::raw(1, 0), Point::new(r64::new((c as i64)-(a as i64), c as i64), r64::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(r64::new(1,1), r64::new(0,1)), "x".to_string())) // 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(r64::new(0,1), r64::new(1,1)), "y".to_string())) // 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(r64::new((a as i64), (c as i64)), r64::new((b as i64), (c as i64) ) ), "z".to_string())) // new point is one unit below in the z-direction, and e3 is projected to (-a/c, -b/c)
        } else {
            Err(format!("Point passed to `advance`, {}, is inside the projected cube but falls on a projected constraint plane", pt))
        }
    }
}

lib.rs

// scaletheory.rs
// Inthar's higher-rank scale theory library
// v0.0.0
use std::collections::BTreeSet;
use std::collections::HashSet;
use num_rational::Rational64 as r64;
mod plane_geometry;
use plane_geometry::*;
use std::cmp::{min, max};

type JIScaleProperty = fn(s: [r64]) -> bool;
type CentsScaleProperty = fn(s: [f64]) -> bool;
type AbstractScaleProperty = fn(s: &str) -> bool;
type EquivalenceRelation<T> = fn(&T, &T) -> bool;

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

// Assumes that equiv : T x T -> bool is an equivalence relation.
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
}

fn gcd(m: i64, n: i64) -> i64 {
	let mut x = m;
	let mut y = n;
	while x != 0 && y != 0 {
		if x > y {
			x = x % y;
		} else {
			y = y % x;
		}
	}
	if y == 0 {x} else {y}
}

fn gcd_u32(m: u32, n: u32) -> u32 {
	let mut x = m;
	let mut y = n;
	while x != 0 && y != 0 {
		if x > y {
			x = x % y;
		} else {
			y = y % x;
		}
	}
	if y == 0 {x} else {y}
}

fn lcm(m: i64, n: i64) -> i64 {
	m*n/gcd(m,n)
}

// Is the scale word s = s(x, y, z) abstractly mv3?
pub fn maximum_variety(s: &str) -> usize {
    let mut result = 0;
	let floor_half: usize = s.len()/2;
	for l in 1..(floor_half+1) {
		let mut sizes: BTreeSet<(usize, usize, usize)> = BTreeSet::new();
		for b in 0..(s.len()) {
			let mut size: (usize, usize, usize) = (0, 0, 0);
			let sl = &slice_cyclic_string(s, b, l);
			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
}

// Is the scale word s = s(x, y, z) pairwise well-formed (pwf)?
pub fn is_pwf(s: &str) -> bool {
	let s1 = s.replace("z", "y");
	let s2 = s.replace("z", "x");
	let s3 = s.replace("y", "x").replace("z", "y");
	return is_sv2(&s1)&&is_sv2(&s2)&&is_sv2(&s3);
}

// Is the scale word s = s(a,b,c) pairwise well-formed (pwf)?
pub fn is_pmos(s: &str) -> bool {
	let s1 = s.replace("z", "y");
	let s2 = s.replace("z", "x");
	let s3 = s.replace("y", "x").replace("z", "y");
	maximum_variety(&s1) == 2 && maximum_variety(&s2) == 2  && maximum_variety(&s3) == 2
}

// Is the scale word s = s(x, y, z), x > y > z > 0, monotone-mos?
pub fn is_mmos(scale: &str) -> bool {
	let s1 = scale.replace("y", "x").replace("z", "y"); // m = L
	let s2 = scale.replace("z", "y"); // m = s
	let s3 = scale.replace("z", "");  // s = 0
	maximum_variety(&s1) == 2 && maximum_variety(&s2) == 2  && maximum_variety(&s3) == 2
}

// Is the scale s = s(x, y) a single period mos (aka well-formed, wf)?
pub fn is_sv2(scale: &str) -> bool {
	if scale.len() <= 1{
		return false;
	}
	let floor_half : usize = scale.len()/2;
	for slice_length in 1..(floor_half+1) {
		let mut sizes: BTreeSet<(usize, usize)> = BTreeSet::new();
		for basepoint in 0..(scale.len()) {
			let mut size: (usize, usize) = (0, 0);
			let sl = &slice_cyclic_string(scale, basepoint, slice_length);
			let chars_in_sl = sl.chars();
			for c in chars_in_sl {
				match c {
					'x' => {size.0 += 1;},
					'y' => {size.1 += 1;},
					_ => {panic!()},
				}
			}
			sizes.insert(size);
		}
		if sizes.len() != 2 {
			return false;
		}
	}
	return true;
}

// 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. 
pub fn max_variety_and_balance(s: &str) -> (usize, usize) {
    if s.len() <= 1 {
        (s.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 = s.len()/2;
	    for l in 1..(floor_half+1) {
		    let mut sizes: BTreeSet<(usize, usize, usize)> = BTreeSet::new();
		    for b in 0..s.len() {
			    let mut size: (usize, usize, usize) = (0, 0, 0);
			    let sl = &slice_cyclic_string(s, b, l);
			    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)
    }
}


// Return a Christoffel word with `a` x's and `b` y's.
// Algorithm from Bulgakova et al, 2023, "On balanced and abelian properties of circular words over a ternary alphabet".
pub fn christoffel_word(a: u32, b: u32) -> String {
    let d = gcd_u32(a, b);
    if d == 1 {
        let mut result : String = String::from("");
        let (mut current_x, mut current_y) = (0u32, 0u32); // Start from the (0,0) vector.
        while current_x < a || current_y < b {
            if (current_y+1)*a <= b * (current_x) { // if making the (0,1) step doesn't lead to going above the line
                current_y += 1; // append the b step and reflect that in the plane vector.
                result.push('y');
            } else {
                current_x += 1;
                result.push('x');
            }
        }
        result
    } else {
        std::iter::repeat(christoffel_word(a/d, b/d)).take(d.try_into().unwrap()).collect::<String>()
    }
}

/// Billiard scales with the given signature up to rotation.
pub fn billiard_scales(a: u32, b: u32, c: u32) -> Vec<String> {
    let mut result = Vec::<String>::new();
    let regions: Vec<ConvexPolygon> = plane_geometry::project_and_partition(a, b, c);
    let n = a + b + c;
    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(word);
        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
}

main.rs

use scale_theory;
use std::env;
use std::fs::File;
use std::io::Write;
use std::time::Instant;

/// 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 main() -> std::io::Result<()> {	
    env::set_var("RUST_BACKTRACE", "1");
    let mut f = File::create("billiard_scales_mediawiki")?;
    let start_time = Instant::now();
    for a in 1..=29 {
        for b in 1..=a {
            for c in 1..=b {
                if a + b + c <= 31 && gcd(gcd(a as i64, b as i64), c as i64) == 1 {
                    writeln!(&f, "== {}x{}y{}z ==", a, b, c);
                    let billiard_scales = scale_theory::billiard_scales(a, b, c);
                    for scale in &billiard_scales{
                        writeln!(&f, "# {}\n#* maximum variety {}, balance {}", scale, scale_theory::maximum_variety(scale), scale_theory::balance(scale));
                    }
                }
            }
        }
    }
    let elapsed_time = start_time.elapsed();
    println!("Done in {:?}", elapsed_time);

    Ok(())
}

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