1

I would like to know the epsilon of a floating point number around a given value.

std::numeric_limits<floating_point_type>::epsilon() provides that only for the number 1.0, while I would like a function to work on any number.

Is there any standard library solution to this? If not - how should I implement the function?

Eric Postpischil
  • 195,579
  • 13
  • 168
  • 312
Tim Kuipers
  • 1,705
  • 2
  • 16
  • 26
  • 1
    Mostly, the epsilon at `x` is `scalb(epsilon, ilogb(x))` or may be calculated as `nexttoward(x, INFINITY) - x`. But there are issues involving whether `x` is subnormal, zero, NaN, infinity, negative (for the latter expression), or a power of the base of the floating-point format. In the latter case, the distance to the next representable value is different in one direction than the other. Do you really need a complete solution? What are you trying to do? – Eric Postpischil Sep 13 '18 at 15:39
  • Take the number twice as large, divide the difference in half until it stops changing. – stark Sep 13 '18 at 15:39
  • @stark: That description is not clear. `double x = 1; double TwiceAsLarge = 2*x; double Difference = TwiceAsLarge - x; while (Difference != Difference/2) Difference /= 2;` produces zero. – Eric Postpischil Sep 13 '18 at 15:45
  • @EricPostpischil I am trying to see whether two values are 'equal'. The `std::nexttoward` can satisfy my needs. Thanks. – Tim Kuipers Sep 13 '18 at 15:45
  • 1
    @TimKuipers: Ascertaining whether two computers floating-point numbers would be equal if they had been computed using exact mathematics is an entirely different problem. That floating-point resolution at the numbers is not a particularly useful value. Correct code requires an analysis of the numerical error of the preceding computations and an evaluation of the effects of false positives. – Eric Postpischil Sep 13 '18 at 15:47
  • @Eric Anyone can deliberately misinterpret someone to make them look wrong. – stark Sep 13 '18 at 15:51
  • @EricPostpischil Such analysis would cost computation time, but in my application it's more important to be fast than to be accurate. The equation you provided gives me a quick and dirty upper bounds to the error I can expect. – Tim Kuipers Sep 13 '18 at 15:54
  • @TimKuipers: Initial analysis should be done at design time. – Eric Postpischil Sep 13 '18 at 15:58
  • Static analysis of FP error bounds does *not* cost computation time. `std::nextafter` is the answer to the question you asked, but I would never use it for approximate comparisons. – Sneftel Sep 13 '18 at 15:58
  • It's not possible to do static error bounds analysis if the amount of error depends on the magnitude of the input. – Tim Kuipers Sep 13 '18 at 16:00
  • @stark: The goal is not to make you look wrong but to point out that your text is unclear. Certainly people familiar with floating-point will be able to figure out what you mean, but we cannot expect that of people learning the subject. – Eric Postpischil Sep 13 '18 at 16:00
  • @TimKuipers It is! It generally comes down to multiplying some multiple of `epsilon` by the maximum magnitude of the inputs. That's likely to be cheaper than `nextafter`, since it doesn't involve a function call. – Sneftel Sep 13 '18 at 16:04
  • @Sneftel So instead of e.g. `6 * getResolutionAt(input_number)` I might just as well do `6 * input_number * std::numeric_limits::epsilon()`. That makes sense. – Tim Kuipers Sep 13 '18 at 16:15
  • 1
    Yep, though make it `input_number * (6 * epsilon)` to encourage the compiler to do constant folding and save you a multiply. – Sneftel Sep 13 '18 at 16:16

1 Answers1

5

Well, the easiest solution to find the epsilon immediately above the value (that is, the distance from that value to the next representable value) would just be 

std::nextafter(x, std::numeric_limits<floating_point_type>::infinity()) - x

Similarly to find the epsilon below the value, you could do 

x - std::nextafter(x, -std::numeric_limits<floating_point_type>::infinity())

Note that those two won't be the same if x is an exact power of two.

Now, there is one slight caveat there: the calculated epsilon above FLT_MAX will be infinity (arguably that's kind of the correct answer, but it doesn't quite match IEEE-754's rounding rules) and the epsilon above infinity will be NaN (which, well, I don't know how I feel about that). In all other cases, the result will be exact.

Tim Kuipers
  • 1,705
  • 2
  • 16
  • 26
Sneftel
  • 40,271
  • 12
  • 71
  • 104