# Implementing std::nextafter: Should denormals-are-zero mode affect it? If so, how?

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.

C's 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.5 should 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 +0.0 or -0.0 was -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 when x equals y, as required by POSIX, and by the C11 standard (7.12.11.3 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.

The 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 nextafterf(float,float) / 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).

How does 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 floats and doubles (but not 80bit long double for nextafterl, or for the direction arg of nexttowards).

• I can't speak to how often nextafter is used generally, but personally I've found it essential for stress testing geometric predicates (orientation, incircle). For example the supporting code for this paper uses their own hand-rolled implementation of nextafter; I'm not sure why they didn't use the libc version. – Daniel Shapero Feb 23 '16 at 15:55
• I had to ponder the same question in the context of the CUDA standard library. The NVIDIA GPUs that CUDA runs on provide a FTZ (flush-to-zero) mode, for single-precision operations only, that essentially combines the DAZ/FTZ modifiers found on x86. Since the standard says that nextafter should return the next "representable" value, in FTZ mode CUDA's nextafterf() jumps from +/-0 to +/- FLT_MIN. The motivation for my choice: it should be possible to step, one-by-one, through all representable floating-point numbers with nextafterf(x,pos_infinity) without getting stuck somewhere. – njuffa May 28 '16 at 23:51
• @njuffa: That's the behaviour I think makes the most sense in DAZ mode. Current glibc doesn't do that, though. Also, detecting DAZ mode and behaving differently would actually require something clever. Like maybe check input == output in the x != y branch. ( My WIP implementation is only affected by DAZ if x and y are both denormals: they compare equal so we return y). Otherwise it's all integer operations. We don't get stuck, just "pause at zero" for a long time. – Peter Cordes May 29 '16 at 1:43
• Any idea if changing the behaviour of glibc nextafter would be a problem for anyone? Probably it's important to have nextafterf and nextafterl behave the same, and also for the behaviour to be the same across platforms. (Currently I only have an x86-64 asm proof-of-concept implementation; will probably have to work backwards from that to get portable C that compiles to something not bad, if I go with mostly scalar rather than doing integer operations in SSE float/vector registers.) – Peter Cordes May 29 '16 at 1:46
• @PeterCordes I would not expect much input from users, all flavors of nextafter are rarely used. Occasional use to bump results up and down one ulp, e.g. for the proper handling of intervals, is all I would expect, and uses near 0 are probably rarer still. I never got any feedback on the behavior of CUDA's nextafterf, other than a question from the OpenCL folks about the "jump" behavior in FTZ mode; my explanation (same as above) wasn't challenged, which I took as tacit agreement with my reasoning. FTZ mode on GPUs is (like the rounding modes) something known at compile time, not dynamic. – njuffa May 29 '16 at 2:58