This might be the wrong stackexchange site for this question. math.SE, cs.SE, programmers.SE, and of course stackoverflow are all possibilities. I'm hoping to reach an audience that might actually use this function, to get some feedback on how they'd like it to work.
double nextafter(double x, double y) function returns the next representable value from x, in the direction of y. (commonly + or - infinity, but can be anything).
I'm tidying up GNU libc/libm's implementation to compile to nicer code. (It currently extracts each double to two 32bit integers. It's pretty clunky and uses lots of branches, and doesn't take advantage of
int64_t.). If I'm going to change it at all, I might as well make as much improvement as possible, but using an FP compare-equal would change the results in DAZ mode.
How should it behave when the denormals-are-zero and/or flush-to-zero are active? On x86, there are two separate flags which sacrifice gradual underflow for performance. By default, neither is set. One or both can be set (they're often used together, but don't have to be):
DAZ is an input filter: When an FP math / compare instruction reads its inputs, denormals are considered +/-0.0. So a compare between two denormals finds they're equal. Arithmetic can easily produce denormal results, though.
FTZ is an output filter: Denormal inputs work normally, with no effect on compare instructions. However, ordinary math instructions can't generate denormals. Denormal results are flushed to +/-0.0.
DBL_MIN/2 + DBL_MIN/1.5should give the same result as without FTZ (e.g. if those input denormals are loaded as constants).
In both cases (esp DAZ), one might argue that the next representable value below
DBL_MIN (the smallest normal number) was
+0.0, and that the next number below
-DBL_MIN. Implementing this would require checking the FTZ/DAZ flags.
One might argue instead that in DAZ mode, two different denormals should be considered equal, because they do compare equal. Thus, it's correct to return
y, as required by POSIX, and by the C11 standard (126.96.36.199 pg 256 and Annex F.10.8.3 pg 529). (The motivation for this specification is for
nextafter(0.0, -0.0) to return
-0.0, and vice versa.)
The 2nd behaviour allows an extremely fast & compact implementation (no taken branches and 3 not-taken branches polluting the branch-predictor, and ~9 clock cycles latency on Intel Haswell (for comparison, an FP multiply takes 5 cycles) for the common case (x!=y, and x!=0, and the result doesn't overflow(infinity) or underflow (denormal)). Less than half the total code size compared to glibc's current version.
The C11 standard says:
"Even though underflow or overflow can occur, the returned value is independent of the current rounding direction mode."
But DAZ isn't a rounding direction mode.
boost::math::ulp function's documentation observes that "our experience is that
std::nextafter often breaks depending which optimizations and hardware flags are in effect".
Consistency with other implementations is desirable for glibc's
nextafter function, but consistency with glibc's existing behaviour may take precedence. The current implementation is not influenced at all by any FP settings, because it uses only integer comparisons on the IEEE float bit patterns. Consistency with
nexttoward, and with
nextafterl(long double, long double) is also important. I haven't looked into the long double versions yet. Integer compares are clunkier there, because
int64_t can't hold 80 bits.
Using integer-only compares is only slightly slower / bulkier to implement (using a 64bit integer compare, and the sign of x and y, instead of just an FP compare and the sign of x). IEEE floats are cleverly designed so their bits can be compared as sign-magnitude integers (and thus as 2's complement integers with a check for the both-negative case).
nextafter behave on other systems (OS X, MSVC++, anything else)? If people are interested in testing it, I might write a little test program.
How do you wish
nextafter behaved? (
nextafter is pretty much fully specified except for DAZ/FTZ, so wishes other than that might be better fulfilled by different functions.)
Did I miss anything?
I know the current implementation is "good enough" for such a rarely-used function. It's ugly, but it's used so rarely that it doesn't have much impact. (Actually, is it rarely used? Has anyone ever seen code that used it a lot?)
Part of the motivation here is that I enjoy optimizing things in assembly language, and integer ops on FP data is interesting, especially on x86-64 where integer instructions can be used on the vector registers that are also used for
doubles (but not 80bit
long double for
nextafterl, or for the direction arg of