Floating-point accuracy visualization example: the ULP error for each math function coming with Julia

nsajko · February 1, 2026, 9:02pm

Introduction

The implementation of a math function, such as cos(::Float64) in Julia, should be fast, but also needs to be at least somewhat accurate. Here I propose a method for visualizing empirical accuracy data.

Consider the discrepancy between the exact value of a univariate real function, such as the sine, and the approximate value that is returned by a call such as sin(0.3). By putting an interval that subsets the domain of the function on the x-axis of a plot, and the value of this discrepancy, the error, onto the y-axis of the same plot, it might be possible to tell on which part of the function’s domain is the implementation less accurate than on other parts. That is, it might be possible to read off the plot the worst-case regions for the approximation.

I choose the error measured in units in the last place/units of least precision (ULPs) as the type of error to visualize here. Other types of error that are often used are absolute error, which is not usually appropriate for floating-point numbers, and relative error, which has these drawbacks to measuring the error in ULPs:

it is usually less immediately useful and less intuitive
plotting software tends to have trouble with the tiny values

One advantage of the error in ULPs is that it has some convenient and easy interpretations in the context of floating-point numbers. For example:

If the error between the exact value and the approximate (floating-point) value is less than half, the approximate value is the number closest to the exact value among the numbers belonging to that floating-point representation. The technical term is correctly rounded.
If the error is merely less than one, the approximate value is one of the two floating-point numbers closest to the exact value. The technical term is faithfully rounded (although faithful rounding is not technically a rounding).

If we tried to, say, evaluate the error in ULPs on an evenly-spaced grid on the chosen interval, the plot would just look like noise if each evaluation was plotted as a data point. However, given that the worst cases are what is most interesting, it might be possible to achieve a comprehensible visualization by giving up on visualizing the best cases. In other words, we’re interested in the upward spikes of the error, and by eliding the downward spikes it might be possible to make the plot a lot less noisy.

To accomplish this, the approach I choose here is vaguely similar to downsampling/decimation, from signal analysis: take n values, where n is some large positive integer, and represent them on the plot by an aggregate value: the maximum value among the n points (in this case). In a signal analysis context, it is also common to apply a lowpass filter to reduce noise/high-frequency components of the signal, before the decimation itself. This helps prevent artifacts in the resulting output. In this case, a sliding window maximum seems like an appropriate lowpass filter to smooth out the data before decimation, to prevent artifacts.

Julia app on Github

The app used for creating these visualizations is published on Github:

GitHub - JuliaMath/VisualizeULPError.jl: Visualize the numerical error of a real univariate function

Neither the package, nor the app it exposes, are registered as of writing this.

The plots

Might be useful to both contributors and users of Julia to better understand where there’s room for improvement in the current implementations. In the cases of some functions, though, the worst case error spikes are difficult or impossible to fix efficiently, without reaching for some arbitrary-precision implementation like MPFR (BigFloat) or Arb.

`acos(::Float64)`

`acosd(::Float64)`

`acosh(::Float64)`

`acot(::Float64)`

`acotd(::Float64)`

`acoth(::Float64)`

`acsc(::Float64)`

`acscd(::Float64)`

`acsch(::Float64)`

`asec(::Float64)`

`asecd(::Float64)`

`asech(::Float64)`

`asin(::Float64)`

`asind(::Float64)`

`asinh(::Float64)`

`atan(::Float64)`

`atand(::Float64)`

`atanh(::Float64)`

`cos(::Float64)`

`cosc(::Float64)`

`cosd(::Float64)`

`cosh(::Float64)`

`cospi(::Float64)`

`cot(::Float64)`

`cotd(::Float64)`

`coth(::Float64)`

`csc(::Float64)`

`cscd(::Float64)`

`csch(::Float64)`

`deg2rad(::Float64)`

`exp(::Float64)`

`exp10(::Float64)`

`exp2(::Float64)`

`expm1(::Float64)`

`log(::Float64)`

`log10(::Float64)`

`log1p(::Float64)`

`log2(::Float64)`

`rad2deg(::Float64)`

`sec(::Float64)`

`secd(::Float64)`

`sech(::Float64)`

`sin(::Float64)`

`sinc(::Float64)`

`sind(::Float64)`

`sinh(::Float64)`

`sinpi(::Float64)`

`tan(::Float64)`

`tand(::Float64)`

`tanh(::Float64)`

`tanpi(::Float64)`

Miscellaneous

