Add partial cycle-finding functions that skip minimal μ computation

src/directed/cycle_detection.rs

@@ -1,57 +1,140 @@
 //! Identify a cycle in an infinite sequence.
 
-/// Identify a cycle in an infinite sequence using Floyd's algorithm.
-/// Return the cycle size, the first element, and the index of first element.
+/// Identify a cycle in an infinite sequence using Floyd's algorithm (partial version).
+/// Return the cycle size, an element in the cycle, and an upper bound on the index of
+/// the first element.
+///
+/// This function computes the cycle length λ and returns an element within the cycle,
+/// along with an upper bound μ̃ on the index of the first cycle element. The upper bound
+/// μ̃ satisfies μ ≤ μ̃ < μ + λ, where μ is the minimal index.
+///
+/// This is faster than [`floyd`] as it skips the computation of the minimal μ.
+/// The upper bound μ̃ is sufficient for many applications, such as computing f^n(x) for
+/// large n, where knowing the exact starting point of the cycle is not necessary.
+///
+/// # Example
+///
+/// ```
+/// use pathfinding::prelude::floyd_partial;
+///
+/// let (lam, _elem, mu_tilde) = floyd_partial(1, |x| (x * 2) % 7);
+/// assert_eq!(lam, 3); // Cycle length
+/// assert!(mu_tilde == 0); // Upper bound on mu (start is in cycle, so mu = 0)
+/// ```
 ///
 /// # Warning
 ///
 /// If no cycle exist, this function loops forever.
-pub fn floyd<T, FS>(start: T, successor: FS) -> (usize, T, usize)
+#[expect(clippy::needless_pass_by_value)]
+pub fn floyd_partial<T, FS>(start: T, successor: FS) -> (usize, T, usize)
 where
     T: Clone + PartialEq,
     FS: Fn(T) -> T,
 {
     let mut tortoise = successor(start.clone());
     let mut hare = successor(successor(start.clone()));
+    let mut tortoise_steps = 1;
     while tortoise != hare {
         (tortoise, hare) = (successor(tortoise), successor(successor(hare)));
+        tortoise_steps += 1;
     }
-    let mut mu = 0;
-    tortoise = start;
-    while tortoise != hare {
-        (tortoise, hare, mu) = (successor(tortoise), successor(hare), mu + 1);
-    }
+    // tortoise and hare met at position tortoise_steps
     let mut lam = 1;
     hare = successor(tortoise.clone());
     while tortoise != hare {
         (hare, lam) = (successor(hare), lam + 1);
     }
-    (lam, tortoise, mu)
+    // Handle edge case where they meet at the start position (pure cycle, mu = 0)
+    // In this case, tortoise_steps equals lam, and to satisfy mu_tilde < mu + lam,
+    // we must return 0.
+    let mu_tilde = if tortoise == start { 0 } else { tortoise_steps };
+    (lam, tortoise, mu_tilde)
 }
 
-/// Identify a cycle in an infinite sequence using Brent's algorithm.
+/// Identify a cycle in an infinite sequence using Floyd's algorithm.
 /// Return the cycle size, the first element, and the index of first element.
 ///
 /// # Warning
 ///
 /// If no cycle exist, this function loops forever.
