Add extension for Arblib.jl

[Joel-Dahne]

As discussed in #721, it would be nice to make the Arblib.jl and IntervalArithmetic.jl packages talk to each other. This is an attempt at making basic conversion between Interval and Arb work.

julia> using Arblib, IntervalArithmetic

julia> x = interval(BigFloat, π)
[3.14159, 3.1416]₂₅₆_com

julia> y = Arb(x)
[3.1415926535897932384626433832795028841971693993751058209749445923078164062862 +/- 5.38e-77]

julia> cos(x)
[-1.0, -0.999999]₂₅₆_com

julia> cos(y)
[-1.000000000000000000000000000000000000000000000000000000000000000000000000000 +/- 1.78e-77]

It is also possible to construct matrices and do basic linear algebra

julia> A = interval.(rand(BigFloat, 2, 2))
2×2 Matrix{Interval{BigFloat}}:
 [0.15383, 0.153831]₂₅₆_com  [0.655025, 0.655026]₂₅₆_com
 [0.165259, 0.16526]₂₅₆_com  [0.828643, 0.828644]₂₅₆_com

julia> B = ArbMatrix(A)
2×2 ArbMatrix:
 [0.15{…73 digits…}85 ± 3.64e-78]  [0.65{…73 digits…}67 ± 4.69e-78]
 [0.16{…73 digits…}30 ± 4.82e-78]  [0.82{…73 digits…}91 ± 1.48e-78]

julia> Arblib.overlaps(B * B, ArbMatrix(A * A))
true

Implementation

Construction of Arb from Interval is straightforward. It only requires overloading Arblib.set! to handle inputs of type Interval. To make some of the internal precision handling correct, I also updated Arblib._precision to be aware of the Interval type.

I realized my first attempt at the construction of Interval from Arb was incorrect. The second attempt was super complicated and used a lot of internals of IntervalArithmetic. The final attempt ended up being fairly simple. For construction, we basically want to treat Arb as equivalent to Interval{BigFloat}. For this, we define

IA._inf(x::Arblib.ArbOrRef) = BigFloat(Arblib.lbound(x))
IA._sup(x::Arblib.ArbOrRef) = BigFloat(Arblib.ubound(x))

and make sure that IntervalArithmetic.promote_numtype handles promotion with Arb in the way we want. Note that this generates quite a lot of allocation, to compute _inf(x) it first computes the lower bound as and Arf, then converts it to a BigFloat and then further converts it to whatever type the interval should be. I could however note figure out a better way to do it at this point.

I also made the following updates:

• Made Base.promote_rule explicitly throw an error for promotion between Arb and Interval.
• Handle some ambiguities related to ExactReal as well as Arf.

I added a lot of tests. This was partially because my first few attempts of implementing the constructor were fairly complicated and I wanted to make sure that I checked all edge cases. The final implementation is a lot simpler and maybe doesn't warrant quite as many tests. But testing more never hurts!

Note that a few of the implemented tests are marked as broken. While doing the implementation, I noticed some issues with the Arblib.getinterval function for infinite intervals with infinite endpoints. This is due to an issue in Flint that I reported.

[OlivierHnt]

This looks great!
So presumably BigFloat(Arblib.lbound(x)) is equivalent to BigFloat(Arblib.lbound(x); precision = precision(Arblib.lbound(x)))?

In #721, we also mentioned special functions; do you think doing something like

for f ∈ list_of_special_functions
    @eval Arblib.$f(x::Interval) = interval($f(Arblib.Arb(x)))
end

is worth it?

In the same spirit, maybe we can address #717 as well, which is to add the Arblib routines for multiprecision matrix product. I am not sure what is the best way of doing so, maybe

function IA._mul!(::IntervalArithmetic.MatMulMode{:fast}, C::AbstractVecOrMat{<:RealOrComplexI{<:BigFloat}}, A, B, α, β)
    # perform the fast Arblib algorithm in extended precision
end

is ok? One thing that could be puzzling is that loading Arblib will change the result (though still a rigorous enclosure) of the matrix product when matmul is set to :fast... or we disable the :fast algorithm on BigFloat unless Arblib is loaded...

[Joel-Dahne]

So presumably BigFloat(Arblib.lbound(x)) is equivalent to BigFloat(Arblib.lbound(x); precision = precision(Arblib.lbound(x)))?

Yes, so the conversion to BigFloat is exact in this case (as long as there is no overflow in the exponent, in which case you would get Inf or -Inf).

I think directly wrapping the special functions come with some issues. Ball arithmetic (like Arb uses) and infsup interval arithmetic (which IntervalArithmetic uses) tend to have fairly different assumptions about how in particular wide inputs are handled. As an example consider evaluating the gamma function on the interval $[2, 20]$ (where it is increasing). With direction evaluation with Arb we get

