ULP Difference of Float Numbers

The Unit in the Last Place (ULP) is the smallest difference between two adjacent floating-point numbers at a specific precision. It measures the granularity of floating-point representations and is essential for understanding rounding errors and numerical stability. For example, a ULP near 1.0 is typically smaller than a ULP near 10.0, as the step size depends on the exponent of the number.

ULP differences are often used to evaluate the accuracy of floating-point operations. When comparing two floating-point numbers, a difference of 1 ULP indicates that the numbers are as close as they can be while still being distinct. Smaller ULP differences imply better precision and numerical accuracy. When implementing floating-point operations or algorithms (e.g., the tanhf function), the smallest possible ULP difference is 0.5 indicating that the exact result is correctly rounded to the nearest floating-point number.

Calculating ULP

To calculate the ULP of a floating-point number, you can use the following C function:

#include <math.h>

float flt_ulp_size(float x) {
    x = fabs(x);

    return nextafterf(x, INFINITY) - x;
}

This function calculates the ULP size of a given floating-point number x. It first takes the absolute value of x using fabs(x), then uses nextafterf(x, INFINITY) to find the next representable float value after x in the direction of positive infinity. The difference between these two values gives the ULP size for x.

ULP Example Values 32-bit Float

Floating-Point Number	ULP (Macro)	ULP (Float)	ULP (Float Hex Representation)	Description
1.0	FLT_EPSILON	1.1920929e-07	0x1.0p-23	Smallest step near 1.0
0.0	FLT_TRUE_MIN	1.4012985e-45	0x1.0p-149	Smallest positive float
0.5	FLT_EPSILON/2	5.9604645e-08	0x1.0p-24	Half the step size near 1.0
10.0	FLT_EPSILON*8	9.536743e-07	0x1.0p-20	Larger step size near 10.0

Calculating ULP Difference

Simple ULP Difference

It is relative hard to correctly calculate the ULP difference between two floating-point numbers. The following C function is a simple way to calculate the ULP difference between a reference value and a given value:

float flt_ulp_diff(float ref, float val) {
    return fabs(ref - val) / flt_ulp_size(ref);
}

This function simple divides the absolute difference by the ULP. The problem with this approach is which ULP to use. You can use the ULP of the reference value, calculated or maximum/minimum/average or both values. However, whatever ULP you choose, it will not be perfect. E.g. in the provided function flt_ulp_diff(0.9999, 2.0) will be around factor 2 greater than flt_ulp_diff(1.0, 2.0) which is not very intuitive.

float flt_ulp_diff(float ref, float val) {
    return fabs(ref - val) / flt_ulp_size(ref);
}

Correct ULP Difference

To correctly calculate the ULP difference between two floating-point numbers, you can use the following C function:

uint32_t ulp_intdiff_float(float f1, float f2) {
    if (signbit(f1) != signbit(f2)) {
        return ulp_intdiff_float(0.0f, fabs(f1)) + ulp_intdiff_float(0.0f, fabs(f2)) + 1;
    }

    uint32_t i1 = *(uint32_t*)&f1;
    uint32_t i2 = *(uint32_t*)&f2;

    return i1 > i2 ? i1 - i2 : i2 - i1;
}

This function calculates the integer difference in ULP between two floating-point numbers f1 and f2. The trick is to reinterpret the floating-point numbers as 32-bit unsigned integers and then calculate the difference between them. The floating-point IEEE 754 representation is designed in such a way that the integer difference between two floating-point numbers (with the same sign) is the ULP difference between them.

To also consider floating-point numbers with different signs, the function calculates the ULP difference towards zero for both numbers and adds the differences together. Additionally, it adds 1 to the result that ULP difference between negative and positive zero (ulp_intdiff_float(-0.0f, 0.0f)) is 1. Adding 1 to the result is not necessary if you do not care about the difference between negative and positive zero.

This implementation satisfies both the properties of triangular additivity (ulp_intdiff_float(x, y) + ulp_intdiff_float(y, z) == ulp_intdiff_float(x, z) for x < y < z) and commutativity (ulp_intdiff_float(x, y) == ulp_intdiff_float(y, x)) of ULP differences, ensuring that the ULP difference between any two floating-point numbers is consistent and order-independent when summed across intermediates.

Examples

Example Results for `ulp_intdiff_float` Function

