5

I have a Python float, and I want to have the floats which are 1 ULP greater and lesser.

In Java, I would do this with Math.nextUp(x) and Math.nextAfter(x, Double.NEGATIVE_INFINITY).

Is there a way to do this in Python? I thought about implementing it myself with math.frexp and math.ldexp but as far as I know Python doesn't specify the size of floating point types.

rlibby
  • 5,931
  • 20
  • 25
  • 1
    Same question was asked here: http://bytes.com/topic/python/answers/739926-next-float, with good answers... – Benjamin Apr 21 '11 at 20:05

2 Answers2

5

Update: In Python 3.9+ there is math.nextafter():

>>> import math
>>> x = 4
>>> math.nextafter(x, math.inf)
4.000000000000001

Old answer:

You could look at how Decimal.next_plus()/Decimal.next_minus() are implemented:

>>> from decimal import Decimal as D
>>> d = D.from_float(123456.78901234567890)
>>> d
Decimal('123456.789012345674564130604267120361328125')
>>> d.next_plus()
Decimal('123456.7890123456745641306043')
>>> d.next_minus()
Decimal('123456.7890123456745641306042')
>>> d.next_toward(D('-inf'))
Decimal('123456.7890123456745641306042')

Make sure that decimal context has values that you need:

>>> from decimal import getcontext
>>> getcontext()
Context(prec=28, rounding=ROUND_HALF_EVEN, Emin=-999999999, Emax=999999999,
capitals=1, flags=[], traps=[InvalidOperation, DivisionByZero, Overflow])

Alternatives:

  • Call C99 nextafter() using ctypes:

      >>> import ctypes
      >>> nextafter = ctypes.CDLL(None).nextafter
      >>> nextafter.argtypes = ctypes.c_double, ctypes.c_double
      >>> nextafter.restype = ctypes.c_double
      >>> nextafter(4, float('+inf'))
      4.000000000000001
      >>> _.as_integer_ratio()
      (4503599627370497, 1125899906842624)
    

    Using numpy:

      >>> import numpy
      >>> numpy.nextafter(4, float('+inf'))
      4.0000000000000009
      >>> _.as_integer_ratio()
      (4503599627370497, 1125899906842624)
    

    Despite different repr(), the result is the same.

  • If we ignore edge cases then a simple frexp/ldexp solution from @S.Lott answer works:

      >>> import math, sys
      >>> m, e = math.frexp(4.0)
      >>> math.ldexp(2 * m + sys.float_info.epsilon, e - 1)
      4.000000000000001
      >>> _.as_integer_ratio()
      (4503599627370497, 1125899906842624)
    
  • pure Python next_after(x, y) implementation by @Mark Dickinson that takes into account edge cases. The result is the same in this case.

jfs
  • 399,953
  • 195
  • 994
  • 1,670
  • It looks good by need some polish - since float in python has 53 bits for mantissa http://docs.python.org/2/tutorial/floatingpoint.html. I think better is use frexp ldexp - can be faster. – Chameleon Nov 27 '13 at 15:10
  • @Chameleon: yes. You can use `frexp()`/`ldexpr()` to find the "next after". – jfs Jan 14 '17 at 14:58
1

I am not sure if this is what you want, but sys.float_info.epsilon is the "difference between 1 and the least value greater than 1 that is representable as a float", and you could do x * (1 + sys.float_info.epsilon).

http://docs.python.org/library/sys.html#sys.float_info

Artur Gaspar
  • 4,407
  • 1
  • 26
  • 28
  • 1
    It's a related idea, but not exactly what I want. Essentially it is correct for an exponent of 0, but for smaller exponents it is too large, and for larger it is too small. Thank you for pointing out `sys.float_info`. At least with this structure I can reliably determine the size of the underlying floating point type, even if I still have to guess that it's IEEE-754. With this I can work out a solution along the lines of @Benjamin's suggestions, specifically with the struct module. – rlibby Apr 22 '11 at 02:38
  • `x + epsilon`? wouldn't that be rather `x * (1 + epsilon)`? – vartec Apr 22 '11 at 14:34
  • @vartec I don't know, but I corrected it in the answer. – Artur Gaspar Apr 23 '11 at 23:37
  • @rlibby: to take into account the exponent: [`next_after_x = (m+sys.float_info.epsilon)*2**e` where `m, e = math.frexp(x)` from @S.Lott's answer](http://stackoverflow.com/a/6064066/4279) – jfs Jan 14 '17 at 14:21
  • or closer to `nextafter(x, float('+inf'))` value: `((2*m+sys.float_info.epsilon)*2**(e-1))` because `.5 <= m < 1`. – jfs Jan 14 '17 at 14:32