IEEE 754 nextup()/nextdown() of a floating point number

Users
1

Does anybody need this IEEE 754 functionality?

2

yes! I use them often when doing interval arithmetic stuff.

I have my own implementation which I use (hardcoded for float64, but generalizing should be simple). It's based on algorithm 5 from this paper by Rump: https://www.tuhh.de/ti3/paper/rump/RuOgMoOi16.pdf

 const FMAX       = 0x1.FFFFFFFFFFFFFp+1023;
    const EPS        = 0x1p-53;
    const PHI        = 0x1.0000000000001p-53;
    const EPSINV     = 0x1p53;
    const ETA        = 0x1p-1074;
    const FMIN2      = 0x1p-1021;
    const EPSINVFMIN = 0x1p-969;

    /*
    :arg c: ``real``

    :returns: The next floating point number after ``c``. ``NAN`` and ``INFINITY`` are left unchanged.
    */
    proc nextup(c : real) {
        if abs(c) >= EPSINVFMIN then
            return if c == -INFINITY then -FMAX else c + PHI * abs(c);
        else if abs(c) < FMIN2 then return c + ETA;
        else {
        var C = EPSINV * c,
            e = PHI * abs(C);
        return EPS * (C + e);
        }
    }

    /*
    :arg c: ``real``

    :returns: The previous floating point number before ``c``. ``NAN`` and ``-INFINITY`` are left unchanged.
    */
    proc nextdown(c : real) {
        if abs(c) >= EPSINVFMIN then
            return if c == INFINITY then FMAX else c - PHI * abs(c);
        else if abs(c) < FMIN2 then return c - ETA;
        else {
        var C = EPSINV * c,
            e = PHI * abs(C);
        return EPS * (C - e);
        }
    }

here's the link to the code: Chimera/src/utils.chpl at main · lucaferranti/Chimera · GitHub (word of warning though, that's a very early prototype of interval arithmetic which I've never found the time/energy to come back too, so if you are looking at it, set your expectations very low and then lower them a bit more :slight_smile: I do hope to revamp it this year at some point)

I thought it would be nice to have those in chapel (not necessarily my implementation, but some) but never found the time to properly open an issue about it.

3

If I get any other takers, I will put a generic implementation on Github and then somebody can decide what they want to do. Otherwise, you are more than welcome to mine which, because it exploits Brad's .transmute() method, is a lot leaner than Rump's.

4

where can I find your implementation? (haven't been in discourse for a while, so maybe I missed some old post of yours)

5

I never posted it because I did not want to be the source of useless noise. You are the only other person who has expressed any interest.

6

In the last post of mine in #26933, there is a link to ATO for a test program which contains my nextup() routine. Feel free to play with it. This routine requires the precision() routines which appear early in the listing. The Unit Test cases need some of the other miscellaneous routines so you can probably ignore these.

Notice that the routine likely() inside nextup() is a NO-OP at this stage in that code. So, instead of the current line which says:

return if likely((_t & _E) < _E) then

You probably only need to have:

return if (_t & _E) < _E then

for your purposes.

Internally the nextup() routine uses some limiting cases of the use of signed and unsigned integer arithmetic which I needed to exploit to achieve optimal code.

The purpose of #26933 is to see whether we can use the concept of likely and unlikely attributes in Chapel without having to actually implement them as a language change. And from Brad's testing, it appears so. Whether we see the same optimal code as you see in a C program using such an approach is still being investigated (but I am leaving that to those much smarter than I am).

7

This is far less complicated than my earlier version and is less than 1% slower.

/*
 * Return the smallest floating-point number of the
 * identical type as x that compares greater than 'x'
 */
inline proc nextup(x : real(?w))
{
    /*
     * map any input of
     * a) a -0 to +0 or
     * b) a signaling NaN to a quiet NaN
     * C) anything else to itself
     */
    const t = x + 0;

    if t == t && t < inf then
    {
        /*
         * compute the IEEE 754 encoding of 't'
         */
        const _t = t.transmute(uint(w));
        /*
         * exploiting sign-carry with a right-shift on int(w),
         * compute (-1) if the hard negative field of t is 1,
         * or +1 otherwise; and then add this to the encoding
         */
        const _dt = ((_t:int(w) >> (w - 1)) | 1:int(w)):uint(w);

        return (_t + _dt).transmute(real(w));
    }
    return t;
}

Notice my edit which eliminated the likely and based on advice from @mppf, 2 params