Implement rounding related methods

src/Arblib.jl

@@ -49,7 +49,7 @@ end
 
 include("arb_types.jl")
 include("fmpz.jl")
-include("rounding.jl")
+include("rounding_types.jl")
 include("types.jl")
 include("hash.jl")
 include("serialize.jl")
@@ -73,6 +73,7 @@ include("rand.jl")
 include("float.jl")
 include("interval.jl")
 include("multi-argument.jl")
+include("rounding.jl")
 
 include("ref.jl")
 include("vector.jl")

src/conversion.jl

@@ -99,6 +99,19 @@ Base.BigInt(x::ArbOrRef) =
 Base.BigInt(x::AcbOrRef) =
     is_int(x) ? BigInt(midref(realref(x))) : throw(InexactError(:BigInt, BigInt, x))
 
+function (::Type{T})(x::Union{ArfOrRef,AcfOrRef,ArbOrRef,AcbOrRef}) where {T<:Integer}
+    if typemax(T) <= typemax(Int)
+        convert(T, Int(x))
+    else
+        convert(T, BigInt(x))
+    end
+end
+(::Type{Integer})(x::Union{ArfOrRef,AcfOrRef,ArbOrRef,AcbOrRef}) = BigInt(x)
+
+# Ambiguity
+Base.Bool(x::Union{ArfOrRef,AcfOrRef,ArbOrRef,AcbOrRef}) =
+    iszero(x) ? false : isone(x) ? true : throw(InexactError(:Bool, Bool, x))
+
 ## Conversion to Complex
 # TODO: This currently allows construction of Complex{ArfRef} and
 # Complex{ArbRef}, which we probably don't want.

src/rounding.jl

@@ -1,44 +1,72 @@
-#=
-    Specifies the rounding mode for the result of an approximate operation.
-ARF_RND_NEAR   (4)
-    Round to the nearest representable number, rounding to even if there is a tie.
-ARF_RND_FLOOR  (2)
-    Round to the nearest representable number in the direction towards minus infinity.
-ARF_RND_CEIL   (3)
-    Round to the nearest representable number in the direction towards plus infinity.
-ARF_RND_DOWN   (0)
-    Round to the nearest representable number in the direction towards zero.
-ARF_RND_UP     (1)
-    Round to the nearest representable number in the direction away from zero.
-=#
-
-@enum(
-    arb_rnd::Cint,
-    ArbRoundToZero, # ARF_RND_DOWN
-    ArbRoundFromZero, # ARF_RND_UP
-    ArbRoundDown, # ARF_RND_FLOOR
-    ArbRoundUp, # ARF_RND_CEIL
-    ArbRoundNearest, # ARF_RND_NEAR
-)
+# IMPROVE: In principle we could try to implement the currently
+# unsupported rounding modes. One would however have to be a bit
+# careful to get them correct, so for now we just throw an error.
+
+# TODO: This is implemented in
+# https://github.com/flintlib/flint/pull/2294 and we can add it once
+# that version is released.
+_round!(res::ArfOrRef, x::ArfOrRef, r::typeof(RoundNearest)) =
+    throw(ArgumentError("rounding mode $r not supported for Arf"))
+_round!(res::ArfOrRef, x::ArfOrRef, r::typeof(RoundNearestTiesAway)) =
+    throw(ArgumentError("rounding mode $r not supported for Arf"))
+_round!(res::ArfOrRef, x::ArfOrRef, r::typeof(RoundNearestTiesUp)) =
+    throw(ArgumentError("rounding mode $r not supported for Arf"))
+_round!(res::ArfOrRef, x::ArfOrRef, r::typeof(RoundToZero)) =
+    signbit(x) ? _round!(res, x, RoundUp) : _round!(res, x, RoundDown)
+_round!(res::ArfOrRef, x::ArfOrRef, r::typeof(RoundFromZero)) =
+    signbit(x) ? _round!(res, x, RoundDown) : _round!(res, x, RoundUp)
+_round!(res::ArfOrRef, x::ArfOrRef, ::typeof(RoundUp)) = ceil!(res, x)
+_round!(res::ArfOrRef, x::ArfOrRef, ::typeof(RoundDown)) = floor!(res, x)
+
+# TODO: There is currently a bug in arb_nint, see
+# https://github.com/flintlib/flint/issues/2293. The implementation
+# here is a correct implementation of the version in Flint. Once the
+# fixed version in Flint is released we should however switch to that.
+function _round!(res::ArbOrRef, x::ArbOrRef, r::typeof(RoundNearest))
+    _round!(res, x + 0.5, RoundDown)
 
