Approximate equality for floating-point

Hi all --

I'm planning to propose the following additions to FloatingPoint:

extension FloatingPoint {
  /// Test approximate equality with relative tolerance.
  ///
  /// Do not use this function to check if a number is approximately
  /// zero; no reasoned relative tolerance can do what you want for
  /// that case. Use `isAlmostZero` instead for that case.
  ///
  /// The relation defined by this predicate is symmetric and reflexive
  /// (except for NaN), but *is not* transitive. Because of this, it is
  /// often unsuitable for use for key comparisons, but it can be used
  /// successfully in many other contexts.
  ///
  /// The internet is full advice about what not to do when comparing
  /// floating-point values:
  ///
  /// - "Never compare floats for equality."
  /// - "Always use an epsilon."
  /// - "Floating-point values are always inexact."
  ///
  /// Much of this advice is false, and most of the rest is technically
  /// correct but misleading. Almost none of it provides specific and
  /// correct recommendations for what you *should* do if you need to
  /// compare floating-point numbers.
  ///
  /// There is no uniformly correct notion of "approximate equality", and
  /// there is no uniformly correct tolerance that can be applied without
  /// careful analysis. This function considers two values to be almost
  /// equal if the relative difference between them is smaller than the
  /// specified `tolerance`.
  ///
  /// The default value of `tolerance` is `sqrt(.ulpOfOne)`; this value
  /// comes from the common numerical analysis wisdom that if you don't
  /// know anything about a computation, you should assume that roughly
  /// half the bits may have been lost to rounding. This is generally a
  /// pretty safe choice of tolerance--if two values that agree to half
  /// their bits but are not meaningfully almost equal, the computation
  /// is likely ill-conditioned and should be reformulated.
  ///
  /// For more complete guidance on an appropriate choice of tolerance,
  /// consult with a friendly numerical analyst.
  ///
  /// - Parameters:
  ///   - other: the value to compare with `self`
  ///   - tolerance: the relative tolerance to use for the comparison.
  ///     Should be in the range (.ulpOfOne, 1).
  ///
  /// - Returns: `true` if `self` is almost equal to `other`; otherwise
  ///   `false`.
  @inlinable
  public func isAlmostEqual(
    to other: Self,
    tolerance: Self = Self.ulpOfOne.squareRoot()
  ) -> Bool {
    // tolerances outside of [.ulpOfOne,1) yield well-defined but useless results,
    // so this is enforced by an assert rathern than a precondition.
    assert(tolerance >= .ulpOfOne && tolerance < 1, "tolerance should be in [.ulpOfOne, 1).")
    // The simple computation below does not necessarily give sensible
    // results if one of self or other is infinite; we need to rescale
    // the computation in that case.
    guard self.isFinite && other.isFinite else {
      return rescaledAlmostEqual(to: other, tolerance: tolerance)
    }
    // This should eventually be rewritten to use a scaling facility to be
    // defined on FloatingPoint suitable for hypot and scaled sums, but the
    // following is good enough to be useful for now.
    let scale = max(abs(self), abs(other), .leastNormalMagnitude)
    return abs(self - other) < scale*tolerance
  }
  
  /// Test if this value is nearly zero with a specified `absoluteTolerance`.
  ///
  /// This test uses an *absolute*, rather than *relative*, tolerance,
  /// because no number should be equal to zero when a relative tolerance
  /// is used.
  ///
  /// Some very rough guidelines for selecting a non-default tolerance for
  /// your computation can be provided:
  ///
  /// - If this value is the result of floating-point additions or
  ///   subtractions, use a tolerance of `.ulpOfOne * n * scale`, where
  ///   `n` is the number of terms that were summed and `scale` is the
  ///   magnitude of the largest term in the sum.
  ///
  /// - If this value is the result of floating-point multiplications,
  ///   consider each term of the product: what is the smallest value that
  ///   should be meaningfully distinguished from zero? Multiply those terms
  ///   together to get a tolerance.
  ///
  /// - More generally, use half of the smallest value that should be
  ///   meaningfully distinct from zero for the purposes of your computation.
  ///
  /// For more complete guidance on an appropriate choice of tolerance,
  /// consult with a friendly numerical analyst.
  ///
  /// - Parameter absoluteTolerance: values with magnitude smaller than
  ///   this value will be considered to be zero. Must be greater than
  ///   zero.
  ///
  /// - Returns: `true` if `abs(self)` is less than `absoluteTolerance`.
  ///            `false` otherwise.
  @inlinable
  public func isAlmostZero(
    absoluteTolerance tolerance: Self = Self.ulpOfOne.squareRoot()
  ) -> Bool {
    assert(tolerance > 0)
    return abs(self) < tolerance
  }
  