julia> using Arblib, IntervalArithmetic

julia> interval(Arblib.SpecialFunctions.gamma(Arb((2, 20))))
[-1.67773e+09, 1.21646e+17]₂₅₆_com

whereas most user would probably expect something more like

julia> interval(Arblib.SpecialFunctions.gamma(Arb(2)), Arblib.SpecialFunctions.gamma(Arb(20)))
[1.0, 1.21646e+17]₂₅₆_com

If the user explicitly converts to Arb, calls the special functions and then converts back it would hopefully be more clear to them where this discrepancy is coming from.

For #717 I'm not sure what the best approach is, it is certainly a bit awkward to have to different package extensions interact in this way. That LinearAlgebra is always implicitly loaded by Julia at least makes it easier. But I could definitely see people being confused about the result of matrix multiplication changing if they load another package. In particular since they don't even have to load it explicitly, but it is enough that they load another package that depend on it.

I was thinking that I could add some page to the documentation, something like "Interfacing with Arblib.jl", and write something about the discussions here. Then users would at least be aware of the existence of these tools. Then we could leave any direct implementations to the future. Does that sound reasonable to you? In that case I should hopefully be able to write some sort of documentation by the beginning of next week.

[Kolaru]

One thing that could be puzzling is that loading Arblib will change the result (though still a rigorous enclosure) of the matrix product when matmul is set to :fast...

I think that this should never happen. I would prefer to disable the :fast matmul without Arblib, as you mentionned.

Ball arithmetic (like Arb uses) and infsup interval arithmetic (which IntervalArithmetic uses) tend to have fairly different assumptions about how in particular wide inputs are handled. As an example consider evaluating the gamma function on the interval [ 2 , 20 ] (where it is increasing). With direction evaluation with Arb we get

Improving the tightness of our results is not considering breaking, so it not too much of a problem: we can have loose support right now for convenience's sake, and improve it latter. In particular, in the future, we could change the definition of monotonous special functions to do the smarter thing, and (now that we have piecewise functions) we could even do it for piecewise monotonous ones :D

That being said, I don't have a strong opinion either way. In fact, it may be more fitting for a SpecialFunctions.jl extension.

[OlivierHnt]

I think adding a documentation page with short descriptions of all our package extensions is a great idea. Since we haven’t done this before, @Joel-Dahne, feel free to get it started in whatever format you think works best.
After that, let us merge this PR.

[Joel-Dahne]

I have now added some documentation for the interface. The goal was to give users some ideas of what type of computations could be reasonable to use Arblib.jl for. In the future we might add a similar page of the documentation to Arblib.jl to discuss what type of computations it could be reasonable to use IntervalArithmetic.jl for.

I have also added conversion between complex intervals and Acb and that I have fixed conversion to so that it doesn't add the NG flag for implicit conversions from Arb and Acb.

Let me know if you have any comments, in particular regarding the documentation! Otherwise I think this should in principle be ready to merge.

I had some issues building the documentation locally, it didn't load the package extension as I expected. It seems to work in the CI though.

[codecov-commenter]

⚠️ Please install the to ensure uploads and comments are reliably processed by Codecov.

Codecov Report

Attention: Patch coverage is 72.72727% with 12 lines in your changes missing coverage. Please review.

Project coverage is 78.16%. Comparing base (31bfb51) to head (b3611b2).
Report is 13 commits behind head on master.

Files with missing lines	Patch %	Lines
ext/IntervalArithmeticArblibExt.jl	72.72%	12 Missing ⚠️

❗ Your organization needs to install the Codecov GitHub app to enable full functionality.

Additional details and impacted files

@@            Coverage Diff             @@
##           master     #723      +/-   ##
==========================================
+ Coverage   78.02%   78.16%   +0.14%     
==========================================
  Files          30       31       +1     
  Lines        2930     2977      +47     
==========================================
+ Hits         2286     2327      +41     
- Misses        644      650       +6     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[Joel-Dahne]

Added some more tests, hopefully Codecov is happier now!

[OlivierHnt]

Is a default value for the decoration really needed? And if so, shouldn't that be IA.com?

[Joel-Dahne]

Good point, it should be IA.com. Indeed, it is probably not needed. I copied it from the here and I see now that the these have default values for d but the versions above don't. I'll remove the default value.

[Joel-Dahne]

I removed the default argument for d as mentioned. I also updated the handling of some tests that were broken due to a bug in Flint. The bug is fixed in the latest release of Flint (from yesterday) and I thought I would write the tests in a way that they check if the bug is there or not.