From 0b2ab21956fbab8eab6d064060d4544499730316 Mon Sep 17 00:00:00 2001 From: Victor Stinner Date: Mon, 13 Jan 2020 12:44:35 +0100 Subject: bpo-39310: Add math.ulp(x) (GH-17965) Add math.ulp(): return the value of the least significant bit of a float. --- Doc/library/math.rst | 26 ++++++++++ Doc/library/sys.rst | 22 +++++---- Doc/whatsnew/3.9.rst | 4 ++ Lib/test/test_math.py | 55 +++++++++++----------- .../2020-01-12-13-34-42.bpo-39310.YMRdcj.rst | 1 + Modules/clinic/mathmodule.c.h | 41 +++++++++++++++- Modules/mathmodule.c | 32 +++++++++++++ 7 files changed, 144 insertions(+), 37 deletions(-) create mode 100644 Misc/NEWS.d/next/Library/2020-01-12-13-34-42.bpo-39310.YMRdcj.rst diff --git a/Doc/library/math.rst b/Doc/library/math.rst index c9f2a38..c4c1800 100644 --- a/Doc/library/math.rst +++ b/Doc/library/math.rst @@ -226,6 +226,8 @@ Number-theoretic and representation functions * ``math.nextafter(x, 0.0)`` goes towards zero. * ``math.nextafter(x, math.copysign(math.inf, x))`` goes away from zero. + See also :func:`math.ulp`. + .. versionadded:: 3.9 .. function:: perm(n, k=None) @@ -284,6 +286,30 @@ Number-theoretic and representation functions :class:`~numbers.Integral` (usually an integer). Delegates to :meth:`x.__trunc__() `. +.. function:: ulp(x) + + Return the value of the least significant bit of the float *x*: + + * If *x* is a NaN (not a number), return *x*. + * If *x* is negative, return ``ulp(-x)``. + * If *x* is a positive infinity, return *x*. + * If *x* is equal to zero, return the smallest positive + *denormalized* representable float (smaller than the minimum positive + *normalized* float, :data:`sys.float_info.min `). + * If *x* is equal to the largest positive representable float, + return the value of the least significant bit of *x*, such that the first + float smaller than *x* is ``x - ulp(x)``. + * Otherwise (*x* is a positive finite number), return the value of the least + significant bit of *x*, such that the first float bigger than *x* + is ``x + ulp(x)``. + + ULP stands for "Unit in the Last Place". + + See also :func:`math.nextafter` and :data:`sys.float_info.epsilon + `. + + .. versionadded:: 3.9 + Note that :func:`frexp` and :func:`modf` have a different call/return pattern than their C equivalents: they take a single argument and return a pair of diff --git a/Doc/library/sys.rst b/Doc/library/sys.rst index 0aae263..351a8e4 100644 --- a/Doc/library/sys.rst +++ b/Doc/library/sys.rst @@ -479,8 +479,10 @@ always available. +---------------------+----------------+--------------------------------------------------+ | attribute | float.h macro | explanation | +=====================+================+==================================================+ - | :const:`epsilon` | DBL_EPSILON | difference between 1 and the least value greater | - | | | than 1 that is representable as a float | + | :const:`epsilon` | DBL_EPSILON | difference between 1.0 and the least value | + | | | greater than 1.0 that is representable as a float| + | | | | + | | | See also :func:`math.ulp`. | +---------------------+----------------+--------------------------------------------------+ | :const:`dig` | DBL_DIG | maximum number of decimal digits that can be | | | | faithfully represented in a float; see below | @@ -488,20 +490,24 @@ always available. | :const:`mant_dig` | DBL_MANT_DIG | float precision: the number of base-``radix`` | | | | digits in the significand of a float | +---------------------+----------------+--------------------------------------------------+ - | :const:`max` | DBL_MAX | maximum representable finite float | + | :const:`max` | DBL_MAX | maximum representable positive finite float | +---------------------+----------------+--------------------------------------------------+ - | :const:`max_exp` | DBL_MAX_EXP | maximum integer e such that ``radix**(e-1)`` is | + | :const:`max_exp` | DBL_MAX_EXP | maximum integer *e* such that ``radix**(e-1)`` is| | | | a representable finite float | +---------------------+----------------+--------------------------------------------------+ - | :const:`max_10_exp` | DBL_MAX_10_EXP | maximum integer e such that ``10**e`` is in the | + | :const:`max_10_exp` | DBL_MAX_10_EXP | maximum integer *e* such that ``10**e`` is in the| | | | range of representable finite floats | +---------------------+----------------+--------------------------------------------------+ - | :const:`min` | DBL_MIN | minimum positive normalized float | + | :const:`min` | DBL_MIN | minimum representable positive *normalized* float| + | | | | + | | | Use :func:`math.ulp(0.0) ` to get the | + | | | smallest positive *denormalized* representable | + | | | float. | +---------------------+----------------+--------------------------------------------------+ - | :const:`min_exp` | DBL_MIN_EXP | minimum integer e such that ``radix**(e-1)`` is | + | :const:`min_exp` | DBL_MIN_EXP | minimum integer *e* such that ``radix**(e-1)`` is| | | | a normalized float | +---------------------+----------------+--------------------------------------------------+ - | :const:`min_10_exp` | DBL_MIN_10_EXP | minimum integer e such that ``10**e`` is a | + | :const:`min_10_exp` | DBL_MIN_10_EXP | minimum integer *e* such that ``10**e`` is a | | | | normalized float | +---------------------+----------------+--------------------------------------------------+ | :const:`radix` | FLT_RADIX | radix of exponent representation | diff --git a/Doc/whatsnew/3.9.rst b/Doc/whatsnew/3.9.rst index a686d64..340079c 100644 --- a/Doc/whatsnew/3.9.rst +++ b/Doc/whatsnew/3.9.rst @@ -184,6 +184,10 @@ Add :func:`math.nextafter`: return the next floating-point value after *x* towards *y*. (Contributed by Victor Stinner in :issue:`39288`.) +Add :func:`math.ulp`: return the value of the least significant bit +of a float. +(Contributed by Victor Stinner in :issue:`39310`.) + nntplib ------- diff --git a/Lib/test/test_math.py b/Lib/test/test_math.py index b64fd41..6d10227 100644 --- a/Lib/test/test_math.py +++ b/Lib/test/test_math.py @@ -53,30 +53,6 @@ def to_ulps(x): return n -def ulp(x): - """Return the value of the least significant bit of a - float x, such that the first float bigger than x is x+ulp(x). - Then, given an expected result x and a tolerance of n ulps, - the result y should be such that abs(y-x) <= n * ulp(x). - The results from this function will only make sense on platforms - where native doubles are represented in IEEE 754 binary64 format. - """ - x = abs(float(x)) - if math.isnan(x) or math.isinf(x): - return x - - # Find next float up from x. - n = struct.unpack(' double + + x: double + / + +Return the value of the least significant bit of the float x. +[clinic start generated code]*/ + +static double +math_ulp_impl(PyObject *module, double x) +/*[clinic end generated code: output=f5207867a9384dd4 input=31f9bfbbe373fcaa]*/ +{ + if (Py_IS_NAN(x)) { + return x; + } + x = fabs(x); + if (Py_IS_INFINITY(x)) { + return x; + } + double inf = m_inf(); + double x2 = nextafter(x, inf); + if (Py_IS_INFINITY(x2)) { + /* special case: x is the largest positive representable float */ + x2 = nextafter(x, -inf); + return x - x2; + } + return x2 - x; +} + + static PyMethodDef math_methods[] = { {"acos", math_acos, METH_O, math_acos_doc}, {"acosh", math_acosh, METH_O, math_acosh_doc}, @@ -3366,6 +3397,7 @@ static PyMethodDef math_methods[] = { MATH_PERM_METHODDEF MATH_COMB_METHODDEF MATH_NEXTAFTER_METHODDEF + MATH_ULP_METHODDEF {NULL, NULL} /* sentinel */ }; -- cgit v0.12