MSVC is experimenting with correctly-rounded math (LLVM libc)

Offtopic
1

Since 2015, the math functions in cmath and cstdlib are supplied by the C runtime via math.h and stdlib.h, which Microsoft ships as part of the Universal C Runtime (UCRT).
…

In the modern day, the UCRT math functions have known mathematical inaccuracies. See Accuracy of Mathematical Functions in Single, Double, Double Extended, and Quadruple Precision by Gladman, et al. for more info (it’s hot off the press). These inaccuracies stem from implementations that are quite old, dating to times when accuracy was often computationally infeasible.

We found that the LLVM C Library had already made incredible progress tackling the math portions of the C runtime. The project upholds accuracy as the primary goal, aiming for correct rounding in all rounding modes, with Gladman, et al. confirming the project’s success. By targeting mathematical accuracy, the library effectively codifies a stable interface that allows for implementation flexibility and optimization. The library is actively maintained and has a healthy community supporting it.

MSVC C++23: constexpr cmath with LLVM Libc - C++ Team Blog

Also note that:

  • LLVM-libc provides some correctly rounded functions (see links below): claimed accuracy of LLVM functions
  • seven binary64 functions are integrated into the GNU libc up from release 2.43 (…), and the following binary32 functions are integrated into the GNU libc up from release 2.42 (…)
  • AMD libm (since 4.2) integrates tanhf from CORE-MATH, and uses some CORE-MATH test cases for its erf function
  • the Intel Math Library or IML (checked with 2026.1.0) includes some (apparently undocumented) cr_xxx functions

– The CORE-MATH project

And Rust libm (compiler-builtins)

2

It seems more and more projects are moving toward correctly-rounded libm.
Does Julia need to follow this wave?

Right now Julia’s math is ported from openlibm (and FDLIBM),
which is not correctly rounded.

And maybe now is a good time to remove openlibm as a dependency?

3

This is certainly something to consider, because of e.g. https://core-math.gitlabpages.inria.fr/cbrt64.pdf Not even OpenLibm 0.8.7 is accurate for sin, cos (or tan, or any other trigonometry seemingly like atanh), even for Float32 (I doubt 0.8.8 from last week changes anything).

While Julia doesn’t actually use OpenLibm for much if anything anymore, I believe the Julia standard library was ported from it. I see cbrt, sin and cos are basically unchanged from 9 years ago.

Intriguing: “pown has strictly mandated IEEE 754 compliance” (unlike powi), and was only adopted in C23, was in IEEE since 2008, optional in C. Julia needs neither, uses multiple dispatch, but likely does something closer to the less accurate powi.

Julia does have cbrt accurate, but only for Float32. LLVM libc has it for Float64 too.

cbrt calls _approx_cbrt and _improve_cbrt neither in OpenLibm. simplify cbrt codepaths

rsqrt for Float32 is accurate (0.500 in table) in LLVM and IML 2026.1.0. The latter has rsqrt for Float64 most accurate (0.501), and atan2pi (0.528). I think Julia has no version fast or slow or rsqrt (I’ve suggested that fast approximation before).

LLVM (“LLVM libc extension”) has f16sqrt correctly rounded, something to consider for Julia.

GNU libc 2.44 has lgamma, tgamma, erf and erfc accurate also for Float64. Where do you get these and compoundn?

I believe all univarate functions are already correctly rounded in some math library for Float32 (reading the table i.e. where 0.500), since can be exhaustively checked, but such is not possible for Float64. So I’m intrigued to see there the proofs of rounding.

The expections or problematic I see are atan2pi, atanpi, compoundn, erfc, lgamma, log10p1, log2p1, logp1, rootn, tgamma.

j0, and j1 and y0 and y1, are hard also apparently, IML 2026.1.0 does them best.

