User:Inthar/Code: Difference between revisions

Inthar (talk | contribs)
mNo edit summary
Inthar (talk | contribs)
No edit summary
Line 1: Line 1:
== planegeometry.rs ==
== 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>.
=== planegeometry.rs ===
<syntaxhighlight lang="rs">// plane_geometry.rs
<syntaxhighlight lang="rs">// plane_geometry.rs
// v0.0 by inthar
// v0.0.0 by inthar
// 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.
Line 637: Line 643:
     }
     }
}
}
</syntaxhighlight>
=== lib.rs ===
<syntaxhighlight lang="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;


#[cfg(test)]
pub fn are_rotationally_equivalent(scale1: &String, scale2: &String)-> bool{
mod tests {
     if scale1.len() != scale2.len() {
     use num_rational::Rational64 as r64;
        return false;
    use crate::plane_geometry::Slope;
     } else {
     use crate::plane_geometry::Point;
        let length = scale1.len();
    use crate::plane_geometry::Line;
        for i in 0..length {
    use crate::plane_geometry::intersection;
            if slice_cyclic_string(scale1, i, length) == *scale2 {
    use crate::plane_geometry::ConvexPolygon;
                return true;
    use crate::plane_geometry::PointLineConfiguration;
             }
    use crate::plane_geometry::projected_cube;
        }
   
        false
    // Tests for the struct `Slope`.
    #[test]
    fn slope_bad_slope() {
        let bad_slope : Result<Slope, String> = Slope::new(0,0);
             assert_eq!(bad_slope,
                    Err("Attempted to create a Slope with both `numer` and `denom` equal to 0".to_string()));
     }
     }
   
}
    #[test]
 
    fn slope_infinite_slope() {
// Assumes that equiv : T x T -> bool is an equivalence relation.
        let infinite_slope : Result<Slope, String> = Slope::new(1,0);
pub fn equivalence_class_representatives<T>(set: Vec<T>, equiv: EquivalenceRelation<T>) -> Vec<T>
        assert_eq!(infinite_slope, Ok(Slope::infinity()));
where T : Clone + PartialEq {
    }
     let mut result : Vec<T> = Vec::new();
   
     for i in 0..set.len() {
    #[test]
         let mut found_class_equivalent = false;
     fn slope_zero_slope() {
         'a: for j in 0..result.len(){ // All vecs in result are nonempty.
        let zero_slope : Result<Slope, String> = Slope::new(0,1);
            if equiv(&result[j], &set[i]){
        assert_eq!(zero_slope, Ok(Slope::zero()));
                found_class_equivalent = true;
     }
                continue 'a;
   
            }
    #[test]
         }
    fn slope_integer_slope() {
        if found_class_equivalent == false {
         let slope_four : Slope = Slope::integer(4i64);
            result.push(set[i].clone());
         assert_eq!(slope_four, Slope::raw(4,1));
         }
    }
   
    #[test]
    fn slope_reduction() {
        let slope_two : Result<Slope, String> = Slope::new(8,4);
        assert_eq!(slope_two, Ok(Slope::raw(2,1)));
        let slope_half : Result<Slope, String> = Slope::new(4,8);
         assert_eq!(slope_half, Ok(Slope::raw(1,2)));
    }
   
    #[test]
    fn slope_as_rational() {
        let good_rational : Result<r64, String> = Slope::new(3,2).expect("Bad slope").as_rational();
        assert_eq!(good_rational, Ok(r64::new(3,2)));
          
        let bad_rational : Result<r64, String> = Slope::new(3,0).expect("Bad slope").as_rational();
        assert_eq!(bad_rational,
            Err("Attempted to convert a Slope with `denom` == 0 into a Rational64".to_string()));
     }
     }
    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}
}


    // Tests for the module `Point`.
