FIX: compare to epsilon during rref

[developing-human]

I found a matrix that was numerically unstable while calculating reduced row echelon form. I've updated the comparisons to zero so they instead compare to epsilon, which has been more stable for my use case.

I added a test case for the unstable scenario, as well as a more basic test since one wasn't present yet for rref.

Updating numerical crates is outside my usual wheelhouse, so I'm very open to feedback here!

[djmaxus]

This particular value of epsilon feels a bit arbitrary when applied universally here.

Would you be interested in considering floating-point constants like machine epsilon or smallest positive normal number? I also assume the choice of a particular threshold constant should be justified formally if possible

[Axect]

I agree with @djmaxus, but if we use f64::EPSILON, then since it's too strict, so it may not be practical. How's your opinion? @developing-human @djmaxus.

Furthermore, I think it's better to use relative tolerance as follows:

let max_abs = self.data.iter().fold(0f64, |acc, &x| acc.max(x.abs()));                  
let epsilon = (max_abs * 1e-12).max(1e-15);

[developing-human]

This is good feedback, thank you.

I had been thinking of having the function accept a value for epsilon (sympy does similar) but this would be a breaking change.

I think computing epsilon based on the values in the matrix is a good solution. I confirmed it works for my use case and pushed an update.

[djmaxus]

Thanks @developing-human and @Axect

[developing-human]

Happy new year @Axect. Do you need anything else from me to get this merged?