proposal: math: add an IsAlmostEqual function

[johnrs]

Proposal

Add an IsAlmostEqual function to the math package.

Justification

Although equating floating-point numbers is generally a bad idea,
sometimes it is necessary. For example, an iterative square root
algorithm might need an equate to determine when to quit.

Possible Solution

Add something like this function to the math package.

// IsAlmostEqual returns the near equality of 2 float64's based on
// a variable sized window, "nrSteps" big.
// A step's size is the amount by which adjacent numbers differ.
// The step size is based on the larger of the two numbers being compared.
//
// 0 steps requires an exact (generally not what you want) match.
// 1 or 2 steps is a common setting for simple calculations.
// Complicated calculations with many accumulated errors might require more steps.
//
//   Special cases are handled as follows (regardless of parameter order):
//       +Inf, +Inf  ==>  true
//       -Inf, -Inf  ==>  true
//       +Inf, x     ==>  false  (x != +Inf)
//       -Inf, x     ==>  false  (x != -Inf)
//        NaN, x     ==>  false
//
func IsAlmostEqual(a, b float64, nrSteps int) (equal bool) {
    if nrSteps < 0 {
        return false
    }

    // This allows matching infinities to pass.
    if a == b {
        return true
    }

    max := math.Max(math.Abs(a), math.Abs(b))
    stepSize := math.Nextafter(max, math.Inf(+1)) - max // max of either +Inf or NaN ==> NaN
    return math.Abs(a-b) <= float64(nrSteps)*stepSize   // NaN fails
}

[ianlancetaylor]

Personally I think this this requires too much knowledge of the floating point format to use correctly, and experts who need this function can easily write it themselves. My opinion could change if the math library for some other language provides a similar function.

In general I feel that this kind of function can be written in several different way, and which way to use requires understanding the algorithm you are trying to implement. That means that adding one such function to the standard library is an attractive nuisance for people who are not experts.

[johnrs]

While it requires some knowledge of floating-point to understand, I don't see it as hard to use. For example, to see if A and B (both float64's) are almost equal you would write:

    if math.IsAlmostEqual(A, B, 1) { ... }

I did look around about before settling on the implementation I proposed. Perhaps I didn't look hard enough, but this seemed to be what was used in the overwhelming number of cases I found.

I agree that it can be misused, but that's true of lots of things. :) I would argue that the name of the function serves as a fairly strong warning that it isn't the same as "==".

[griesemer]

The suspect IsAlmostEqual is is almost never what you want. If you know what you are doing in a numerical algorithm and understand float math, the function most likely does too much or too little and you end up writing your own specialized code anyway. On the other hand, for people unsure about floating point properties, IsAlmostEqual may lure them into a false sense correctness.

I'm not convinced this function carries its weight.

[johnrs]

When you say it might do too much or too little, I don't understand. Could you expand on this or perhaps give examples?

I would argue that == is almost never what you want. Someone who doesn't understand is probably going to start with ==, then use a fixed delta (ex: 0.000001). By definition, he isn't going to understand enough to do it correctly. So he's stuck.

Reference: Knuth, sec 4.2.2, page 217: A programmer using floating point arithmetic almost never wants to test if two computed values are exactly equal to each other (or at least he hardly ever should try to do so), because this is an extremely improbably occurrence.

[griesemer]

Agreed that == is usually wrong. But the proposed function does quite a bit of work, and I suspect ( I don't have proof) that more often than not it does too much (no need to take care of Inf or NaN), and as a consequence may be slower than necessary. It's also harder to analyze the numerical code's error behavior because one now has to understand the detailed properties of IsAlmostEqual. Testing that abs(x - y) < eps is in fact often enough for the right eps.

We have been reluctant to add functionality to math that's not essential, and this function doesn't seem to fall in that category either.

[johnrs]

I don't understand some of the points you are making.

Does too much work: There is but a single if needed to fully handle infinities, and nothing at all for NaN's. This doesn't seem like too much of a sacrifice to achieve correctness.

Analyzing the properties: How does someone writing their own routine make this easier?

Testing that abs(x - y) < eps is in fact often enough for the right eps.

That's exactly what this function does. It's the last line. The two prior lines determine the right (relative) eps. Are you suggesting that there is a better of doing this?

Not essential: I would argue that it's more essential than the floating point ==. Indeed, an extremest (not me!) might argue that == should not be allowed at all (or at least flagged by vet) for floating point numbers. Or replace == with this function. This would prevent/fix many bugs. For those rare times you actually needed an exact equate you could simply use a<=b && a>=b. :)

I would argue that it is an occassionally needed function. Since it's not obvious to many programmers how to implement it, they typically use a kludge involving a fixed eps - which is unusable in a library and is a landmine for refactoring.

[griesemer]

What I mean by "too much" is that floating-point algorithms that care about speed will want to do the absolute minimum - no need to deal with infinities or NaNs, and possible compute the epsilon beforehand (e.g. a loop) rather than each time. In my experience with floating-point arithmetic, conditions like this one are often required to determine termination conditions of loops.

I am not disagreeing with you that this is an "occasionally needed" function - but so are many functions. That's not a sufficient criteria for inclusion in a standard library.

Let's have other people chime in here; ideally folks with experience in numerical algorithms. I've simply stated my opinion, others may feel differently and may have concrete evidence to substantiate their claims.

[johnrs]

I see your point about calculating the eps once when there's a loop.

Looking forward to seeing usage feedback.

[ALTree]

Expanding over this:

floating-point algorithms that care about speed will want to do the absolute minimum

one should also note that IsAlmostEqual as it was proposed here won't be inlined by the current compiler (because is non-leaf); so good luck using it as the termination condition of a tight loop..

[as]

Regardless of whether adding this to stdlib is a good idea or not, what is the reason for using the word Almost in IsAlmostEqual? I would expect the function to be called Equal and the signature of the function to define the expected conditions for equality.

[johnrs]

Almost is what differentiates it from strict equality. It allows for limited precision of floating point.

[as]

@johnrs

func IsAlmostEqual(a, b float64, nrSteps int) (equal bool) {
func IsEqual(a, b float64, nrSteps int) (equal bool) {

For me the "Almost" adverb removes clarity; the third parameter to the func should clarify the conditions under which a is equal to b. The purpose of nrSteps isn't clear because it requires me to know something floating point representation. In other languages, the IsAlmostEqual functions either have a relative difference parameter (instead of nrSteps) or carry an implementation-defined relative difference . The func name is the same, but semantics differ, and my first impression was "what is nrSteps and why is it an int?".

[johnrs]

nrSteps: OK, I agree that "steps" perhaps isn't a great name, especially because perhaps I didn't define it well enough. Thus I would change

// A step's size is the amount by which adjacent numbers differ.

to

// A step's size (relative epsilon) is the amount by which adjacent numbers differ.

If that doesn't seem clear enough, then "nrSteps" could be changed to something like "nrRelativeEpsilons".

And I just noticed that I should have said that the window was 2 x "nrSteps" big, since it's double-sided.

Understanding this setting does require some understanding. I don't think that the function header is a good place for a tutorial on floating point, however. To make it easier for the average user I did show the common setting (1 or 2 steps).

This setting could be eliminated by locking it down (to 1 or 2), but that would limit the function's use. What might be best would be to make the setting optional. Perhaps 2 functions (one with and one without it), or kludge it by making it a variant (expecting just 0 or 1 instances).

The setting is an positive integer because the only window sizes that make sense are integral multiples of the relative epsilon. You are dealing with mantissa differences of a few bits, and you can't have a fraction of a bit.