Calculating ULP distance between two floating point numbers (quickly)

General Usage
1

I need to compare floating point numbers for near-equality. I would like to compare them based on the ULP distance (see here if you aren’t familiar/need a refresher).

The most straightforward way to do this, as far as I know, is to convert the bits that make up the floating point numbers into their equivalent integers, and then take the difference of the two resulting integers as the ULP distance.

The very naive solution I have come up with is to do something like this in Julia is:

function ULPsBetween(a::Float64, b::Float64)
    a = bitstring(a)
    b = bitstring(b)
    
    aintegerrepresentation = parse(Int, a; base=2)
    bintegerrepresentation = parse(Int, b; base=2)
    
    return aintegerrepresentation - bintegerrepresentation
end

This code works correctly, but it’s very slow since it calls bitstring() and parse(). I used BenchmarkTools to evaluate the performance, and its on the order of 1 microsecond. For the number of times this function will be called in my program, that is very slow.

In C, we can do something like this instead of the parse() function:

int64_t floatToInt(double f)
{
    int64_t r;
    memcpy(&r, &f, sizeof(double));
    return r;
}

This code executes on the order of a few nanoseconds.

Obviously, I’m not using Julia correctly. But what would a more elegant/faster solution look like? Any ideas?

For now, I will probably write a Julia function that calls the C code above with ccall(). But I would really prefer to have everything in pure Julia. Finally, I searched quite a bit and couldn’t find anything about calculating the ULP distance in Julia on here or StackOverflow, so this post may be useful to others who find themselves in the same situation in the future.

Thank you!

2

Either use reinterpret to convert your floats to integers and extract the bit patterns from that or rescale the output of significand or frexp.

3

Is using isapprox an option? You can give it an absolute or relative allowed difference.

In case it isn’t, the code from your link would look like this in julia:

function ulpsBetween(a::Float64, b::Float64)
    a == b && return 0
    isnan(a) || isnan(b) && return typemax(Int)
    isinf(a) || isinf(b) && return typemax(Int)
    
    a_int = reinterpret(Int64, a)
    b_int = reinterpret(Int64, b)

    (a_int < 0) != (b_int < 0) && return typemax(Int)

    return abs(a_int - b_int)
end

Still, the link discusses precisely why this sort of comparison usually isn’t a good idea (and isapprox is usually the best choice anyway). It’s implemented like this:

function isapprox(x::Number, y::Number;                                                                                          
                  atol::Real=0, rtol::Real=rtoldefault(x,y,atol),                                                                
                  nans::Bool=false, norm::Function=abs)                                                                          
    x == y || (isfinite(x) && isfinite(y) && norm(x-y) <= max(atol, rtol*max(norm(x), norm(y)))) || (nans && isnan(x) && isnan(y)
end

which you can find via @edit isapprox(0.5, 1.0) at base/floatfuncs.jl:360.

4

Like @GunnarFarneback suggested, something like this?

ulpdistance(a, b) = -(reinterpret.(Int, (b, a))...)

Could make it a bit less terse, maybe… Runtime is a couple of ns.

Hm, I see that @Sukera wrote something similar, it looked totally different at first glance.

5

That approach ignores special handling NaN and Inf though.

6

Yeah, yeah, boring :wink:

7

Wow, I wasn’t expecting such excellent answers so quickly! What a great community. I’m glad I adopted Julia for this project.

reinterpret() was exactly what I needed.

I will look into using isapprox() instead of the ULP method.

8

If what you actually want is approximate comparison of floats, then that’s exactly what isapprox is for.

9

It seems to me that a more reliable way than comparing integer representations (which won’t work if the two numbers are more than a factor of two apart) would be just to use (b - a) / eps(a):

julia> a, b = 1.1, nextfloat(1.1)
(1.1, 1.1000000000000003)

julia> (b - a) / eps(a)
1.0