Xenharmonic Wiki

Code for billiard scales: Difference between revisions

← Older edit
VisualWikitext
BudjarnLambeth (talk | contribs)
autopatrolled, emailconfirmed
21,264 edits
Categorised uncategorised page
Inthar (talk | contribs)
autopatrolled, emailconfirmed
15,668 edits
 
(5 intermediate revisions by the same user not shown)
Line 10: Line 10:
// Geometry of points and lines in the plane, implemented with rational numbers. The implementation is janky, though.
// 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.
// 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::Add;
use std::ops::Neg;
use std::ops::Sub;
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;
use num_rational::Rational32 as r32;
 
use crate::helpers::gcd;


#[derive(Copy, Clone, Debug, Hash, PartialEq, Eq)]
#[derive(Copy, Clone, Debug, Hash, PartialEq, Eq)]
/// The possible configurations of a point and an oriented line on the plane.
/// The possible configurations of a point and an oriented line on the plane.
pub enum PointLineConfiguration {
pub enum PointLineConfiguration {
     Left, // The point is on the left half-plane.
     Left,     // The point is on the left half-plane.
     Right, // The point is on the right half-plane.
     Right,     // The point is on the right half-plane.
     OnTheLine, // The point is on the line.
     OnTheLine, // The point is on the line.
}
}


