fredrikj.net / blog /

# Floating-point wrappers in Arb

*September 25, 2021*

The biggest internal change in the just-released Arb version 2.21 (changelog) is a faster, completely rewritten gamma function. I will not write about the gamma function here since it is already discussed in a very detailed paper; instead, I will highlight a new convenience feature: the floating-point wrapper module (arb_fpwrap.h).

Most users don't need mathematical functions
with 1000-digit precision or sophisticated error tracking;
an implementation of a Bessel function
will usually be more than adequate
if it takes a `double` as input
and outputs a `double` accurate to 16 digits.
It is "easy" to create such a function by wrapping Arb.
However, in practice, this is not quite trivial unless you are
familiar with the library: you have to create Arb variables,
find the right conversion functions, write a loop that checks
the accuracy, and so on.
This evidently stumps a lot of people;
a common mistake is to omit the precision refinement (leading to unexpected inaccurate results),
and I get occasional emails asking about how to convert Arb numbers to floating-point numbers.
Some people probably give up before they even consider asking for help (I'm often that person when I'm checking out other software).
Some time ago, I created an external header file (arbcmath.h)
demonstrating how to wrap Arb,
but the interface was lacking, it was not kept up to date,
and as a header file, it was useless for interfacing from non-C/C++ languages.

So why not save everyone some trouble and provide floating-point wrapper functions in Arb itself? Such functions now exist, and they are this easy to use:

#include "arb_fpwrap.h" int main() { double x, y; complex_double cx, cy; int flags = 0; /* default options */ x = 2.0; cx.real = 0.5; cx.imag = 123.0; arb_fpwrap_double_zeta(&y, x, flags); arb_fpwrap_cdouble_zeta(&cy, cx, flags); printf("zeta(%g) = %.16g\n", x, y); printf("zeta(%g + %gi) = %.16g + %.16gi\n", cx.real, cx.imag, cy.real, cy.imag); flint_cleanup(); return 0; }

> build/examples/fpwrap zeta(2) = 1.644934066848226 zeta(0.5 + 123i) = 0.006252861175594465 + 0.08206030514520983i

The documentation includes sample code for calling wrapper functions from Python (using `ctypes`) and Julia (using `ccall`).
It should be just as easy in any other language with a C foreign function interface.
(The only tricky detail is that
the output value is passed by reference; this is necessary to allow functions to return
an error code and so that functions with multiple outputs can have a consistent interface.)

A motivating example: Julia users might be disappointed when they learn that the SpecialFunctions.jl package doesn't support Bessel functions of complex order:

julia> import SpecialFunctions: besselj julia> besselj(2.0 + 3.0im, 3.0 + 4.0im) ERROR: MethodError: no method matching besselj(::Complex{Float64}, ::Complex{Float64}) Closest candidates are: besselj(::Float64, ::Complex{Float64}) at /home/fredrik/.julia/packages/SpecialFunctions/f769j/src/bessel.jl:380 besselj(::T, ::Complex{T}) where T<:AbstractFloat at /home/fredrik/.julia/packages/SpecialFunctions/f769j/src/bessel.jl:590 besselj(::Real, ::Complex) at /home/fredrik/.julia/packages/SpecialFunctions/f769j/src/bessel.jl:584 Stacktrace: [1] top-level scope at REPL[2]:1

Arb to the rescue! We can just `ccall` the `arb_fpwrap_cdouble_bessel_j` function:

julia> function besselj(nu::Complex{Float64}, x::Complex{Float64}) cy = Ref{Complex{Float64}}() if Bool(ccall((:arb_fpwrap_cdouble_bessel_j, :libarb), Cint, (Ptr{Complex{Float64}}, Complex{Float64}, Complex{Float64}, Cint), cy, nu, x, 0)) error("unable to evaluate accurately at ", x) end return cy[] end besselj (generic function with 11 methods) julia> besselj(2.0 + 3.0im, 3.0 + 4.0im) 0.44014236912314847 - 0.11722919001237514im

Granted, we still needed some boilerplate, but it only involves native Julia types.

Floating-point wrappers around Arb will usually be much (10-1000 times) slower than optimized hardware-precision implementations of special functions. The advantage of the Arb functions is that they guarantee high accuracy (optionally, guaranteed correct rounding) on a large domain. If Arb is unable to compute a certified accurate result, the floating-point wrappers output NaN and return a status flag instead of silently returning plausible-looking but wrong numerical values.

In rare cases, Arb can be faster than hardware-precision numerical libraries when those are not particularly optimized:

julia> using BenchmarkTools julia> import SpecialFunctions:zeta julia> @btime zeta(0.5 + 1e8im) 11.415 s (0 allocations: 0 bytes) -3.3628392091065953 + 1.4072345711319583im julia> function myzeta(x::Complex{Float64}) cy = Ref{Complex{Float64}}() if Bool(ccall((:arb_fpwrap_cdouble_zeta, :libarb), Cint, (Ptr{Complex{Float64}}, Complex{Float64}, Int), cy, x, 0)) error("unable to evaluate accurately at ", x) end return cy[] end myzeta (generic function with 1 method) julia> @btime myzeta(0.5 + 1e8im) 1.653 ms (0 allocations: 0 bytes) -3.362839487530728 + 1.407234559646448im

The code is still a bit experimental and has not been tested thoroughly.
Only `double` precision is supported for now, but
the wrapper module could just as well be expanded to provide
accurate math functions for other formats like IEEE 754 quadruple precision,
decimal types,
arbitrary-precision types (arf, MPFR, MPC), etc.

Other future improvements might include more wrapped functions, support for different rounding modes, support for negative zero (currently handled as positive zero), evaluation at infinities (currently considered an error and giving NaN), and providing correct rounding at special points where the correct rounding test currently fails.

fredrikj.net | Blog index | RSS feed | Follow me on Mastodon | Become a sponsor