-Base.convert(::Type{arb_rnd}, ::RoundingMode{:ToZero}) = ArbRoundToZero
-Base.convert(::Type{arb_rnd}, ::RoundingMode{:FromZero}) = ArbRoundFromZero
-Base.convert(::Type{arb_rnd}, ::RoundingMode{:Down}) = ArbRoundDown
-Base.convert(::Type{arb_rnd}, ::RoundingMode{:Up}) = ArbRoundUp
-Base.convert(::Type{arb_rnd}, ::RoundingMode{:Nearest}) = ArbRoundNearest
-
-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
-    else
-        throw(ArgumentError("invalid Arb rounding mode code: $r"))
+    u = (2x - 1) / 4
+    if isinteger(u)
+        sub!(res, res, 1)
+    elseif contains_int(u)
+        sub!(res, res, unit_interval!(zero(res)))
     end
+
+    return res
+end
+
+_round!(res::ArbOrRef, x::ArbOrRef, r::typeof(RoundNearestTiesAway)) =
+    throw(ArgumentError("rounding mode $r not supported for Arb"))
+_round!(res::ArbOrRef, x::ArbOrRef, r::typeof(RoundNearestTiesUp)) =
+    throw(ArgumentError("rounding mode $r not supported for Arb"))
+_round!(res::ArbOrRef, x::ArbOrRef, r::typeof(RoundToZero)) = trunc!(res, x)
+_round!(res::ArbOrRef, x::ArbOrRef, r::typeof(RoundFromZero)) =
+    throw(ArgumentError("rounding mode $r not supported for Arb"))
+_round!(res::ArbOrRef, x::ArbOrRef, ::typeof(RoundUp)) = ceil!(res, x)
+_round!(res::ArbOrRef, x::ArbOrRef, ::typeof(RoundDown)) = floor!(res, x)
+
+Base.round(x::Union{ArfOrRef,ArbOrRef}, r::RoundingMode) = _round!(zero(x), x, r)
+# Handle ambiguities
+Base.round(x::Union{ArfOrRef,ArbOrRef}, r::typeof(RoundNearestTiesAway)) =
+    _round!(zero(x), x, r)
+Base.round(x::Union{ArfOrRef,ArbOrRef}, r::typeof(RoundNearestTiesUp)) =
+    _round!(zero(x), x, r)
+Base.round(x::Union{ArfOrRef,ArbOrRef}, r::typeof(RoundFromZero)) = _round!(zero(x), x, r)
+
+function Base.round(
+    z::Union{AcfOrRef,AcbOrRef},
+    rr::RoundingMode = RoundNearest,
+    ri::RoundingMode = RoundNearest,
+)
+    res = zero(z)
+    _round!(realref(res), realref(z), rr)
+    _round!(imagref(res), imagref(z), ri)
+    return res
+end
+
+# There is no Arf version of this since getting the correct rounding
+# for that would require some more work.
+function Base.div(x::ArbOrRef, y::ArbOrRef, r::RoundingMode)
+    res = x / y
+    return _round!(res, res, r)
 end

src/rounding_types.jl

