Add RelativeEq trait

[ILyoan]

Fix for #5316
Current Implementation of ApproxEq trait uses absolute error method, and this method is working well for small float point numbers but not suitable for large numbers.
This patch provides RelativeEq trait and implementation of it uses relative error method that will be better for comparing big numbers.

[brendanzab]

@Erik-S alerted me to this problem today on IRC. If this is an inherent problem with ApproxEq, you are welcome to change its implementation to use this method, thereby removing the need for RelativeEq.

[brson]

Yeah, let's try to stop proliferating eq traits and merge these two.

[ILyoan]

All right I'll try merge these two.
Close this PR.

[ILyoan]

I have two options in my mind.

i) changing approx_eq_eps() like this

fn approx_eq_eps(&self, other: &f32, approx_epsilon: &f32) -> bool {
    // check absolute error.
    if (*self - *other).abs() < *approx_epsilon { true }    
    else {
        // check relative error.
        (*self - *other).abs() < fmax(self.abs(), other.abs()) * (*approx_epsilon)
    }
}

ii) add methods for relative error checking something like this

pub trait ApproxEq<Eps> {
    fn approx_epsilon() -> Eps;
    fn approx_eq(&self, other: &Self) -> bool;
    fn approx_eq_eps(&self, other: &Self, approx_epsilon: &Eps) -> bool;
+   fn relative_epsilon() -> Eps;
+   fn relative_eq(&self, other: &Self) -> bool;
+   fn relative_eq_eps(&self, other: &Self, approx_epsilon: &Eps) -> bool;
}

impl ApproxEq<f32> for f32 {
     #[inline(always)]
     fn approx_epsilon() -> f32 { 1.0e-6 }

     #[inline(always)]
     fn approx_eq(&self, other: &f32) -> bool {
         self.approx_eq_eps(other, &ApproxEq::approx_epsilon::<f32, f32>())
     }

     #[inline(always)]
     fn approx_eq_eps(&self, other: &f32, approx_epsilon: &f32) -> bool {
         (*self - *other).abs() < *approx_epsilon
     }

+   #[inline(always)]
+   fn relative_epsilon() -> f32 { 1.0e-6 }
+
+   #[inline(always)]
+   fn relative_eq(&self, other: &f32) -> bool {
+       self.relative_eq_eps(other, &ApproxEq::relative_epsilon::<f32, f32>())
+   }
+
+   #[inline(always)]
+   fn relative_eq_eps(&self, other: &f32, relative_epsilon: &f32) -> bool {
+       (*self - *other).abs() < fmax(self.abs(), other.abs()) * (*relative_epsilon)
+   }
}

Thoughts?

[brendanzab]

I would like to see what @Erik-S has to say on this (hint hint!). The folks on #4819 might be interested too. Ideally it's be nice to have a simpler trait. Could you explain the reasoning behind having a separate relative_eq method?

[Erik-S]

In my opinion, these two properties should hold:

approx_eq(1e-12,2e-12) = false

approx_eq(1e15,1e15+1) = true

The first is false because those numbers are "far apart" (one is twice the other). The second is true because they are "close" - 2 units in the least significant bit, if I did my math right (and assuming f64).

The second property shouldn't cause much trouble. The first could cause trouble if there is subtractive cancellation, which can lead to large relative errors. But if an algorithm has catastrophic subtractive cancellation, the algorithm is broken (usually).

Like most of floating point math, there is no one right answer. We need to provide a reasonable solution, and then people can code their own algorithms if they need more.

[ILyoan]

Because I think just one method couldn't be a solution for all kind of problem.
I thought there might be a situation that is not suitable with relative error method and providing many options for programmer could relax this kind of conflicts because they can choose proper method.
But if you prefer simpler trait, I don't insist providing many methods. I'm okay for any approach.

[brendanzab]

What are your thoughts @pcwalton, seeing as you use ApproxEq in rust_geom?

[Erik-S]

Part of my dislike for the relative_eq() - style math (not scaling eps by the input magnitude) is that there is literally no value of eps that will be reasonable across even a small part of the range of floating point numbers (roughly 1e-300 to 1e300 - 600 orders of magnitude). The approx_eq() (scaling eps by magnitude) still doesn't have one right value for eps, but you can define a reasonable value over the whole range.

We could provide a version of relative_eq() that takes eps as an argument, but that starts getting complicated. This would be a great place to use default values on arguments, if we could.

[ILyoan]

@bjz @Erik-S Thanks for your comments.
On second thought, I prefer the first option. How about this style?
https://gist.github.com/ILyoan/5538590
This separates the absolute and relative epsilons, and check the equality based on these values. so you can even cancel relative checking by giving absolute_epsilon as proper value and relative_epsilon 0.
By default value, approx_eq() will hold Erik-S's properties.

[brendanzab]

I have concerns about the duplicated computations in your gist. @pcwalton and I are going to be using these for matrix and vector comparisons, so efficiency is important.

Would this be slightly better?

fn approx_eq_eps(&self, other: &f32, absolute_epsilon: &f32, relative_epsilon: &f32) -> bool {
        match (*self - *other).abs() {
            // check absolute error
            n if n < *absolute_epsilon => true,
            // check relative error
            n => n < self.abs().max(&other.abs()) * (*relative_epsilon),
        }
    }
}

You'd also have to rename approx_eq_eps, as unfortunately (or fortunately?) Rust doesn't support method overloading. I reiterate @Erik-S when he says that this is a good use case for default arguments.

[dobkeratops]

i like the function that takes epsilon as a parameter - you'd want to rely on the compiler inlining constants

[Erik-S]

@ILyoan That looks pretty good. I'd recommend deleting absolute_epsilon() (since the default isn't used, and will rarely be correct).

I'd recommend setting eps to:
f32: 7.62939453125e-6 (2^-17)
f64: 1.4210854715202003717422485351563e-14 (2^-46)
Both are 64 * the machine epsilon.

[brendanzab]

@Erik-S could we just do this to be self-documenting?:

#[inline(always)] fn approx_epsilon() -> f32 { 64f32 * Float::epsilon() }
#[inline(always)] fn approx_epsilon() -> f64 { 64f64 * Float::epsilon() }

It would hopefully be inlined and evaluated at compile time by the compiler (wishing for associated constants here...).

Also, why the choice of 64? Is that arbitrary?