diff --git a/src/Arblib.jl b/src/Arblib.jl index 124ff0d..6b47e34 100644 --- a/src/Arblib.jl +++ b/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") diff --git a/src/conversion.jl b/src/conversion.jl index c6037ae..b84564e 100644 --- a/src/conversion.jl +++ b/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. diff --git a/src/rounding.jl b/src/rounding.jl index 4bd9fa1..1e615d8 100644 --- a/src/rounding.jl +++ b/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 diff --git a/src/rounding_types.jl b/src/rounding_types.jl new file mode 100644 index 0000000..4bd9fa1 --- /dev/null +++ b/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 diff --git a/test/conversion.jl b/test/conversion.jl index ddd8ef4..527993e 100644 --- a/test/conversion.jl +++ b/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 diff --git a/test/rounding.jl b/test/rounding.jl index f769b82..42470df 100644 --- a/test/rounding.jl +++ b/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 diff --git a/test/rounding_types.jl b/test/rounding_types.jl new file mode 100644 index 0000000..7d23bdc --- /dev/null +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 311f064..86a127d 100644 --- a/test/runtests.jl +++ b/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")