traverseda
/
cpython
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

bpo-39310: Add math.ulp(x) (GH-17965) 

Add math.ulp(): return the value of the least significant bit
of a float.
This commit is contained in:
Victor Stinner 2020-01-13 12:44:35 +01:00 committed by GitHub
parent 7ba6f18de2
commit 0b2ab21956
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
7 changed files with 144 additions and 37 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
Doc/library/math.rst
View File

@ -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

22
Doc/library/sys.rst
View File

@ -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 |

4
Doc/whatsnew/3.9.rst
View File

@ -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
-------

55
Lib/test/test_math.py
View File

@ -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

1
Misc/NEWS.d/next/Library/2020-01-12-13-34-42.bpo-39310.YMRdcj.rst Normal file
View File

@ -0,0 +1 @@
Add :func:`math.ulp`: return the value of the least significant bit of a float.

41
Modules/clinic/mathmodule.c.h generated
View File

@ -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]*/

32
Modules/mathmodule.c
View File

@ -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 */
};