@@ -0,0 +1,44 @@
+#=
+    Specifies the rounding mode for the result of an approximate operation.
+ARF_RND_NEAR   (4)
+    Round to the nearest representable number, rounding to even if there is a tie.
+ARF_RND_FLOOR  (2)
+    Round to the nearest representable number in the direction towards minus infinity.
+ARF_RND_CEIL   (3)
+    Round to the nearest representable number in the direction towards plus infinity.
+ARF_RND_DOWN   (0)
+    Round to the nearest representable number in the direction towards zero.
+ARF_RND_UP     (1)
+    Round to the nearest representable number in the direction away from zero.
+=#
+
+@enum(
+    arb_rnd::Cint,
+    ArbRoundToZero, # ARF_RND_DOWN
+    ArbRoundFromZero, # ARF_RND_UP
+    ArbRoundDown, # ARF_RND_FLOOR
+    ArbRoundUp, # ARF_RND_CEIL
+    ArbRoundNearest, # ARF_RND_NEAR
+)
+
+Base.convert(::Type{arb_rnd}, ::RoundingMode{:ToZero}) = ArbRoundToZero
+Base.convert(::Type{arb_rnd}, ::RoundingMode{:FromZero}) = ArbRoundFromZero
+Base.convert(::Type{arb_rnd}, ::RoundingMode{:Down}) = ArbRoundDown
+Base.convert(::Type{arb_rnd}, ::RoundingMode{:Up}) = ArbRoundUp
+Base.convert(::Type{arb_rnd}, ::RoundingMode{:Nearest}) = ArbRoundNearest
+
+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
+    else
+        throw(ArgumentError("invalid Arb rounding mode code: $r"))
+    end
+end

test/conversion.jl

@@ -222,6 +222,53 @@
             @test_throws InexactError BigInt(Acb(setball(Arb, 1, 1)))
             @test_throws InexactError BigInt(Acb(2, 1))
         end
+
+        @testset "Other" begin
+            @test !Bool(Arf(0))
+            @test Bool(Arf(1))
+            @test_throws InexactError Bool(Arf(2))
+            @test Integer(Arf(1)) isa BigInt
+            @test Integer(Arf(1)) == 1
+            @test UInt(Arf(1)) isa UInt
+            @test UInt(Arf(1)) == 1
+            @test Int16(Arf(2)) isa Int16
+            @test Int16(Arf(2)) == 2
+            @test Int128(Arf(typemax(Int)) + 1) == Int128(typemax(Int)) + 1
+
+            @test !Bool(Acf(0))
+            @test Bool(Acf(1))
+            @test_throws InexactError Bool(Acf(1, 1))
+            @test Integer(Acf(1)) isa BigInt
+            @test Integer(Acf(1)) == 1
+            @test UInt(Acf(1)) isa UInt
+            @test UInt(Acf(1)) == 1
+            @test Int16(Acf(2)) isa Int16
+            @test Int16(Acf(2)) == 2
+            @test Int128(Acf(typemax(Int)) + 1) == Int128(typemax(Int)) + 1
+
+
+            @test !Bool(Arb(0))
+            @test Bool(Arb(1))
+            @test_throws InexactError Bool(Arb((-1, 1)))
+            @test Integer(Arb(1)) isa BigInt
+            @test Integer(Arb(1)) == 1
+            @test UInt(Arb(1)) isa UInt
+            @test UInt(Arb(1)) == 1
+            @test Int16(Arb(2)) isa Int16
+            @test Int16(Arb(2)) == 2
+            @test Int128(Arb(typemax(Int)) + 1) == Int128(typemax(Int)) + 1
+
+            @test !Bool(Acb(0))
+            @test Bool(Acb(1))
+            @test_throws InexactError Bool(Acb((-1, 1)))
+            @test Integer(Acb(1)) isa BigInt
+            @test Integer(Acb(1)) == 1
+            @test UInt(Acb(1)) isa UInt
+            @test UInt(Acb(1)) == 1
+            @test Int16(Acb(2)) isa Int16
+            @test Int16(Acb(2)) == 2
+            @test Int128(Acb(typemax(Int)) + 1) == Int128(typemax(Int)) + 1
+        end
     end
 
     @testset "Complex" begin

test/rounding.jl

@@ -1,13 +1,105 @@
 @testset "rounding" begin
