[Python-ideas] PEP 485: A Function for testing approximate equality
Andrew Barnert
abarnert at yahoo.com
Sat Jan 24 00:09:32 CET 2015
On Friday, January 23, 2015 1:31 PM, Nathaniel Smith <njs at pobox.com> wrote:
> I'd strongly consider expanding the scope of this PEP a bit so that
> it's proposing both a relative/absolute-error-based function *and* a
> ULP-difference function. There was a plausible-looking one using
> struct posted in the other thread
If you're thinking about the post I think you are (mine), I wouldn't suggest using that.
The main issue is that was a bits-difference function, not an ulps-difference function--in other words, bits_difference(x, y) is the number of times you have to do y = nexttoward(y, x) to get y == x. In the case where the difference is <= 2, or where x and y are finite numbers with the same sign and exponent, they happen to be the same, but otherwise, they don't. For example, consider x as the 5th largest number with one exponent, and y as the 5th smallest number with the next. They're 10 bits away, but 12.5 ulp(x) away and 7.5 ulp(y) away. Most algorithms that you want to test for ulp difference are specified to be within 0.5, 1, or 2 ulp, and for C lib functions it's always 0.5 or 1 (except pow in certain cases), or to no more than double the ulp difference--but definitely not _all_, so it would be misleading to offer a bits-difference function as an ulps-difference function.
Secondarily, even as a bit-difference function, what I posted isn't complete (but I think the version in https://github.com/abarnert/floatextras is), and makes various decisions and assumptions that aren't necessarily the only option. Also, there's nothing else in the stdlib that directly accesses the bits of a float in Python, which seems a little weird. Finally, neither Python nor the C89 standard that CPython implies require that float actually be an IEEE 754-1985 double (much less an IEEE 754-2005 binary64, the later standard I actually have a copy of...). In particular, sys.float_info doesn't assume it.
I think if we wanted this, we'd want to implement nexttoward in C (by calling the C99/POSIX2001 function if present, and maybe our own bit-twiddling-IEEE-doubles-in-C implementation for Windows, but it's not there otherwise), then define ulp (in C or Python) in terms of nexttoward, then define ulp_difference(x, y) (ditto) in terms of ulp(y). This does require a bit of care to make sure that, e.g., ulp_difference(float_info.max, inf) comes out as 1 or as an error, whichever one you want, and so on. (That means it also requires deciding what to do for each edge case, since they're not standardized by IEEE 754-1985, IEEE 754-2008, C99, or POSIX2001.) This would work correctly and consistently on almost every *nix platform (even some that don't use IEEE double) and on Windows, and wouldn't exist on platforms where it won't work correctly.
Of course other implementations would have to come up with some other compatible implementation, but at least Java has an ulp function, and if .NET doesn't, it can probably make assumptions about the underlying platform.
If we also want a bits_difference function in the stdlib (and I'm not sure we do), I'd suggest also writing that in C, by pointer-casting from double to int64_t and using the information in C99 math.h/limits.h (and again maybe special-casing Windows), rather than twiddling IEEE bits in Python.
More information about the Python-ideas
mailing list