Repository navigation
Tight CORE-MATH bounds for every elementary function, from C, without MPFR: useful for :correct? #784
Description
IntervalArithmetic's default :correct rounding uses CRlibm for 18 functions and MPFR for the rest: exp2, exp10, cbrt, tanh, the inverse hyperbolics, and pow's and atan2's corners. CORE-MATH could give those bounds, correctly rounded down and up, but Julia can't set the rounding mode, so the directed results have to come from C.
ival (MIT, part of crnumerics) is a C library that does this.
- Functions: 32 CORE-MATH functions as tightest intervals, the mode set inside the call.
- Calls: arrays of intervals per call, or one interval at a time at 17 to 54 ns a call (add to exp, on an x86-64 laptop), so a scalar path through
ccallworks too. - The rest of IEEE 1788.1: arithmetic, reverse operations, the constructors.
- Checks: against independent MPFR references, and against ITF1788, where all 7,451 bare tests it runs give tight results. Ten of ITF1788's expected results turned out not to be the tightest; each correction comes with a proof (two are in pow_rev.itl: two powRev2 expected results are not the tightest ITF1788.jl#8, with a fix in pow_rev.itl: the tightest results for two powRev2 tests ITF1788.jl#9).
- An accurate mode, as
:ulpdoes but over arrays, through crmvec's vector functions: exp 12 ns, sin 29 ns per interval. - Installing: released (crnumerics 0.2.0), with
make install, a versionedlibival.so.0and pkg-config, which a BinaryBuilder recipe can use as it is.
Against IntervalArithmetic.jl's :correct itself, on the same 4,096 narrow intervals, every bound is the same, bit for bit. Times are ns per interval on one EPYC 7773X core, with IntervalArithmetic 1.0.12 on Julia 1.11.7:
| exp (CRlibm) | exp2 (MPFR) | tanh (MPFR) | |
|---|---|---|---|
IntervalArithmetic.jl, :correct |
70 | 3,130 | 3,016 |
| ival, tight | 30 | 31 | 60 |
Two questions:
- Would a JLL of ival be welcome as
:correct's path for the functions CRlibm lacks? - Or would directed entry points (
cr_exp_rdand the like) in CoreMath.jl's own build fit better?
The second is smaller, and I'd be glad to help with either.
anun333. Written with an AI assistant (Claude); every figure is reproducible from the repository (make ival-julia-compare, and make ival-compare for the C libraries).
Activity
OlivierHnt commentedon Oct 11, 2026
I am not sure this is correct. As far as I understand, the problem is mainly on the side of LLM which does not respect directed rounding.
So I am not sure how your suggestion fixes this.
anun333 commentedon Oct 11, 2026
If I understand you right, the problem is that LLVM assumes round-to-nearest in the code it compiles, so Julia can't safely change the rounding mode itself.
ival doesn't ask it to. Each call changes the mode only inside ival's own C code (compiled with
-frounding-math) and restores the caller's floating-point environment before returning, so no Julia code runs in another mode.I checked that from Julia: through
ccall, on 4,096 intervals each of exp, exp2, exp10, tanh, cbrt, asinh, acosh and atanh, every bound equals:correct's, and Julia is in round-to-nearest before and after every call (IntervalArithmetic 1.0.12, Julia 1.11.7). The script is ival/julia/ccall-check.jl (make ival-julia-ccall).Written with an AI assistant (Claude);
make ival-julia-ccallreproduces the check.OlivierHnt commentedon Oct 11, 2026
It is possible that compiling small C kernels for the functions in CORE-MATH (eg
exp_rd,exp_ru) work with the proper floating-point rounding flag. Though I am not entirely sure, the compiler works in mysterious ways.Currently we emulate, which is not ideal, but I think at least is thread safe.
anun333 commentedon Oct 11, 2026
On thread safety: the rounding mode is part of each thread's floating-point environment (C11 7.6: it "has thread storage duration"; MXCSR on x86-64, FPCR on AArch64), so a change inside a call is seen only by the calling thread, and ival restores it before returning. I measured that as well: with 3 and with 7 Julia threads calling ival at once (25,600 calls of 4,096 intervals), every bound equals
:correct's, while a task on another thread checked its own rounding mode 86 to 152 million times and always found round-to-nearest (threads-check.jl,make ival-julia-threads).On the compiler, you're right to be wary: ival's C is built with
-frounding-math -ffp-contract=off, which GCC documents as not guaranteed to stop every transformation the mode affects. That is why it is checked against ITF1788 and MPFR references on every platform; on macOS those checks caught Apple Clang moving inlined square roots past a mode change, now prevented with avolatile.If small kernels are the way you'd like to go, I could add
_rd/_ruentry points for the functions:correcttakes from MPFR to the CoreMath recipe on Yggdrasil (each sets the mode, callscr_*and restores it), with the matching calls in CoreMath.jl, for you to review.Written with an AI assistant (Claude);
make ival-julia-threadsreproduces the check.