-    @test convert(Arblib.arb_rnd, RoundToZero) == Arblib.ArbRoundToZero
-    @test convert(Arblib.arb_rnd, RoundFromZero) == Arblib.ArbRoundFromZero
-    @test convert(Arblib.arb_rnd, RoundDown) == Arblib.ArbRoundDown
-    @test convert(Arblib.arb_rnd, RoundUp) == Arblib.ArbRoundUp
-    @test convert(Arblib.arb_rnd, RoundNearest) == Arblib.ArbRoundNearest
-
-    @test convert(RoundingMode, Arblib.ArbRoundToZero) == RoundToZero
-    @test convert(RoundingMode, Arblib.ArbRoundFromZero) == RoundFromZero
-    @test convert(RoundingMode, Arblib.ArbRoundDown) == RoundDown
-    @test convert(RoundingMode, Arblib.ArbRoundUp) == RoundUp
-    @test convert(RoundingMode, Arblib.ArbRoundNearest) == RoundNearest
+    @testset "round" begin
+        xs = range(-3, 3, length = 25)
+        # RoundNearest
+        @test_throws ArgumentError round(Arf(0))
+        @test_throws ArgumentError round(Arf(0), RoundNearest)
+        @test round.(xs) == round.(Arb.(xs)) == round.(Arb.(xs), RoundNearest)
+        # RoundNearestTiesAway
+        @test_throws ArgumentError round(Arf(0), RoundNearestTiesAway)
+        @test_throws ArgumentError round(Arb(0), RoundNearestTiesAway)
+        # RoundNearestTiesUp
+        @test_throws ArgumentError round(Arf(0), RoundNearestTiesUp)
+        @test_throws ArgumentError round(Arb(0), RoundNearestTiesUp)
+        # RoundToZero
+        @test round.(xs, RoundToZero) ==
+              round.(Arf.(xs), RoundToZero) ==
+              round.(Arb.(xs), RoundToZero)
+        # RoundFromZero
+        @test round.(xs, RoundFromZero) == round.(Arf.(xs), RoundFromZero)
+        @test_throws ArgumentError round(Arb(0), RoundFromZero)
+        # RoundUp
+        @test round.(xs, RoundUp) == round.(Arf.(xs), RoundUp) == round.(Arb.(xs), RoundUp)
+        # RoundDown
+        @test round.(xs, RoundDown) ==
+              round.(Arf.(xs), RoundDown) ==
+              round.(Arb.(xs), RoundDown)
+
+        # Make sure it works correctly for intervals
+        @test Arblib.contains(round(Arb((0.4, 0.6))), 0)
+        @test Arblib.contains(round(Arb((0.4, 0.6))), 1) # Bug in Flint
+
+        @test Arblib.contains(round(Arb((0.9, 1.1)), RoundToZero), 0)
+        @test Arblib.contains(round(Arb((0.9, 1.1)), RoundToZero), 1)
+
+        @test Arblib.contains(round(Arb((0.9, 1.1)), RoundDown), 0)
+        @test Arblib.contains(round(Arb((0.9, 1.1)), RoundDown), 1)
+
+        @test Arblib.contains(round(Arb((0.9, 1.1)), RoundUp), 1)
+        @test Arblib.contains(round(Arb((0.9, 1.1)), RoundUp), 2)
+
+        # Complex
+        ys = transpose(xs)
+        for rr in (RoundToZero, RoundFromZero, RoundUp, RoundDown)
+            for ri in (RoundToZero, RoundFromZero, RoundUp, RoundDown)
+                @test round.(complex.(xs, ys), rr, ri) == round.(Acf.(xs, ys), rr, ri)
+            end
+        end
+        for rr in (RoundNearest, RoundToZero, RoundUp, RoundDown)
+            for ri in (RoundNearest, RoundToZero, RoundUp, RoundDown)
+                @test round.(complex.(xs, ys), rr, ri) == round.(Acb.(xs, ys), rr, ri)
+            end
+        end
+        @test_throws ArgumentError round(Acf(1, 1), RoundNearest, RoundDown)
+        @test_throws ArgumentError round(Acf(1, 1), RoundDown, RoundNearest)
+        @test_throws ArgumentError round(Acb(1, 1), RoundNearestTiesAway, RoundDown)
+        @test_throws ArgumentError round(Acb(1, 1), RoundDown, RoundNearestTiesAway)
+    end
+
+    @testset "div" begin
+        xs = range(-10, 10, length = 41)
+        ys = transpose(xs)
+
+        # Checks if x and y are equal, treating zeros as equal
+        # independent of sign and NaN as equal independent of
+        # representation.
+        _isequal_or_both_zero_or_nan(x, y) =
+            isequal(x, y) || (iszero(x) && iszero(y)) || (isnan(x) && isnan(y))
+
+        # RoundNearest
+        @test all(_isequal_or_both_zero_or_nan.(div.(xs, ys), div.(Arb.(xs), Arb.(ys))))
+        # RoundNearestTiesAway
+        @test_throws ArgumentError div(Arb(1), Arb(1), RoundNearestTiesAway)
+        # RoundNearestTiesUp
+        @test_throws ArgumentError div(Arf(1), Arb(1), RoundNearestTiesUp)
+        @test_throws ArgumentError div(Arb(1), Arb(1), RoundNearestTiesUp)
+        # RoundToZero
+        @test all(
+            _isequal_or_both_zero_or_nan.(
+                div.(xs, ys, RoundToZero),
+                div.(Arb.(xs), Arb.(ys), RoundToZero),
+            ),
+        )
+        # RoundFromZero
+        @test_throws ArgumentError div(Arb(1), Arb(1), RoundFromZero)
+        # RoundUp
+        @test all(
+            _isequal_or_both_zero_or_nan.(
+                div.(xs, ys, RoundUp),
+                div.(Arb.(xs), Arb.(ys), RoundUp),
+            ),
+        )
+        # RoundDown
+        @test all(
+            _isequal_or_both_zero_or_nan.(
+                div.(xs, ys, RoundDown),
+                div.(Arb.(xs), Arb.(ys), RoundDown),
+            ),
+        )
+
+        # Test that creating ranges work. This uses div internally and
+        # was the initial motivation for implementing div
+        @test (Arb(1):Arb(10)) == (Arb.(1:10))
+        @test_throws InexactError Arb(π):(Arb(π)+1)
+    end
 end