PRs featuring this visualization approach

Version, platform info

julia> versioninfo()
Julia Version 1.14.0-DEV.1670
Commit 6e64d0c3442 (2026-02-02 03:21 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × AMD Ryzen 3 5300U with Radeon Graphics
  WORD_SIZE: 64
  LLVM: libLLVM-20.1.8 (ORCJIT, znver2)
  GC: Built with stock GC
Threads: 5 default, 1 interactive, 5 GC (on 8 virtual cores)

jebej · February 1, 2026, 11:54pm

Nice! But the plots are hard to see due to the low contrast, low resolution and small text, if it’s not too much trouble, I’m sure an svg (or higher res) version with better sizing and colors would be appreciated!

OlivierHnt · February 2, 2026, 1:36am

Interesting, thank you for sharing.
We recently had a conversation about ULPs over at IntervalArithmetic.
We used to have a rounding mode based on prevfloat/nextfloat (fast, not tight, but GPU compatible), though we ultimately pulled the plug because determining the ULPs turned out to be tricky and architecture-dependent.

I am not an export on this subject, though Julia itself provides implementations of some functions:

Where can we find the list of implemented functions?
For these functions, this means that the ULP is architecture independent, correct?
Do any of these implementations provide proven bounds on the ULP error?
I would not expect correct rounding in general, but even a pessimistic, certified upper bound on the ULP would already be very useful for our purposes.
Maybe one day we can even get an implementation of CORE-Math functions.

PatrickHaecker · February 3, 2026, 4:56am

This sounds interesting. Am I the only one who can’t see the plots? I only see an “X” below each function name and I thought “X” never, ever marks the spot.

nsajko · February 3, 2026, 7:30am

Yeah, the post is ~~not done yet~~

For context, this was originally a blog post on the Julia Forem. However Forem does not accept SVG, so I provided the plots as PNG. However, Forem downscaled all the images, which is something I only realized after publishing the Forem post and linking to it from here. Since then I got sidetracked by a possible Julia bug, so I haven’t gotten around to regenerating the plots.

EDIT: the post is finished now, excuse the confusion.

nsajko · February 3, 2026, 4:39pm

The images are now in SVG. Hope it is OK now. Could also have titled the plots and labeled the axes, now that I think about it

nsajko · February 3, 2026, 4:52pm

Can you link some relevant discussion? I’m not completely sure what you are referring to.

A possible explanation is how the result of muladd(::Float64, ::Float64, ::Float64) is architecture-dependent, because modern platforms have a fused multiply-add (FMA) instruction, so muladd may be implemented with fma, but on older architectures muladd may suffer from double rounding.

The docs: Mathematics · The Julia Language

Are you asking about the approximation error? If so, yeah, as mentioned just now above. An arch that lacks FMA may suffer from worse accuracy, but modern platforms do have an FMA. The Julia test suite adapts like this:

github.com/JuliaLang/julia

test/math.jl

8f91f510c


      
          # has_fma has no runtime support.
          # So we need function wrappers to make this work.
          has_fma_Int() = Core.Compiler.have_fma(Int)
          has_fma_Float32() = Core.Compiler.have_fma(Float32)
          has_fma_Float64() = Core.Compiler.have_fma(Float64)
          
          has_fma = Dict(
              Int => has_fma_Int(),
              Rational{Int} => has_fma_Int(),
              Float32 => has_fma_Float32(),
              Float64 => has_fma_Float64(),
              BigFloat => true,
          )

If you’re, on the other hand, just asking about how to calculate the ULP error between a and b, the formula is basically (b - a) / eps(a). See the entire implementation here:

https://github.com/nsajko/VisualizeULPError.jl/blob/368d4291b0fa37241e3af1c7b340ee6f9dbf7e37/src/ULPError.jl

Pretty sure the answer is no.

GeorgeGkountouras · February 3, 2026, 7:53pm

Cool! Does it work for arbitrary Julia code e.g. my_cos() ? What about C functions?

woclass · February 4, 2026, 4:58am

I did some similar things for Bessels.jl, but just use relative error.

Relative error plot for Besseli(nu, x) · JuliaMath/Bessels.jl · Discussion #126

edit: remove images. Maybe there is something wrong with img upload.

OlivierHnt · February 4, 2026, 6:42am

Thx! Huh so all of these listed functions are implemented in pure Julia.

I was thinking of how, e.g., the gnu lib is architecture dependent.
So the only safe cases are the functions which are part of the IEEE 754. Anything beyond this, I suppose, would have to be a separate Julia package.

nsajko · February 4, 2026, 6:52am

Thank you!

Sure, just define the function inline when calling the app. I had to make slight adjustments (invokelatest) to make that work. The readme on Github now includes the following example:

To measure the error of an arbitrary function, define it when providing the command options:

visualize_ulp_error plot 'let sin_fast
    function sin_fast(x::BigFloat)  # accurate reference
        sin(x)
    end
    function sin_fast(x::Float64)   # fast approximation of `sin(x)` for `abs(x) ≤ 0.125`
        # Found with Sollya, using this command: `fpminimax(sin(x), [|1, 3, 5|], [|D...|], [0; 0.125]);`.
        p = (0.9999999999763268, -0.16666663941291426, 0.008328683245851542)
        x * evalpoly(x * x, p)  # odd polynomial composed from simpler polynomials
    end
    (;
        parent_dir = "/home/nsajko/ulp_error_plots",
        downsampled_length = 650,
        factor = 6000,
        window_size = 10,
        bf_precision = 140,
        no_scoped_values = true,
        width = 1100px,
        height = 500px,
        func = sin_fast,
        itv = (-0.125, 0.125),
    )
end'

The result:

A catch is that a Julia app can not(?) load more Julia packages. So, if you want to measure the error of a function defined in some package, you would need to either:

Redefine the function inline when calling the app, without referring to any package.

Use the package directly from the REPL, instead of using the Julia app. There is no supported API yet, but if you just want to hack around, something like this works at the moment:

using LogExpFunctions: logit            # Example: measure the error of `logit` from LogExpFunctions.jl
using VisualizeULPError: do_plot_impl
using Gadfly: px
options = (;
    parent_dir = "/tmp",
    downsampled_length = 650,
    factor = 5000,
    window_size = 10,
    bf_precision = 140,
    no_scoped_values = true,
    width = 1200px,
    height = 500px,
    func = logit,
    itv = (0.0, 1.0),
)
do_plot_impl(options)

Result for latest release of LogExpFunctions.jl, v0.3.29:

NB, that error spike is fixed on the master branch of LogExpFunctions.jl.

I suppose one could wrap a C function using @ccall.

Tamas_Papp · February 4, 2026, 7:36am

The 10^6 error for asech(1.0) looks bogus.

It would be better to make the plots

black, or at least dark blue,
with a larger axis labelling,
using thicker lines.

nsajko · February 4, 2026, 9:31am

I doubt that. In general, the error may be underestimated, but I can’t think of how it could be overestimated. A closeup, with parameters chosen to disable all smart functionality, showing the error for the rightmost hundred Float64 numbers in the domain of asech:

visualize_ulp_error -t5 -- plot '(;
    parent_dir = "/home/nsajko/ulp_error_plots",
    downsampled_length = 1000,        # this is excessive but there is no harm in that
    factor = 1,                       # disable downsampling
    window_size = 1,                  # disable smoothing
    bf_precision = 300,               # extra high precision, just to be safe
    no_scoped_values = true,
    width = 500px,
    height = 300px,
    func = asech,
    itv = (prevfloat(1.0, 99), 1.0),  # only consider the last hundred `Float64` numbers in the domain of `asech`
)'

