2 files changed, 72 insertions, 1 deletions
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 */
 };