Arblib.jl ❤️ IntervalArithmetic.jl

[Joel-Dahne]

As discussed in #717, there is potential in making Arblib.jl and IntervalArithmetic.jl work together. I see two potential benefits:

1. Users can easily switch between doing different parts of their computations in different packages depending on which package suits that particular part of the computation best.
2. The IntervalArithmetic.jl library could use some of the Arblib.jl functionality to speed up certain computations (in particular high precision computations).

For the first point it would suffice to create a package extension, which in principle could live in either of the packages (though as I extend upon below it would likely be more natural for it to live in the IntervalArithmetic.jl package). The second point would mean IntervalArithmetic.jl taking on Arblib.jl as a full dependency.

For users

To make the packages work well together for users the main thing to implement would be conversions between intervals and the types in Arblib.jl (primarily Arb and Acb, representing real and complex balls respectively). The user could then easily convert between the two packages to suit their needs.

Some things that Arblib.jl excel at:

1. Fast high precision computations
2. A large library of special functions
3. Support for mutable arithmetic
4. Taylor series expansions (in particular at high degrees)

Some things that IntervalArithmetic.jl excel at:

1. Fast computations at Float32 and Float64 precision
2. Computations with wide intervals (Arblib.jl is not optimized for this, though the situation is improving)
3. Built in safety features such as decorations (it is much easier to shoot yourself in the foot with Arblib.jl)
4. Plotting (this is actually my primary use case for using it)

In principle the conversion functions could be very simple. I have been using

Arblib.Arb(x::Union{Interval,BareInterval}) =
    isempty_interval(x) ? indeterminate(Arb) : Arb((inf(x), sup(x)))

Arblib.Acb(z::Complex{<:Interval}) = Acb(Arb(real(z)), Arb(imag(z)))

IntervalArithmetic.interval(::Type{T}, x::Arb) where {T} =
    isnan(x) ? nai(T) : interval(T, getinterval(BigFloat, x)...)

IntervalArithmetic.interval(::Type{T}, (x, y)::NTuple{2,Arb}) where {T} =
    (isnan(Arblib.midref(x)) || isnan(Arblib.midref(y))) ? nai(T) :
    interval(T, BigFloat(lbound(x)), BigFloat(ubound(y)))

There are a number of design decisions to make though. Some that come to mind are:

1. Should implicit conversion from Arb to Interval set the NG flag? There is no decorations for Arb so in principle there are no guarantees given.
2. What should the default T type be (should there be a default?) when converting Arb to Interval{T}?
3. Should we allow intervals of type Interval{Arb}?

Most of these design decisions are on the IntervalArithmetic.jl side. This is the reason I believe (as I mentioned above) that a potential extension should probably live on the IntervalArithmetic.jl side.

For IntervalArithmetic.jl

I believe the primary use case for Arblib.jl in IntervalArithmetic.jl would be to speed up certain computations with BigFloat intervals.

One use case is for linear algebra, where Arblib.jl has very fast implementations for certain things. This includes for example matrix multiplication, computation of inverses and solving linear systems. I would expect the Arblib.jl implementations to beat the BigFloat implementations by about 2 orders of magnitude in many cases. See the Flint documentation for real matrices and complex matrices for what type of things are implemented (not all of these have (convenient) wrappers in Arblib.jl at the moment though). Here one could argue that if we implement convenient conversion between Arblib.jl types and intervals then the user could be responsible for using the fast algorithms in Arblib.jl.

Another use case is for interval extensions of special functions (anything not currently implemented by BigFloat). Note that Arblib.jl usually doesn't handle wide intervals well. The approach to take would therefore likely be to use Arblib.jl to compute correctly point values and then implement the logic for which points the function needs to be computed at in IntervaArithmetic.jl. Though handling special functions might be outside the scope of IntervalArithmetic.jl?

There might be more use cases that I haven't thought about as well!

[OlivierHnt]

Thank you @Joel-Dahne for opening such a comprehensive issue. The way forward is to add Arblib.jl as a weak dependency (package extension); this goes for special functions, and the matrix multiplication algorithm as well.

1. Should implicit conversion from Arb to Interval set the NG flag? There is no decorations for Arb so in principle there are no guarantees given.

The "NG" flag is not tied to decorations; it's a marker for conversions between numerical types and Interval. Since Arb computations are rigorous (i.e., take into account numerical uncertainties), constructing an Interval from Arb does not warrant an "NG" flag in my opinion. It should probably just match how interval(::BareInterval) is handled.

2. What should the default T type be (should there be a default?) when converting Arb to Interval{T}?

I would probably suggest BigFloat, since Arblib operates with multiple precision. We'll have to make sure to match the precision properly.

3. Should we allow intervals of type Interval{Arb}?

I don't think so.