Add extension for Arblib.jl

Project.toml

@@ -11,6 +11,7 @@ Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 RoundingEmulator = "5eaf0fd0-dfba-4ccb-bf02-d820a40db705"
 
 [weakdeps]
+Arblib = "fb37089c-8514-4489-9461-98f9c8763369"
 DiffRules = "b552c78f-8df3-52c6-915a-8e097449b14b"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
@@ -19,6 +20,7 @@ RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [extensions]
+IntervalArithmeticArblibExt = "Arblib"
 IntervalArithmeticDiffRulesExt = "DiffRules"
 IntervalArithmeticForwardDiffExt = "ForwardDiff"
 IntervalArithmeticIntervalSetsExt = "IntervalSets"
@@ -27,6 +29,7 @@ IntervalArithmeticRecipesBaseExt = "RecipesBase"
 IntervalArithmeticSparseArraysExt = "SparseArrays"
 
 [compat]
+Arblib = "1.3.0"
 CRlibm = "1.0.2"
 DiffRules = "1"
 ForwardDiff = "0.10, 1"

docs/Project.toml

@@ -1,3 +1,6 @@
 [deps]
+Arblib = "fb37089c-8514-4489-9461-98f9c8763369"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"

docs/make.jl

@@ -22,6 +22,9 @@ makedocs(;
         "Manual" => [
             "Constructing intervals" => "manual/construction.md",
             "Usage" => "manual/usage.md",
+            "Interfaces" => [
+              "Arblib.jl" => "manual/interfaces/Arblib.md"
+            ],
             "API" => "manual/api.md"
         ],
         "Examples" => ["Rigorous computation of ``\\pi``" => "examples/pi.md"],

docs/src/manual/interfaces/Arblib.md

@@ -0,0 +1,116 @@
+# Interfacing with Arblib.jl
+
+[Arblib.jl](https://github.com/kalmarek/Arblib.jl/) is a Julia package
+for rigorous numerics based on ball arithmetic. It provides a thin and
+efficient wrapper around the parts of the C library
+[FLINT](https://flintlib.org/) that are concerned with real and
+complex numbers.
+
+The purpose of this interface is to make conversions between the types
+used by IntervalArithmetic.jl and Arblib.jl straightforward. This
+allows for easy switching between doing different parts of
+computations in different packages, depending on which package suits
+that particular part of the computation best.
+
+Some of the things that Arblib.jl excels at are:
+1. Fast high precision computations, including linear algebra routines.
+2. A large library of special functions.
+3. [Support for mutable arithmetic](https://kalmarek.github.io/Arblib.jl/stable/interface-mutable/).
+4. [Taylor series expansions](https://kalmarek.github.io/Arblib.jl/stable/interface-series/).
+
+Some things that IntervalArithmetic.jl excels at:
+1. Fast computations at `Float32` and `Float64` precision.
+2. Computations with wide intervals (FLINT is in general not optimized
+   for this, though the situation is improving).
+3. Built-in safety features, such as decorations and the `NG` flag.
+   Arblib.jl provides fewer built-in safeguards, requiring more
+   careful handling.
+
+## Conversion between `Interval` and Arblib.jl types
+The fundamental types in Arblib.jl are `Arb` and `Acb`, corresponding
+to real and complex balls respectively. Conversion between `Interval`
+and `Arb` as well as between `Complex{Interval}` and `Acb` is done
+through the standard constructors. If no type is specified it defaults
+to `Interval{BigFloat}` for `Arb` and `Complex{Interval{BigFloat}}`
+for `Acb`.
+
+``` @repl 1
+using Arblib
+
+x = interval(π)
+Arb(x)
+
+y = Arb(π)
+interval(y) # Defaults to BigFloat
+interval(Float64, y) # Other type can be explicitly specified
+
+z = complex(interval(2), interval(1 // 3))
+Acb(z)
+
+w = Acb(2, 1 // 3)
+interval(w)
+interval(Float64, w)
+```
+
+## Linear algebra
+To use the optimized linear algebra routines in FLINT the matrices
+should be converted to `ArbMatrix` or `AcbMatrix` (depending on
+whether they are real or complex). Basic methods can then be used
+directly:
+
+``` @repl
+using Arblib, LinearAlgebra
+
+A = interval.(BigFloat, [1 2; 3 4])
+
+B = ArbMatrix(A)
+
+B * B
+
+inv(B)
+
+B \ B
+
+eigvals(AcbMatrix(B)) # Eigenvalues are only supported for AcbMatrix
+```
+
+Full documentation about supported methods can be found in the FLINT
+documentation for [arb_mat](https://flintlib.org/doc/arb_mat.html) and
+[acb_mat](https://flintlib.org/doc/acb_mat.html). Note that many of
+these methods are not wrapped in Arblib.jl, but most of them can be
+used through the [low level
+wrapper](https://kalmarek.github.io/Arblib.jl/stable/wrapper-methods/).
+
+## Special functions
+Arblib.jl wraps a large number of the special functions in
+[SpecialFunctions.jl](https://github.com/JuliaMath/SpecialFunctions.jl).
+
+``` @repl
+using Arblib, SpecialFunctions
+
+x = interval(BigFloat, 2)
+
+interval(besselj0(Arb(x))) # Convert to Arb and then back
+
+interval(gamma(Arb(x)))
+
+interval(zeta(Arb(x)))
+```
+
+Note that FLINT implements several special functions that are not in
+SpecialFunctions.jl, such as the [confluent hypergeometric
+functions](https://flintlib.org/doc/acb_hypgeom.html#confluent-hypergeometric-functions)
+and the [Lerch
+transcendent](https://flintlib.org/doc/acb_dirichlet.html#lerch-transcendent);
+they can be used through the [low-level
+wrapper](https://kalmarek.github.io/Arblib.jl/stable/wrapper-methods/):
+
+``` @repl
+using Arblib
+
+z = interval(BigFloat, 2 + 3im)
+s = interval(BigFloat, 2)
+a = interval(BigFloat, 4)
+
+interval(Arblib.dirichlet_lerch_phi!(Acb(), Acb(z), Acb(s), Acb(a)))
+```

ext/IntervalArithmeticArblibExt.jl

@@ -0,0 +1,153 @@
+module IntervalArithmeticArblibExt
+
+import IntervalArithmetic as IA
+import Arblib
+
+# Promotion rules
+
+# These are needed for ambiguity reasons
+Base.promote_rule(
+    ::Type{IA.Interval{T}},
+    ::Type{S},
+) where {T<:IA.NumTypes,S<:Arblib.ArfOrRef} = IA.Interval{promote_type(T, S)}
+Base.promote_rule(
+    ::Type{T},
+    ::Type{IA.Interval{S}},
+) where {T<:Arblib.ArfOrRef,S<:IA.NumTypes} = IA.Interval{promote_type(T, S)}
+
+Base.promote_rule(
+    ::Type{IA.Interval{T}},
+    ::Type{S},
+) where {T<:IA.NumTypes,S<:Arblib.ArbOrRef} =
+    throw(ArgumentError("promotion between Interval and Arb is purposely not supported"))
+Base.promote_rule(
+    ::Type{T},
+    ::Type{IA.Interval{S}},
+) where {T<:Arblib.ArbOrRef,S<:IA.NumTypes} =
+    throw(ArgumentError("promotion between Interval and Arb is purposely not supported"))
+
+# Construction of Arb from Interval
+# Note that this automatically gives support for construction of Acb.
+
+Arblib._precision(x::Union{IA.Interval,IA.BareInterval}) =
+    Arblib._precision(IA.inf(x), IA.sup(x))
+
+function Arblib.set!(
+    res::Arblib.ArbLike,
+    x::Union{IA.Interval,IA.BareInterval};
+    prec::Integer = precision(res),
+)
+    if IA.isempty_interval(x)
+        Arblib.indeterminate!(res)
+    else
+        Arblib.set!(res, (IA.inf(x), IA.sup(x)); prec)
+    end
+    return res
+end
+
+# Construction of Interval from Arb
+
+# We essentially treat Arb as Interval{BigFloat} for the purposes of
+# constructing Intervals. This means that we only have to change the
+# behavior of _inf and _sup since these are called early on in the
+# construction process to handle endpoints which are Intervals.
+
+IA._inf(x::Arblib.ArbOrRef) = BigFloat(Arblib.lbound(x))
+IA._sup(x::Arblib.ArbOrRef) = BigFloat(Arblib.ubound(x))
+
+# Make promote_numtype for Arb behave like promotion with BigFloat
+IA.promote_numtype(::Type{<:Arblib.ArbOrRef}, ::Type{<:Arblib.ArbOrRef}) = BigFloat
+IA.promote_numtype(::Type{<:Arblib.ArbOrRef}, ::Type{S}) where {S} =
+    IA.promote_numtype(BigFloat, S)
+IA.promote_numtype(::Type{<:Arblib.ArbOrRef}, ::Type{S}) where {S<:IA.NumTypes} =
+    IA.promote_numtype(BigFloat, S)
+IA.promote_numtype(::Type{T}, ::Type{<:Arblib.ArbOrRef}) where {T} =
+    IA.promote_numtype(T, BigFloat)
+IA.promote_numtype(::Type{T}, ::Type{<:Arblib.ArbOrRef}) where {T<:IA.NumTypes} =
+    IA.promote_numtype(T, BigFloat)
+
+# Construction of Complex{Interval} from Acb
+IA.numtype(::Type{<:Arblib.AcbOrRef}) = Arblib.Arb
+
+# Handle Acb in the same way as Complex.
+IA._interval_infsup(
+    ::Type{T},
+    a::Arblib.AcbOrRef,
+    b::Arblib.AcbOrRef,
+    d::IA.Decoration,
+) where {T<:IA.NumTypes} = complex(
+    IA._interval_infsup(T, real(a), real(b), d),
+    IA._interval_infsup(T, imag(a), imag(b), d),
+)
+IA._interval_infsup(
+    ::Type{T},
+    a::Arblib.AcbOrRef,
+    b,
+    d::IA.Decoration,
+) where {T<:IA.NumTypes} = complex(
+    IA._interval_infsup(T, real(a), real(b), d),
+    IA._interval_infsup(T, imag(a), imag(b), d),
+)
+IA._interval_infsup(
+    ::Type{T},
+    a,
+    b::Arblib.AcbOrRef,
+    d::IA.Decoration,
+) where {T<:IA.NumTypes} = complex(
+    IA._interval_infsup(T, real(a), real(b), d),
+    IA._interval_infsup(T, imag(a), imag(b), d),
+)
+
+# Handle ambiguities
+for S in [Union{IA.BareInterval,IA.Interval}, Complex]
+    @eval IA._interval_infsup(
+        ::Type{T},
+        a::Arblib.AcbOrRef,
+        b::$S,
+        d::IA.Decoration,
+    ) where {T<:IA.NumTypes} = complex(
+        IA._interval_infsup(T, real(a), real(b), d),
+        IA._interval_infsup(T, imag(a), imag(b), d),
+    )
+    @eval IA._interval_infsup(
+        ::Type{T},
+        a::$S,
+        b::Arblib.AcbOrRef,
+        d::IA.Decoration,
+    ) where {T<:IA.NumTypes} = complex(
+        IA._interval_infsup(T, real(a), real(b), d),
+        IA._interval_infsup(T, imag(a), imag(b), d),
+    )
+end
+
+# Make it so that implicit conversion from Arb to Interval doesn't set
+# the NG flag.
+
+Base.convert(::Type{IA.Interval{T}}, x::Arblib.ArbOrRef) where {T<:IA.NumTypes} =
+    IA.interval(T, x)
+
+function Base.convert(::Type{IA.Interval{T}}, x::Arblib.AcbOrRef) where {T<:IA.NumTypes}
+    isreal(x) || return throw(DomainError(x, "imaginary part must be zero"))
+    return convert(IA.Interval{T}, real(x))
+end
+
+Base.convert(::Type{Complex{IA.Interval{T}}}, x::Arblib.ArbOrRef) where {T<:IA.NumTypes} =
+    complex(IA.interval(T, x))
+
+# Handle ambiguities related to ExactReal
+
+Arblib.Mag(x::IA.ExactReal) = convert(Arblib.Mag, x)
+Arblib.Arf(x::IA.ExactReal) = convert(Arblib.Arf, x)
+Arblib.Arb(x::IA.ExactReal) = convert(Arblib.Arb, x)
+
+Base.promote_rule(::Type{T}, ::Type{IA.ExactReal{S}}) where {T<:Arblib.ArfOrRef,S<:Real} =
+    promote_type(T, S)
+Base.promote_rule(::Type{IA.ExactReal{T}}, ::Type{S}) where {T<:Real,S<:Arblib.ArfOrRef} =
+    promote_type(T, S)
+
+Base.promote_rule(::Type{T}, ::Type{IA.ExactReal{S}}) where {T<:Arblib.ArbOrRef,S<:Real} =
+    promote_type(T, S)
+Base.promote_rule(::Type{IA.ExactReal{T}}, ::Type{S}) where {T<:Real,S<:Arblib.ArbOrRef} =
+    promote_type(T, S)
+
+end

test/Project.toml

@@ -1,5 +1,6 @@
 [deps]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+Arblib = "fb37089c-8514-4489-9461-98f9c8763369"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"