fn gcd_u32(m: u32, n: u32) -> u32 {
   
let mut x = m;
    #[test]
let mut y = n;
    fn point_create_point() {
while x != 0 && y != 0 {
        let point : Point = Point::new(r64::new(2,1),r64::new(2,1));
if x > y {
        assert_eq!(point, Point{ x: r64::new(2,1), y: r64::new(2,1)})
x = x % y;
    }
} else {
   
y = y % x;
    #[test]
}
    fn point_midpoint_vertical() {
}
        let point1 : Point = Point::new(r64::new(2,1),r64::new(2,1));
if y == 0 {x} else {y}
        let point2 : Point = Point::new(r64::new(2,1),r64::new(1,1));
}
        assert_eq!(Point::midpoint(&point1, &point2), Point{ x: r64::new(2,1), y: r64::new(3,2)})
    }
   
    #[test]
    fn point_midpoint_nonvertical() {
        let point1 : Point = Point::new(r64::new(1,1),r64::new(2,1));
        let point2 : Point = Point::new(r64::new(2,1),r64::new(1,1));
        assert_eq!(Point::midpoint(&point1, &point2), Point{ x: r64::new(3,2), y: r64::new(3,2)})
    }


    #[test]
fn lcm(m: i64, n: i64) -> i64 {
    fn point_midpoint_identical() {
m*n/gcd(m,n)
        let point1 : Point = Point::new(r64::new(1,1),r64::new(1,1));
}
        let point2 : Point = Point::new(r64::new(1,1),r64::new(1,1));
        assert_eq!(Point::midpoint(&point1, &point2), Point{ x: r64::new(1,1), y: r64::new(1,1)})
    }
   
    #[test]
    fn point_collinear() {
        let point1 : Point = Point::new(r64::new(1,1),r64::new(1,1));
        let point2 : Point = Point::new(r64::new(3,1),r64::new(4,1));
        let point3 : Point = Point::new(r64::new(1,1),r64::new(4,1));
        let point4 : Point = Point::new(r64::new(2,1),r64::new(5,2));
       
        assert!(Point::are_collinear(point1, point1, point1));
        assert!(Point::are_collinear(point1, point1, point2));
        assert!(Point::are_collinear(point1, point2, point1));
        assert!(Point::are_collinear(point2, point1, point1));
        assert!(Point::are_collinear(point1, point2, point4));
        assert!(!Point::are_collinear(point1, point2, point3));
    }
   
    #[test]
    fn line_parallel_implies_equal_slope() {
        let line1 = Line::from_points(Point::new(r64::new(0,1),r64::new(0,1)), Point::new(r64::new(1,1),r64::new(1,1))).unwrap();
        let line2 = Line::from_points(Point::new(r64::new(0,1),r64::new(1,1)), Point::new(r64::new(1,1),r64::new(2,1))).unwrap();
        assert_eq!(line1.slope(), line2.slope());
    }


    // Tests for the module `Line`.
// Is the scale word s = s(x, y, z) abstractly mv3?
    #[test]
pub fn maximum_variety(s: &str) -> usize {
    fn line_equality() {
    let mut result = 0;
        let line1 = Line::from_points(Point::new(r64::new(0,1),r64::new(0,1)), Point::new(r64::new(1,1),r64::new(1,1))).unwrap();
let floor_half: usize = s.len()/2;
        let line2 = Line::from_slope_and_point(Slope::integer(1), Point::new(r64::new(2,1),r64::new(2,1))); // equal to `line1`
for l in 1..(floor_half+1) {
        assert_eq!(line1, line2);
let mut sizes: BTreeSet<(usize, usize, usize)> = BTreeSet::new();
        let line3 = Line::from_slope_and_point(Slope::integer(1), Point::new(r64::new(1,2),r64::new(1,2)));
for b in 0..(s.len()) {
        assert_eq!(line3, line2);
let mut size: (usize, usize, usize) = (0, 0, 0);
       
let sl = &slice_cyclic_string(s, b, l);
        let line4 = Line::from_points(Point::new(r64::new(0,1),r64::new(0,1)), Point::new(r64::new(2,1),r64::new(1,1))).unwrap();
let chars_in_sl = sl.chars();
        assert_ne!(line1, line4);
for ch in chars_in_sl {
        let line5 = Line::from_points(Point::new(r64::new(1,1),r64::new(2,1)), Point::new(r64::new(0,1),r64::new(1,1))).unwrap();
match ch {
        assert_ne!(line1, line5);
'x' => {size.0 += 1;},
    }
'y' => {size.1 += 1;},
   
'z' => {size.2 += 1;},
    #[test]
_ => {panic!()},
    fn line_intersection_parallel_equal() {
}
        let line1 = Line::from_points(Point::new(r64::new(0,1),r64::new(0,1)), Point::new(r64::new(1,1),r64::new(1,1))).unwrap();
}
         let line2 = Line::from_slope_and_point(Slope::integer(1), Point::new(r64::new(2,1),r64::new(2,1))); // equal to `line1`
sizes.insert(size);
        assert_eq!(intersection(line1, line2), None);
}
    }
         result = max(result, sizes.len()); // update result
}
result
}


   
// Is the scale word s = s(x, y, z) pairwise well-formed (pwf)?
    #[test]