fn is_between(a:r64, b:r64, c:r64) -> bool {
fn is_between_r32(a: r32, b: r32, c: r32) -> bool {
     if b == c {
     if b == c {
         a == b
         a == b
Line 35: Line 37:
         c < a && a < b
         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.
/// A rational number or unsigned infinity; represents a valid slope for a line in Q^2.
#[derive(Copy, Clone, Debug, Hash)]
#[derive(Copy, Clone, Debug)]
pub struct Slope {
pub struct Slope {
     numer: i64,
     numer: i32,
     denom: i64,
     denom: i32,
}
}


Line 95: Line 60:
impl Slope {
impl Slope {
     /// Converts a `Slope` to a rational; returns Err if `denom` is 0.
     /// Converts a `Slope` to a rational; returns Err if `denom` is 0.
     pub fn as_rational(&self) -> Result<r64, String> {
     pub fn as_rational(&self) -> Result<r32, String> {
         if self.denom != 0 {
         if self.denom != 0 {
             Ok(r64::new(self.numer, self.denom))
             Ok(r32::new(self.numer, self.denom))
         } else {
         } else {
             Err("Attempted to convert a Slope with `denom` == 0 into a Rational64".to_string())
             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.
     /// Converts the projective rational to a rational; panics if the denominator is zero.
     pub fn as_rational_raw(&self) -> r64 {
     pub fn as_rational_raw(&self) -> r32 {
         r64::new(self.numer, self.denom)
         r32::new(self.numer, self.denom)
     }
     }
   
 
     /// Create a new reduced Slope; panics if both numer and denom are 0.
     /// Create a new reduced Slope; panics if both numer and denom are 0.
     pub fn new(n: i64, d: i64) -> Result<Self, String> {
     pub fn new(n: i32, d: i32) -> Result<Self, String> {
         let (mut numer, mut denom) = (n, d);
         let (mut numer, mut denom) = (n, d);
         if numer == 0 && denom == 0 {
         if numer == 0 && denom == 0 {
Line 118: Line 83:
                 denom = -denom;
                 denom = -denom;
             }
             }
             let d = gcd(numer, denom);
             let d = gcd(numer as u32, denom as u32) as i32;
             if d > 1 {
             if d > 1 {
                 numer = numer/d;
                 numer /= d;
                 denom = denom/d;
                 denom /= d;
             }
             }
             Ok(Self {
             Ok(Self { numer, denom })
                numer,
                denom,
            })
         }
         }
     }
     }
   
 
     /// Create a new Slope without checking that `numer` == `denom` == 0, or reducing the rational afterwards.
     /// Create a new Slope without checking that `numer` == `denom` == 0, or reducing the rational afterwards.
     pub fn raw(numer: i64, denom: i64) -> Self {
     pub fn raw(numer: i32, denom: i32) -> Self {
         Self {
         Self { numer, denom }
            numer,
            denom,
        }
     }
     }
   
 
     /// Create a new `Slope` from a rational.
     /// Create a new `Slope` from a rational.
     pub fn rational(rat: r64) -> Self {
     pub fn rational(rat: r32) -> Self {
         Self {
         Self {
             numer: *rat.numer(),
             numer: *rat.numer(),
Line 145: Line 104:
         }
         }
     }
     }
   
 
     /// Create a new `Slope` from an integer.
     /// Create a new `Slope` from an integer.
     pub fn integer(n: i64) -> Self {
     pub fn integer(n: i32) -> Self {
         // use 'raw' since an integer is always reduced
         // use 'raw' since an integer is always reduced
         Self::raw(n, 1)
         Self::raw(n, 1)
     }
     }
   
 
     /// Shorthand for slope 0.
     /// Shorthand for slope 0.
     pub fn zero() -> Self {
     pub fn zero() -> Self {
Line 161: Line 120:
         Self::raw(1, 1)
         Self::raw(1, 1)
     }
     }
   
 
     /// Shorthand for infinite slope.
     /// Shorthand for infinite slope.
     pub fn infinity() -> Self {
     pub fn infinity() -> Self {
Line 167: Line 126:
     }
     }
}
}


/// A finite point; a member of Q^2.
/// A finite point; a member of Q^2.
#[derive(Copy, Clone, Debug, Hash, PartialEq)]
#[derive(Copy, Clone, Debug, Hash, PartialEq, Eq)]
pub struct Point {
pub struct Point {
     // The x-coordinate.
     // The x-coordinate.
     x: r64,
     x: r32,
     // The y-coordinate.
     // The y-coordinate.
     y: r64,
     y: r32,
}
}


Line 189: Line 147:
}
}


impl Sub for Point{
impl Sub for Point {
     type Output = Self;
     type Output = Self;


Line 200: Line 158:
}
}


impl Neg for Point{
impl Neg for Point {
     type Output = Self;
     type Output = Self;


Line 212: Line 170:


impl Point {
impl Point {
     pub fn new(x: r64, y: r64) -> Self {
     pub fn new(x: r32, y: r32) -> Self {
         Self {
         Self { x, y }
            x,
            y,
        }
     }
     }
   
 
     pub fn zero() -> Self {
     pub fn zero() -> Self {
         Self {
         Self {
             x: r64::new(0,1),
             x: r32::new(0, 1),
             y: r64::new(0,1),
             y: r32::new(0, 1),
         }
         }
     }
     }
   
 
     pub fn scalar_mult(lambda: r64, v: Point) -> Self {
     pub fn scalar_mult(lambda: r32, v: Point) -> Self {
         Self {
         Self {
             x: lambda*v.x,
             x: lambda * v.x,
             y: lambda*v.y,
             y: lambda * v.y,
         }
         }
     }
     }


     pub fn average(points: &Vec<Self>) -> Self {
     pub fn average(points: &Vec<Self>) -> Self {
Line 239: Line 193:
             sum = sum + *p;
             sum = sum + *p;
         }
         }
         Self::scalar_mult(r64::new(1,1)/r64::new(points.len() as i64,1), sum)
         Self::scalar_mult(r32::new(1, 1) / r32::new(points.len() as i32, 1), sum)
     }
     }
   
 
     pub fn theta(v: Point) -> f64 {
     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 )
         (*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].
     /// 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> {
     fn ordered_ccw(mut points: Vec<Self>) -> Vec<Self> {
         let centroid = Self::average(&points);
         let centroid = Self::average(&points);
         points.sort_by(|a,b| Self::theta(*a-centroid).partial_cmp(&Self::theta(*b-centroid)).unwrap() );
         points.sort_by(|a, b| {
            Self::theta(*a - centroid)
                .partial_cmp(&Self::theta(*b - centroid))
                .unwrap()
        });
         points
         points
     }
     }
   
 
     pub fn midpoint(p1: &Self, p2: &Self) -> Self {
     pub fn midpoint(p1: &Self, p2: &Self) -> Self {
         Self {
         Self {
             x: (p1.x+p2.x)/r64::new(2,1),
             x: (p1.x + p2.x) / r32::new(2, 1),
             y: (p1.y+p2.y)/r64::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.
     // Use the slope-intercept formula including the case where the slope is infinite.
     pub fn are_collinear(p1: Self, p2: Self, p3: Self) -> bool {
     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)
         (p1.y - p3.y) * (p1.x - p2.x) == (p1.y - p2.y) * (p1.x - p3.x)
     }
     }
}
}
Line 272: Line 230:


/// An oriented line in Q^2 of rational or infinite slope.
/// An oriented line in Q^2 of rational or infinite slope.
#[derive(Copy, Clone, Debug, Hash)]
#[derive(Copy, Clone, Debug)]
pub struct Line {
pub struct Line {
     point1: Point,
     point1: Point,
Line 285: Line 243:


impl Line {
impl Line {
     pub fn slope(self: Line) -> Slope {
     pub fn slope(&self) -> Slope {
         Slope::new(
         Slope::new(
             (*(self.point2.y-self.point1.y).numer()) * (*(self.point2.x-self.point1.x).denom()),  
             (*(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()
             (*(self.point2.x - self.point1.x).numer()) * (*(self.point2.y - self.point1.y).denom()),
        )
        .unwrap()
     }
     }
   
 
     pub fn point(self: Line) -> Point {
     pub fn point(&self) -> Point {
         self.point1
         self.point1
     }
     }
   
 
     pub fn from_slope_and_point(sl: Slope, point: Point) -> Self {
     pub fn from_slope_and_point(sl: &Slope, point: &Point) -> Self {
         Self {
         Self {
             point1: point,
             point1: *point,
             point2: point + Point::new(r64::new_raw(sl.denom, 1), r64::new_raw(sl.numer, 1)),
             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> {
     pub fn from_points(point1: &Point, point2: &Point) -> Result<Self, String> {
         if point1 == point2 {
         if point1 == point2 {
             Err("The points {point1} and {point2} are equal".to_string())
             Err("The points {point1} and {point2} are equal".to_string())
         } else {
         } else {
             Ok(Line{point1, point2})
             Ok(Line {
                point1: *point1,
                point2: *point2,
            })
         }
         }
     }
     }
   
 
     /// Test whether a line is vertical.
     /// Test whether a line is vertical.
     fn is_vertical(&self) -> bool {
     fn is_vertical(&self) -> bool {
Line 314: Line 277:
     }
     }


    /// 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.
     /// Test whether `p` satisfies the equation of the line.
     fn has_point(&self, p: Point) -> bool {
     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)
         (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.
     /// Determine whether `p1` and `p2` are on the same side of the line.
     fn on_same_side(&self, p1: Point, p2: Point) -> bool {
     fn on_same_side(&self, p1: Point, p2: Point) -> bool {
         if self.is_vertical() {
         if self.is_vertical() {
             (p1.x < self.point1.x && p2.x < self.point1.x ) || (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.x > self.point1.x && p2.x > self.point1.x)
             ((p1.y - self.point1.y)*r64::new_raw(self.slope().denom,1) > r64::new_raw(self.slope().numer,1)*(p1.x - self.point1.x) &&
         } else {
                 (p2.y - self.point1.y)*r64::new_raw(self.slope().denom,1) > r64::new_raw(self.slope().numer,1)*(p2.x - self.point1.x))
            // 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`
             || // ... 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) &&
             ((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)*r64::new_raw(self.slope().denom,1) < r64::new_raw(self.slope().numer,1)*(p2.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.
     /// 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 {
     pub fn point_line_config(&self, p: Point) -> PointLineConfiguration {
         if self.point2.x != self.point1.x { // For non-vertical lines
         if self.point2.x != self.point1.x {
                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) {
            // For non-vertical lines
                    PointLineConfiguration::Left
            if ((p.y - self.point().y) * r32::new_raw(self.slope().denom, 1)
                } 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) {
                - r32::new_raw(self.slope().numer, 1) * (p.x - self.point().x))
                    PointLineConfiguration::OnTheLine
                * (self.point2.x - self.point1.x).signum()
                } else {
                > r32::new_raw(0, 1)
                    PointLineConfiguration::Right
            {
                }
                PointLineConfiguration::Left
         } else { // For vertical lines
            } else if ((p.y - self.point().y) * r32::new_raw(self.slope().denom, 1)
                // If point2.y < point1.y, then the line is oriented downwards, so multiply by -1 to compensate.  
                - r32::new_raw(self.slope().numer, 1) * (p.x - self.point().x))
             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.
                * (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
                 PointLineConfiguration::Left
             } else if (p.x-self.point().x) * (self.point2.y-self.point1.y) == r64::new_raw(0,1) {
             } else if (p.x - self.point().x) * (self.point2.y - self.point1.y) == r32::new_raw(0, 1)
            {
                 PointLineConfiguration::OnTheLine
                 PointLineConfiguration::OnTheLine
             } else {
             } else {
Line 364: Line 332:
     }
     }
}
}


impl fmt::Display for Line {
impl fmt::Display for Line {
Line 373: Line 340:


/// For two non-parallel ones, find the intersection; for two parallel lines, return `None`.
/// For two non-parallel ones, find the intersection; for two parallel lines, return `None`.
pub fn intersection(l1: Line, l2: Line) -> Option<Point> {
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);
     let (x1, y1, x2, y2) = (l1.point().x, l1.point().y, l2.point().x, l2.point().y);
     match (l1.slope(), l2.slope()) {
     match (l1.slope(), l2.slope()) {
         (s1, s2) if s1 == s2 => None, // Same 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 {
         (s1, s2) if s1 == Slope::infinity() && s2 != Slope::infinity() => Some(Point {
             x: x1,
             x: x1,
             y: s2.as_rational_raw() * (x1 - x2) + y2,
             y: s2.as_rational_raw() * (x1 - x2) + y2,
         }),
         }),
        // Unequal slopes, s2 is infinite
         (s1, s2) if s1 != Slope::infinity() && s2 == Slope::infinity() => Some(Point {
         (s1, s2) if s1 != Slope::infinity() && s2 == Slope::infinity() => Some(Point {
             x: x2,
             x: x2,
             y: s1.as_rational_raw() * (x2 - x1) + y1,
             y: s1.as_rational_raw() * (x2 - x1) + y1,
         }),
         }),
        // Unequal slopes, neither infinite
         (s1, s2) => Some(Point {
         (s1, s2) => Some(Point {
             x: (y1 - s1.as_rational_raw() * x1 + s2.as_rational_raw() * x2 - y2)
             x: (y1 - s1.as_rational_raw() * x1 + s2.as_rational_raw() * x2 - y2)
                 / (s2.as_rational_raw() - s1.as_rational_raw()),
                 / (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)
             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()),
                 / (s2.as_rational_raw() - s1.as_rational_raw()),
         }),
         }),
Line 394: Line 368:
}
}


#[derive(Clone, Debug, Hash)]
#[derive(Clone, Debug)]
pub struct ConvexPolygon {
pub struct ConvexPolygon {
     vertices: Vec<Point>, // Assumes that the vertices are ordered in clockwise or counterclockwise order.
     vertices: Vec<Point>, // Assumes that the vertices are ordered in clockwise or counterclockwise order.
Line 414: Line 388:
}
}


impl ConvexPolygon { // Assumes that the vertices form a convex polygon. Does not check if the `Vec` of vertices consists of distinct elements.
impl ConvexPolygon {
     pub fn new(vertices: Vec<Point>) -> Option<Self> {
    // 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 {
         if vertices.len() >= 3 {
             Some(Self {
             Some(Self {
                 vertices: Point::ordered_ccw(vertices), // This will work for vertices of a convex polygon.
                 vertices: Point::ordered_ccw(vertices.to_vec()), // This will work for vertices of a convex polygon.
             })
             })
         } else {
         } else {
Line 424: Line 399:
         }
         }
     }
     }
   
 
    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 {
     pub fn centroid(&self) -> Point {
         Point::average(&self.vertices)
         Point::average(&self.vertices)
     }
     }
   
 
     pub fn has_point_inside(&self, point: Point) -> bool {
     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?
         // 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 {
         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]) {
             if !Line::from_points(&self.vertices[i], &self.vertices[i + 1])
                .unwrap()
                .on_same_side(point, self.vertices[i + 2])
            {
                 return false;
                 return false;
             }
             }
         }
         }
         // Check last two triples.
         // 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]) {
         if !Line::from_points(
            return false;
            &self.vertices[self.vertices.len() - 2],
         } else if !Line::from_points(self.vertices[self.vertices.len()-1], self.vertices[0]).unwrap().on_same_side(point, self.vertices[1]) {
            &self.vertices[self.vertices.len() - 1],
             return false;
        )
        .unwrap()
        .on_same_side(point, self.vertices[0])
         {
             false
         } else {
         } else {
             return true;
             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.
     /// A vector storing the edges of the polygon in order.
     fn lines_for_edges(&self) -> Vec<Line> {
     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();
         let mut result: Vec<Line> = (0..self.vertices.len() - 1)
         result.push(Line::from_points(self.vertices[self.vertices.len()-1], self.vertices[0]).unwrap());
            .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
         result
     }
     }
   
 
     /// Subdivide a polygon with the given line and return the two pieces.
     /// Subdivide a polygon with the given line and return the two pieces.
     pub fn subdivide(&self, line: Line) -> (Option<ConvexPolygon>, Option<ConvexPolygon>) {
     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.
         // Collect vertices that fall to the left, on the line, and to the right, respectively.
         let mut vertices_left: Vec<Point> = self.vertices
        // The left and right collections are missing two vertices; we need to append them to make them `ConvexPolygon`s.
                .iter()
         let mut vertices_left: Vec<Point> = self
                .filter(|&x| line.point_line_config(*x) == PointLineConfiguration::Left).cloned().collect();
            .vertices
         let mut vertices_right: Vec<Point> = self.vertices
            .iter()
                .iter()
            .filter(|&x| line.point_line_config(*x) == PointLineConfiguration::Left)
                .filter(|&x| line.point_line_config(*x) == PointLineConfiguration::Right).cloned().collect();
            .cloned()
         let vertices_on_line: Vec<Point> = self.vertices
            .collect();
                .iter()
         let mut vertices_right: Vec<Point> = self
                .filter(|&x| line.point_line_config(*x) == PointLineConfiguration::OnTheLine).cloned().collect();
            .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() {
         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.
             // Add any vertices on the line both to the new left polygon and the new right polygon.
Line 478: Line 469:
                 vertices_right.push(*point);
                 vertices_right.push(*point);
             }
             }
             // Look for intersection points on edges with the line.
             // Look for intersection points on edges with the line.
             for edge in &self.lines_for_edges() {
             for edge in &self.lines_for_edges() {
                 if let Some(intsn) = intersection(*edge, line) {
                 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) {
                     && is_between_r32(intsn.x, edge.point1.x, edge.point2.x)
                        vertices_left.push(intsn);
                    && is_between_r32(intsn.y, edge.point1.y, edge.point2.y)
                        vertices_right.push(intsn);
                {
                    }
                    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.
             // 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))
         (
            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
//! # 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};
 
pub type Letter = usize;
 
// ERRORS


    fn triangle_from_lines(l1: Line, l2: Line, l3: Line) -> Option<Self> {
#[derive(Copy, Clone, PartialEq, Eq, Debug)]
        if !l1.is_parallel(l2) && !l2.is_parallel(l3) && !l3.is_parallel(l1) {
/// Error type for illegal billiard trajectory states.
             let p1 = intersection(l1, l2).unwrap();
pub enum BadBilliardState {
            let p2 = intersection(l2, l3).unwrap();
    // Projected cube is invalid.
            let p3 = intersection(l3, l1).unwrap();
    ProjectedCubeInvalid((u32, u32, u32)),
             Self::new(vec![p1, p2, p3])
    // Point is not inside the projected cube,
        } else {
    PointNotInCube(Point),
             None
    // 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
             ),
         }
         }
     }
     }
   
}
    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 {
impl std::error::Error for BadBilliardState {}
             let p1 = intersection(l1, m1).unwrap();
 
             let p2 = intersection(l1, m2).unwrap();
/// Rotate a scale word left by `degree` positions (change mode).
             let p3 = intersection(l2, m1).unwrap();
pub fn rotate<T: std::clone::Clone>(slice: &[T], degree: usize) -> Vec<T> {
             let p4 = intersection(l2, m2).unwrap();
    let degree = degree % slice.len();
             Self::new(vec![p1, p2, p3, p4])
    if degree == 0 {
         slice.to_vec()
    } else {
        [&slice[degree..slice.len()], &slice[0..degree]].concat()
    }
}
 
/// The lexicographically least mode of a word (where the letters are in their usual order).
pub fn least_mode(scale: &[Letter]) -> Vec<Letter> {
    rotate(scale, booth(scale))
}
 
/// The rotation required from the current word to the
/// lexicographically least mode of a word.
/// Booth's algorithm requires at most 3*n* comparisons and *n* storage locations where *n* is the input word's length.
/// See Booth, K. S. (1980). Lexicographically least circular substrings.
/// Information Processing Letters, 10(4-5), 240–242. doi:10.1016/0020-0190(80)90149-0
pub fn booth(scale: &[Letter]) -> usize {
    let scale_len = scale.len();
    // `failure_func` is the failure function of the least rotation; `usize::MAX` is used as a null value.
    // null indicates that the failure function does not point backwards in the string.
    // `usize::MAX` will behave the same way as -1 does, assuming wrapping unsigned addition
    let mut failure_func = vec![usize::MAX; 2 * scale_len];
    let mut least_rotation: usize = 0;
    // `scan_pos` loops over `scale` twice.
    for scan_pos in 1..2 * scale_len {
        let mut match_len = failure_func[scan_pos - least_rotation - 1];
        while match_len != usize::MAX
            && scale[scan_pos % scale_len]
                != scale[least_rotation.wrapping_add(match_len).wrapping_add(1) % scale_len]
        {
             // (1) If the scan_pos-th letter is less than s[(least_rotation + match_len + 1) % scale_len] then change least_rotation to scan_pos - match_len - 1,
            // in effect left-shifting the failure function and the input string.
            // This appropriately compensates for the new, shorter least substring.
            if scale[scan_pos % scale_len]
                < scale[least_rotation.wrapping_add(match_len).wrapping_add(1) % scale_len]
             {
                least_rotation = scan_pos.wrapping_sub(match_len).wrapping_sub(1);
             }
            match_len = failure_func[match_len];
        }
        if match_len == usize::MAX
            && scale[scan_pos % scale_len]
                != scale[least_rotation.wrapping_add(match_len).wrapping_add(1) % scale_len]
        {
            // See note (1) above.
             if scale[scan_pos % scale_len]
                < scale[least_rotation.wrapping_add(match_len).wrapping_add(1) % scale_len]
            {
                least_rotation = scan_pos;
             }
            failure_func[scan_pos - least_rotation] = usize::MAX;
         } else {
         } else {
             None
             failure_func[scan_pos - least_rotation] = match_len.wrapping_add(1);
         }
         }
        // The induction hypothesis is that
        // at this point `failure_func[0..scan_pos - least_rotation]` is the failure function of `s[least_rotation..(least_rotation+scan_pos)%scale_len]`,
        // and `least_rotation` is the lexicographically least subword of the letters scanned so far.
     }
     }
    least_rotation
}
}


fn projected_constraint_planes(a: u32, b: u32, c: u32) -> (Vec<Line>, Vec<Line>, Vec<Line>) {
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 first = (0..b + c + 1)
     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)
        .map(|i| {
     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)
            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)
     (first, second, third)
}
}


/// Project the unit cube in R^3 via a linear map that has (a,b,c) in its kernel.
/// 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> {
fn projected_cube(a: u32, b: u32, c: u32) -> Option<ConvexPolygon> {
     if a > 0 && b > 0 && c > 0 {
     if a > 0 && b > 0 && c > 0 {
         let projected_cube_vertices = vec![Point::new(r64::new(0,1), r64::new(1,1)),
         let projected_cube_vertices = vec![
                    Point::new(r64::new(1,1), r64::new(1,1)),  
            Point::new(r32::new(0, 1), r32::new(1, 1)),
                    Point::new(r64::new(1,1), r64::new(0,1)),
            Point::new(r32::new(1, 1), r32::new(1, 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(r32::new(1, 1), r32::new(0, 1)),
                    Point::new(r64::new(-(a as i64), c as i64), r64::new(-(b as i64), c as i64)),  
            Point::new(
                    Point::new(r64::new(-(a as i64), c as i64) + r64::new(1,1), r64::new(-(b as i64), c as i64))];
                r32::new(-(a as i32), c as i32),
         ConvexPolygon::new(projected_cube_vertices)
                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 {
     } else {
         None
         None
Line 540: Line 669:


/// Partition the projected unit cube with the projected constraint planes.
/// Partition the projected unit cube with the projected constraint planes.
pub fn project_and_partition(a: u32, b: u32, c: u32) -> Vec<ConvexPolygon> {
fn project_and_partition(a: u32, b: u32, c: u32) -> Option<Vec<ConvexPolygon>> {
     let projected_cube = projected_cube(a, b, c).unwrap();
     projected_cube(a, b, c).map(|projected_cube| {
    let lineses = projected_constraint_planes(a, b, c); // it's a vector of vector of `Line`s
        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 first_shapes = Vec::<ConvexPolygon>::new();
    let mut second_shapes = Vec::<ConvexPolygon>::new();
        let mut second_shapes = Vec::<ConvexPolygon>::new();
    let mut third_shapes = Vec::<ConvexPolygon>::new();
        let mut third_shapes = Vec::<ConvexPolygon>::new();
    let mut i: usize = 0;
        let mut i: usize = 0;
    let mut remaining: ConvexPolygon = projected_cube;
        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 {
         loop {
             let (maybe_remaining, maybe_next) = remaining2.subdivide((lineses.1)[j]); // The remaining portion is to the Left.
             let (maybe_remaining, maybe_next) = remaining.subdivide(&(lineses.0)[i]);
            // The remaining portion is to the Left.
             match (maybe_remaining, maybe_next) {
             match (maybe_remaining, maybe_next) {
                 (Some(r), None) => {
                 (Some(r), None) => {
                     remaining2 = r;
                     remaining = r;
                 },
                 }
                 (Some(r), Some(next)) => {
                 (Some(r), Some(next)) => {
                     remaining2 = r;
                     remaining = r;
                     second_shapes.push(next);
                     first_shapes.push(next);
                 },
                 }
                 (None, Some(next)) => {
                 (None, Some(next)) => {
                     second_shapes.push(next);
                     // This is the last `next` you need to look at.
                    first_shapes.push(next);
                     break;
                     break;
                 },
                 }
                 (_, _) => {
                 (_, _) => {
                     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;
             }
             }
            j = j+1;
         }
         }
     }
        third_shapes
     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.
/// Advance a projected point, `pt`, by one letter according to what integer coordinate it next hits as it traverses in the given velocity,
             match (maybe_remaining, maybe_next) {
/// and return the next point and the corresponding letter, one of "0", "1", and "2".
                 (Some(r), None) => {
fn advance(a: u32, b: u32, c: u32, pt: Point) -> Result<(Point, Letter), BadBilliardState> {
                    remaining3 = r;
     if let Some(projected_cube) = projected_cube(a, b, c) {
                 },
         // Consider: _ /
                 (Some(r), Some(next)) => {
        //            |
                    remaining3 = r;
         if !projected_cube.has_point_inside(pt) {
                    third_shapes.push(next);
            Err(BadBilliardState::PointNotInCube(pt))
                 },
         } else {
                 (None, Some(next)) => {
            let projected_x_y_equals_1_1 = Line::from_slope_and_point(
                    third_shapes.push(next);
                &Slope::raw(b as i32, a as i32),
                    break;
                &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(
                     break;
                &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))
             }
             }
            j = j+1;
         }
         }
    } else {
        Err(BadBilliardState::ProjectedCubeInvalid((a, b, c)))
     }
     }
    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,
/// Generate all billiard scales with the given step signature.
/// 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> {
/// Returns scales in canonical form (lexicographically first rotation),
     let projected_cube = projected_cube(a, b, c).unwrap();
/// sorted and deduplicated.
     // Consider: _ /
///
    //            |
/// # Arguments
    if !projected_cube.has_point_inside(pt) {
///
        Err(format!("Point passed to `advance`, {}, is not inside the projected cube", pt))
/// * `a` - Number of L (large) steps
     } else {  
/// * `b` - Number of m (medium) steps
         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 /
/// * `c` - Number of s (small) steps
        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
/// # Returns
        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
/// Empty vector if the signature is invalid (any component is 0).
                 && 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
pub fn billiard_scales(a: usize, b: usize, c: usize) -> Vec<Vec<Letter>> {
            Ok((pt - Point::new(r64::new(1,1), r64::new(0,1)), "x".to_string())) // new point is one unit to the left
     let (a, b, c) = (a as u32, b as u32, c as u32);
        } else if projected_x_z_equals_1_1.point_line_config(pt) == PointLineConfiguration::Left
     let d = gcd(gcd(a, b), c);
                 && 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
     if d != 0 {
            Ok((pt - Point::new(r64::new(0,1), r64::new(1,1)), "y".to_string())) // new point is one unit below
         let mut result = Vec::<Vec<Letter>>::new();
        } else if projected_x_z_equals_1_1.point_line_config(pt) == PointLineConfiguration::Right
        let n = a / d + b / d + c / d;
                 && 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
        if let Some(regions) = project_and_partition(a / d, b / d, c / d) {
            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)
            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 {
         } else {
             Err(format!("Point passed to `advance`, {}, is inside the projected cube but falls on a projected constraint plane", pt))
             vec![]
         }
         }
    } else {
        vec![]
     }
     }
}
}
Line 650: Line 863:
mod plane_geometry;
mod plane_geometry;
use plane_geometry::*;
use plane_geometry::*;
mod billiard;
use billiard::*;
use std::cmp::{min, max};
use std::cmp::{min, max};


Line 954: Line 1,169:
name = "billiard_scales"
name = "billiard_scales"
version = "0.0.0"
version = "0.0.0"
edition = "2021"
edition = "2024"


# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
Line 963: Line 1,178:
</syntaxhighlight>
</syntaxhighlight>


[[Category:Software]]
[[Category:Code]]
Retrieved from "https://en.xen.wiki/w/Code_for_billiard_scales"