Merge pull request #179 from kalmarek/float-interface-improvements

Float interface improvements - and some other things

Project.toml

@@ -1,7 +1,7 @@
 name = "Arblib"
 uuid = "fb37089c-8514-4489-9461-98f9c8763369"
 authors = ["Marek Kaluba <[email protected]>", "Sascha Timme <Sascha Timme <[email protected]>", "Joel Dahne <[email protected]>"]
-version = "1.0"
+version = "1.1"
 
 [deps]
 FLINT_jll = "e134572f-a0d5-539d-bddf-3cad8db41a82"

src/Arblib.jl

@@ -60,6 +60,7 @@ include("precision.jl")
 
 include("setters.jl")
 include("constructors.jl")
+include("conversion.jl")
 include("predicates.jl")
 include("show.jl")
 include("promotion.jl")

src/constructors.jl

@@ -66,22 +66,6 @@ function Acb(re::AbstractString, im::AbstractString; prec::Integer = DEFAULT_PRE
     return res
 end
 
-Base.Int(x::ArfOrRef; rnd::Union{arb_rnd,RoundingMode} = RoundNearest) =
-    is_int(x) ? get_si(x, rnd) : throw(InexactError(:Int64, Int64, x))
-
-Base.Float64(x::MagOrRef) = get(x)
-Base.Float64(x::ArfOrRef; rnd::Union{arb_rnd,RoundingMode} = RoundNearest) = get_d(x, rnd)
-Base.Float64(x::ArbOrRef) = Float64(midref(x))
-
-Base.ComplexF64(z::AcbOrRef) = Complex(Float64(realref(z)), Float64(imagref(z)))
-
-function Base.BigFloat(x::ArfOrRef)
-    y = BigFloat(; precision = precision(x))
-    get!(y, x)
-    return y
-end
-Base.BigFloat(x::ArbOrRef) = BigFloat(midref(x))
-
 Base.zero(::Union{Mag,Type{Mag}}) = Mag(UInt64(0))
 Base.one(::Union{Mag,Type{Mag}}) = Mag(UInt64(1))
 Base.zero(x::T) where {T<:Union{Arf,Arb,Acb}} = T(0, prec = precision(x))

src/conversion.jl

@@ -0,0 +1,102 @@
+# Conversion to float types in Base
+
+## Float64
+Base.Float64(x::MagOrRef, r::RoundingMode = RoundUp) = Float64(x, convert(arb_rnd, r))
+Base.Float64(x::ArfOrRef, r::RoundingMode) = Float64(x, convert(arb_rnd, r))
+Base.Float64(x::ArbOrRef, r::RoundingMode = RoundNearest) = Float64(x, convert(arb_rnd, r))
+
+function Base.Float64(x::MagOrRef, r::arb_rnd)
+    r == ArbRoundUp ||
+        throw(ArgumentError("only supports rounding up when converting Mag to Float64"))
+    return get(x)
+end
+Base.Float64(x::ArfOrRef, r::arb_rnd) = get_d(x, r)
+Base.Float64(x::ArbOrRef, r::arb_rnd) = Float64(midref(x), r)
+
+# Deprecated
+Base.Float64(x::ArfOrRef; rnd::arb_rnd = ArbRoundNearest) = Float64(x, rnd)
+
+## Float16 and Float32
+Base.Float16(x::MagOrRef, r::RoundingMode = RoundUp) = Float16(Float64(x, r), r)
+Base.Float16(x::MagOrRef, r::arb_rnd) = Float16(x, convert(RoundingMode, r))
+Base.Float32(x::MagOrRef, r::RoundingMode = RoundUp) = Float32(Float64(x, r), r)
+Base.Float32(x::MagOrRef, r::arb_rnd) = Float32(x, convert(RoundingMode, r))
+Base.Float16(x::Union{ArfOrRef,ArbOrRef}, r::RoundingMode = RoundNearest) =
+    Float16(Float64(x, r), r)
+Base.Float16(x::Union{ArfOrRef,ArbOrRef}, r::arb_rnd) = Float16(x, convert(RoundingMode, r))
+Base.Float32(x::Union{ArfOrRef,ArbOrRef}, r::RoundingMode = RoundNearest) =
+    Float32(Float64(x, r), r)
+Base.Float32(x::Union{ArfOrRef,ArbOrRef}, r::arb_rnd) = Float32(x, convert(RoundingMode, r))
+
+## BigFloat
+
+# Note that this uses different default values than the constructors
+# in Base. It always defaults to RoundNearest and uses the precision
+# given by the input argument. This means that it doesn't depend on
+# the global rounding mode and precision.
+
+Base.BigFloat(x::Union{ArfOrRef,ArbOrRef}, r::RoundingMode; precision = Base.precision(x)) =
+    BigFloat(x, convert(Base.MPFR.MPFRRoundingMode, r); precision)
+
+Base.BigFloat(x::Union{ArfOrRef,ArbOrRef}, r::arb_rnd; precision = Base.precision(x)) =
+    BigFloat(x, convert(RoundingMode, r); precision)
+
+function Base.BigFloat(
+    x::ArfOrRef,
+    r::Base.MPFR.MPFRRoundingMode = Base.MPFR.RoundNearest;
+    precision = Base.precision(x),
+)
+    y = BigFloat(; precision)
+    get!(y, x, r)
+    return y
+end
+Base.BigFloat(
+    x::ArbOrRef,
+    r::Base.MPFR.MPFRRoundingMode = Base.MPFR.RoundNearest;
+    precision = Base.precision(x),
+) = BigFloat(midref(x), r; precision)
+
+## Conversion to integers
+
+# IMPROVE: We currently don't support any rounding to integers. This
+# would maybe be nice to do.
+# IMPROVE: Is there any point in supporting conversion to other
+# integers?
+
+function Base.Int(x::ArfOrRef)
+    is_int(x) || throw(InexactError(:Int, Int, x))
+    typemin(Int) <= x <= typemax(Int) || throw(InexactError(:Int, Int, x))
+    return get_si(x, ArbRoundNearest)
+end
+
+Base.Int(x::ArbOrRef) = is_int(x) ? Int(midref(x)) : throw(InexactError(:Int, Int, x))
+
+Base.Int(x::AcbOrRef) =
+    is_int(x) ? Int(midref(realref(x))) : throw(InexactError(:Int, Int, x))
+
+function Base.BigInt(x::ArfOrRef)
+    is_int(x) || throw(InexactError(:BigInt, BigInt, x))
+    n = fmpz_struct()
+    ccall(
+        @libflint(arf_get_fmpz),
+        Cint,
+        (Ref{fmpz_struct}, Ref{arf_struct}, Ref{arb_rnd}),
+        n,
+        x,
+        ArbRoundNearest,
+    )
+    return BigInt(n)
+end
+
+Base.BigInt(x::ArbOrRef) =
+    is_int(x) ? BigInt(midref(x)) : throw(InexactError(:BigInt, BigInt, x))
+
+Base.BigInt(x::AcbOrRef) =
+    is_int(x) ? BigInt(midref(realref(x))) : throw(InexactError(:BigInt, BigInt, x))
+
+## Conversion to Complex
+
+# TODO: This currently allows construction of Complex{ArbRef}, which
+# we probably don't want.
+Base.Complex{T}(z::AcbOrRef) where {T} = Complex{T}(realref(z), imagref(z))
+Base.Complex(z::AcbOrRef) = Complex{Arb}(z)

src/elementary.jl

@@ -66,6 +66,17 @@ sqrtpos(x::ArbOrRef) = sqrtpos!(zero(x), x)
 sqrt1pm1(x::ArbOrRef) = sqrt1pm1!(zero(x), x)
 rsqrt(x::Union{ArbOrRef,AcbOrRef}) = rsqrt!(zero(x), x)
 
+function Base.log2(x::Union{ArbOrRef,AcbOrRef})
+    res = log(x)
+    log_2 = const_log2!(Arb(prec = precision(x)))
+    return div!(res, res, log_2)
+end
+function Base.log10(x::Union{ArbOrRef,AcbOrRef})
+    res = log(x)
+    log_10 = const_log10!(Arb(prec = precision(x)))
+    return div!(res, res, log_10)
+end
+
 Base.sinpi(x::Union{ArbOrRef,AcbOrRef}) = sin_pi!(zero(x), x)
 Base.cospi(x::Union{ArbOrRef,AcbOrRef}) = cos_pi!(zero(x), x)
 tanpi(x::Union{ArbOrRef,AcbOrRef}) = tan_pi!(zero(x), x)

src/float.jl

@@ -16,3 +16,34 @@ function Base.eps(T::Type{<:Union{ArfOrRef,ArbOrRef}})
     return eps!(res, res)
 end
 Base.eps(x::Union{ArfOrRef,ArbOrRef}) = eps!(zero(x), x)
+
+Base.typemin(::Type{<:MagOrRef}) = zero!(Mag())
+Base.typemin(x::Union{ArfOrRef,ArbOrRef}) = neg_inf!(zero(x))
+Base.typemin(T::Type{<:Union{ArfOrRef,ArbOrRef}}) = neg_inf!(zero(T))
+
+Base.typemax(::Type{<:MagOrRef}) = inf!(Mag())
+Base.typemax(x::Union{ArfOrRef,ArbOrRef}) = pos_inf!(zero(x))
+Base.typemax(T::Type{<:Union{ArfOrRef,ArbOrRef}}) = pos_inf!(zero(T))
+
+function Base.frexp(x::ArfOrRef)
+    m = zero(x)
+    e = fmpz_struct()
+    ccall(
+        @libflint(arf_frexp),
+        Nothing,
+        (Ref{arf_struct}, Ref{fmpz_struct}, Ref{arf_struct}),
+        m,
+        e,
+        x,
+    )
+
+    return m, BigInt(e)
+end
+
+function Base.frexp(x::ArbOrRef)
+    # Compute for midpoint and just scale radius
+    _, e = frexp(midref(x))
+    return ldexp(x, -e), e
+end
+
+Base.ldexp(x::Union{ArfOrRef,ArbOrRef}, n::Integer) = mul_2exp!(zero(x), x, n)

src/interval.jl

@@ -5,11 +5,12 @@ export radius,
     radius([T, ] x::ArbOrRef)
 
 Returns the radius of `x` as a `Mag`. If `T` is given convert to this
-type, supports `Mag`, `Arf` and `Arb`.
+type, supports `Mag`, `Arf`, `Arb` and `Float64`.
 """
 radius(::Type{Mag}, x::ArbOrRef) = Mag(radref(x))
 radius(::Type{Arf}, x::ArbOrRef) = Arf(radref(x), prec = precision(x))
 radius(::Type{Arb}, x::ArbOrRef) = Arb(radref(x), prec = precision(x))
+radius(::Type{Float64}, x::ArbOrRef) = Float64(radref(x))
 radius(x::ArbOrRef) = radius(Mag, x)
 
 """
@@ -133,11 +134,12 @@ 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`.
+to this type, supports `Arf` and `Arb`.
 
 See also [`setball`](@ref) and [`getinterval`](@ref).
 """
 getball(x::ArbOrRef) = (Arf(midref(x)), Mag(radref(x)))
+getball(::Type{Arf}, x::ArbOrRef) = (Arf(midref(x)), Arf(radref(x), prec = precision(x)))
 getball(::Type{Arb}, x::ArbOrRef) = (Arb(midref(x)), Arb(radref(x), prec = precision(x)))
 
 """

src/precision.jl

@@ -17,16 +17,32 @@ else
     # overloading that we automatically support giving base argument.
 
     # Types
-    Base._precision(::Type{<:Union{ArbTypes,ArbStructTypes,Ptr{<:ArbStructTypes}}}) =
-        DEFAULT_PRECISION[]
-    # Types not storing their precision
-    Base._precision(::Union{ArbStructTypes,Ptr{<:ArbStructTypes}}) = DEFAULT_PRECISION[]
-    # Types storing their precision
-    Base._precision(x::ArbTypes) = x.prec
-    # Mag doesn't store a precision
-    Base._precision(::MagOrRef) = DEFAULT_PRECISION[]
-    # ArbSeries and AcbSeries don't store their precision directly
-    Base._precision(x::Union{ArbSeries,AcbSeries}) = Base._precision(x.poly)
+    if VERSION < v"1.11.0-DEV"
+        Base._precision(::Type{<:Union{ArbTypes,ArbStructTypes,Ptr{<:ArbStructTypes}}}) =
+            DEFAULT_PRECISION[]
+        # Types not storing their precision
+        Base._precision(::Union{ArbStructTypes,Ptr{<:ArbStructTypes}}) = DEFAULT_PRECISION[]
+        # Types storing their precision
+        Base._precision(x::ArbTypes) = x.prec
+        # Mag doesn't store a precision
+        Base._precision(::MagOrRef) = DEFAULT_PRECISION[]
+        # ArbSeries and AcbSeries don't store their precision directly
+        Base._precision(x::Union{ArbSeries,AcbSeries}) = Base._precision(x.poly)
+    else
+        Base._precision_with_base_2(
+            ::Type{<:Union{ArbTypes,ArbStructTypes,Ptr{<:ArbStructTypes}}},
+        ) = DEFAULT_PRECISION[]
+        # Types not storing their precision
+        Base._precision_with_base_2(::Union{ArbStructTypes,Ptr{<:ArbStructTypes}}) =
+            DEFAULT_PRECISION[]
+        # Types storing their precision
+        Base._precision_with_base_2(x::ArbTypes) = x.prec
+        # Mag doesn't store a precision
+        Base._precision_with_base_2(::MagOrRef) = DEFAULT_PRECISION[]
+        # ArbSeries and AcbSeries don't store their precision directly
+        Base._precision_with_base_2(x::Union{ArbSeries,AcbSeries}) =
+            Base._precision_with_base_2(x.poly)
+    end
 
     # Base.precision only allows AbstractFloat, we want to be able to use
     # all ArbLib types.

src/rounding.jl

@@ -28,3 +28,21 @@ Base.convert(::Type{arb_rnd}, ::RoundingMode{:Down}) = ArbRoundDown
 Base.convert(::Type{arb_rnd}, ::RoundingMode{:Up}) = ArbRoundUp
 Base.convert(::Type{arb_rnd}, ::RoundingMode{:Nearest}) = ArbRoundNearest
 Base.convert(::Type{arb_rnd}, ::RoundingMode{:Exact}) = ArbRoundExact
+
+function Base.convert(::Type{RoundingMode}, r::arb_rnd)
+    if r == ArbRoundToZero
+        return RoundToZero
+    elseif r == ArbRoundFromZero
+        return RoundFromZero
+    elseif r == ArbRoundDown
+        return RoundDown
+    elseif r == ArbRoundUp
+        return RoundUp
+    elseif r == ArbRoundNearest
+        return RoundNearest
+    elseif r == ArbRoundExact
+        return RoundingMode{:Exact}()
+    else
+        throw(ArgumentError("invalid Arb rounding mode code: $r"))
+    end
+end

src/setters.jl

@@ -69,9 +69,9 @@ end
 
 function set!(res::ArbLike, ::Irrational{:φ}; prec::Integer = precision(res))
     set!(res, 5)
-    sqrt!(res, res, prec)
-    add!(res, res, 1, prec)
-    return div!(res, res, 2, prec)
+    sqrt!(res, res; prec)
+    add!(res, res, 1; prec)
+    return mul_2exp!(res, res, -1)
 end
 
 function set!(

src/show.jl

@@ -112,6 +112,8 @@ function load_string!(x::Union{MagLike,ArfLike,ArbLike}, str::AbstractString)
     return x
 end
 
+load_string(T::Type{<:Union{Mag,Arf,Arb}}, str::AbstractString) = load_string!(zero(T), str)
+
 function dump_string(x::Union{MagLike,ArfLike,ArbLike})
     char_ptr = dump(x)
     str = unsafe_string(char_ptr)

test/constructors.jl

@@ -128,28 +128,6 @@
         @test isequal(Acb((1, 2), (3, 4)), Acb(Arb((1, 2)), Arb((3, 4))))
     end
 
-    @testset "Conversion" begin
-        @test Int(Arf(2.0)) isa Int
-        @test Int(Arf(2.0)) == 2
-        @test_throws InexactError Int(Arf(2.5))
-
-        @test Float64(Mag(2)) isa Float64
-        @test Float64(Mag(2)) == 2.0
-        @test Float64(Arf(2)) isa Float64
-        @test Float64(Arf(2)) == 2.0
-        @test Float64(Arb(2)) isa Float64
-        @test Float64(Arb(2)) == 2.0
-
-        @test ComplexF64(Acb(2 + 3im)) isa ComplexF64
-        @test ComplexF64(Acb(2 + 3im)) == 2.0 + 3.0im
-
-        @test Arblib.get_si(Arf(0.5)) == 0
-        @test Arblib.get_si(Arf(0.5); rnd = Arblib.ArbRoundFromZero) == 1
-
-        @test BigFloat(Arf(2.5)) == 2.5
-        @test BigFloat(Arb(2.5)) == 2.5
-    end
-
     @testset "zeros/ones" begin
         for T in [Arf, Arb, Acb]
             @test zeros(T, 2) == [zero(T), zero(T)]