From b5d298901373da84b4fda21a6629b90997756b14 Mon Sep 17 00:00:00 2001
From: Developing Human <developinghuman@protonmail.com>
Date: Fri, 19 Dec 2025 22:13:13 -0500
Subject: [PATCH 1/3] FIX: compare to epsilon during rref

---
 src/structure/matrix.rs |  5 ++--
 tests/matrix.rs         | 51 +++++++++++++++++++++++++++++++++++++++++
 2 files changed, 54 insertions(+), 2 deletions(-)

diff --git a/src/structure/matrix.rs b/src/structure/matrix.rs
index 72220e29..be937174 100644
--- a/src/structure/matrix.rs
+++ b/src/structure/matrix.rs
@@ -3293,6 +3293,7 @@ impl LinearAlgebra<Matrix> for Matrix {
     ///
     /// Implementation of [RosettaCode](https://rosettacode.org/wiki/Reduced_row_echelon_form)
     fn rref(&self) -> Matrix {
+        let epsilon = 1e-10;
         let mut lead = 0usize;
         let mut result = self.clone();
         'outer: for r in 0..self.row {
@@ -3300,7 +3301,7 @@ impl LinearAlgebra<Matrix> for Matrix {
                 break;
             }
             let mut i = r;
-            while result[(i, lead)] == 0f64 {
+            while result[(i, lead)].abs() < epsilon {
                 i += 1;
                 if self.row == i {
                     i = r;
@@ -3314,7 +3315,7 @@ impl LinearAlgebra<Matrix> for Matrix {
                 result.swap(i, r, Row);
             }
             let tmp = result[(r, lead)];
-            if tmp != 0f64 {
+            if tmp.abs() > epsilon {
                 unsafe {
                     result.row_mut(r).iter_mut().for_each(|t| *(*t) /= tmp);
                 }
diff --git a/tests/matrix.rs b/tests/matrix.rs
index defb9149..d3a21107 100644
--- a/tests/matrix.rs
+++ b/tests/matrix.rs
@@ -111,3 +111,54 @@ fn test_kronecker() {
     let c1 = a1.kronecker(&b1);
     assert_eq!(c1, ml_matrix("0 5 0 10;6 7 12 14;0 15 0 20;18 21 24 28"));
 }
+
+#[test]
+fn test_rref() {
+    let a = ml_matrix(
+        r#"
+        -3 2 -1 -1;
+        6 -6 7 -7;
+        3 -4 4 -6"#,
+    );
+    let b = a.rref();
+
+    assert_eq!(
+        b,
+        ml_matrix(
+            r#"
+    1 0 0 2;
+    0 1 0 2;
+    0 0 1 -1"#
+        )
+    );
+}
+
+#[test]
+fn test_rref_unstable() {
+    let epsilon = 1e-10;
+
+    // this matrix has a tendency to become unstable during rref,
+    let a = ml_matrix(
+        r#"
+    1 1 0 0 0 1 0 1 31;
+    1 1 1 1 0 0 1 1 185;
+    0 0 1 0 0 1 1 1 165;
+    1 0 1 0 1 1 0 1 32;
+    1 0 1 0 0 0 1 1 174;
+    0 0 1 0 1 1 1 1 171;
+    0 1 1 0 1 1 0 1 27;
+    1 0 0 1 0 1 0 0 20;
+    1 0 1 1 0 1 0 0 23"#,
+    );
+
+    let b = a.rref();
+
+    // creating a row like "0 0 0 0 0 0 0 0 1" will "prove" 0 == 1
+    // which is a tell of numeric instability
+    for row in 0..b.row {
+        let ends_in_1 = (b[(row, b.col - 1)] - 1.0).abs() < epsilon;
+        let rest_zeroes = (0..b.col - 1).all(|col| b[(row, col)].abs() < epsilon);
+
+        assert!(!(ends_in_1 && rest_zeroes));
+    }
+}

From e4518f4e0332b746a27b0bb1e63875526fb00b0f Mon Sep 17 00:00:00 2001
From: Developing Human <developinghuman@protonmail.com>
Date: Fri, 19 Dec 2025 22:26:11 -0500
Subject: [PATCH 2/3] cleanup comment

---
 tests/matrix.rs | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/tests/matrix.rs b/tests/matrix.rs
index d3a21107..3ce5d3c0 100644
--- a/tests/matrix.rs
+++ b/tests/matrix.rs
@@ -137,7 +137,7 @@ fn test_rref() {
 fn test_rref_unstable() {
     let epsilon = 1e-10;
 
-    // this matrix has a tendency to become unstable during rref,
+    // this matrix used to become unstable during rref
     let a = ml_matrix(
         r#"
     1 1 0 0 0 1 0 1 31;

From f987e345a1808ecebd5514addb15f475848affbc Mon Sep 17 00:00:00 2001
From: Developing Human <developinghuman@protonmail.com>
Date: Sat, 20 Dec 2025 12:43:04 -0500
Subject: [PATCH 3/3] compute epsilon based on the scale of matrix data

---
 src/structure/matrix.rs | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/src/structure/matrix.rs b/src/structure/matrix.rs
index be937174..7e0241db 100644
--- a/src/structure/matrix.rs
+++ b/src/structure/matrix.rs
@@ -3293,7 +3293,9 @@ impl LinearAlgebra<Matrix> for Matrix {
     ///
     /// Implementation of [RosettaCode](https://rosettacode.org/wiki/Reduced_row_echelon_form)
     fn rref(&self) -> Matrix {
-        let epsilon = 1e-10;
+        let max_abs = self.data.iter().fold(0f64, |acc, &x| acc.max(x.abs()));
+        let epsilon = (max_abs * 1e-12).max(1e-15);
+
         let mut lead = 0usize;
         let mut result = self.clone();
         'outer: for r in 0..self.row {