diff options
-rw-r--r-- | Doc/library/math.rst | 26 | ||||
-rw-r--r-- | Doc/library/sys.rst | 22 | ||||
-rw-r--r-- | Doc/whatsnew/3.9.rst | 4 | ||||
-rw-r--r-- | Lib/test/test_math.py | 55 | ||||
-rw-r--r-- | Misc/NEWS.d/next/Library/2020-01-12-13-34-42.bpo-39310.YMRdcj.rst | 1 | ||||
-rw-r--r-- | Modules/clinic/mathmodule.c.h | 41 | ||||
-rw-r--r-- | Modules/mathmodule.c | 32 |
7 files changed, 144 insertions, 37 deletions
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__() <object.__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 <sys.float_info>`). + * 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 + <sys.float_info>`. + + .. 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) <math.ulp>` 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('<q', struct.pack('<d', x))[0] - x_next = struct.unpack('<d', struct.pack('<q', n + 1))[0] - if math.isinf(x_next): - # Corner case: x was the largest finite float. Then it's - # not an exact power of two, so we can take the difference - # between x and the previous float. - x_prev = struct.unpack('<d', struct.pack('<q', n - 1))[0] - return x - x_prev - else: - return x_next - x - # Here's a pure Python version of the math.factorial algorithm, for # documentation and comparison purposes. # @@ -470,9 +446,9 @@ class MathTests(unittest.TestCase): def testCos(self): self.assertRaises(TypeError, math.cos) - self.ftest('cos(-pi/2)', math.cos(-math.pi/2), 0, abs_tol=ulp(1)) + self.ftest('cos(-pi/2)', math.cos(-math.pi/2), 0, abs_tol=math.ulp(1)) self.ftest('cos(0)', math.cos(0), 1) - self.ftest('cos(pi/2)', math.cos(math.pi/2), 0, abs_tol=ulp(1)) + self.ftest('cos(pi/2)', math.cos(math.pi/2), 0, abs_tol=math.ulp(1)) self.ftest('cos(pi)', math.cos(math.pi), -1) try: self.assertTrue(math.isnan(math.cos(INF))) @@ -1445,7 +1421,7 @@ class MathTests(unittest.TestCase): self.assertRaises(TypeError, math.tanh) self.ftest('tanh(0)', math.tanh(0), 0) self.ftest('tanh(1)+tanh(-1)', math.tanh(1)+math.tanh(-1), 0, - abs_tol=ulp(1)) + abs_tol=math.ulp(1)) self.ftest('tanh(inf)', math.tanh(INF), 1) self.ftest('tanh(-inf)', math.tanh(NINF), -1) self.assertTrue(math.isnan(math.tanh(NAN))) @@ -2036,7 +2012,7 @@ class IsCloseTests(unittest.TestCase): def assertEqualSign(self, x, y): """Similar to assertEqual(), but compare also the sign. - Function useful to check to signed zero. + Function useful to compare signed zeros. """ self.assertEqual(x, y) self.assertEqual(math.copysign(1.0, x), math.copysign(1.0, y)) @@ -2087,6 +2063,29 @@ class IsCloseTests(unittest.TestCase): self.assertTrue(math.isnan(math.nextafter(1.0, NAN))) self.assertTrue(math.isnan(math.nextafter(NAN, NAN))) + @requires_IEEE_754 + def test_ulp(self): + self.assertEqual(math.ulp(1.0), sys.float_info.epsilon) + # use int ** int rather than float ** int to not rely on pow() accuracy + self.assertEqual(math.ulp(2 ** 52), 1.0) + self.assertEqual(math.ulp(2 ** 53), 2.0) + self.assertEqual(math.ulp(2 ** 64), 4096.0) + + # min and max + self.assertEqual(math.ulp(0.0), + sys.float_info.min * sys.float_info.epsilon) + self.assertEqual(math.ulp(FLOAT_MAX), + FLOAT_MAX - math.nextafter(FLOAT_MAX, -INF)) + + # special cases + self.assertEqual(math.ulp(INF), INF) + self.assertTrue(math.isnan(math.ulp(math.nan))) + + # 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)) + def test_main(): from doctest import DocFileSuite diff --git a/Misc/NEWS.d/next/Library/2020-01-12-13-34-42.bpo-39310.YMRdcj.rst b/Misc/NEWS.d/next/Library/2020-01-12-13-34-42.bpo-39310.YMRdcj.rst new file mode 100644 index 0000000..a787f69 --- /dev/null +++ b/Misc/NEWS.d/next/Library/2020-01-12-13-34-42.bpo-39310.YMRdcj.rst @@ -0,0 +1 @@ +Add :func:`math.ulp`: return the value of the least significant bit of a float. diff --git a/Modules/clinic/mathmodule.c.h b/Modules/clinic/mathmodule.c.h index f34633c..f95d291 100644 --- a/Modules/clinic/mathmodule.c.h +++ b/Modules/clinic/mathmodule.c.h @@ -856,4 +856,43 @@ math_nextafter(PyObject *module, PyObject *const *args, Py_ssize_t nargs) exit: return return_value; } -/*[clinic end generated code: output=e4ed1a800e4b2eae input=a9049054013a1b77]*/ + +PyDoc_STRVAR(math_ulp__doc__, +"ulp($module, x, /)\n" +"--\n" +"\n" +"Return the value of the least significant bit of the float x."); + +#define MATH_ULP_METHODDEF \ + {"ulp", (PyCFunction)math_ulp, METH_O, math_ulp__doc__}, + +static double +math_ulp_impl(PyObject *module, double x); + +static PyObject * +math_ulp(PyObject *module, PyObject *arg) +{ + PyObject *return_value = NULL; + double x; + double _return_value; + + if (PyFloat_CheckExact(arg)) { + x = PyFloat_AS_DOUBLE(arg); + } + else + { + x = PyFloat_AsDouble(arg); + if (x == -1.0 && PyErr_Occurred()) { + goto exit; + } + } + _return_value = math_ulp_impl(module, x); + if ((_return_value == -1.0) && PyErr_Occurred()) { + goto exit; + } + return_value = PyFloat_FromDouble(_return_value); + +exit: + return return_value; +} +/*[clinic end generated code: output=9b51d215dbcac060 input=a9049054013a1b77]*/ diff --git a/Modules/mathmodule.c b/Modules/mathmodule.c index 632a421..5e8e485 100644 --- a/Modules/mathmodule.c +++ b/Modules/mathmodule.c @@ -3314,6 +3314,37 @@ math_nextafter_impl(PyObject *module, double x, double y) } +/*[clinic input] +math.ulp -> 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 */ }; |