diff --git a/Cargo.lock b/Cargo.lock
index 32bcd38..5182d99 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -29,6 +29,12 @@ version = "1.0.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
 checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd"
 
+[[package]]
+name = "either"
+version = "1.15.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "48c757948c5ede0e46177b7add2e67155f70e33c07fea8284df6576da70b3719"
+
 [[package]]
 name = "equivalent"
 version = "1.0.2"
@@ -111,11 +117,21 @@ dependencies = [
  "hashbrown",
 ]
 
+[[package]]
+name = "itertools"
+version = "0.14.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "2b192c782037fadd9cfa75548310488aabdbf3d2da73885b31bd0abd03351285"
+dependencies = [
+ "either",
+]
+
 [[package]]
 name = "kasuari"
 version = "0.4.7"
 dependencies = [
  "hashbrown",
+ "itertools",
  "rstest",
  "thiserror",
 ]
diff --git a/Cargo.toml b/Cargo.toml
index 2f40cd2..504c958 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -34,3 +34,4 @@ thiserror = { version = "2.0", default-features = false }
 
 [dev-dependencies]
 rstest = "0.25"
+itertools = "0.14.0"
diff --git a/src/row.rs b/src/row.rs
index fb9bf6c..5154102 100644
--- a/src/row.rs
+++ b/src/row.rs
@@ -30,6 +30,10 @@ impl Symbol {
     pub fn kind(&self) -> SymbolKind {
         self.1
     }
+
+    pub fn id(&self) -> usize {
+        self.0
+    }
 }
 
 pub fn near_zero(value: f64) -> bool {
diff --git a/src/solver.rs b/src/solver.rs
index 8528e80..d1b52c3 100644
--- a/src/solver.rs
+++ b/src/solver.rs
@@ -647,18 +647,24 @@ impl Solver {
 
     /// Compute the entering variable for a pivot operation.
     ///
-    /// This method will return first symbol in the objective function which
+    /// This method will return a symbol in the objective function which
     /// is non-dummy and has a coefficient less than zero. If no symbol meets
     /// the criteria, it means the objective function is at a minimum, and an
     /// invalid symbol is returned.
     /// Could return an External symbol
     fn get_entering_symbol(objective: &Row) -> Symbol {
+        let mut entering = Symbol::invalid();
+        let mut min_id = usize::MAX;
+
         for (symbol, value) in &objective.cells {
             if symbol.kind() != SymbolKind::Dummy && *value < 0.0 {
-                return *symbol;
+                if symbol.id() < min_id {
+                    min_id = symbol.id();
+                    entering = *symbol;
+                }
             }
         }
-        Symbol::invalid()
+        entering
     }
 
     /// Compute the entering symbol for the dual optimize operation.
@@ -709,13 +715,17 @@ impl Solver {
     fn get_leaving_row(&mut self, entering: Symbol) -> Option<(Symbol, Box<Row>)> {
         let mut ratio = f64::INFINITY;
         let mut found = None;
+        let mut min_id = usize::MAX;
+
         for (symbol, row) in &self.rows {
             if symbol.kind() != SymbolKind::External {
                 let temp = row.coefficient_for(entering);
                 if temp < 0.0 {
                     let temp_ratio = -row.constant / temp;
-                    if temp_ratio < ratio {
+                    // in case of duplicates, choose symbol with the smallest id to prevent cycling
+                    if temp_ratio < ratio || (temp_ratio == ratio && symbol.id() < min_id) {
                         ratio = temp_ratio;
+                        min_id = symbol.id();
                         found = Some(*symbol);
                     }
                 }
@@ -816,3 +826,38 @@ impl Solver {
             .unwrap_or(0.0)
     }
 }
+
+#[cfg(test)]
+mod tests {
+    use itertools::Itertools;
+
+    use super::*;
+    use crate::WeightedRelation::EQ;
+
+    #[test]
+    fn regression_18() {
+        const N: usize = 90;
+
+        let mut solver = Solver::new();
+
+        let variable_count = N * 2 + 2;
+        let variables = core::iter::repeat_with(Variable::new)
+            .take(variable_count)
+            .collect_vec();
+        let segments: Vec<(&Variable, &Variable)> = variables.iter().skip(1).tuples().collect_vec();
+
+        let constraints = alloc::vec![1.0; N];
+
+        for ((&left_constraint, &left_segment), (&right_constraint, &right_segment)) in
+            constraints.iter().zip(segments.iter()).tuple_combinations()
+        {
+            solver
+                .add_constraint(
+                    (right_constraint * (*left_segment.1 - *left_segment.0))
+                        | EQ(Strength::MEDIUM.div_f64(10.0))
+                        | (left_constraint * (*right_segment.1 - *right_segment.0)),
+                )
+                .unwrap();
+        }
+    }
+}