@@ Line 619: / Line 619: @@
 /// 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".
-/// Assumes gcd(a, b, c) == 1.
 pub fn advance(a: u32, b: u32, c: u32, pt: Point) -> Result<(Point, String), String> {
      let projected_cube = projected_cube(a, b, c).unwrap();
@@ Line 643: / Line 642: @@
          }
      }
+}
+</syntaxhighlight>
+=== lib.rs ===
+<syntaxhighlight lang="rs">use std::collections::BTreeSet;
+use num_rational::Rational64 as r64;
+mod plane_geometry;
+use plane_geometry::*;
+use std::cmp::{min, max};
+type EquivalenceRelation<T> = fn(&T, &T) -> bool;
+/// Treating `scale` as a cyclic string (that is, "scale[i] == scale[i % scale.len()]"),
+/// take a slice of length `slice_length` from `basepoint`; assumes `slice_length` <= `scale`.len()
+fn slice_cyclic_string(scale: &str, basepoint: usize, slice_length: usize) -> String {
+	let arr = scale;
+	if basepoint + slice_length < scale.len() {
+		return String::from(&arr[basepoint..(basepoint+slice_length)]);
+	} else {
+		let mut ret = String::new();
+		ret.push_str(&arr[basepoint..]);
+		ret.push_str(&arr[..(basepoint + slice_length - scale.len())]);
+		return ret;
+	}
+}
+pub fn are_rotationally_equivalent(scale1: &String, scale2: &String)-> bool{
+    if scale1.len() != scale2.len() {
+        return false;
+    } else {
+        let length = scale1.len();
+        for i in 0..length {
+            if slice_cyclic_string(scale1, i, length) == *scale2 {
+                return true;
+            }
+        }
+        false
+    }
+}
+/// Assuming that "`equiv` : T x T -> bool" is an equivalence relation,
+/// return the equivalence class representatives `set` under `equiv`.
+pub fn equivalence_class_representatives<T>(set: Vec<T>, equiv: EquivalenceRelation<T>) -> Vec<T>
+where T : Clone + PartialEq {
+    let mut result : Vec<T> = Vec::new();
+    for i in 0..set.len() {
+        let mut found_class_equivalent = false;
+        'a: for j in 0..result.len(){ // All vecs in result are nonempty.
+            if equiv(&result[j], &set[i]){
+                found_class_equivalent = true;
+                continue 'a;
+            }
+        }
+        if found_class_equivalent == false {
+            result.push(set[i].clone());
+        }
+    }
+    result
+}
+/// Gives the greatest common denominator of the two inputs, unless that's 2^63.
+/// 2^63 doesn't fit in an `i64`, so it returns -2^63, which does.
+pub fn gcd(u: i64, v: i64) -> i64 {
+    // `wrapping_abs` gives a number's absolute value, unless that's 2^63. 2^63
+    // won't fit in `i64`, so it gives -2^63 instead.
+    let mut v = v.wrapping_abs() as u64;
+    if u == 0 {
+        return v as i64;
+    }
+    let mut u = u.wrapping_abs() as u64;
+    if v == 0 {
+        return u as i64;
+    }
+    // `|` is bitwise OR. `trailing_zeros` quickly counts a binary number's
+    // trailing zeros, giving its prime factorization's exponent on two.
+    let gcd_exponent_on_two = (u | v).trailing_zeros();
+    // `>>=` divides the left by two to the power of the right, storing that in
+    // the left variable. `u` divided by its prime factorization's power of two
+    // turns it odd.
+    u >>= u.trailing_zeros();
+    v >>= v.trailing_zeros();
+    while u != v {
+        if u < v {
+            // Swap the variables' values with each other.
+            core::mem::swap(&mut u, &mut v);
+        }
+        u -= v;
+        u >>= u.trailing_zeros();
+    }
+    // `<<` multiplies the left by two to the power of the right.
+    (u << gcd_exponent_on_two) as i64
+}
+fn lcm(m: i64, n: i64) -> i64 {
+	m*n/gcd(m,n)
+}
+/// Computes the maximum variety of a ternary scale word.
+pub fn maximum_variety(scale: &str) -> usize {
+    let mut result = 0;
+	let floor_half: usize = scale.len()/2;
+	for slice_length in 1..(floor_half+1) {
+		let mut sizes: BTreeSet<(usize, usize, usize)> = BTreeSet::new();
+		for basepoint in 0..(scale.len()) {
+			let mut size: (usize, usize, usize) = (0, 0, 0);
+			let sl = &slice_cyclic_string(scale, basepoint, slice_length);
+			let chars_in_sl = sl.chars();
+			for ch in chars_in_sl {
+				match ch {
+					'x' => {size.0 += 1;},
+					'y' => {size.1 += 1;},
+					'z' => {size.2 += 1;},
+					_ => {panic!()},
+				}
+			}
+			sizes.insert(size);
+		}
+        result = max(result, sizes.len()); // update result
+	}
+	result
+}
+/// Return the balance of `s` = `s`(x, y, z), defined as max { | |w|_{x_i} - |w'|_{x_i} | : x_i is a letter of `s` and k = len(w) = len(w') }.
+pub fn balance(scale: &str) -> usize {
+    if scale.len() <= 1 {
+        scale.len()
+    } else {
+        let mut result = 0;
+	    let floor_half: usize = scale.len()/2;
+	    for slice_length in 1..(floor_half+1) {
+            let (mut min_x, mut max_x, mut min_y, mut max_y, mut min_z, mut max_z) = (usize::MAX, 0, usize::MAX, 0, usize::MAX, 0); // min and max for each length
+		    let mut sizes: BTreeSet<(usize, usize, usize)> = BTreeSet::new();
+		    for basepoint in 0..(scale.len()) {
+			    let mut size: (usize, usize, usize) = (0, 0, 0);
+			    let sl = &slice_cyclic_string(scale, basepoint, slice_length);
+			    let chars_in_sl = sl.chars();
+			    for ch in chars_in_sl {
+				    match ch {
+					    'x' => {size.0 += 1;},
+					    'y' => {size.1 += 1;},
+					    'z' => {size.2 += 1;},
+					    _ => {panic!()},
+				    }
+			    }
+			    sizes.insert(size);
+		    }
+            for vector in &sizes {
+                min_x = min(min_x, vector.0);
+                max_x = max(max_x, vector.0);
+                min_y = min(min_y, vector.1);
+                max_y = max(max_y, vector.1);
+                min_z = min(min_z, vector.2);
+                max_z = max(max_z, vector.2);
+            }
+            result = max(result, max(max(max_x.abs_diff(min_x), max_y.abs_diff(min_y)), max_z.abs_diff(min_z)));
+	    }
+        result
+    }
+}
+/// Function for computing two properties, maximum variety and balance, while executing the subroutine to extract words only once at a time for both properties.
+pub fn max_variety_and_balance(scale: &str) -> (usize, usize) {
+    if scale.len() <= 1 {
+        (scale.len(), 0)
+    } else {
+        let mut mv_result = 0;
+        let mut balance_result = 0;
+        let (mut min_x, mut max_x, mut min_y, mut max_y, mut min_z, mut max_z) = (usize::MAX, 0, usize::MAX, 0, usize::MAX, 0);
+	    let floor_half: usize = scale.len()/2;
+	    for slice_length in 1..(floor_half+1) {
+		    let mut sizes: BTreeSet<(usize, usize, usize)> = BTreeSet::new();
+		    for basepoint in 0..scale.len() {
+			    let mut size: (usize, usize, usize) = (0, 0, 0);
+			    let sl = &slice_cyclic_string(scale, basepoint, slice_length);
+			    let chars_in_sl = sl.chars();
+			    for ch in chars_in_sl {
+				    match ch {
+					    'x' => {size.0 += 1;},
+					    'y' => {size.1 += 1;},
+					    'z' => {size.2 += 1;},
+					    _ => {panic!()},
+				    }
+			    }
+			    sizes.insert(size);
+		    }
+            for vector in &sizes {
+                min_x = min(min_x, vector.0);
+                max_x = max(max_x, vector.0);
+                min_y = min(min_y, vector.1);
+                max_y = max(max_y, vector.1);
+                min_z = min(min_z, vector.2);
+                max_z = max(max_z, vector.2);
+            }
+            mv_result = max(mv_result, sizes.len()); // update result
+            balance_result = max(balance_result, max(max(max_x.abs_diff(min_x), max_y.abs_diff(min_y)), max_z.abs_diff(min_z)));
+	    }
+	    (mv_result, balance_result)
+    }
+}
+/// 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(a as i64, b as i64) as u32;
+    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 d = gcd(gcd(a as i64, b as i64), c as i64) as u32; // take gcd(a, b, c)
+    let mut result = Vec::<String>::new();
+    let regions: Vec<ConvexPolygon> = plane_geometry::project_and_partition(a/d, b/d, c/d);
+    let n = a/d + b/d + c/d;
+    for region in regions {
+        let mut word = String::new();
+        let mut next_point = region.centroid();
+        let first_point = region.centroid().clone();
+        for i in 0..n {
+            let res = plane_geometry::advance(a, b, c, next_point).unwrap();
+            next_point = res.0;
+            debug_assert!(i == n-1 || first_point != next_point); // There should not be a subperiod in the orbit of the "billiard ball".
+            word.push_str(&(res.1));
+        }
+        result.push(std::iter::repeat(word).take(d.try_into().unwrap()).collect::<String>()); // Each word should be repeated d times.
+        debug_assert_eq!(first_point, next_point); // The last point should come back to the first point. When i = n - 1, next_point gets set to point number n.
+    }
+    result = equivalence_class_representatives::<String>(result, are_rotationally_equivalent);
+    result
 }
 </syntaxhighlight>
 === main.rs ===
-<syntaxhighlight lang="rust>
+<syntaxhighlight lang="rust">
 use billiard_scales;
 use std::env;

User:Inthar/Code: Difference between revisions