fredrikj.net / blog /

# Arb 2.0 preview: intervals at the speed of points

*May 11, 2014*

The forthcoming Arb 2.0 will feature a new floating-point type `arf_t` and a new ball interval type `arb_t` with better performance than the present `fmpr_t` and `fmprb_t` types.

Arb 1.0 did achieve the goal of providing interval arithmetic with negligible overhead compared to floating-point arithmetic... at precisions in the thousands of digits. Calculations done at more typical precisions of perhaps 10-100 digits could be a factor 2-4 slower in the worst case . Arb 2.0 is intended to bring this down to a factor 1.0-1.2.

I'm pleased to report that work on the
new types is going very well. The new `arb_t` ball type is
already complete except for transcendental functions. I will probably
implement transcendental functions initially by wrapping the
`fmprb_t` versions, and only speeding them up later. Complex numbers, polynomials and matrices should be fairly easy to port to the new types at this point.

Below, I've benchmarked various arithmetic operations with the old and new types. The plots show time normalized against MPFR (lower is better). Red is old types, blue is new. The thick solid lines show the cost of ball arithmetic (`fmprb_t`, `arb_t`), and the thin dashed lines show the cost of floating-point arithmetic (`fmpr_t`, `arf_t`) to give an idea of how much time is wasted on propagating error bounds.

Addition is now typically no more than 20% slower than MPFR. Note that `fmpr_t` special-cases addition up to 64 bits, which is why it is significantly faster than all other types in this range. I will also implement this optimization for the `arf_t` type in the future.

The improvement for multiplication is huge. This could be as much as three times slower than MPFR with the `fmprb_t` type, but like addition, it is now within 20% even at low precision, where propagating error bounds is a lot of extra work (if the input balls are $[x \pm r_x], [y \pm r_y]$, one has to compute $[xy \pm (|x|r_y + |y|r_x + r_x r_y+\varepsilon_{\text{rounding}})]$ rather than just $xy$). Note also that from about 2048 bits and above, `arf_mul` now also piggybacks on MPFR's mulhigh code, giving a 20% improvement at precisions in the thousands of digits.

The most important operation is fused multiply-add or addmul, i.e.
computing $z \gets z + xy$. This is the building block of matrix
multiplication, classical polynomial multiplication, and many other
things. In this plot, I have used the cost of `mpfr_fma` as the
baseline, but this function is actually much slower than doing an
inexact `mpfr_mul` followed by an `mpfr_add` (assuming
that you already have a temporary variable allocated), illustrated by
the lower black curve. This problem was discussed at the MPFR-MPC
developers meeting back in January, and should hopefully be fixed in
a future version of MPFR. Anyway, the new `arb_addmul` is much
faster than `fmprb_addmul`. Unlike `arf_mul/arb_mul`, it
does not
implement the mulhigh trick yet (to be added in the future).

An important variation is to compute $z \gets z + nx$ where $n$ is a fixed word-size integer. Here I used `mpfr_mul_ui`+`mpfr_add` with a preallocated temporary as the baseline. Once again, there's a significant improvement for the `arb_t` type.

Division has been improved a lot as well. The new `arb_div` is at most two times slower than the raw floating-point division and comes close around 100 digits of precision (compared to "six times slower" and "1000 digits" for `fmprb_div`). Most of the speedup comes from computing error bounds faster, and removing general overhead. I don't see an easy way to get rid of the remaining overhead, but it might be possible.

Finally, the improvement for square roots is similar to that for division.

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