  /// Rescales self and other to give meaningful results when one of them
  /// is infinite. We also handle NaN here so that the fast path doesn't
  /// need to worry about it.
  @usableFromInline
  internal func rescaledAlmostEqual(to other: Self, tolerance: Self) -> Bool {
    // NaN is considered to be not approximately equal to anything, not even
    // itself.
    if self.isNaN || other.isNaN { return false }
    if self.isInfinite {
      if other.isInfinite { return self == other }
      let scaledSelf = Self(signOf: self, magnitudeOf: 1)
      let scaledOther = Self(sign: other.sign, exponent: -1,
                          significand: other.significand)
      return scaledSelf.isAlmostEqual(to: scaledOther, tolerance: tolerance)
    }
    // If self is finite and other is infinite, flip order and use scaling
    // defined above, since this relation is symmetric.
    return other.rescaledAlmostEqual(to: self, tolerance: tolerance)
  }
}

Some quick obvious responses to questions not directly addressed in the comments above:

Q: Why not change the behavior of ==?
A: Approximate equality is useful to have, but is a terrible default because it isn't transitive, and does not imply substitutability under any useful model of floating-point values. This breaks all sorts of implicit invariants that programs depend on.

Q. Why not introduce a new operator? a ≅ b? a == b ± t?
A. There's a bunch of clever things that one might do in this direction, but it's not at all obvious which (if any) of them is actually a good idea. Defining good semantics for approximate equality via a method in the standard library makes it easy for people to experiment with clever operators if they want to, and provides very useful functionality now.

Q: Why is a relative tolerance used?
A: Because--as a library--we do not know the scaling of user's inputs a priori, a relative tolerance is the only sensible option. Any absolute tolerance will be wrong more often than it is right.

Q: Why do we need a special comparison for zero?
A: Because zero is the one and only value with which a relative comparison tolerance is essentially never useful. For well-scaled values, x.isAlmostEqual(to: 0, tolerance: t) will be false for any valid t.

Q: Why not a single function using a hybrid relative/absolute tolerance?
A: Because--as a library--we do not know the scaling of user's inputs a priori, this will be wrong more often than right. The scale invariance (up to underflow) of the main isAlmostEqual function is a highly-desirable property, because it mirrors the way the rest of the floating-point number system works.

21 Likes

The comments, before considering the implementations, make this a worthwhile addition, in my opinion.

2 Likes

Big +1 from me. I’ve got similar behaviours in my code which would be more elegant written like this, and even more so if they were in the Standard Library.

Btw I’m stealing them until they turn up :shushing_face::rofl:

1 Like

It's so very easy to get floating point "equality" wrong (I usually refer to this page), so this would be a welcome addition. Plus, we have already implemented this in our own codebase (with potential subtle bugs?), so it would be good to be able to rely on a standard implementation.

Just a question purely out of curiosity, since I'm not well-versed in numerical analysis enough, what does this code do in the case of subnormal numbers? My understanding was that these introduce additional headaches because the standard error margins don't apply.

The only semi-sensible thing is to apply a relative tolerance as though the scaling were that of the the smallest normal number (because the rounding point for all subnormals is in the same location as it is for that number). That's what this is doing:

let scale = max(abs(self), abs(other), .leastNormalMagnitude)

The scale relative to which the relative tolerance is applied is the greater in magnitude of self and other, or .leastNormalMagnitude if both are subnormal or zero. This handles both rolling off the tolerance for subnormals and also makes the relation symmetric.

2 Likes

It would be useful if there is a way to use approximately equal on switch case.

Could this be done by providing an API on FloatingPoint that returns a Range<Self> whose bounds are ± the tolerance around a particular value? Then you could do something like:

switch value {
case 1.0.approximateRange():
  ...
case 2.0.approximateRange(tolerance: whatever):
  ...
default:
  ...
}
1 Like

Caveat: I haven't thought very carefully about this but...

I'd probably do it with a wrapper type i.e.

struct Approximately<T: FloatingPoint> {
  let value: T
  init(_ v: T) { value = v }
  
  // I can never remember which is l and which is r here....
  static func ~= (lhs: Approximately, rhs: T) -> Bool {
    return lhs.value.isAlmostEqual(to: rhs)
  }
}

switch 0.1 + 0.2 {
case Approximately(0.3): print("yay!")
default: print("boo!")
}
5 Likes