Input 1	Input 2	ULP Difference	Explanation
1.0	1.0 + FLT_EPSILON	1	Adjacent floating-point values
1.0	1.0 - FLT_EPSILON	2	One step away in the opposite direction
FLT_TRUE_MIN	0.0	1	Smallest positive float to zero
-FLT_TRUE_MIN	FLT_TRUE_MIN	3	Crossing zero, including sign change
-1.0	-1.0 - FLT_EPSILON	1	Adjacent floating-point values on the negative side
-1.0	-1.0 + FLT_EPSILON	2	One step away in the opposite direction on the negative side
-0.0	0.0	1	Difference between signed zero values
0.0	0.0	0	No difference between identical values
FLT_MAX	INFINITY	1	Transition from the largest finite float to infinity
-FLT_MAX	-INFINITY	1	Transition from the largest negative finite float to negative infinity
1.0	0.5	8388608	Number of representable floats between 1.0 and 0.5 (2^23)
1.0	2.0	8388608	Number of representable floats between 1.0 and 2.0 (2^23)

Calculating Rational ULP Differences

When comparing a reference value provided as a double and a test value as a float, rational ULP differences can provide a finer granularity of comparison. The function below implements this computation by combining integer ULP differences with a fractional component derived from the difference between the double reference value and its float representation.

The following function calculates rational ULP differences:

double ulp_diff(double ref, float b) {
    // Handle NaN cases: returns 0 if both are NaN, INFINITY otherwise
    if (isnan(b) || isnan(ref)) {
        return isnan(b) == isnan(ref) ? 0.0 : INFINITY;
    }

    // Convert reference value to float for integer ULP difference calculation
    float fref = (float)ref;

    // Calculate integer ULP difference between the float reference and test value
    double iulp = (double) ulp_intdiff_float(fref, b);

    // Ensure positive values for further calculations
    ref = fabs(ref);
    fref = fabsf(fref);
    b = fabsf(b);

    // Calculate the ULP size of the rounded-down reference value
    float ulpRef = flt_ulp_size(double_to_float_rounddown(ref));

    // Calculate the fractional difference in ULPs
    double diff = (ref - (double)fref);
    double fulp = diff / ulpRef;
    if (b > ref)
        fulp = -fulp;

    // Return the sum of integer ULPs and fractional ULPs
    return iulp + fulp;
}

The function first handles special cases like NaN values, then calculates the integer ULP difference between the float-converted reference value and the test value. It then computes the fractional ULP difference based on the difference between the original double reference value and its float representation. We can use the method from the simple ULP difference to calculate the fractional ULP difference. That is working here because the ULP size for the fractional part is clearly defined by the location of the reference value.

One problem is the rounding when converting the double to float. E.g. 1.0 - FLT_EPSILON/4 require the ULP size of the range [0.5, 1.0) but would be converted to 1.0 in a normal double to float conversion. This is why the function uses double_to_float_rounddown to convert the double to float. This function converts the double to float and rounds down to the next smaller float value:

float double_to_float_rounddown(double value) {
    float float_value = (float)value;

    if (float_value > value) {
        return nextafterf(float_value, -INFINITY);
    }

    return float_value;
}

Examples

The following table provides examples of the results obtained using the ulp_diff function:

Example Results for ulp_diff Function

Reference Value	Test Value	ULP Difference	Explanation
1.0 + FLT_EPSILON/2	1.0	0.5	The test value is halfway between two floats.
1.0 + FLT_EPSILON/4	1.0	0.25	The test value is a quarter ULP away.
1.0 - FLT_EPSILON/2	1.0	1.0	The test value is 1 ULP lower.
1.0 - FLT_EPSILON/4	1.0	0.5	The test value is halfway to the next float.
1.0 + FLT_EPSILON/2	1.0 + FLT_EPSILON*10	9.5	Integer difference with fractional part.
1.0 + FLT_EPSILON/4	1.0 + FLT_EPSILON*10	9.75	Includes both integer and fractional differences.
1.0 - FLT_EPSILON/2	1.0 - FLT_EPSILON*10	19.0	The difference spans multiple ULPs.
1.0 - FLT_EPSILON/4	1.0 - FLT_EPSILON*10	19.5	Combined fractional and integer ULP difference.
0.0 + FLT_TRUE_MIN/100	0.0	0.01	A very small fractional ULP difference.