correct round of reals?

[Tim Peters]

[Peter Schneider-Kamp, wants an IEEE nearest/even rint]

[Fredrik Lundh]
> def rint(v, m=4503599627370496.0):
>     if abs(v) > m:
>         return v
>     if v > 0:
>         return (v + m) - m
>     else:
>         return (v - m) + m

That's pretty!  I like it.

> (this assumes IEEE doubles)

And inherits the rounding mode in effect.  Unclear whether Peter wants that,
or wants to force nearest/even.

[back to Peter]
> Thanks. Do you know if e.g. VC++ uses IEEE doubles?
>
> The platforms this code has to run on are those that
> do not have an rint(3) function. So a connected
> question is if these platforms use IEEE doubles?

Almost all current platforms use IEEE doubles.  /F's method can be
generalized by preceding it with, e.g.,

# Find largest positive integet _m (in float format) that's
# a power of 2 and for which _m + 1. is exactly representable.

_m = 1.0
i = 0
while (_m + 1.0) - _m == 1.0 and i < 500:
    # _m's not too big yet -- boost it.
    _m = _m * 2.0
    i = i + 1
# Now _m is twice too big, or we never found a limit.
_m = _m * 0.5

if not 10 < i < 120:
    raise SystemError("The floating-point on this machine "
                      "looks bizarre")
del i

and replacing the "def" line via:

def rint(v, m=_m):

The new block of code is executed only upon first module import, so this is
just as efficient in use, but avoids /F's 53-bit mantissa assumption.
However, this rint inherits the rounding discipline used by the platform fp,
and, as above, it's unclear (to me) whether that's what you want.

If you want to force the issue, and be platform-independent, life is more
complicated, in your choice of two ways:  you can write a hairier algorithm
(as I've done below), or you can write a whole bunch of simple (but obscure)
algorithms #ifdef'ed on platform-specific symbols (they're "simple but
obscure" by virtue of exploiting platform quirks, just as /F's simple but
obscure routine exploits 754 defaults).

best-of-all-is-to-learn-to-love-truncation<0.9-wink>-ly y'rs  - tim


from math import modf

def rint(x):
    """x -> nearest integer, as a float.  Halfway to nearest even.

    Caution:  If x is NaN or an infinity, this is as unpredictable as
    the ISO/ANSI C routines it reiies on.  If abs(x) <= 0.5, a zero is
    returned, but whether the sign bit is preserved is also inherited
    from the platform C.  Ditto for the effect on the 754 "inexact"
    flag.

    In return for those minor insecurites, this computes 754 nearest/
    even rounding on non-754 platforms too, provided C's double is
    a binary floating-point format.
    """

    f, i = modf(x)
    absf = abs(f)
    if absf >= 0.5:
        i = i + (f < 0 and -1. or 1.)
        # Now i is the "add a half and chop" result.  This is the
        # final result unless it's a halfway case.
        # Then:  If i was odd before, this is the correct nearest/even
        # result, and its last bit is clear now.  But if i was even
        # before, its last bit is set now, and needs to be cleared --
        # & since i was even, adding 1 didn't propagate a carry.  So,
        # in either case, clearing the last bit "works" (either leaves
        # a good result unchanged, or repairs a bad result).
        if absf == 0.5:
            i = modf(i*0.5)[1]*2. # i.e., clear last bit
    return i

assert rint(0.) == rint(.4) == rint(.5) == rint(-.4) == rint(-.5) == 0.
assert rint(.6) == rint(.9) == rint(1.) == rint(1.4) == 1.
assert rint(1.5) == rint(1.9) == rint(2.) == rint(2.4) == rint(2.5) == 2.
assert rint(-.6) == rint(-1.) == rint(-1.4) == -1.
assert rint(-1.5) == rint(-1.9) == rint(-2.4) == rint(-2.5) == -2.
assert rint(1.23e300) == 1.23e300
assert rint(-1.23e300) == -1.23e300