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