Units in the last place (ULP)

What is ULP?

Units in the last place (ULP) is the distance between two adjacent floating-point numbers at a given value. It represents the smallest possible change in a number’s representation at that magnitude.

ULP is a function of the number’s exponent and can be calculated with the following expression:

    

        
        ULP(x) = nextafter(x, +inf) - x

        
    

ULP example: binary model

Building on the example from the previous section:

Fixed B = 2, p = 3, e^max = 2, e^min = -1

Based on the above definition, you can compute the ULP values for the numbers in this set as follows:

    

        
        
ULP(0.625) = nextafter(0.625, +inf) - 0.625 = 0.75-0.625 = 0.125

    

    

        
        
ULP(4.0) = 1.0

    

As the exponent of x increases, ULP(x) increases exponentially. That is, the spacing between floating-point values grows with the magnitude of x.

Numbers with the same exponent have the same ULP.

ULP in IEEE-754

For normalized IEEE-754 floating-point numbers, a similar behavior is observed: the distance between two adjacent representable values — that is, ULP(x) — is a power of two that depends only on the exponent of x.

Optimized expression

A faster, commonly used expression for ULP is:

    

        
        
ULP(x) = 2^(e-p+1)

    

Where:

• e is the unbiased exponent (in the IEEE-754 definition of single precision this is E-127)
• p is the precision (23 for IEEE-754 single-precision)

When computing the ULP of IEEE-754 floats, this expression becomes:

    

        
        
ULP(x) = 2^(e-23)

    

This expression is often used in mathematical computations of ULP since it offers performance benefits.

    

        
        
min_pos_denormal = 2⁻²³ × 2⁻¹²⁶ = 2⁻¹⁴⁹

    

    

        
        
second_min_pos_denormal = 2⁻²² × 2⁻¹²⁶ = 2⁻¹⁴⁸ = 2 × 2⁻¹⁴⁹

    

ULP implementation in C

Below is an example of an implementation of the ULP function of a number.

Use a text editor to save the code below in a file named ulp.h.

    

        
        
#include <stdint.h>
#include <string.h>
#include <math.h>

// Bit cast float to uint32_t
static inline uint32_t asuint(float x) {
    uint32_t u;
    memcpy(&u, &x, sizeof(u));
    return u;
}

// Compute exponent of ULP spacing at x
static inline int ulpscale(float x) {
    // Recover the biased exponent E
    int e = asuint(x) >> 23 & 0xff;
    if (e == 0)
        e++;  // handle subnormals

    // Compute the ULP exponent: e - p = E - 127 - 23
    return e - 127 - 23;
}

// Compute ULP spacing at x using ulpscale and scalbnf
static float ulp(float x) {
    return scalbnf(1.0f, ulpscale(x));
}

    

There are three key functions in this implementation:

• the asuint(x) function reinterprets the bit pattern of a float as a 32-bit unsigned integer, allowing the extraction of specific bit fields such as the exponent.
• the ulpscale(x) function returns the base-2 exponent of the ULP spacing at a given float value x, which is the result of log2(ULP(x)). The e variable in this function corresponds to the quantity E previously mentioned (the bitwise value of the exponent).
• the scalbnf(m, n) function (a standard function declared in math.h) efficiently evaluates m × 2^n.

Here’s an example program that calls ulp() to compute the spacing near a float value.

Use a text editor to save the code below in a file named ulp.c.

    

        
        
#include <stdio.h>
#include "ulp.h"

int main() {
    float x = 1.00000001f;
    float spacing = ulp(x);

    printf("ULP of %.8f is %.a\n", x, spacing);
    return 0;
}

    

Compile the program with GCC.

    

        
        
gcc -O2 ulp.c -o ulp

    

Run the program:

    

        
        
./ulp

    

On most systems, the output will print:

    

        
        ULP of 1.00000000 is 0x1p-23

        
    

This is the correct ULP spacing for values near 1.0f in IEEE-754 single-precision format.