pub fn is_pwf(s: &str) -> bool {
    fn line_intersection_parallel_unequal() {
let s1 = s.replace("z", "y");
        let line1 = Line::from_points(Point::new(r64::new(0,1),r64::new(0,1)), Point::new(r64::new(1,1),r64::new(1,1))).unwrap();
let s2 = s.replace("z", "x");
        let line2 = Line::from_points(Point::new(r64::new(0,1),r64::new(1,2)), Point::new(r64::new(1,1),r64::new(3,2))).unwrap(); // not equal but parallel to `line1`
let s3 = s.replace("y", "x").replace("z", "y");
        assert_eq!(intersection(line1, line2), None);
return is_sv2(&s1)&&is_sv2(&s2)&&is_sv2(&s3);
}


        let line3 = Line::from_points(Point::new(r64::new(0,1),r64::new(1,1)), Point::new(r64::new(1,1),r64::new(2,1))).unwrap(); // not equal but parallel to `line1`
// Is the scale word s = s(a,b,c) pairwise well-formed (pwf)?
        assert_eq!(intersection(line1, line3), None);
pub fn is_pmos(s: &str) -> bool {
    }
let s1 = s.replace("z", "y");
   
let s2 = s.replace("z", "x");
    #[test]
let s3 = s.replace("y", "x").replace("z", "y");
    fn line_intersection_nonparallel() {
maximum_variety(&s1) == 2 && maximum_variety(&s2) == 2  && maximum_variety(&s3) == 2
        let line1 = Line::from_points(Point::new(r64::new(0,1),r64::new(0,1)), Point::new(r64::new(1,1),r64::new(1,1))).unwrap();
}
        let line2 = Line::from_points(Point::new(r64::new(0,1),r64::new(0,1)), Point::new(r64::new(1,1),r64::new(2,1))).unwrap();
        assert_eq!(intersection(line1, line2), Some(Point::new(r64::new(0,1), r64::new(0,1))));
    }
   
    #[test]
    fn point_order_coordinate_vecs() {
        let e1 = Point::new(r64::new(1,1),r64::new(0,1));
        let e2 = Point::new(r64::new(0,1),r64::new(1,1));


        let mut points = vec![e1, -e1, e2, -e2, e1+e2, e1-e2, -e1-e2, -e1+e2];
// Is the scale word s = s(x, y, z), x > y > z > 0, monotone-mos?
        points = Point::ordered_ccw(points);
pub fn is_mmos(scale: &str) -> bool {
        assert!(points[0] == -e1-e2 && points[1] == -e2 && points[2] == e1-e2 && points[3] == e1 && points[4] == e1+e2 && points[5] == e2 && points[6] == -e1+e2 && points[7] == -e1 );
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
    #[test]
maximum_variety(&s1) == 2 && maximum_variety(&s2) == 2  && maximum_variety(&s3) == 2
    fn convexpolygon_disallow_polygon_with_fewer_than_3_vertices() {
}
        let p1 = Point::new(r64::new(1,1),r64::new(0,1));
        let p2 = Point::new(r64::new(0,1),r64::new(1,1));
        let result0 = ConvexPolygon::new(vec![]);
        assert_eq!(result0, None);
        let result1 = ConvexPolygon::new(vec![p1]);
        assert_eq!(result1, None);
        let result2 = ConvexPolygon::new(vec![p1, p2]);
        assert_eq!(result2, None);
    }
   
    #[test]
    fn convexpolygon_vertices_of_triangle_are_sorted() {
        let p1 = Point::new(r64::new(1,1),r64::new(0,1));
        let p2 = Point::new(r64::new(0,1),r64::new(1,1));
        let p3 = Point::new(r64::new(0,1),r64::new(0,1));


        let triangle = ConvexPolygon::new(vec![p1, p2, p3]);
// Is the scale s = s(x, y) a single period mos (aka well-formed, wf)?
        assert!(triangle.is_some());
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;
}


         let vertices = triangle.unwrap().vertices;
// 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') }.
        assert_eq!(vertices, vec![p3, p1, p2]);
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
     }
     }
   
}
    #[test]
    fn convexpolygon_vertices_of_square_are_sorted() {
        let p1 = Point::new(r64::new(1,1),r64::new(0,1));
        let p2 = Point::new(r64::new(0,1),r64::new(1,1));
        let p3 = Point::new(r64::new(0,1),r64::new(0,1));
        let p4 = Point::new(r64::new(1,1),r64::new(1,1));
       
        let square = ConvexPolygon::new(vec![p1, p2, p3, p4]);
        assert!(square.is_some());


         let vertices = square.unwrap().vertices;
// Function for computing two properties, maximum variety and balance, while executing the subroutine to extract words only once.
        assert_eq!(vertices, vec![p3, p1, p4, p2]);
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)
     }
     }
}


    #[test]
    fn convexpolygon_vertices_of_5gon_are_sorted() {
        let p1 = Point::new(r64::new(1,1),r64::new(0,1));
        let p2 = Point::new(r64::new(0,1),r64::new(1,1));
        let p3 = Point::new(r64::new(0,1),r64::new(0,1));
        let p4 = Point::new(r64::new(1,1),r64::new(1,1));
        let p5 = Point::new(r64::new(1,2),r64::new(2,1));
       
        let pentagon  = ConvexPolygon::new(vec![p1, p2, p3, p4, p5]);
        assert!(pentagon.is_some());


         let vertices = pentagon.unwrap().vertices;
// Return a Christoffel word with `a` x's and `b` y's.
         assert_eq!(vertices, vec![p3, p1, p4, p5, p2]);
// 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>()
     }
     }
   
}
    #[test]
    fn convexpolygon_triangle_has_centroid_inside() {
        let p1 = Point::new(r64::new(1,1),r64::new(0,1));
        let p2 = Point::new(r64::new(0,1),r64::new(1,1));
        let p3 = Point::new(r64::new(0,1),r64::new(0,1));


        let triangle = ConvexPolygon::new(vec![p1, p2, p3]);
/// Billiard scales with the given signature up to rotation.
         let centroid = ConvexPolygon::centroid(&triangle.clone().unwrap());
pub fn billiard_scales(a: u32, b: u32, c: u32) -> Vec<String> {
        assert!(triangle.clone().unwrap().has_point_inside(centroid));
    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);
     #[test]
     result
    fn convexpolygon_vertex_is_not_inside() {
}
        let p1 = Point::new(r64::new(1,1),r64::new(0,1));
</syntaxhighlight>
        let p2 = Point::new(r64::new(0,1),r64::new(1,1));
=== main.rs ===
        let p3 = Point::new(r64::new(0,1),r64::new(0,1));
<syntaxhighlight lang="rust>
use scale_theory;
use std::env;
use std::fs::File;
use std::io::Write;
use std::time::Instant;


        let triangle = ConvexPolygon::new(vec![p1, p2, p3]);
/// Gives the greatest common denominator of the two inputs, unless that's 2^63.
         assert!(!triangle.unwrap().has_point_inside(p1));
/// 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;
     }
     }
   
    #[test]
    fn convexpolygon_point_on_edge_is_not_inside() {
        let p1 = Point::new(r64::new(1,1),r64::new(0,1));
        let p2 = Point::new(r64::new(0,1),r64::new(1,1));
        let p3 = Point::new(r64::new(0,1),r64::new(0,1));


        let triangle = ConvexPolygon::new(vec![p1, p2, p3]);
    // `|` is bitwise OR. `trailing_zeros` quickly counts a binary number's
        let p_edge = Point::midpoint(&p1, &p2);
    // trailing zeros, giving its prime factorization's exponent on two.
        assert!(!triangle.unwrap().has_point_inside(p_edge)); 
     let gcd_exponent_on_two = (u | v).trailing_zeros();
    }
      
    #[test]
    fn convexpolygon_point_outside_is_not_inside() {
        let p1 = Point::new(r64::new(1,1),r64::new(0,1));
        let p2 = Point::new(r64::new(0,1),r64::new(1,1));
        let p3 = Point::new(r64::new(0,1),r64::new(0,1));


        let triangle = ConvexPolygon::new(vec![p1, p2, p3]);
    // `>>=` divides the left by two to the power of the right, storing that in
        let p_outside = Point::new(r64::new(1,1),r64::new(1,1));
    // the left variable. `u` divided by its prime factorization's power of two
        assert!(!triangle.unwrap().has_point_inside(p_outside));  
    // turns it odd.
    }
    u >>= u.trailing_zeros();
    v >>= v.trailing_zeros();


     #[test]
     while u != v {
    fn line_point_line_config_test_vertical_line() {
        if u < v {
        let p_left = Point::new(r64::new(-1,1),r64::new(0,1));
            // Swap the variables' values with each other.
        let p_on_line = Point::new(r64::new(0,1),r64::new(1,1));
            core::mem::swap(&mut u, &mut v);
         let p_right = Point::new(r64::new(1,1),r64::new(1,1));
         }
         let vertical_line_up = Line::from_slope_and_point(Slope::infinity(), Point::zero());
         u -= v;
         assert_eq!(vertical_line_up.point_line_config(p_left), PointLineConfiguration::Left);
         u >>= u.trailing_zeros();
        assert_eq!(vertical_line_up.point_line_config(p_on_line), PointLineConfiguration::OnTheLine);
        assert_eq!(vertical_line_up.point_line_config(p_right), PointLineConfiguration::Right);
        let vertical_line_down = Line::from_points(Point::zero(), Point::new(r64::new(0,1), r64::new(-1,1)));
        assert_eq!(vertical_line_down.clone().unwrap().point_line_config(p_left), PointLineConfiguration::Right);
        assert_eq!(vertical_line_down.clone().unwrap().point_line_config(p_on_line), PointLineConfiguration::OnTheLine);
        assert_eq!(vertical_line_down.clone().unwrap().point_line_config(p_right), PointLineConfiguration::Left);
     }
     }


     #[test]
     // `<<` multiplies the left by two to the power of the right.
    fn line_point_line_config_test_pos_slope() {
     (u << gcd_exponent_on_two) as i64
        let p_above = Point::new(r64::new(0,1),r64::new(1,1));
}
        let p_on_line = Point::new(r64::new(1,1),r64::new(1,1));
        let p_below = Point::new(r64::new(1,1),r64::new(0,1));
        let line_up_right = Line::from_slope_and_point(Slope::integer(1), Point::zero());
        assert_eq!(line_up_right.point_line_config(p_above), PointLineConfiguration::Left);
        assert_eq!(line_up_right.point_line_config(p_on_line), PointLineConfiguration::OnTheLine);
        assert_eq!(line_up_right.point_line_config(p_below), PointLineConfiguration::Right);
        let line_down_left = Line::from_points(Point::zero(), Point::new(r64::new(-1,1), r64::new(-1,1)));
        assert_eq!(line_down_left.clone().unwrap().point_line_config(p_above), PointLineConfiguration::Right);
        assert_eq!(line_down_left.clone().unwrap().point_line_config(p_on_line), PointLineConfiguration::OnTheLine);
        assert_eq!(line_down_left.clone().unwrap().point_line_config(p_below), PointLineConfiguration::Left);
    }
   
     #[test]
    fn line_point_line_config_test_neg_slope() {
        let p_below = Point::new(r64::new(-1,1),r64::new(0,1));
        let p_on_line = Point::new(r64::new(-1,1),r64::new(1,1));
        let p_above = Point::new(r64::new(0,1),r64::new(1,1));
        let line_down_right = Line::from_slope_and_point(Slope::integer(-1), Point::zero());
        assert_eq!(line_down_right.point_line_config(p_below), PointLineConfiguration::Right);
        assert_eq!(line_down_right.point_line_config(p_on_line), PointLineConfiguration::OnTheLine);
        assert_eq!(line_down_right.point_line_config(p_above), PointLineConfiguration::Left);
        let line_up_left = Line::from_points(Point::zero(), Point::new(r64::new(-1,1), r64::new(1,1)));
        assert_eq!(line_up_left.clone().unwrap().point_line_config(p_below), PointLineConfiguration::Left);
        assert_eq!(line_up_left.clone().unwrap().point_line_config(p_on_line), PointLineConfiguration::OnTheLine);
        assert_eq!(line_up_left.clone().unwrap().point_line_config(p_above), PointLineConfiguration::Right);
    }
   
    #[test]
    fn convexpolygon_if_all_vertices_are_to_one_side() {
        let p1 = Point::new(r64::new(1,1),r64::new(0,1));
        let p2 = Point::new(r64::new(0,1),r64::new(1,1));
        let p3 = Point::new(r64::new(0,1),r64::new(0,1));
        let p4 = Point::new(r64::new(1,1),r64::new(1,1));
        let square = ConvexPolygon::raw(vec![p3, p1, p4, p2]);
       
        let line_to_the_right = Line::from_slope_and_point(Slope::infinity(), Point::new(r64::new(2,1), r64::new(0,1)));
        let line_to_the_left = Line::from_slope_and_point(Slope::infinity(), Point::new(r64::new(-1,1), r64::new(0,1)));
        let line_above = Line::from_slope_and_point(Slope::zero(), Point::new(r64::new(0,1), r64::new(2,1)));
        let line_below = Line::from_slope_and_point(Slope::zero(), Point::new(r64::new(0,1), r64::new(-1,1)));
       
        assert_eq!(square.subdivide(line_to_the_right), (Some(ConvexPolygon::raw(vec![p3, p1, p4, p2])), Option::<ConvexPolygon>::None) );
        assert_eq!(square.subdivide(line_to_the_left), (Option::<ConvexPolygon>::None, Some(ConvexPolygon::raw(vec![p3, p1, p4, p2]))) );
        assert_eq!(square.subdivide(line_above), (Option::<ConvexPolygon>::None, Some(ConvexPolygon::raw(vec![p3, p1, p4, p2])))  );
        assert_eq!(square.subdivide(line_below), (Some(ConvexPolygon::raw(vec![p3, p1, p4, p2])), Option::<ConvexPolygon>::None) );
    }
   
    #[test]
    fn convexpolygon_new_line_cuts_two_edges() {
        let p1 = Point::new(r64::new(1,1),r64::new(0,1));
        let p2 = Point::new(r64::new(0,1),r64::new(1,1));
        let p3 = Point::new(r64::new(0,1),r64::new(0,1));
        let p4 = Point::new(r64::new(1,1),r64::new(1,1));
        let square = ConvexPolygon::raw(vec![p3, p1, p4, p2]);
       
        let line = Line::from_slope_and_point(Slope::infinity(), Point::new(r64::new(1,2), r64::new(0,1)));
        let (left_polygon, right_polygon) = square.subdivide(line);
        assert_eq!(left_polygon, Some(ConvexPolygon::raw(vec![Point::new(r64::new(0,1),r64::new(0,1)),
                Point::new(r64::new(1,2),r64::new(0,1)),
                Point::new(r64::new(1,2),r64::new(1,1)),
                Point::new(r64::new(0,1),r64::new(1,1))
                ])), "`left_polygon` is not correct");
        assert_eq!(right_polygon, Some(ConvexPolygon::raw(vec![Point::new(r64::new(1,2),r64::new(0,1)),
                Point::new(r64::new(1,1),r64::new(0,1)),
                Point::new(r64::new(1,1),r64::new(1,1)),
                Point::new(r64::new(1,2),r64::new(1,1))
                ])), "`right_polygon` is not correct");
    }
   
    #[test]
    fn convexpolygon_new_line_grazes_one_vertex() {
        let p1 = Point::new(r64::new(1,1),r64::new(0,1));
        let p2 = Point::new(r64::new(0,1),r64::new(1,1));
        let p3 = Point::new(r64::new(0,1),r64::new(0,1));
        let p4 = Point::new(r64::new(1,1),r64::new(1,1));
        let square = ConvexPolygon::raw(vec![p3, p1, p4, p2]);
       
        let line = Line::from_slope_and_point(Slope::integer(1), Point::new(r64::new(0,1), r64::new(1,1)));
        let (left_polygon, right_polygon) = square.subdivide(line);
       
        assert_eq!(left_polygon, None, "`left_polygon` is not correct");
        assert_eq!(right_polygon, Some(ConvexPolygon::raw(vec![p3, p1, p4, p2])), "`right_polygon` is not correct");


        let line2 = Line::from_slope_and_point(Slope::integer(1), Point::new(r64::new(1,1), r64::new(0,1)));
fn main() -> std::io::Result<()> {
        let (left_polygon, right_polygon) = square.subdivide(line2);
    env::set_var("RUST_BACKTRACE", "1");
       
    let mut f = File::create("billiard_scales_mediawiki")?;
        assert_eq!(left_polygon, Some(ConvexPolygon::raw(vec![p3, p1, p4, p2])), "`left_polygon` is not correct");
    let start_time = Instant::now();
        assert_eq!(right_polygon, None, "`right_polygon` is not correct");
    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();
    #[test]
    println!("Done in {:?}", elapsed_time);
    fn convexpolygon_new_line_coincides_with_one_edge() {
        let p1 = Point::new(r64::new(1,1),r64::new(0,1));
        let p2 = Point::new(r64::new(0,1),r64::new(1,1));
        let p3 = Point::new(r64::new(0,1),r64::new(0,1));
        let p4 = Point::new(r64::new(1,1),r64::new(1,1));
        let square = ConvexPolygon::raw(vec![p3, p1, p4, p2]);
       
        let line = Line::from_slope_and_point(Slope::infinity(), Point::new(r64::new(0,1), r64::new(0,1)));
        let (left_polygon, right_polygon) = square.subdivide(line);
       
        assert_eq!(left_polygon, None, "`left_polygon` is not correct");
        assert_eq!(right_polygon, Some(ConvexPolygon::raw(vec![p3, p1, p4, p2])), "`right_polygon` is not correct");


        let line = Line::from_slope_and_point(Slope::integer(0), Point::new(r64::new(1,1), r64::new(0,1)));
     Ok(())
        let (left_polygon, right_polygon) = square.subdivide(line);
       
        assert_eq!(left_polygon, Some(ConvexPolygon::raw(vec![p3, p1, p4, p2])), "`left_polygon` is not correct");
        assert_eq!(right_polygon, None, "`right_polygon` is not correct");
    }
   
    #[test]
    fn convexpolygon_new_line_cuts_one_edge() {
        let p1 = Point::new(r64::new(1,1),r64::new(0,1));
        let p2 = Point::new(r64::new(0,1),r64::new(1,1));
        let p3 = Point::new(r64::new(0,1),r64::new(0,1));
        let p4 = Point::new(r64::new(1,1),r64::new(1,1));
        let square = ConvexPolygon::raw(vec![p3, p1, p4, p2]);
       
        let line = Line::from_slope_and_point(Slope::integer(2), Point::new(r64::new(0,1), r64::new(0,1)));
        let (left_polygon, right_polygon) = square.subdivide(line);
       
        assert_eq!(left_polygon, Some(ConvexPolygon::raw(vec![Point::new(r64::new(0,1),r64::new(0,1)),
                Point::new(r64::new(1,2),r64::new(1,1)),
                Point::new(r64::new(0,1),r64::new(1,1)),
                ])), "`left_polygon is not correct");
        assert_eq!(right_polygon, Some(ConvexPolygon::raw(vec![Point::new(r64::new(0,1),r64::new(0,1)),
                Point::new(r64::new(1,1),r64::new(0,1)),
                Point::new(r64::new(1,1),r64::new(1,1)),
                Point::new(r64::new(1,2),r64::new(1,1))
                ])), "`right_polygon` is not correct");
    }
   
    #[test]
    fn convexpolygon_new_line_cuts_no_edges() {
        let p1 = Point::new(r64::new(1,1),r64::new(0,1));
        let p2 = Point::new(r64::new(0,1),r64::new(1,1));
        let p3 = Point::new(r64::new(0,1),r64::new(0,1));
        let p4 = Point::new(r64::new(1,1),r64::new(1,1));
        let square = ConvexPolygon::raw(vec![p3, p1, p4, p2]);
       
        let line = Line::from_slope_and_point(Slope::integer(1), Point::new(r64::new(0,1), r64::new(0,1)));
        let (left_polygon, right_polygon) = square.subdivide(line);
       
        assert_eq!(left_polygon, Some(ConvexPolygon::raw(vec![Point::new(r64::new(0,1),r64::new(0,1)),
                Point::new(r64::new(1,1),r64::new(1,1)),
                Point::new(r64::new(0,1),r64::new(1,1))
                ])), "`left_polygon` is not correct");
        assert_eq!(right_polygon, Some(ConvexPolygon::raw(vec![Point::new(r64::new(0,1),r64::new(0,1)),
                Point::new(r64::new(1,1),r64::new(0,1)),
                Point::new(r64::new(1,1),r64::new(1,1))
                ])), "`right_polygon` is not correct");
    }
   
    #[test]
     fn convexpolygon_subdivide_test() {
        let hexagon = projected_cube(5, 2, 3);
        assert_eq!(hexagon, ConvexPolygon::new(vec![Point::new(r64::new(1,1),r64::new(0,1)), Point::new(r64::new(0,1),r64::new(1,1)), Point::new(r64::new(1,1),r64::new(1,1)),
                Point::new(r64::new(-5,3),r64::new(-2,3)),  Point::new(r64::new(-2,3),r64::new(-2,3)), Point::new(r64::new(-5,3),r64::new(1,3))]) );
        let line = Line::from_slope_and_point(Slope::integer(0), Point::new(r64::new(0,1), r64::new(-1,3)));
        let divided = hexagon.unwrap().subdivide(line);
        assert_eq!(divided.0, ConvexPolygon::new(vec![Point::new(r64::new(1,1),r64::new(0,1)), Point::new(r64::new(0,1),r64::new(1,1)), Point::new(r64::new(1,1),r64::new(1,1)),
                Point::new(r64::new(-5,3),r64::new(-1,3)),  Point::new(r64::new(1,6),r64::new(-1,3)), Point::new(r64::new(-5,3),r64::new(1,3))]) );
        assert_eq!(divided.1, ConvexPolygon::new(vec![ Point::new(r64::new(-5,3),r64::new(-1,3)),  Point::new(r64::new(1,6),r64::new(-1,3)),
                Point::new(r64::new(-5,3),r64::new(-2,3)), Point::new(r64::new(-2,3),r64::new(-2,3)) ]) );
        let line2 = Line::from_slope_and_point(Slope::integer(0), Point::new(r64::new(0,1), r64::new(0,1)));
        let divided_again = (divided.0).unwrap().subdivide(line2);
       
        assert_eq!(divided_again.0, ConvexPolygon::new(vec![Point::new(r64::new(1,1),r64::new(0,1)), Point::new(r64::new(0,1),r64::new(1,1)), Point::new(r64::new(1,1),r64::new(1,1)),
                Point::new(r64::new(-5,3),r64::new(0,1)), Point::new(r64::new(-5,3),r64::new(1,3))]) );
        assert_eq!(divided_again.1, ConvexPolygon::new(vec![Point::new(r64::new(1,1),r64::new(0,1)), Point::new(r64::new(1,6),r64::new(-1,3)),
                Point::new(r64::new(-5,3),r64::new(0,1)), Point::new(r64::new(-5,3),r64::new(-1,3))]) );
    }
}
}
</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>
</syntaxhighlight>