-pub fn brent<T, FS>(start: T, successor: FS) -> (usize, T, usize)
+pub fn floyd<T, FS>(start: T, successor: FS) -> (usize, T, usize)
+where
+    T: Clone + PartialEq,
+    FS: Fn(T) -> T,
+{
+    let (lam, mut hare, _) = floyd_partial(start.clone(), &successor);
+    // Find the exact mu
+    let (mut mu, mut tortoise) = (0, start);
+    while tortoise != hare {
+        (tortoise, hare, mu) = (successor(tortoise), successor(hare), mu + 1);
+    }
+    (lam, tortoise, mu)
+}
+
+/// Identify a cycle in an infinite sequence using Brent's algorithm (partial version).
+/// Return the cycle size, an element in the cycle, and an upper bound on the index of
+/// the first element.
+///
+/// This function computes the cycle length λ and returns an element within the cycle,
+/// along with an upper bound μ̃ on the index of the first cycle element. The upper bound
+/// satisfies μ ≤ μ̃. Due to the nature of Brent's algorithm with its power-of-2 stepping,
+/// the bound may be looser than `μ + λ` in some cases, but is still reasonable for
+/// practical applications.
+///
+/// This is faster than [`brent`] as it skips the computation of the minimal μ.
+/// The upper bound μ̃ is sufficient for many applications, such as computing f^n(x) for
+/// large n, where knowing the exact starting point of the cycle is not necessary.
+///
+/// # Example
+///
+/// ```
+/// use pathfinding::prelude::brent_partial;
+///
+/// let (lam, _elem, mu_tilde) = brent_partial(1, |x| (x * 2) % 7);
+/// assert_eq!(lam, 3); // Cycle length
+/// assert!(mu_tilde >= 1); // Upper bound on mu
+/// ```
+///
+/// # Warning
+///
+/// If no cycle exist, this function loops forever.
+pub fn brent_partial<T, FS>(start: T, successor: FS) -> (usize, T, usize)
 where
     T: Clone + PartialEq,
     FS: Fn(T) -> T,
 {
     let mut power = 1;
     let mut lam = 1;
     let mut tortoise = start.clone();
-    let mut hare = successor(start.clone());
+    let mut hare = successor(start);
+    let mut hare_steps = 1;
     while tortoise != hare {
         if power == lam {
             (tortoise, power, lam) = (hare.clone(), power * 2, 0);
         }
         (hare, lam) = (successor(hare), lam + 1);
+        hare_steps += 1;
     }
+    // Use hare_steps as the upper bound, as it represents where we detected the cycle.
+    (lam, hare, hare_steps)
+}
+
+/// Identify a cycle in an infinite sequence using Brent's algorithm.
+/// Return the cycle size, the first element, and the index of first element.
+///
+/// # Warning
+///
+/// If no cycle exist, this function loops forever.
+pub fn brent<T, FS>(start: T, successor: FS) -> (usize, T, usize)
+where
+    T: Clone + PartialEq,
+    FS: Fn(T) -> T,
+{
+    let (lam, _hare, _hare_steps) = brent_partial(start.clone(), &successor);
+    // Find the exact mu
     let mut mu = 0;
-    (tortoise, hare) = (start.clone(), (0..lam).fold(start, |x, _| successor(x)));
+    let mut tortoise = start.clone();
+    let mut hare = (0..lam).fold(start, |x, _| successor(x));
     while tortoise != hare {
         (tortoise, hare, mu) = (successor(tortoise), successor(hare), mu + 1);
     }

tests/cycle_detection.rs

@@ -9,3 +9,128 @@ fn floyd_works() {
 fn brent_works() {
     assert_eq!(brent(-10, |x| (x + 5) % 6 + 3), (3, 6, 2));
 }
+
+#[test]
+fn floyd_partial_works() {
+    let (lam, elem, mu_tilde) = floyd_partial(-10, |x| (x + 5) % 6 + 3);
+    // Check that we get the correct cycle length
+    assert_eq!(lam, 3);
+    // Check that elem is in the cycle (cycle is 6, 8, 4)
+    assert!([4, 6, 8].contains(&elem));
+    // Check that mu_tilde is an upper bound: mu <= mu_tilde < mu + lam
+    // We know mu = 2 from the full algorithm
+    assert!(2 <= mu_tilde);
+    assert!(mu_tilde < 2 + 3);
+}
+
+#[test]
+fn brent_partial_works() {
+    let (lam, elem, mu_tilde) = brent_partial(-10, |x| (x + 5) % 6 + 3);
+    // Check that we get the correct cycle length
+    assert_eq!(lam, 3);
+    // Check that elem is in the cycle (cycle is 6, 8, 4)
+    assert!([4, 6, 8].contains(&elem));
+    // Check that mu_tilde is an upper bound: mu <= mu_tilde
+    // We know mu = 2 from the full algorithm
+    assert!(2 <= mu_tilde);
+}
+
+#[test]
+fn partial_functions_match_full_functions() {
+    // Test that partial functions return the same lambda as full functions
+    let f = |x| (x + 5) % 6 + 3;
+
+    let (lam_floyd, elem_floyd, mu_floyd) = floyd(-10, f);
+    let (lam_floyd_partial, elem_floyd_partial, mu_tilde_floyd) = floyd_partial(-10, f);
+
+    let (lam_brent, elem_brent, mu_brent) = brent(-10, f);
+    let (lam_brent_partial, elem_brent_partial, mu_tilde_brent) = brent_partial(-10, f);
+
+    // Lambda should be the same
+    assert_eq!(lam_floyd, lam_floyd_partial);
+    assert_eq!(lam_brent, lam_brent_partial);
+    assert_eq!(lam_floyd, lam_brent);
+
+    // Elements should be in the cycle
+    assert!([4, 6, 8].contains(&elem_floyd));
+    assert!([4, 6, 8].contains(&elem_floyd_partial));
+    assert!([4, 6, 8].contains(&elem_brent));
+    assert!([4, 6, 8].contains(&elem_brent_partial));
+
+    // Mu values from full algorithms should be the same
+    assert_eq!(mu_floyd, mu_brent);
+
+    // Mu_tilde should be valid upper bounds
+    assert!(mu_floyd <= mu_tilde_floyd);
+    assert!(mu_tilde_floyd < mu_floyd + lam_floyd);
+    assert!(mu_brent <= mu_tilde_brent);
+}
+
+#[test]
+fn test_longer_cycle() {
+    // Test with a longer cycle: sequence from 0 to 99, then cycles
+    let f = |x: i32| (x + 1) % 100;
+
+    let (lam_floyd, elem_floyd, mu_floyd) = floyd(0, f);
+    let (lam_floyd_partial, elem_floyd_partial, mu_tilde_floyd) = floyd_partial(0, f);
+
+    let (lam_brent, elem_brent, mu_brent) = brent(0, f);
+    let (lam_brent_partial, elem_brent_partial, mu_tilde_brent) = brent_partial(0, f);
+
+    // Cycle length should be 100 (0, 1, 2, ..., 99, 0, ...)
+    assert_eq!(lam_floyd, 100);
+    assert_eq!(lam_floyd_partial, 100);
+    assert_eq!(lam_brent, 100);
+    assert_eq!(lam_brent_partial, 100);
+
+    // Mu should be 0 (cycle starts immediately)
+    assert_eq!(mu_floyd, 0);
+    assert_eq!(mu_brent, 0);
+
+    // Elements should be in the cycle
+    assert!((0..100).contains(&elem_floyd));
+    assert!((0..100).contains(&elem_floyd_partial));
+    assert!((0..100).contains(&elem_brent));
+    assert!((0..100).contains(&elem_brent_partial));
+
+    // Mu_tilde should be valid upper bounds
+    assert!(mu_floyd <= mu_tilde_floyd);
+    assert!(mu_tilde_floyd < mu_floyd + lam_floyd);
+    assert!(mu_brent <= mu_tilde_brent);
+}
+
+#[test]
+fn test_short_cycle_large_mu() {
+    // Sequence starting from -100, adds 1 each step,
+    // but when value reaches 10, it resets to 0
+    // This creates: -100, -99, ..., -1, 0, 1, ..., 9, 10, 0, 1, ..., 9, 10, 0, ...
+    // mu = 100 (steps to reach 0, which is the start of the cycle), lambda = 11 (0 to 10 inclusive)
+    let f = |x: i32| if x == 10 { 0 } else { x + 1 };
+
+    let (lam_floyd, elem_floyd, mu_floyd) = floyd(-100, f);
+    let (lam_floyd_partial, elem_floyd_partial, mu_tilde_floyd) = floyd_partial(-100, f);
+
+    let (lam_brent, elem_brent, mu_brent) = brent(-100, f);
+    let (lam_brent_partial, elem_brent_partial, mu_tilde_brent) = brent_partial(-100, f);
+
+    // Cycle length should be 11 (0, 1, 2, ..., 9, 10, 0, ...)
+    assert_eq!(lam_floyd, 11);
+    assert_eq!(lam_floyd_partial, 11);
+    assert_eq!(lam_brent, 11);
+    assert_eq!(lam_brent_partial, 11);
+
+    // Mu should be 100 (it takes 100 steps to get from -100 to 0, then cycles)
+    assert_eq!(mu_floyd, 100);
+    assert_eq!(mu_brent, 100);
+
+    // Elements should be in the cycle (0 to 10)
+    assert!((0..=10).contains(&elem_floyd));
+    assert!((0..=10).contains(&elem_floyd_partial));
+    assert!((0..=10).contains(&elem_brent));
+    assert!((0..=10).contains(&elem_brent_partial));
+
+    // Mu_tilde should be valid upper bounds
+    assert!(mu_floyd <= mu_tilde_floyd);
+    assert!(mu_tilde_floyd < mu_floyd + lam_floyd);
+    assert!(mu_brent <= mu_tilde_brent);
+}