bpo-39310: Add math.ulp(x)

[vstinner]

Add math.ulp(): return the unit in the last place of x.

https://bugs.python.org/issue39310

[vstinner]

The documentation should be enhanced 😊

[vstinner]

The documentation should be enhanced blush

Done. I copied the test_math.ulp() docstring.

[tim-one]

I think Java does as well as can be done here:

https://www.geeksforgeeks.org/java-math-ulp-method-examples/

Short course:

• Return a NaN for a NaN input.
• Return +inf for an infinite input (regardless of sign).
• Return the smallest positive representable double for +0 or -0 input.

[vstinner]

@mdickinson, @tim-one, @serhiy-storchaka: I updated the PR. Would you mind to review it again?

I rewrote the documentation to clarify the behavior on corner cases: zero, min/max, inf, NaN.

[mdickinson]

One request for an explicit test for ulp(-0.0). Apart from that, this LGTM.

[serhiy-storchaka]

Agree with @mdickinson. The rest LGTM.

[vstinner]

Interesting case, math.ulp(x) == x is for one value, the smallest denormalized positive float:

>>> import math, sys
>>> x = math.ulp(0)
>>> math.ulp(x) == x
True
>>> x
5e-324

[vstinner]

I'm not sure that math.ulp(sys.float_info.max) is correct: it looks 2x larger than what we should get, since (max - math.ulp(max)/2) != max and (max + math.ulp(max)/2) != max.

vstinner@apu$ ./python
Python 3.9.0a2+ (heads/nextafter:6bbcab5fa2, Jan 11 2020, 16:24:52) 
>>> import math, sys
>>> max = sys.float_info.max
>>> inf = math.inf

# what do we have between max and inf?

>>> math.nextafter(max, inf)
inf
>>> math.nextafter(inf, -inf) == max
True
>>> inf - max
inf

# from max to first float before max (prev)

>>> prev=math.nextafter(max, -inf)
>>> (max - math.ulp(max)) == prev
True
>>> (max - math.ulp(max)/2) == prev  # oh!!!
True
>>> (max - math.ulp(max)/4) == prev
False
>>> (max - math.ulp(max)/4) == max
True

# from max to inf

>>> max + math.ulp(max)
inf
>>> max + math.ulp(max)/2  # oh!!!
inf
>>> max + math.ulp(max)/4
1.7976931348623157e+308
>>> max + math.ulp(max)/4 == max
True

# max ULP value
>>> math.ulp(max)
1.99584030953472e+292
>>> math.ulp(max)/2
9.9792015476736e+291

[mdickinson]

>>> (max - math.ulp(max)/2) == prev  # oh!!!
True

This is fine; it's just round-ties-to-even at work. max - math.ulp(max) / 2 is exactly halfway between two representable floats: max and prev. Of those two, prev has its last bit clear, so the tie is rounded to prev rather than max.

[vstinner]

One request for an explicit test for ulp(-0.0). Apart from that, this LGTM.

I wrote tests differently to better highlight that we expect ulp(-x) = ulp(x), and not copy/paste results from tests on positive numbers to tests on negative numbers:

        # negative number: ulp(-x) == ulp(x)
        for x in (0.0, 1.0, 2 ** 52, 2 ** 64, INF):
            with self.subTest(x=x):
                self.assertEqual(math.ulp(-x), math.ulp(x))

[vstinner]

Oh, nntplib test fail on the Documentation job of Travis CI:

Warning, treated as error:

**********************************************************************

File "library/nntplib.rst", line ?, in default

Failed example:

    s = NNTP('news.gmane.io')

Exception raised:

    Traceback (most recent call last):

      File "/home/travis/build/python/cpython/Lib/doctest.py", line 1329, in __run

        exec(compile(example.source, filename, "single",

      File "<doctest default[0]>", line 1, in <module>

        s = NNTP('news.gmane.io')

      File "/home/travis/build/python/cpython/Lib/nntplib.py", line 1049, in __init__

        super().__init__(file, host, readermode, timeout)

      File "/home/travis/build/python/cpython/Lib/nntplib.py", line 331, in __init__

        self.welcome = self._getresp()

      File "/home/travis/build/python/cpython/Lib/nntplib.py", line 456, in _getresp

        raise NNTPTemporaryError(resp)

    nntplib.NNTPTemporaryError: 400 load at 18.62, try later

I retry the job.

[vstinner]

This is fine; it's just round-ties-to-even at work. max - math.ulp(max) / 2 is exactly halfway between two representable floats: max and prev. Of those two, prev has its last bit clear, so the tie is rounded to prev rather than max.

Floating points are so complex :-(

[vstinner]

Oops, I forgot to update the documentation in the docstring and commit message: updated.

[vstinner]

Thanks a lot @tim-one, @serhiy-storchaka and @mdickinson: the review was required and really helpful here ;-)