Is the function math.ulp more precise than the explicit formula in Python?

In floating-point arithmetic, the unit in the last place (ULP) of a floating-point number is the spacing between that number and the consecutive one, i.e. the value of its least significant digit (rightmost digit) if it is 1. It is given by this formula:

ULP(x) = b^−(p−1) • |x|

where b is the base (2 for binary numerals) and p (53 for double-precision significands) is the precision.

Python 3.9 introduced a new function math.ulp to compute the ULP of a floating-point number.

With this function, the previous formula is verified as expected for the ULP of 1:

>>> math.ulp(1)
2.220446049250313e-16
>>> 2**(-(53 - 1)) * abs(1)
2.220446049250313e-16

but it is not verified for the ULP of 10⁻¹⁰ for instance:

>>> math.ulp(1e-10)
1.2924697071141057e-26
>>> 2**(-(53 - 1)) * abs(1e-10)
2.2204460492503132e-26

Is math.ulp(x) more precise than 2**(-(53 - 1)) * abs(x)? Why?

The CPython implementation is in Modules/mathmodule.c#L3408-L3427 but I cannot find the implementation of the called function nextafter to understand:

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;
}

2^−(53−1) • |x| (or 2**(-(53 - 1)) * abs(x)) is not a formula for ULP(x) (or math.ulp(x)) because it does not give the value of a 1 in the position of the lowest bit of x but rather the value of the significand of x (1.something) scaled to the position of the lowest bit of x. When x is not a power of two, its significand exceeds 1, and the formula is too high.

The correct formula is 2^−(53−1) • 2^{max(e, −1022)} where e is the IEEE 754 normalized exponent of x, i.e. 2^e ≤ |x| < 2^e+1 (or 2**(-(53 - 1)) * 2**max(math.floor(math.log2(x)), -1022)).