Plot:

Reading off the plot: the error alternates between huge and small for successive floating-point numbers at the right edge of the domain. Adding this as a usage example to the readme.

Indeed, it is easy to check this in the REPL starting from first principles:

julia> const asech_accurate = Float64 ∘ asech ∘ big
Float64 ∘ asech ∘ big

julia> f(x) = (asech_accurate(x), asech(x))
f (generic function with 1 method)

julia> ulp_err(a, b) = abs(b - a)/eps(a)  # simplified definition, does not handle edge cases
ulp_err (generic function with 1 method)

julia> x = 1.0
1.0

julia> y = f(x)
(0.0, 0.0)

julia> x = prevfloat(x)
0.9999999999999999

julia> y = f(x)
(1.4901161193847656e-8, 2.1073424255447017e-8)

julia> ulp_err(y...)
1.865452045155277e15

Thank you. I intend to fix these issues without getting bogged down in such details of the presentation, by exporting the final data as Vega-friendly JSON, and/or relying on Vega.jl or Deneb.jl. Then the ecosystems around Vega should offer sane defaults for things like coloring, labeling, sizing, interactive zooming, etc. Issue: JSON export? Vega export? Use Vega.jl? Use Deneb.jl? · Issue #7 · JuliaMath/VisualizeULPError.jl · GitHub

