Next floating point number

[Steven D'Aprano]

I'm looking for some way to get the next floating point number after any
particular float.

(Any mathematicians out there: I am aware that there is no "next real
number". But floats are not real numbers, they only have a finite
precision.)

According to the IEEE standard, there should be a routine like next(x,y)
which returns the next float starting from x in the direction of y.

Unless I have missed something, Python doesn't appear to give an interface
to the C library's next float function. I assume that most C floating
point libraries will have such a function.

So I came up with a pure Python implementation, and hoped that somebody
who had experience with numerical programming in Python would comment.



def nextfloat(x, y=None):
"""Return the next floating point number starting at x in the
direction of y. If y is not given, the direction is away from zero.
"""
from operator import add, sub
from math import frexp # Returns a tuple (mantissa, exponent)
x = float(x)
if y is None:
if x < 0.0:
op = sub
else:
op = add
elif y > x:
op = add
elif y < x:
op = sub
else:
raise ValueError("direction float is equal to starting float")
# We find the next float by finding the smallest float epsilon which
# is distinguishable from x. We know epsilon will be a power of two,
# because floats are implemented internally in binary.
# We get our initial estimate for epsilon from the floating point
# exponent of x.
epsilon = 2.0**(frexp(x)[1]) # epsilon > abs(x)
lasteps = epsilon # Save a copy of our last epsilon for later use.
while op(x, epsilon) != x:
lasteps = epsilon
epsilon /= 2.0
# When we get here, we've halved epsilon one time too many.
# We can't just double it, because that fails for the case where x is
# zero - epsilon will underflow to zero. So we need to save and use
# the previous iteration of epsilon.
return op(x, lasteps)



Thanks,

 

[Robert Kern]

Unless I have missed something, Python doesn't appear to give an interface
to the C library's next float function. I assume that most C floating
point libraries will have such a function.

http://www.mkssoftware.com/docs/man3/nextafter.3.asp

--
Robert Kern
(e-mail address removed)

"In the fields of hell where the grass grows high
Are the graves of dreams allowed to die."
-- Richard Harter

 

[Tim Peters]

[Steven D'Aprano]

I'm looking for some way to get the next floating point number after any
particular float. ....
According to the IEEE standard, there should be a routine like next(x,y)
which returns the next float starting from x in the direction of y.

Unless I have missed something, Python doesn't appear to give an interface
to the C library's next float function. I assume that most C floating
point libraries will have such a function.

While the C99 standard defines such a function (several, actually),
the C89 standard does not, so Python can't rely on one being
available. In general, Python's `math` module exposes only standard
C89 libm functions, plus a few extras it can reliably and portably
build itself on top of those. It does not expose platform-specific
libm functions. You can argue with that policy, but not successfully

So I came up with a pure Python implementation, and hoped that somebody
who had experience with numerical programming in Python would comment.

If you're happy with what you wrote, who needs comments ;-) Here's a
careful, "kinda portable" implementation in C:

http://www.netlib.org/toms/722

If you ignore all the end cases (NaNs, infinities, signaling underflow
and overflow, ...), the heart of it is just adding/subtracting 1
to/from the 64-bit double representation, viewing it as an 8-byte
integer. That works fine for the IEEE-754 float representations (but
does not work for all float representations).

I've used this simple routine based on that observation, which ignores
all endcases, and only works if both input and result are >= 0:

"""
from struct import pack, unpack

def next(x, direction=+1):
bits = unpack(">Q", pack(">d", x))[0]
return unpack(">d", pack(">Q", bits + direction))[0]
"""

For example,
9.9999999999999982e+099

 

[Bengt Richter]

I'm looking for some way to get the next floating point number after any
particular float.

(Any mathematicians out there: I am aware that there is no "next real
number". But floats are not real numbers, they only have a finite
precision.)

According to the IEEE standard, there should be a routine like next(x,y)
which returns the next float starting from x in the direction of y.

Unless I have missed something, Python doesn't appear to give an interface
to the C library's next float function. I assume that most C floating
point libraries will have such a function.

So I came up with a pure Python implementation, and hoped that somebody
who had experience with numerical programming in Python would comment.



def nextfloat(x, y=None):
"""Return the next floating point number starting at x in the
direction of y. If y is not given, the direction is away from zero.
"""
from operator import add, sub
from math import frexp # Returns a tuple (mantissa, exponent)
x = float(x)
if y is None:
if x < 0.0:
op = sub
else:
op = add
elif y > x:
op = add
elif y < x:
op = sub
else:
raise ValueError("direction float is equal to starting float")
# We find the next float by finding the smallest float epsilon which
# is distinguishable from x. We know epsilon will be a power of two,
# because floats are implemented internally in binary.
# We get our initial estimate for epsilon from the floating point
# exponent of x.
epsilon = 2.0**(frexp(x)[1]) # epsilon > abs(x)
lasteps = epsilon # Save a copy of our last epsilon for later use.
while op(x, epsilon) != x:
lasteps = epsilon
epsilon /= 2.0
# When we get here, we've halved epsilon one time too many.
# We can't just double it, because that fails for the case where x is
# zero - epsilon will underflow to zero. So we need to save and use
# the previous iteration of epsilon.
return op(x, lasteps)

I wonder if this won't work (for IEEE 754 double that is)

from math import frexp
def nextf(x, y):
f,e = frexp(x)
if (f==0.5 or f==-0.5) and x>=y: eps = 2.0**-54
else: eps = 2.0**-53
if x<y: return (f+eps)*2.0**e
else: return (f-eps)*2.0**e


Regards,
Bengt Richter

 

[Steven D'Aprano]

I wonder if this won't work (for IEEE 754 double that is)

from math import frexp
def nextf(x, y):
f,e = frexp(x)
if (f==0.5 or f==-0.5) and x>=y: eps = 2.0**-54
else: eps = 2.0**-53
if x<y: return (f+eps)*2.0**e
else: return (f-eps)*2.0**e

Doesn't that assume floats are 64 bits? Is it possible that they might
not be on some platform or another?

 

[Steven D'Aprano]

While the C99 standard defines such a function (several, actually),
the C89 standard does not, so Python can't rely on one being
available. In general, Python's `math` module exposes only standard
C89 libm functions, plus a few extras it can reliably and portably
build itself on top of those. It does not expose platform-specific
libm functions. You can argue with that policy, but not successfully
unless you take over maintenance of mathmodule.c <0.5 wink>.

*rueful smile*

If you're happy with what you wrote, who needs comments ;-)

Because I don't know _quite_ enough numeric programming to tell whether my
happiness is justified or just the happiness of somebody who hasn't
noticed the great big enormous bug staring him right in the face

Here's a
careful, "kinda portable" implementation in C:

http://www.netlib.org/toms/722

If you ignore all the end cases (NaNs, infinities, signaling underflow
and overflow, ...), the heart of it is just adding/subtracting 1 to/from
the 64-bit double representation, viewing it as an 8-byte integer. That
works fine for the IEEE-754 float representations (but does not work for
all float representations).

Will Python always use 64-bit floats?

I've used this simple routine based on that observation, which ignores
all endcases, and only works if both input and result are >= 0:

"""
from struct import pack, unpack

def next(x, direction=+1):
bits = unpack(">Q", pack(">d", x))[0] return unpack(">d", pack(">Q",
bits + direction))[0]
"""

Nice! I played around with struct.pack and unpack using a format string of
"f", but didn't get anywhere, so thanks for this.

 

[Bengt Richter]

I wonder if this won't work (for IEEE 754 double that is) ^^^^^^^^^^^^^^^^^^^[1]

from math import frexp
def nextf(x, y):
f,e = frexp(x)
if (f==0.5 or f==-0.5) and x>=y: eps = 2.0**-54
else: eps = 2.0**-53
if x<y: return (f+eps)*2.0**e
else: return (f-eps)*2.0**e

Doesn't that assume floats are 64 bits? Is it possible that they might
not be on some platform or another?

yes [1], and yes ;-)

I wonder if frexp is always available, and if the above would generalize
to something that could be self-tuning via a decorator or conditional defs.
ISTM I recall a FP format based on shifting 4 bits at a time to normalize,
and keeping an exponent as a multiple of 16, which might need a bunch
of epsilons. Obviously too, you wouldn't want to calculate things like 2.**-53
more than once. Also for a limited range of e integers, 2.**e is probably
faster looked up from a precalculated list p2[e]. The math module could
also expose an efficient multiply by a power of two using FSCALE if
the pentium FPU is there. And there's probably other stuff to take advantage
of, like doing it mostly without the FPU, a la Tim's struct thing. Specialized
C code to do the function could be pretty fast.
I'm a little curious what you are using this for.

Regards,
Bengt Richter

 

[Steven D'Aprano]

^^^^^^^^^^^^^^^^^^^[1]

Ah! I read that, and blanked it out of my mind :-(

I wonder if frexp is always available,

I believe so, since it is part of the math module.

[snip]

I'm a little curious what you are using this for.

Mostly curiosity -- I'm fascinated by numeric programming, but haven't got
the time to really sit down and learn it properly. But periodically I find
time to dabble for a few hours.

At a more practical level, for many numeric algorithms it is useful to
know the machine accuracy (which of course varies over the range of
floats). The machine accuracy at float x is nextfloat(x) - x. That gives
you an idea of how far apart floats are -- differences smaller than that
number are impossible for the algorithm to resolve.

 

[Tim Peters]

[Bengt Richter]

I wonder if frexp is always available,

Yes, `frexp` is a standard C89 libm function. Python's `math` doesn't
contain any platform-specific functions.

....

The math module could also expose an efficient multiply by a power
of two using FSCALE if the pentium FPU is there.

`ldexp` is also a standard C89 libm function, so Python also exposes
that. Whether it uses FSCALE is up to the platform C libm.

 

[Tim Peters]

[Steven D'Aprano]

...
Will Python always use 64-bit floats?

A CPython "float" is whatever the platform C compiler means by
"double". The C standard doesn't define the size of a double, so
neither does Python define the size of a float.

That said, I don't know of any Python platform to date where a C
double was not 64 bits. It's even possible that all _current_ Python
platforms use exactly the same format for C double.(the IEEE-754
double format), modulo endianness.

There's a subtlety related to that in my pack/unpack code, BTW: using
a ">d" format forces `struct` to use a "standard" big-endian double
encoding, which is 64 bits, regardless of how the platform C stores
doubles and regardless of how many bytes a native C double occupies.
So, oddly enough, the pack/unpack code would work even if the platform
C double used some, e.g., non-IEEE 4-byte VAX float encoding. The
only real docs about this "standard" encoding are in Python's
pickletools module:

"""
float8 = ArgumentDescriptor(
name='float8',
n=8,
reader=read_float8,
doc="""An 8-byte binary representation of a float, big-endian.

The format is unique to Python, and shared with the struct
module (format string '>d') "in theory" (the struct and cPickle
implementations don't share the code -- they should). It's
strongly related to the IEEE-754 double format, and, in normal
cases, is in fact identical to the big-endian 754 double format.
On other boxes the dynamic range is limited to that of a 754
double, and "add a half and chop" rounding is used to reduce
the precision to 53 bits. However, even on a 754 box,
infinities, NaNs, and minus zero may not be handled correctly
(may not survive roundtrip pickling intact).
""")
"""