test/rounding_types.jl

@@ -0,0 +1,13 @@
+@testset "rounding_types" begin
+    @test convert(Arblib.arb_rnd, RoundToZero) == Arblib.ArbRoundToZero
+    @test convert(Arblib.arb_rnd, RoundFromZero) == Arblib.ArbRoundFromZero
+    @test convert(Arblib.arb_rnd, RoundDown) == Arblib.ArbRoundDown
+    @test convert(Arblib.arb_rnd, RoundUp) == Arblib.ArbRoundUp
+    @test convert(Arblib.arb_rnd, RoundNearest) == Arblib.ArbRoundNearest
+
+    @test convert(RoundingMode, Arblib.ArbRoundToZero) == RoundToZero
+    @test convert(RoundingMode, Arblib.ArbRoundFromZero) == RoundFromZero
+    @test convert(RoundingMode, Arblib.ArbRoundDown) == RoundDown
+    @test convert(RoundingMode, Arblib.ArbRoundUp) == RoundUp
+    @test convert(RoundingMode, Arblib.ArbRoundNearest) == RoundNearest
+end

test/runtests.jl

@@ -30,7 +30,7 @@ DocMeta.setdocmeta!(Arblib, :DocTestSetup, :(using Arblib); recursive = true)
     include("ArbCall/runtests.jl")
 
     include("arb_types.jl")
-    include("rounding.jl")
+    include("rounding_types.jl")
     include("types.jl")
     include("hash.jl")
     include("serialize.jl")
@@ -50,6 +50,7 @@ DocMeta.setdocmeta!(Arblib, :DocTestSetup, :(using Arblib); recursive = true)
     include("float.jl")
     include("interval.jl")
     include("multi-argument.jl")
+    include("rounding.jl")
     include("vector.jl")
     include("matrix.jl")
     include("eigen.jl")