foobar_lv2 · February 4, 2026, 10:18am

The large error for asech comes from x->(inv(x)-1).

But I don’t think this really matters? People who call asech/acosh for values close to the singularity at 1 deserve what they get, same as people who call log for values close to 1.

This is a good indication that we’re missing an acosh1p / asech1m in Base/stdlib though, analogue to log1p. (we could also consider inv1p)

PS. The fix would be to implement acosh1p and inv1p and then use asech(x)=acosh1p(inv1p) for values close to 1.

An more relaxed definition of relative error would be: err(f_approx, f, x) = abs(f(x)-f_approx(x)) / max(eps(x), eps(f(x))). I would call this the “garbage in, garbage out” principle.

GunnarFarneback · February 4, 2026, 12:08pm

I don’t think users of asech deserve any worse than users of acos. As it happens acos gives a more accurate result for asech than the current Float64 asech implementation itself does for 0.999999991838298 < x < 1.

foobar_lv2 · February 4, 2026, 12:41pm

Fair enough. Then we should add invm1(t) = inv(t) - one(typeof(t)) to Base, with appropriate specializations for floating point (probably use geometric series when t is close to 1, but magic float tricks are not my forte), and use it internally where appropriate. And then we also need to add acosh1p to Base (such that asech(t) = acosh1p(invm1(t)) close to 1). And then there is no good reason not to have asech1m as well.

GunnarFarneback · February 4, 2026, 2:38pm

Or we could implement asech analogously to acosh, e.g. something like

@noinline asech_domain_error(x) = throw(DomainError(x, "asech(x) is only defined for 0 ≤ x ≤ 1."))
function new_asech(x::T) where T <: Union{Float32, Float64}
    isnan(x) && return x

    if x < T(0) || x > T(1)
        return asech_domain_error(x)
    elseif x == T(1)
        return T(0)
    elseif x > T(0.5)
        t = T(1) - x
        return log1p((t + sqrt(T(2) * t - t * t)) / x)
    elseif x > T(0.5)^28
        y = 1 / x
        t = y * y
        return log(T(2) * y - 1 / (y + sqrt(t - T(1))))
    else
        return -log(x) + Base.Math.AH_LN2(T)
    end
end

As a bonus we get a better domain error and reasonable results for denormalized input.

Plot doesn’t look amazing but vastly better than today:

Topic		Replies	Views
Paper on the accuracy of math functions Numerics juliamath , accuracy	1	417	February 1, 2026
Trig functions very slow Performance	67	7170	October 10, 2018
Numerical difference between Julia and (C, Fortran, Python) for a chaotic damped driven pendulum Modelling & Simulations numerics , precision	17	2804	June 8, 2021
How would I use Supposition.jl on computer arithmetic corner cases? General Usage question	18	330	September 19, 2024
Base.runtests("math") failures on chips older than haswell New to Julia	12	747	August 7, 2021

Floating-point accuracy visualization example: the ULP error for each math function coming with Julia

Introduction

Julia app on Github

The plots

acos(::Float64)

acosd(::Float64)

acosh(::Float64)

acot(::Float64)

acotd(::Float64)

acoth(::Float64)

acsc(::Float64)

acscd(::Float64)

acsch(::Float64)

asec(::Float64)

asecd(::Float64)

asech(::Float64)

asin(::Float64)

asind(::Float64)

asinh(::Float64)

atan(::Float64)

atand(::Float64)

atanh(::Float64)

cos(::Float64)

cosc(::Float64)

cosd(::Float64)

cosh(::Float64)

cospi(::Float64)

cot(::Float64)

cotd(::Float64)

coth(::Float64)

csc(::Float64)

cscd(::Float64)

csch(::Float64)

deg2rad(::Float64)

exp(::Float64)

exp10(::Float64)

exp2(::Float64)

expm1(::Float64)

log(::Float64)

log10(::Float64)

log1p(::Float64)

log2(::Float64)

rad2deg(::Float64)

sec(::Float64)

secd(::Float64)

sech(::Float64)

sin(::Float64)

sinc(::Float64)

sind(::Float64)

sinh(::Float64)

sinpi(::Float64)

tan(::Float64)

tand(::Float64)

tanh(::Float64)

tanpi(::Float64)