4 ms·

I think the real problem is that IEEE 754 floating point almost always gets treated as black box that magically provides the behavior of real numbers, which is

by nsajko 6y ago

I think the real problem is that IEEE 754 floating point almost always gets treated as black box that magically provides the behavior of real numbers, which is impossible. Instead of recognizing that the abstraction leaks in cases like these because the floating point representation only has so many bits, people just ignore the mathematical reality. The right thing to do is descending to the level of the representation (bits) and comparing ULP differences when necessary. The most annoying thing is how the bit-level representation of floating point is ignored even when relative error is significant instead of absolute error, because ULP differences are just that - relative error. EDIT: the epsilon fix in the test code in the blog post actually converts the relative error to absolute error so both could be done with the absolute error check. It would be more clear to instead do both checks explicitly. I think it would also be faster, which might matter in some situations.

GuB-42 6y ago
Even though floating point numbers are usually treated like a black box, most programmers know that they are imprecise, and the common advise is usually "never check for equality". Ok, fair enough, but how do you compare then? Standard libraries are rarely helpful here, and more often then not, you have to roll your own. Usually, that's poorly done. I've seen it done wrong way too many times, and that's by competent programmers. Why aren't there ULP-based functions in libmath for instance? I shouldn't have to cast float to int, do bit manipulations, and break portability for that. Do something like an optimized version of int ulpdiff(double a, double b) { int ret = 0; if (a > b) { while (a > b) { b = nextafter(b, a); ++ret; } } else { while (a < b) { a = nextafter(a, b); --ret; } } return ret; }
- nsajko 6y ago
  I used the following C code earlier this year (use at your own peril): https://github.com/nsajko/numericcompfricas/blob/master/check/checker.c#L234 https://github.com/nsajko/numericcompfricas/blob/master/chec... Distance between 0.0 and -0.0 is taken to be 0. If x or y are NaN, the distance is taken to be the greatest positive value of the return type. Replace mfloat_t with double. Replace int64 with long or similar. int64 ud(mfloat_t x, mfloat_t y) { // Handle NaNs. if (x != x || y != y) { return (int64)(~0UL >> 1); } union {mfloat_t X; uint64 a;} u1; union {mfloat_t Y; uint64 b;} u2; u1.X = x; u2.Y = y; uint64 a = u1.a, b = u2.b; if ((a & (1UL << 63)) == (b & (1UL << 63))) { // a and b are of the same sign. a = a - b; if ((a & (1UL << 63))) { return -a; } return a; } return (int64)(a - (1UL << 63) + b); } Using memcpy could be more portable than using unions, but I think it doesn't really matter. EDIT: The Musl libc, for example, also uses unions instead of memcpy: https://git.musl-libc.org/cgit/musl/tree/src/math/nextafterf.c https://git.musl-libc.org/cgit/musl/tree/src/math/nextafterf...
- nestorD 6y ago
  As usual, Rust has a crate for that: https://crates.io/crates/float-cmp https://crates.io/crates/float-cmp
  - vlovich123 6y ago
    That led me to https://crates.io/crates/efloat/0.2.0 https://crates.io/crates/efloat/0.2.0 which sounds even more useful as you can actually track what the error bounds are after you've done all your math. That sounds eminently more useful than just making epsilon math correct, no?
    - nsajko 6y ago
      That just tracks the error bounds on individual "small" operations. I think this is more like what you want, an interval arithmetic library: https://github.com/eseraygun/rust-honestintervals https://github.com/eseraygun/rust-honestintervals The point is that AFAIK you need the whole algorithm (like, e.g., logarithm) analysed to get sensible error bounds.
  - nsajko 6y ago
    That does not provide an ULP distance function, rather it tests for approximate equality while using both relative and absolute error margins. Also, a warning: the code behaves as if numbers of opposite sign are never "approximately equal". Also, the documentation sucks: it doesn't even mention the sign issue, nor does it mention how NaN values are handled, but it nevertheless contains irrelevant anecdotes.
dnautics 6y ago
> I think the real problem is that IEEE 754 floating point almost always gets treated as black box that magically provides the behavior of real numbers, which is impossible. That's true for integers, too, though the edge cases are fewer, lots of problems happen around catastrophic surface areas like overflow.

by nsajko 6y ago

GuB-42 6y ago

Even though floating point numbers are usually treated like a black box, most programmers know that they are imprecise, and the common advise is usually "never check for equality". Ok, fair enough, but how do you compare then? Standard libraries are rarely helpful here, and more often then not, you have to roll your own. Usually, that's poorly done. I've seen it done wrong way too many times, and that's by competent programmers. Why aren't there ULP-based functions in libmath for instance? I shouldn't have to cast float to int, do bit manipulations, and break portability for that. Do something like an optimized version of int ulpdiff(double a, double b) { int ret = 0; if (a > b) { while (a > b) { b = nextafter(b, a); ++ret; } } else { while (a < b) { a = nextafter(a, b); --ret; } } return ret; }

nsajko 6y ago
I used the following C code earlier this year (use at your own peril): https://github.com/nsajko/numericcompfricas/blob/master/check/checker.c#L234 https://github.com/nsajko/numericcompfricas/blob/master/chec... Distance between 0.0 and -0.0 is taken to be 0. If x or y are NaN, the distance is taken to be the greatest positive value of the return type. Replace mfloat_t with double. Replace int64 with long or similar. int64 ud(mfloat_t x, mfloat_t y) { // Handle NaNs. if (x != x || y != y) { return (int64)(~0UL >> 1); } union {mfloat_t X; uint64 a;} u1; union {mfloat_t Y; uint64 b;} u2; u1.X = x; u2.Y = y; uint64 a = u1.a, b = u2.b; if ((a & (1UL << 63)) == (b & (1UL << 63))) { // a and b are of the same sign. a = a - b; if ((a & (1UL << 63))) { return -a; } return a; } return (int64)(a - (1UL << 63) + b); } Using memcpy could be more portable than using unions, but I think it doesn't really matter. EDIT: The Musl libc, for example, also uses unions instead of memcpy: https://git.musl-libc.org/cgit/musl/tree/src/math/nextafterf.c https://git.musl-libc.org/cgit/musl/tree/src/math/nextafterf...
nestorD 6y ago
As usual, Rust has a crate for that: https://crates.io/crates/float-cmp https://crates.io/crates/float-cmp
- vlovich123 6y ago
  That led me to https://crates.io/crates/efloat/0.2.0 https://crates.io/crates/efloat/0.2.0 which sounds even more useful as you can actually track what the error bounds are after you've done all your math. That sounds eminently more useful than just making epsilon math correct, no?
  - nsajko 6y ago
    That just tracks the error bounds on individual "small" operations. I think this is more like what you want, an interval arithmetic library: https://github.com/eseraygun/rust-honestintervals https://github.com/eseraygun/rust-honestintervals The point is that AFAIK you need the whole algorithm (like, e.g., logarithm) analysed to get sensible error bounds.
- nsajko 6y ago
  That does not provide an ULP distance function, rather it tests for approximate equality while using both relative and absolute error margins. Also, a warning: the code behaves as if numbers of opposite sign are never "approximately equal". Also, the documentation sucks: it doesn't even mention the sign issue, nor does it mention how NaN values are handled, but it nevertheless contains irrelevant anecdotes.

dnautics 6y ago

> I think the real problem is that IEEE 754 floating point almost always gets treated as black box that magically provides the behavior of real numbers, which is impossible. That's true for integers, too, though the edge cases are fewer, lots of problems happen around catastrophic surface areas like overflow.