Some new methods

src/constructors.jl

@@ -1,3 +1,5 @@
+export setball
+
 # Mag
 Mag(x) = set!(Mag(), x)
 Mag(x::Union{MagRef,ArfRef}) = Mag(mag_struct(cstruct(x)))
@@ -22,6 +24,24 @@ function Arb(str::AbstractString; prec::Integer = DEFAULT_PRECISION[])
     return res
 end
 
+"""
+    setball(::Type{Arb}, m, r; prec = _precision(m))
+
+Returns an `Arb` with the midpoint and radius set to `m` and `r`
+respectively.
+
+Note that the `m` is converted to an `Arf` and therefore rounded. So
+for example `setball(1 // 3, 0)` will not contain ``1 / 3``.
+
+See also [`getball`](@ref) and [`add_error`](@ref).
+"""
+function setball(::Type{Arb}, m, r; prec = _precision(m))
+    res = Arb(; prec)
+    Arblib.set!(Arblib.midref(res), m)
+    Arblib.set!(Arblib.radref(res), r)
+    return res
+end
+
 ## Acb
 Acb(
     x::Union{Real,arb_struct,arf_struct,Tuple{<:Real,<:Real}};

src/interval.jl

@@ -1,4 +1,5 @@
-export radius, midpoint, lbound, ubound, abs_lbound, abs_ubound, getinterval, getball
+export radius,
+    midpoint, lbound, ubound, abs_lbound, abs_ubound, getinterval, getball, add_error #
 
 """
     radius([T, ] x::ArbOrRef)
@@ -95,6 +96,8 @@ the ball `x`, both of them are of type `Arf`. If `T` is given convert
 to this type, supports `Arf`, `BigFloat` and `Arb`.
 
 If `x` contains `NaN` both `l` and `u` will be `NaN`.
+
+See also [`getball`](@ref).
 """
 function getinterval(::Type{Arf}, x::ArbOrRef)
     l, u = Arf(prec = precision(x)), Arf(prec = precision(x))
@@ -119,6 +122,8 @@ getinterval(x::ArbOrRef) = getinterval(Arf, x)
 Returns a tuple `(m::Arf, r::Mag)` where `m` is the midpoint of the
 ball and `r` is the radius. If `T` is given convert both `m` and `r`
 to this type, supports `Arb`.
+
+See also [`setball`](@ref) and [`getinterval`](@ref). #
 """
 getball(x::ArbOrRef) = (Arf(midref(x)), Mag(radref(x)))
 getball(::Type{Arb}, x::ArbOrRef) = (Arb(midref(x)), Arb(radref(x), prec = precision(x)))
@@ -164,3 +169,17 @@ end
 # TODO: Could be optimized, both for performance and enclosure
 Base.intersect(x::ArbOrRef, y::ArbOrRef, z::ArbOrRef, xs...) =
     foldl(intersect, xs, init = intersect(intersect(x, y), z))
+
+"""
+    add_error(x, err)
+
+Returns a copy of `x` with the absolute value of `err` added to the radius.
+
+For complex `x` it adds the error to both the real and imaginary
+parts. For matrices it adds it elementwise.
+
+See also [`set_ball`](@ref).
+"""
+add_error(x::Union{ArbOrRef,AcbOrRef}, err::Union{MagOrRef,ArfOrRef,ArbOrRef}) =
+    add_error!(copy(x), err)
+add_error(x::Union{ArbMatrixLike,AcbMatrixLike}, err::MagOrRef) = add_error!(copy(x), err)

src/special-functions.jl

@@ -12,7 +12,7 @@ SpecialFunctions.loggamma(x::Series) = lgamma_series!(zero(x), x, length(x))
 # Not implemented by Arb
 
 #SpecialFunctions.logfactorial(x)
-# Only relevant for Arblib
+# Not relevant for Arblib, only implemented for x::Integer
 
 SpecialFunctions.digamma(x::Union{ArbOrRef,AcbOrRef}) = digamma!(zero(x), x)
 SpecialFunctions.digamma(x::Series) = digamma_series!(zero(x), x, length(x))
@@ -152,8 +152,7 @@ SpecialFunctions.erf(x::Series) = hypgeom_erf_series!(zero(x), x, length(x))
 SpecialFunctions.erfc(x::Union{ArbOrRef,AcbOrRef}) = hypgeom_erfc!(zero(x), x)
 SpecialFunctions.erfc(x::Series) = hypgeom_erfc_series!(zero(x), x, length(x))
 
-#SpecialFunctions.erfcinv(x)
-# Not implemented by Arb
+SpecialFunctions.erfcinv(x::ArbOrRef) = hypgeom_erfcinv!(zero(x), x)
 
 #SpecialFunctions.erfcx(x)
 # Not implemented by Arb
@@ -164,12 +163,10 @@ SpecialFunctions.erfc(x::Series) = hypgeom_erfc_series!(zero(x), x, length(x))
 #SpecialFunctions.logerfcx(x)
 # Not implemented by Arb
 
-#SpecialFunctions.erfi(x)
 SpecialFunctions.erfi(x::Union{ArbOrRef,AcbOrRef}) = hypgeom_erfi!(zero(x), x)
 SpecialFunctions.erfi(x::Series) = hypgeom_erfi_series!(zero(x), x, length(x))
 
-#SpecialFunctions.erfinv(x)
-# Not implemented by Arb
+SpecialFunctions.erfinv(x::ArbOrRef) = hypgeom_erfinv!(zero(x), x)
 
 #SpecialFunctions.dawson(x)
 # Not implemented by Arb

test/constructors.jl

@@ -59,6 +59,20 @@
         @test precision(Arb(MathConstants.γ, prec = 80)) == 80
         @test precision(Arb(MathConstants.catalan, prec = 80)) == 80
         @test precision(Arb(MathConstants.φ, prec = 80)) == 80
+
+        # setball
+        @test isone(Arblib.setball(Arb, Arf(1), Mag(0)))
+        @test isone(Arblib.setball(Arb, 1, 0))
+        @test isone(Arblib.setball(Arb, 1.0, 0.0))
+        @test isone(Arblib.setball(Arb, big(1.0), 0))
+        @test isone(Arblib.setball(Arb, 1 // 1, 0))
+
+        @test getinterval(Arblib.setball(Arb, 2, 1)) == (1, 3)
+        @test contains(Arblib.setball(Arb, 0.5, 0.5), 0)
+        @test contains(Arblib.setball(Arb, 0.5, 0.5), 1)
+
+        @test precision(Arblib.setball(Arb, Arf(prec = 80), 0)) == 80
+        @test precision(Arblib.setball(Arb, Arf(prec = 90), 0, prec = 80)) == 80
     end
 
     @testset "Acb" begin

test/interval.jl

@@ -196,4 +196,35 @@
         @test_throws ArgumentError intersect(Arb(1), Arb(2))
         @test_throws ArgumentError intersect([xs; Arb(2)]...)
     end
+
+    @testset "add_error" begin
+        x = Arblib.setball(Arb, 0, 1)
+
+        @test radius(Arblib.add_error(x, Mag(1))) == Mag(2)
+        @test radius(Arblib.add_error(x, Arf(2))) == Mag(1) + Mag(Arf(2))
+        @test radius(Arblib.add_error(x, Arb(3))) == Mag(1) + Mag(Arf(3))
+
+        y = Arblib.add_error(x, Mag(1))
+        Arblib.zero!(y)
+        @test !isequal(x, y)
+
+        x = Acb(Arblib.setball(Arb, 0, 1), Arblib.setball(Arb, 2, 3))
+
+        y = Arblib.add_error(x, Mag(1))
+        @test radius(real(y)) == Mag(2)
+        @test radius(imag(y)) == Mag(1) + Mag(3)
+        y = Arblib.add_error(x, Arf(2))
+        @test radius(real(y)) == Mag(1) + Mag(Arf(2))
+        @test radius(imag(y)) == Mag(3) + Mag(Arf(2))
+        y = Arblib.add_error(x, Arb(3))
+        @test radius(real(y)) == Mag(1) + Mag(Arf(3))
+        @test radius(imag(y)) == Mag(3) + Mag(Arf(3))
+
+        @test precision(Arblib.add_error(Arb(prec = 80), Arb(prec = 90))) == 80
+        @test precision(Arblib.add_error(Acb(prec = 80), Arb(prec = 90))) == 80
+
+        @test all(isone, radius.(Arblib.add_error(ArbMatrix(2, 2), Mag(1))))
+        @test all(isone, radius.(real.(Arblib.add_error(AcbMatrix(2, 2), Mag(1)))))
+        @test all(isone, radius.(imag.(Arblib.add_error(AcbMatrix(2, 2), Mag(1)))))
+    end
 end

test/special-functions.jl

@@ -108,12 +108,18 @@
         @test erfc(ArbSeries([2, 1]))[0] ≈ erfc(2)
         @test erfc(AcbSeries([2 + 2im, 1]))[0] ≈ erfc(2 + 2im)
 
+        @test erfcinv(Arb(0.1)) ≈ erfcinv(0.1)
+        @test erfcinv(Arb(0.2)) ≈ erfcinv(0.2)
+
         @test erfi(Arb(2)) ≈ erfi(2)
         @test erfi(Arb(3)) ≈ erfi(3)
         @test erfi(Acb(2 + 2im)) ≈ erfi(2 + 2im)
         @test erfi(Acb(3 + 3im)) ≈ erfi(3 + 3im)
         @test erfi(ArbSeries([2, 1]))[0] ≈ erfi(2)
         @test erfi(AcbSeries([2 + 2im, 1]))[0] ≈ erfi(2 + 2im)
+
+        @test erfinv(Arb(0.1)) ≈ erfinv(0.1)
+        @test erfinv(Arb(0.2)) ≈ erfinv(0.2)
     end
 
     @testset "Airy Functions" begin