Using this approach, we have created RLIBM-32, a library containing implementations of several elementary functions that produce correctly rounded results for all inputs for the 32-bit float type and the 32-bit posit type.
..
The functions in RLIBM-32 are significantly faster than the state of the art. The above graphs show the speedup of RLIBM-32’s float functions compared to glibc’s and Intel’s libm. On average, RLIBM-32 has 1.1x and 1.2x speedup over glibc’s float and double functions. RLIBM-32 has 1.5x and 1.6x speedup over Intel’s float and double functions. Overall, RLIBM-32 not only produces correctly rounded results for all inputs but it is faster than mainstream math libraries, which have been optimized for decades.

Only LLVM libc claims pow accurate; in Float32 (and 1 ULP for FLoat64) but elsewhere I see only 0.501 elsewhere for LLVM 22.1.8.

I’m a bit confused why Julia has (the OpenLibm dependency any more and) OpenLibm_jll. This was used for 64-bit and 32-bit, and if I recall only not used for 32-bit Windows, why I suggested dropping OpenLibm and 32-bit Windows at the time, now it has been decided to drop 32-bit Windows.

I see Julia developers still maintain OpenLibm, for long double and more that doesn’t apply to Julia:

4

In case you hadn’t seen it, the paper

analyses the accuracy of Julia’s mathematical functions. Apart from the hyperbolic functions which clearly need some love, almost all others have less than 1 ULP of error (the only exceptions being exp10 in single precision with a maximum error of 1.05 ULP, and tan in double precision with worst error of 1.09 ULP, but note that the double precision analysis is non-exhaustive). While almost none of the functions is correctly rounded (only correctly rounded one is sqrt, which just calls llmv.sqrt), I’d say they mostly do a decent job at balancing speed and accuracy.

5

FYI: I’ve learned there is available CoreMath.jl for

Correctly-rounded mathematical functions with CORE-MATH

Intriguingly, it’s also faster for some (Float32) values (but slower for Float64, at least corresponding), for sinh, and it and cosh and all hyperbolic are rather inaccurate in Julia:

julia> a = 10.0f0;

julia> @btime cr_sinh($a)
  14.338 ns (0 allocations: 0 bytes)
11013.232f0

julia> @btime sinh($a)
  18.879 ns (0 allocations: 0 bytes)
11013.233f0

julia> a = 10.0

julia> @btime cr_sinh($a)
  29.177 ns (0 allocations: 0 bytes)
11013.232874703393

julia> @btime sinh($a)
  16.238 ns (0 allocations: 0 bytes)
11013.232874703393

I’m looking into cr_cbrt as we speak.

julia> time_cbrt()

Timing cbrt over all 2^32 Float32 values…

Timing cr_cbrt over all 2^32 Float32 values…

Results

cbrt: 47.377 seconds
cr_cbrt: 46.909 seconds
cbrt: 11.031 ns/call
cr_cbrt: 10.922 ns/call
ratio: 0.990111x

Checksums:
cbrt: 0x00000000
cr_cbrt: 0x00000000
(identical)
(0x0000000b07e441dc, 0x0000000aebf7868c)

Before other way around (might be noisy machine or interference because of rand):

Generating 1000000 random inputs…

======================================================================
Float32

Accuracy:
Julia cbrt: max error = 0 ULP
cr_cbrt: max error = 0 ULP
Julia non-CR results: 0 / 1000000
cr_cbrt non-CR results: 0 / 1000000

Speed:
Julia cbrt: 9.279 ns/call
cr_cbrt: 9.496 ns/call
ratio: 1.023x

[Both numbers lower do not make sense, if I’m also testing rand with it).]

I need to check my code because before mismatch (likely only for NaNs, then that would not be worrying):

julia> mismatches = exhaustive_compare()
Checking all 2^32 Float32 bit patterns…
This is 4294967296 values.

Mismatches: 8388606

First mismatch:
input bits: 0x7f800001
input: NaN
cbrt bits: 0x7f800001
cr_cbrt bits:0x7fc00001
cbrt: NaN
cr_cbrt: NaN
0x00000000007ffffe