diff --git a/src/directed/cycle_detection.rs b/src/directed/cycle_detection.rs index 0c30ee9e..3835e72c 100644 --- a/src/directed/cycle_detection.rs +++ b/src/directed/cycle_detection.rs @@ -1,41 +1,104 @@ //! 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(start: T, successor: FS) -> (usize, T, usize) +#[expect(clippy::needless_pass_by_value)] +pub fn floyd_partial(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(start: T, successor: FS) -> (usize, T, usize) +pub fn floyd(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(start: T, successor: FS) -> (usize, T, usize) where T: Clone + PartialEq, FS: Fn(T) -> T, @@ -43,15 +106,35 @@ where 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(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); } diff --git a/tests/cycle_detection.rs b/tests/cycle_detection.rs index 604b231a..15ec6d75 100644 --- a/tests/cycle_detection.rs +++ b/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); +}