diff --git a/src/Arblib.jl b/src/Arblib.jl index 98161905..2429830b 100644 --- a/src/Arblib.jl +++ b/src/Arblib.jl @@ -71,9 +71,12 @@ include("predicates.jl") include("show.jl") include("promotion.jl") include("arithmetic.jl") +include("elementary.jl") +include("minmax.jl") include("rand.jl") include("float.jl") include("interval.jl") +include("multi-argument.jl") include("ref.jl") include("vector.jl") diff --git a/src/arithmetic.jl b/src/arithmetic.jl index 17848660..930f56c9 100644 --- a/src/arithmetic.jl +++ b/src/arithmetic.jl @@ -2,55 +2,63 @@ const _BitInteger = (Int == Int64 ? Base.BitInteger64 : Base.BitInteger32) const _BitSigned = (Int == Int64 ? Base.BitSigned64 : Base.BitSigned32) const _BitUnsigned = (Int == Int64 ? Base.BitUnsigned64 : Base.BitUnsigned32) -### Mag +### +, -, *, / for (jf, af) in [(:+, :add!), (:-, :sub!), (:*, :mul!), (:/, :div!)] @eval Base.$jf(x::MagOrRef, y::MagOrRef) = $af(zero(x), x, y) -end -Base.:+(x::MagOrRef, y::Integer) = add!(zero(x), x, convert(UInt, y)) -Base.:+(x::Integer, y::MagOrRef) = add!(zero(y), y, convert(UInt, x)) -Base.:*(x::MagOrRef, y::Integer) = mul!(zero(x), x, convert(UInt, y)) -Base.:*(x::Integer, y::MagOrRef) = mul!(zero(y), y, convert(UInt, x)) -Base.:/(x::MagOrRef, y::Integer) = div!(zero(x), x, convert(UInt, y)) - -Base.:(^)(x::MagOrRef, e::Integer) = pow!(zero(x), x, convert(UInt, e)) -rsqrt(x::MagOrRef) = rsqrt!(zero(x), x) -Base.hypot(x::MagOrRef, y::MagOrRef) = hypot!(zero(x), x, y) -root(x::MagOrRef, n::Integer) = root!(zero(x), x, convert(UInt, n)) -neglog(x::MagOrRef) = neg_log!(zero(x), x) -expinv(x::MagOrRef) = expinv!(zero(x), x) -for f in [:inv, :sqrt, :log, :log1p, :exp, :expm1, :atan, :cosh, :sinh] - @eval Base.$f(x::MagOrRef) = $(Symbol(f, :!))(zero(x), x) -end -Base.min(x::MagOrRef, y::MagOrRef) = Arblib.min!(zero(x), x, y) -Base.max(x::MagOrRef, y::MagOrRef) = Arblib.max!(zero(x), x, y) -Base.minmax(x::MagOrRef, y::MagOrRef) = (min(x, y), max(x, y)) - -### Arf -function Base.sign(x::ArfOrRef) - isnan(x) && return Arf(NaN) # Follow Julia and return NaN - return Arf(sgn(x)) -end - -Base.abs(x::ArfOrRef) = abs!(zero(x), x) -Base.:(-)(x::ArfOrRef) = neg!(zero(x), x) -for (jf, af) in [(:+, :add!), (:-, :sub!), (:*, :mul!), (:/, :div!)] @eval function Base.$jf(x::ArfOrRef, y::Union{ArfOrRef,_BitInteger}) z = Arf(prec = _precision(x, y)) $af(z, x, y) return z end + + @eval Base.$jf(x::ArbOrRef, y::Union{ArbOrRef,ArfOrRef,_BitInteger}) = + $af(Arb(prec = _precision(x, y)), x, y) + + @eval Base.$jf(x::AcbOrRef, y::Union{AcbOrRef,ArbOrRef,_BitInteger}) = + $af(Acb(prec = _precision(x, y)), x, y) + + # Avoid one allocation for operations on irrationals or BitInteger + # rationals + @eval function Base.$jf( + x::Union{ArbOrRef,AcbOrRef}, + y::Union{Irrational,Rational{<:_BitInteger}}, + ) + z = zero(x) + z[] = y + return $af(z, x, z) + end + + if jf == :+ || jf == :* + # Addition and multiplication is commutative + @eval Base.$jf(x::_BitInteger, y::ArfOrRef) = $jf(y, x) + + @eval Base.$jf(x::Union{ArfOrRef,_BitInteger,Irrational}, y::ArbOrRef) = $jf(y, x) + + @eval Base.$jf(x::Union{ArbOrRef,_BitInteger,Irrational}, y::AcbOrRef) = $jf(y, x) + else + # Subtraction and division also need optimizations for when + # left argument is an integer + @eval function Base.$jf( + x::Union{Irrational,Rational{<:_BitInteger},_BitInteger}, + y::Union{ArbOrRef,AcbOrRef}, + ) + z = zero(y) + z[] = x + return $af(z, z, y) + end + end end -function Base.:+(x::_BitInteger, y::ArfOrRef) - z = zero(y) - add!(z, y, x) - return z -end -function Base.:*(x::_BitInteger, y::ArfOrRef) - z = zero(y) - mul!(z, y, x) - return z -end + +Base.:-(x::Union{ArfOrRef,ArbOrRef,AcbOrRef}) = neg!(zero(x), x) +Base.inv(x::Union{MagOrRef,ArbOrRef,AcbOrRef}) = inv!(zero(x), x) + +Base.:+(x::MagOrRef, y::Integer) = add!(zero(x), x, convert(UInt, y)) +Base.:+(x::Integer, y::MagOrRef) = y + x +Base.:*(x::MagOrRef, y::Integer) = mul!(zero(x), x, convert(UInt, y)) +Base.:*(x::Integer, y::MagOrRef) = y * x +Base.:/(x::MagOrRef, y::Integer) = div!(zero(x), x, convert(UInt, y)) + function Base.:/(x::_BitUnsigned, y::ArfOrRef) z = zero(y) ui_div!(z, x, y) @@ -62,193 +70,78 @@ function Base.:/(x::_BitSigned, y::ArfOrRef) return z end -function Base.sqrt(x::ArfOrRef) - y = zero(x) - sqrt!(y, x) - return y -end -function rsqrt(x::ArfOrRef) - y = zero(x) - rsqrt!(y, x) - return y +Base.:(/)(x::_BitUnsigned, y::ArbOrRef) = ui_div!(zero(y), x, y) + +function Base.:*(x::AcbOrRef, y::Complex{Bool}) + if real(y) + if imag(y) + z = mul_onei!(zero(x), x) + return add!(z, x, z) + else + return Acb(x) + end + end + imag(y) && return mul_onei!(zero(x), x) + return zero(x) end -function root(x::ArfOrRef, k::Integer) - y = zero(x) - root!(y, x, convert(UInt, k)) - return y +Base.:*(x::Complex{Bool}, y::AcbOrRef) = y * x + +### fma and muladd +function Base.fma(x::ArfOrRef, y::ArfOrRef, z::ArfOrRef) + res = zero(x) + fma!(res, x, y, z) + return res end +Base.fma(x::ArbOrRef, y::ArbOrRef, z::ArbOrRef) = fma!(zero(x), x, y, z) +Base.fma(x::ArbOrRef, y::Union{_BitInteger,ArfOrRef}, z::ArbOrRef) = fma!(zero(x), x, y, z) +Base.fma(x::Union{_BitInteger,ArfOrRef}, y::ArbOrRef, z::ArbOrRef) = fma(y, x, z) -Base.min(x::ArfOrRef, y::ArfOrRef) = Arblib.min!(zero(x), x, y) -Base.max(x::ArfOrRef, y::ArfOrRef) = Arblib.max!(zero(x), x, y) -Base.minmax(x::ArfOrRef, y::ArfOrRef) = (min(x, y), max(x, y)) +Base.muladd(x::ArfOrRef, y::ArfOrRef, z::ArfOrRef) = fma(x, y, z) +Base.muladd(x::ArbOrRef, y::ArbOrRef, z::ArbOrRef) = fma(x, y, z) +Base.muladd(x::ArbOrRef, y::Union{_BitInteger,ArfOrRef}, z::ArbOrRef) = fma(x, y, z) +Base.muladd(x::Union{_BitInteger,ArfOrRef}, y::ArbOrRef, z::ArbOrRef) = fma(x, y, z) -### Arb and Acb -for (jf, af) in [(:+, :add!), (:-, :sub!), (:*, :mul!), (:/, :div!)] - @eval Base.$jf(x::ArbOrRef, y::Union{ArbOrRef,ArfOrRef,_BitInteger}) = - $af(Arb(prec = _precision(x, y)), x, y) - @eval Base.$jf(x::AcbOrRef, y::Union{AcbOrRef,ArbOrRef,_BitInteger}) = - $af(Acb(prec = _precision(x, y)), x, y) - if jf == :(+) || jf == :(*) - @eval Base.$jf(x::Union{ArfOrRef,_BitInteger}, y::ArbOrRef) = - $af(Arb(prec = _precision(x, y)), y, x) - @eval Base.$jf(x::Union{ArbOrRef,_BitInteger}, y::AcbOrRef) = - $af(Acb(prec = _precision(x, y)), y, x) - end -end +### signbit, sign and abs +Base.signbit(x::MagOrRef) = false +Base.signbit(x::ArfOrRef) = !isnan(x) && sgn(x) < 0 +Base.signbit(x::ArbOrRef) = isnegative(x) -Base.:(-)(x::Union{ArbOrRef,AcbOrRef}) = neg!(zero(x), x) -Base.abs(x::ArbOrRef) = abs!(zero(x), x) -Base.:(/)(x::_BitUnsigned, y::ArbOrRef) = ui_div!(zero(y), x, y) +# For Arf sign of NaN is undefined, we follow Julia and return NaN +Base.sign(x::MagOrRef) = iszero(x) ? zero(x) : one(x) +Base.sign(x::ArfOrRef) = isnan(x) ? nan!(zero(x)) : Arf(sgn(x), prec = _precision(x)) +Base.sign(x::Union{ArbOrRef,AcbRef}) = sgn!(zero(x), x) + +Base.abs(x::MagOrRef) = copy(x) +Base.abs(x::Union{ArfOrRef,ArbOrRef}) = abs!(zero(x), x) +Base.abs(z::AcbOrRef) = abs!(Arb(prec = _precision(z)), z) + +### ^ -Base.:(^)(x::ArbOrRef, y::ArbOrRef) = pow!(Arb(prec = _precision(x, y)), x, y) -function Base.:(^)(x::ArbOrRef, y::_BitInteger) +Base.:^(x::MagOrRef, e::Integer) = pow!(zero(x), x, convert(UInt, e)) +Base.:^(x::ArbOrRef, y::ArbOrRef) = pow!(Arb(prec = _precision(x, y)), x, y) +function Base.:^(x::ArbOrRef, y::_BitInteger) z = zero(x) x, n = (y >= 0 ? (x, y) : (inv!(z, x), -y)) return pow!(z, x, convert(UInt, n)) end -Base.:(^)(x::AcbOrRef, y::Union{AcbOrRef,ArbOrRef,_BitInteger}) = +Base.:^(x::AcbOrRef, y::Union{AcbOrRef,ArbOrRef,_BitInteger}) = pow!(Acb(prec = _precision(x, y)), x, y) -# We define the same special cases as Arb does, this avoids some -# overhead +sqr(x::Union{ArbOrRef,AcbOrRef}) = sqr!(zero(x), x) +cube(x::AcbOrRef) = cube!(zero(x), x) + +# We define the same special cases as Arb does, this avoids some overhead +Base.literal_pow(::typeof(^), x::AcbOrRef, ::Val{-3}) = (y = inv(x); cube!(y, y)) Base.literal_pow(::typeof(^), x::Union{ArbOrRef,AcbOrRef}, ::Val{-2}) = (y = inv(x); sqr!(y, y)) #Base.literal_pow(::typeof(^), x::Union{ArbOrRef,AcbOrRef}, ::Val{-1}) - implemented in Base.intfuncs.jl Base.literal_pow(::typeof(^), x::Union{ArbOrRef,AcbOrRef}, ::Val{0}) = one(x) Base.literal_pow(::typeof(^), x::Union{ArbOrRef,AcbOrRef}, ::Val{1}) = copy(x) Base.literal_pow(::typeof(^), x::Union{ArbOrRef,AcbOrRef}, ::Val{2}) = sqr(x) +Base.literal_pow(::typeof(^), x::AcbOrRef, ::Val{3}) = cube(x) -Base.hypot(x::ArbOrRef, y::ArbOrRef) = hypot!(Arb(prec = _precision(x, y)), x, y) - -root(x::Union{ArbOrRef,AcbOrRef}, k::Integer) = root!(zero(x), x, convert(UInt, k)) - -# Unary methods in Base -for f in [ - :inv, - :sqrt, - :log, - :log1p, - :exp, - :expm1, - :sin, - :cos, - :tan, - :cot, - :sec, - :csc, - :atan, - :asin, - :acos, - :sinh, - :cosh, - :tanh, - :coth, - :sech, - :csch, - :atanh, - :asinh, - :acosh, -] - @eval Base.$f(x::Union{ArbOrRef,AcbOrRef}) = $(Symbol(f, :!))(zero(x), x) -end - -sqrtpos(x::ArbOrRef) = sqrtpos!(zero(x), x) -sqrt1pm1(x::ArbOrRef) = sqrt1pm1!(zero(x), x) -rsqrt(x::Union{ArbOrRef,AcbOrRef}) = rsqrt!(zero(x), x) -sqr(x::Union{ArbOrRef,AcbOrRef}) = sqr!(zero(x), x) - -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) -cotpi(x::Union{ArbOrRef,AcbOrRef}) = cot_pi!(zero(x), x) -cscpi(x::Union{ArbOrRef,AcbOrRef}) = csc_pi!(zero(x), x) -# Julias definition of sinc is equivalent to Arbs definition of sincpi -Base.sinc(x::Union{ArbOrRef,AcbOrRef}) = sinc_pi!(zero(x), x) -Base.atan(y::ArbOrRef, x::ArbOrRef) = atan2!(Arb(prec = _precision(y, x)), y, x) - -function Base.sincos(x::Union{ArbOrRef,AcbOrRef}) - s, c = zero(x), zero(x) - sin_cos!(s, c, x) - return (s, c) -end -function Base.sincospi(x::Union{ArbOrRef,AcbOrRef}) - s, c = zero(x), zero(x) - sin_cos_pi!(s, c, x) - return (s, c) -end -function sinhcosh(x::Union{ArbOrRef,AcbOrRef}) - s, c = zero(x), zero(x) - sinh_cosh!(s, c, x) - return (s, c) -end - -Base.min(x::ArbOrRef, y::ArbOrRef) = Arblib.min!(zero(x), x, y) -Base.max(x::ArbOrRef, y::ArbOrRef) = Arblib.max!(zero(x), x, y) -Base.minmax(x::ArbOrRef, y::ArbOrRef) = (min(x, y), max(x, y)) - -### Acb -function Base.:(*)(x::AcbOrRef, y::Complex{Bool}) - if real(y) - if imag(y) - z = mul_onei!(zero(x), x) - return add!(z, x, z) - else - return Acb(x) - end - end - imag(y) && return mul_onei!(zero(x), x) - return zero(x) -end -Base.:(*)(x::Complex{Bool}, y::AcbOrRef) = y * x - -Base.real(z::AcbLike; prec = _precision(z)) = get_real!(Arb(; prec), z) -Base.imag(z::AcbLike; prec = _precision(z)) = get_imag!(Arb(; prec), z) -Base.conj(z::AcbLike) = conj!(Acb(prec = _precision(z)), z) -Base.abs(z::AcbLike) = abs!(Arb(prec = _precision(z)), z) - -### minimum and maximum -# The default implemented in Julia have several issues for Arb types. -# See https://github.com/JuliaLang/julia/issues/45932. -# Note that it works fine for Mag and Arf. - -# Before 1.8.0 there is no way to fix the implementation in Base. -# Instead we define a new method. Note that this doesn't fully fix the -# problem, there is no way to dispatch on for example -# minimum(x -> Arb(x), [1, 2, 3, 4]) or an array with only some Arb. -if VERSION < v"1.8.0-rc3" - function Base.minimum(A::AbstractArray{<:ArbOrRef}) - isempty(A) && - throw(ArgumentError("reducing over an empty collection is not allowed")) - res = copy(first(A)) - for x in Iterators.drop(A, 1) - Arblib.min!(res, res, x) - end - return res - end - - function Base.maximum(A::AbstractArray{<:ArbOrRef}) - isempty(A) && - throw(ArgumentError("reducing over an empty collection is not allowed")) - res = copy(first(A)) - for x in Iterators.drop(A, 1) - Arblib.max!(res, res, x) - end - return res - end -end +### real, imag, conj -# Since 1.8.0 it is possible to fix the Base implementation by -# overloading some internal methods. This also works before 1.8.0 but -# doesn't solve the full problem. - -# The default implementation in Base is not correct for Arb -Base._fast(::typeof(min), x::Arb, y::Arb) = min(x, y) -Base._fast(::typeof(min), x::Arb, y) = min(x, y) -Base._fast(::typeof(min), x, y::Arb) = min(x, y) -Base._fast(::typeof(max), x::Arb, y::Arb) = max(x, y) -Base._fast(::typeof(max), x::Arb, y) = max(x, y) -Base._fast(::typeof(max), x, y::Arb) = max(x, y) - -# Mag, Arf and Arb don't have signed zeros -Base.isbadzero(::typeof(min), x::Union{Mag,Arf,Arb}) = false -Base.isbadzero(::typeof(max), x::Union{Mag,Arf,Arb}) = false +Base.real(z::AcbOrRef) = get_real!(Arb(prec = _precision(z)), z) +Base.imag(z::AcbOrRef) = get_imag!(Arb(prec = _precision(z)), z) +Base.conj(z::AcbOrRef) = conj!(zero(z), z) diff --git a/src/elementary.jl b/src/elementary.jl new file mode 100644 index 00000000..b829807c --- /dev/null +++ b/src/elementary.jl @@ -0,0 +1,92 @@ +### Mag + +rsqrt(x::MagOrRef) = rsqrt!(zero(x), x) +Base.hypot(x::MagOrRef, y::MagOrRef) = hypot!(zero(x), x, y) +root(x::MagOrRef, n::Integer) = root!(zero(x), x, convert(UInt, n)) +neglog(x::MagOrRef) = neg_log!(zero(x), x) +expinv(x::MagOrRef) = expinv!(zero(x), x) +for f in [:sqrt, :log, :log1p, :exp, :expm1, :atan, :cosh, :sinh] + @eval Base.$f(x::MagOrRef) = $(Symbol(f, :!))(zero(x), x) +end + +### Arf + +function Base.sqrt(x::ArfOrRef) + y = zero(x) + sqrt!(y, x) + return y +end +function rsqrt(x::ArfOrRef) + y = zero(x) + rsqrt!(y, x) + return y +end +function root(x::ArfOrRef, k::Integer) + y = zero(x) + root!(y, x, convert(UInt, k)) + return y +end + +### Arb and Acb + +Base.hypot(x::ArbOrRef, y::ArbOrRef) = hypot!(Arb(prec = _precision(x, y)), x, y) + +root(x::Union{ArbOrRef,AcbOrRef}, k::Integer) = root!(zero(x), x, convert(UInt, k)) + +# Unary methods in Base +for f in [ + :sqrt, + :log, + :log1p, + :exp, + :expm1, + :sin, + :cos, + :tan, + :cot, + :sec, + :csc, + :atan, + :asin, + :acos, + :sinh, + :cosh, + :tanh, + :coth, + :sech, + :csch, + :atanh, + :asinh, + :acosh, +] + @eval Base.$f(x::Union{ArbOrRef,AcbOrRef}) = $(Symbol(f, :!))(zero(x), x) +end + +sqrtpos(x::ArbOrRef) = sqrtpos!(zero(x), x) +sqrt1pm1(x::ArbOrRef) = sqrt1pm1!(zero(x), x) +rsqrt(x::Union{ArbOrRef,AcbOrRef}) = rsqrt!(zero(x), x) + +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) +cotpi(x::Union{ArbOrRef,AcbOrRef}) = cot_pi!(zero(x), x) +cscpi(x::Union{ArbOrRef,AcbOrRef}) = csc_pi!(zero(x), x) +# Julias definition of sinc is equivalent to Arbs definition of sincpi +Base.sinc(x::Union{ArbOrRef,AcbOrRef}) = sinc_pi!(zero(x), x) +Base.atan(y::ArbOrRef, x::ArbOrRef) = atan2!(Arb(prec = _precision(y, x)), y, x) + +function Base.sincos(x::Union{ArbOrRef,AcbOrRef}) + s, c = zero(x), zero(x) + sin_cos!(s, c, x) + return (s, c) +end +function Base.sincospi(x::Union{ArbOrRef,AcbOrRef}) + s, c = zero(x), zero(x) + sin_cos_pi!(s, c, x) + return (s, c) +end +function sinhcosh(x::Union{ArbOrRef,AcbOrRef}) + s, c = zero(x), zero(x) + sinh_cosh!(s, c, x) + return (s, c) +end diff --git a/src/minmax.jl b/src/minmax.jl new file mode 100644 index 00000000..aec6e5eb --- /dev/null +++ b/src/minmax.jl @@ -0,0 +1,58 @@ +Base.min(x::MagOrRef, y::MagOrRef) = Arblib.min!(zero(x), x, y) +Base.max(x::MagOrRef, y::MagOrRef) = Arblib.max!(zero(x), x, y) +Base.minmax(x::MagOrRef, y::MagOrRef) = (min(x, y), max(x, y)) + +for T in (ArfOrRef, ArbOrRef) + @eval Base.min(x::$T, y::$T) = + Arblib.min!($(_nonreftype(T))(prec = _precision(x, y)), x, y) + @eval Base.max(x::$T, y::$T) = + Arblib.max!($(_nonreftype(T))(prec = _precision(x, y)), x, y) + @eval Base.minmax(x::$T, y::$T) = (min(x, y), max(x, y)) +end + +### minimum and maximum +# The default implemented in Julia have several issues for the Arb type. +# See https://github.com/JuliaLang/julia/issues/45932. +# Note that it works fine for Mag and Arf. + +# Before 1.8.0 there is no way to fix the implementation in Base. +# Instead we define a new method. Note that this doesn't fully fix the +# problem, there is no way to dispatch on for example +# minimum(x -> Arb(x), [1, 2, 3, 4]) or an array with only some Arb. +if VERSION < v"1.8.0-rc3" + function Base.minimum(A::AbstractArray{<:ArbOrRef}) + isempty(A) && + throw(ArgumentError("reducing over an empty collection is not allowed")) + res = copy(first(A)) + for x in Iterators.drop(A, 1) + Arblib.min!(res, res, x) + end + return res + end + + function Base.maximum(A::AbstractArray{<:ArbOrRef}) + isempty(A) && + throw(ArgumentError("reducing over an empty collection is not allowed")) + res = copy(first(A)) + for x in Iterators.drop(A, 1) + Arblib.max!(res, res, x) + end + return res + end +end + +# Since 1.8.0 it is possible to fix the Base implementation by +# overloading some internal methods. This also works before 1.8.0 but +# doesn't solve the full problem. + +# The default implementation in Base is not correct for Arb +Base._fast(::typeof(min), x::Arb, y::Arb) = min(x, y) +Base._fast(::typeof(min), x::Arb, y) = min(x, y) +Base._fast(::typeof(min), x, y::Arb) = min(x, y) +Base._fast(::typeof(max), x::Arb, y::Arb) = max(x, y) +Base._fast(::typeof(max), x::Arb, y) = max(x, y) +Base._fast(::typeof(max), x, y::Arb) = max(x, y) + +# Mag, Arf and Arb don't have signed zeros +Base.isbadzero(::typeof(min), x::Union{Mag,Arf,Arb}) = false +Base.isbadzero(::typeof(max), x::Union{Mag,Arf,Arb}) = false diff --git a/src/multi-argument.jl b/src/multi-argument.jl new file mode 100644 index 00000000..6b9486e3 --- /dev/null +++ b/src/multi-argument.jl @@ -0,0 +1,64 @@ +# Contains implementation of multi-argument versions of +, *, min and +# max. They are more efficient because the don't need to allocate as +# much. This is similar to the specialised methods for + and * for +# BigFloat. + +for (jf, af) in [(:+, :add!), (:*, :mul!), (:min, :min!), (:max, :max!)] + @eval function Base.$jf(a::MagOrRef, b::MagOrRef, c::MagOrRef) + res = Mag() + $af(res, a, b) + $af(res, res, c) + return res + end + + @eval function Base.$jf(a::MagOrRef, b::MagOrRef, c::MagOrRef, d::MagOrRef) + res = Mag() + $af(res, a, b) + $af(res, res, c) + $af(res, res, d) + return res + end + + @eval function Base.$jf(a::MagOrRef, b::MagOrRef, c::MagOrRef, d::MagOrRef, e::MagOrRef) + res = Mag() + $af(res, a, b) + $af(res, res, c) + $af(res, res, d) + $af(res, res, e) + return res + end +end + +for T in (ArfOrRef, ArbOrRef, AcbOrRef) + @eval @inline _precision(x::$T, y::$T, z::$T, rest::Vararg{S}) where {S<:$T} = + max(precision(x), _precision(y, z, rest...)) + + for (jf, af) in [(:+, :add!), (:*, :mul!), (:min, :min!), (:max, :max!)] + T == AcbOrRef && jf == :min && continue + T == AcbOrRef && jf == :max && continue + + @eval function Base.$jf(a::$T, b::$T, c::$T) + res = $(_nonreftype(T))(prec = _precision(a, b, c)) + $af(res, a, b) + $af(res, res, c) + return res + end + + @eval function Base.$jf(a::$T, b::$T, c::$T, d::$T) + res = $(_nonreftype(T))(prec = _precision(a, b, c, d)) + $af(res, a, b) + $af(res, res, c) + $af(res, res, d) + return res + end + + @eval function Base.$jf(a::$T, b::$T, c::$T, d::$T, e::$T) + res = $(_nonreftype(T))(prec = _precision(a, b, c, d, e)) + $af(res, a, b) + $af(res, res, c) + $af(res, res, d) + $af(res, res, e) + return res + end + end +end diff --git a/test/arithmetic.jl b/test/arithmetic.jl index b9b0d878..3d45bcfc 100644 --- a/test/arithmetic.jl +++ b/test/arithmetic.jl @@ -1,5 +1,6 @@ @testset "Arithmethic" begin @testset "Mag" begin + # +, -, *, / @test Mag(4) <= Mag(1) + Mag(3) <= Mag(5) @test Mag(2) <= Mag(3) - Mag(1) <= Mag(3) @test Mag(0) <= Mag(3) - Mag(5) <= Mag(1) @@ -11,33 +12,24 @@ @test Mag(2) <= Mag(4) / 2 <= Mag(3) @test Mag(1 / 4) <= inv(Mag(4)) <= Mag(1 / 3) - @test Mag(8) <= Mag(2)^3 <= Mag(9) - @test Mag(2) <= sqrt(Mag(4)) <= Mag(3) - @test Mag(1 / 2) <= Arblib.rsqrt(Mag(4)) <= Mag(1) - @test Mag(5) <= hypot(Mag(3), Mag(4)) <= Mag(6) - @test Mag(2) <= Arblib.root(Mag(8), 3) <= Mag(3) - @test Mag(log(8)) <= log(Mag(8)) <= Mag(log(9)) - @test log(Mag(2)) <= Arblib.neglog(inv(Mag(2.0))) <= log(Mag(3)) - @test log(Mag(2)) <= log1p(Mag(1)) <= log(Mag(3)) - @test Mag(exp(2)) <= exp(Mag(2)) <= Mag(exp(3)) - @test Mag(exp(-2)) <= Arblib.expinv(Mag(2)) <= Mag(exp(-1)) - @test Mag(expm1(2)) <= expm1(Mag(2)) <= Mag(expm1(3)) - @test Mag(atan(2)) <= atan(Mag(2)) <= Mag(atan(3)) - @test Mag(cosh(2)) <= cosh(Mag(2)) <= Mag(cosh(3)) - @test Mag(sinh(2)) <= sinh(Mag(2)) <= Mag(sinh(3)) - @test min(Mag(1), Mag(2)) == Mag(1) - @test max(Mag(1), Mag(2)) == Mag(2) - @test minmax(Mag(1), Mag(2)) == minmax(Mag(2), Mag(1)) == (Mag(1), Mag(2)) + # signbit, sign and abs + @test !signbit(Mag(0)) + @test !signbit(Mag(2)) + + @test iszero(sign(Mag(0))) + @test isone(sign(Mag(1))) + @test isone(sign(Mag(Inf))) + + @test abs(Mag(5)) == Mag(5) + + # ^ + @test Mag(8) <= Mag(2)^3 <= Mag(9) end @testset "Arf" begin - @test sign(Arf(2)) == Arf(1) - @test sign(Arf(-2)) == Arf(-1) - @test sign(Arf(0)) == Arf(0) - - @test abs(Arf(-2)) == abs(Arf(2)) == Arf(2) - @test -Arf(2) == Arf(-2) + # +, -, *, / + @test -Arf(2) == -2 @test Arf(1) + Arf(2) == Arf(1) + 2 == @@ -47,14 +39,14 @@ Arf(1) + UInt8(2) == Arf(1) + UInt8(2) == UInt8(2) + Arf(1) == - Arf(3) - @test Arf(1) - Arf(2) == Arf(1) - 2 == Arf(1) - UInt(2) == Arf(-1) + 3 + @test Arf(1) - Arf(2) == Arf(1) - 2 == Arf(1) - UInt(2) == -1 @test Arf(2) * Arf(3) == Arf(2) * 3 == 3 * Arf(2) == Arf(2) * UInt(3) == UInt(3) * Arf(2) == - Arf(6) + 6 @test Arf(6) / Arf(2) == Arf(6) / 2 == 6 / Arf(2) == @@ -62,18 +54,27 @@ UInt(6) / Arf(2) == Arf(6) / UInt8(2) == UInt8(6) / Arf(2) == - Arf(3) + 3 + + # fma and muladd + @test fma(Arf(2), Arf(3), Arf(4)) == muladd(Arf(2), Arf(3), Arf(4)) == 10 - @test sqrt(Arf(4)) == Arf(2) - @test Arblib.rsqrt(Arf(4)) == Arf(1 // 2) - @test Arblib.root(Arf(8), 3) == Arf(2) + # signbit, sign and abs + @test signbit(Arf(-2)) + @test !signbit(Arf(0)) + @test !signbit(Arf(2)) + @test !signbit(Arf(NaN)) - @test min(Arf(1), Arf(2)) == Arf(1) - @test max(Arf(1), Arf(2)) == Arf(2) - @test minmax(Arf(1), Arf(2)) == minmax(Arf(2), Arf(1)) == (Arf(1), Arf(2)) + @test sign(Arf(-2)) == -1 + @test sign(Arf(0)) == 0 + @test sign(Arf(2)) == 1 + @test sign(Arf(NaN)) == NaN + + @test abs(Arf(-2)) == abs(Arf(2)) == 2 end @testset "$T" for T in [Arb, Acb] + # +, -, *, / @test T(1) + T(2) == T(1) + 2 == 2 + T(1) == @@ -81,123 +82,80 @@ UInt(2) + T(1) == T(1) + UInt8(2) == UInt8(2) + T(1) == - T(3) - @test T(1) - T(2) == T(1) - 2 == T(1) - UInt(2) == T(-1) - @test T(2) * T(3) == - T(2) * 3 == - 3 * T(2) == - T(2) * UInt(3) == - UInt(3) * T(2) == - T(6) - @test T(6) / T(2) == T(6) / 2 == T(6) / UInt(2) == T(3) - + 3 + @test T(1) - T(2) == T(1) - 2 == T(1) - UInt(2) == 1 - T(2) == T(-1) + @test T(2) * T(3) == T(2) * 3 == 3 * T(2) == T(2) * UInt(3) == UInt(3) * T(2) == 6 + @test T(6) / T(2) == T(6) / 2 == T(6) / UInt(2) == 6 / T(2) == 3 + + @test isequal(T(1) + π, T(1) + T(π)) + @test isequal(π + T(1), T(1) + T(π)) + @test T(1) + 3 // 2 == 3 // 2 + T(1) == 5 // 2 + @test isequal(T(1) - π, T(1) - T(π)) + @test isequal(π - T(1), T(π) - T(1)) + @test T(1) - 3 // 2 == -(3 // 2 - T(1)) == -1 // 2 + @test isequal(T(2) * π, T(2) * T(π)) + @test isequal(π * T(2), T(π) * T(2)) + @test T(2) * 3 // 2 == 3 // 2 * T(2) == 3 + @test isequal(T(2) / π, T(2) / T(π)) + @test isequal(π / T(2), T(π) / T(2)) + @test T(4) / 4 // 1 == 4 // 1 / T(4) == 1 + + # ^ @test Base.literal_pow(^, T(2), Val(-2)) == T(2)^-2 == Arblib.sqr(inv(T(2))) == - T(1 // 4) - @test Base.literal_pow(^, T(2), Val(0)) == T(2)^0 == one(T) - @test Base.literal_pow(^, T(2), Val(1)) == T(2)^1 == T(2) - @test Base.literal_pow(^, T(2), Val(2)) == T(2)^2 == Arblib.sqr(T(2)) == T(4) - - @test Arblib.root(T(8), 3) == T(2) - - for f in [ - inv, - sqrt, - log, - log1p, - exp, - expm1, - sin, - cos, - tan, - cot, - sec, - csc, - atan, - asin, - acos, - sinh, - cosh, - tanh, - coth, - sech, - csch, - atanh, - asinh, - ] - @test f(T(0.5)) ≈ f(0.5) - end - @test acosh(T(2)) ≈ acosh(2) - - @test Arblib.rsqrt(T(1 // 4)) == T(2) - @test Arblib.sqr(T(3)) == T(9) - - @test sinpi(T(1)) == T(0) - @test cospi(T(1)) == T(-1) - @test Arblib.tanpi(T(1)) == T(0) - @test Arblib.cotpi(T(0.5)) == 0 - @test Arblib.cscpi(T(0.5)) == 1 - @test sinc(T(0.5)) ≈ sinc(0.5) - - @test isequal(sincos(T(1)), (sin(T(1)), cos(T(1)))) - if VERSION >= v"1.6" - @test isequal(sincospi(T(1)), (sinpi(T(1)), cospi(T(1)))) - else - @test isequal(Arblib.sincospi(T(1)), (sinpi(T(1)), cospi(T(1)))) - end - @test isequal(Arblib.sinhcosh(T(1)), (sinh(T(1)), cosh(T(1)))) + 1 // 4 + @test Base.literal_pow(^, T(2), Val(-1)) == T(2)^-1 == inv(T(2)) == 1 // 2 + @test Base.literal_pow(^, T(2), Val(0)) == T(2)^0 == 1 + @test Base.literal_pow(^, T(2), Val(1)) == T(2)^1 == 2 + @test Base.literal_pow(^, T(2), Val(2)) == T(2)^2 == Arblib.sqr(T(2)) == 4 end @testset "Arb - specific" begin - @test Arb(1) + Arf(2) == Arf(2) + Arb(1) == Arb(3) + @test Arb(1) + Arf(2) == Arf(2) + Arb(1) == 3 @test Arb(1) - Arf(2) == Arb(-1) - @test Arb(2) * Arf(3) == Arf(3) * Arb(2) == Arb(6) - @test Arb(6) / Arf(2) == UInt(6) / Arb(2) == UInt8(6) / Arb(2) == Arb(3) - + @test Arb(2) * Arf(3) == Arf(3) * Arb(2) == 6 + @test Arb(6) / Arf(2) == UInt(6) / Arb(2) == UInt8(6) / Arb(2) == 3 + + # fma and muladd + @test fma(Arb(2), Arb(3), Arb(4)) == + fma(Arb(2), 3, Arb(4)) == + fma(2, Arb(3), Arb(4)) == + muladd(Arb(2), Arb(3), Arb(4)) == + muladd(Arb(2), 3, Arb(4)) == + muladd(2, Arb(3), Arb(4)) == + 10 + + # signbit, sign and abs + @test signbit(Arb(-2)) + @test !signbit(Arb(0)) + @test !signbit(Arb(2)) + @test !signbit(Arb((-1, 0))) + @test !signbit(Arb(NaN)) + + @test sign(Arb(-2)) == -1 + @test sign(Arb(0)) == 0 + @test sign(Arb(2)) == 1 + @test isequal(sign(Arb((-3, 1))), Arblib.zero_pm_one!(Arb())) + @test isequal(sign(Arb(NaN)), Arblib.zero_pm_one!(Arb())) + + @test abs(Arb(-2)) == abs(Arb(2)) == 2 + + # ^ @test Arb(2)^Arb(3) == Arb(2)^3 == Arb(2)^Int8(3) == Arb(2)^UInt(3) == Arb(2)^UInt8(3) == (Arb(2)^-3)^-1 == - Arb(8) - @test hypot(Arb(3), Arb(4)) == Arb(5) - - @test Arblib.sqrtpos(Arb(4)) == Arb(2) - @test Arblib.sqrtpos(Arb(-4)) == Arb(0) - @test Arblib.sqrt1pm1(Arb(3)) == Arb(1) - - @test atan(Arb(2), Arb(3)) ≈ atan(2, 3) - - @test min(Arb(1), Arb(2)) == Arb(1) - @test max(Arb(1), Arb(2)) == Arb(2) - @test minmax(Arb(1), Arb(2)) == minmax(Arb(2), Arb(1)) == (Arb(1), Arb(2)) - @test Arblib.contains(min(Arb((0, 2)), Arb((-1, 3))), -1) - @test Arblib.contains(min(Arb((0, 2)), Arb((-1, 3))), 2) - @test !Arblib.contains(min(Arb((0, 2)), Arb((-1, 3))), 3) - @test Arblib.contains(max(Arb((0, 2)), Arb((-1, 3))), 0) - @test Arblib.contains(max(Arb((0, 2)), Arb((-1, 3))), 3) - @test !Arblib.contains(max(Arb((0, 2)), Arb((-1, 3))), -1) - @test all(Arblib.contains.(minmax(Arb((0, 2)), Arb((-1, 3))), (-1, 0))) - @test all(Arblib.contains.(minmax(Arb((0, 2)), Arb((-1, 3))), (2, 3))) - @test all(.!Arblib.contains.(minmax(Arb((0, 2)), Arb((-1, 3))), (3, -1))) + 8 end @testset "Acb - specific" begin - @test Acb(1) + Arb(2) == Arb(2) + Acb(1) == Acb(3) - @test Acb(1) - Arb(2) == Acb(-1) - @test Acb(2) * Arb(3) == Arb(3) * Acb(2) == Acb(6) - @test Acb(6) / Arb(2) == Acb(3) - - @test Acb(2)^Acb(3) == - Acb(2)^Arb(3) == - Acb(2)^3 == - Acb(2)^Int8(3) == - Acb(2)^UInt(3) == - Acb(2)^UInt8(3) == - (Acb(2)^-3)^-1 == - Acb(8) + @test Acb(1) + Arb(2) == Arb(2) + Acb(1) == 3 + @test Acb(1) - Arb(2) == -1 + @test Acb(2) * Arb(3) == Arb(3) * Acb(2) == 6 + @test Acb(6) / Arb(2) == 3 @test Acb(2, 3) * Complex{Bool}(0, 0) == Complex{Bool}(0, 0) * Acb(2, 3) == @@ -212,101 +170,32 @@ Complex{Bool}(1, 1) * Acb(2, 3) == Acb(2, 3) * (1 + 1im) - @test real(Acb(1, 2)) isa Arb - @test real(Acb(1, 2)) == Arb(1) - @test imag(Acb(1, 2)) isa Arb - @test imag(Acb(1, 2)) == Arb(2) - - @test conj(Acb(1, 2)) == Acb(1, -2) - + # abs @test abs(Acb(3, 4)) isa Arb - @test abs(Acb(3, 4)) == Arb(5) - end - - @testset "minimum/maximum/extrema" begin - # See https://github.com/JuliaLang/julia/issues/45932 for - # discussions about the issues with the Base implementation - - # Currently there is no special implementation of extrema, the - # default implementation works well. But to help find future - # issues we test it here as well. - - A = Arb[10:20; 0:9] - @test minimum(A) == 0 - @test maximum(A) == 20 - @test extrema(A) == (0, 20) - - A = [Arb((i, i + 1)) for i = 0:20] - @test contains(minimum(A), Arb((0, 1))) - @test contains(minimum(reverse(A)), Arb((0, 1))) - @test contains(maximum(A), Arb((20, 21))) - @test contains(maximum(reverse(A)), Arb((20, 21))) - @test all(contains.(extrema(A), (Arb((0, 1)), Arb((20, 21))))) - @test all(contains.(extrema(reverse(A)), (Arb((0, 1)), Arb((20, 21))))) - - # Fails in Julia version < 1.8 with default implementation due - # to short circuiting in Base.mapreduce_impl - A = [setball(Arb, 2, 1); zeros(Arb, 257)] - @test iszero(minimum(A)) - @test iszero(maximum(-A)) - @test iszero(extrema(A)[1]) - @test iszero(extrema(-A)[2]) - # Before 1.8.0 these still fails and there is no real way to - # overload them - if VERSION < v"1.8.0-rc3" - @test_broken iszero(minimum(identity, A)) - @test_broken iszero(maximum(identity, -A)) - else - @test iszero(minimum(identity, A)) - @test iszero(maximum(identity, -A)) - end - # These work - @test iszero(extrema(identity, A)[1]) - @test iszero(extrema(identity, -A)[2]) + @test abs(Acb(3, 4)) == 5 - # Fails with default implementation due to Base._fast - #A = [Arb(0); [setball(Arb, 0, i) for i in reverse(0:257)]] - A = [setball(Arb, 0, i) for i = 0:257] - @test contains(minimum(A), -257) - @test contains(maximum(A), 257) - @test contains(extrema(A)[1], -257) - @test contains(extrema(A)[2], 257) - @test contains(minimum(identity, A), -257) - @test contains(maximum(identity, A), 257) - @test contains(extrema(identity, A)[1], -257) - @test contains(extrema(identity, A)[2], 257) - - # Fails with default implementation due to both short circuit - # and Base._fast - A = [setball(Arb, 0, i) for i = 0:1000] - @test contains(minimum(A), -1000) - @test contains(maximum(A), 1000) - @test contains(extrema(A)[1], -1000) - @test contains(extrema(A)[2], 1000) - if VERSION < v"1.8.0-rc3" - @test_broken contains(minimum(identity, A), -1000) - @test_broken contains(maximum(identity, A), 1000) - else - @test contains(minimum(identity, A), -1000) - @test contains(maximum(identity, A), 1000) - end - @test contains(extrema(identity, A)[1], -1000) - @test contains(extrema(identity, A)[2], 1000) - - @test !Base.isbadzero(min, zero(Mag)) - @test !Base.isbadzero(min, zero(Arf)) - @test !Base.isbadzero(min, zero(Arb)) + # ^ + @test Acb(2)^Acb(3) == + Acb(2)^Arb(3) == + Acb(2)^3 == + Acb(2)^Int8(3) == + Acb(2)^UInt(3) == + Acb(2)^UInt8(3) == + (Acb(2)^-3)^-1 == + 8 - @test !Base.isbadzero(max, zero(Mag)) - @test !Base.isbadzero(max, zero(Arf)) - @test !Base.isbadzero(max, zero(Arb)) + @test Base.literal_pow(^, Acb(2), Val(-3)) == + Acb(2)^-3 == + Arblib.cube(inv(Acb(2))) == + 1 // 8 + @test Base.literal_pow(^, Acb(2), Val(3)) == Acb(2)^3 == Arblib.cube(Acb(2)) == 8 - @test !Base.isgoodzero(min, zero(Mag)) - @test !Base.isgoodzero(min, zero(Arf)) - @test !Base.isgoodzero(min, zero(Arb)) + # real, imag, conj + @test real(Acb(1, 2)) isa Arb + @test real(Acb(1, 2)) == 1 + @test imag(Acb(1, 2)) isa Arb + @test imag(Acb(1, 2)) == 2 - @test !Base.isgoodzero(max, zero(Mag)) - @test !Base.isgoodzero(max, zero(Arf)) - @test !Base.isgoodzero(max, zero(Arb)) + @test conj(Acb(1, 2)) == Acb(1, -2) end end diff --git a/test/elementary.jl b/test/elementary.jl new file mode 100644 index 00000000..00ae477c --- /dev/null +++ b/test/elementary.jl @@ -0,0 +1,78 @@ +@testset "Elementary" begin + @testset "Mag" begin + @test Mag(2) <= sqrt(Mag(4)) <= Mag(3) + @test Mag(1 / 2) <= Arblib.rsqrt(Mag(4)) <= Mag(1) + @test Mag(5) <= hypot(Mag(3), Mag(4)) <= Mag(6) + @test Mag(2) <= Arblib.root(Mag(8), 3) <= Mag(3) + @test Mag(log(8)) <= log(Mag(8)) <= Mag(log(9)) + @test log(Mag(2)) <= Arblib.neglog(inv(Mag(2.0))) <= log(Mag(3)) + @test log(Mag(2)) <= log1p(Mag(1)) <= log(Mag(3)) + @test Mag(exp(2)) <= exp(Mag(2)) <= Mag(exp(3)) + @test Mag(exp(-2)) <= Arblib.expinv(Mag(2)) <= Mag(exp(-1)) + @test Mag(expm1(2)) <= expm1(Mag(2)) <= Mag(expm1(3)) + @test Mag(atan(2)) <= atan(Mag(2)) <= Mag(atan(3)) + @test Mag(cosh(2)) <= cosh(Mag(2)) <= Mag(cosh(3)) + @test Mag(sinh(2)) <= sinh(Mag(2)) <= Mag(sinh(3)) + end + + @testset "Arf" begin + @test sqrt(Arf(4)) == 2 + @test Arblib.rsqrt(Arf(4)) == 1 // 2 + @test Arblib.root(Arf(8), 3) == 2 + end + + @testset "$T" for T in [Arb, Acb] + @test Arblib.root(T(8), 3) == T(2) + + for f in [ + sqrt, + log, + log1p, + exp, + expm1, + sin, + cos, + tan, + cot, + sec, + csc, + atan, + asin, + acos, + sinh, + cosh, + tanh, + coth, + sech, + csch, + atanh, + asinh, + ] + @test f(T(0.5)) ≈ f(0.5) + end + @test acosh(T(2)) ≈ acosh(2) + + @test Arblib.rsqrt(T(1 // 4)) == 2 + + @test sinpi(T(1)) == 0 + @test cospi(T(1)) == -1 + @test Arblib.tanpi(T(1)) == 0 + @test Arblib.cotpi(T(0.5)) == 0 + @test Arblib.cscpi(T(0.5)) == 1 + @test sinc(T(0.5)) ≈ sinc(0.5) + + @test isequal(sincos(T(1)), (sin(T(1)), cos(T(1)))) + @test isequal(sincospi(T(1)), (sinpi(T(1)), cospi(T(1)))) + @test isequal(Arblib.sinhcosh(T(1)), (sinh(T(1)), cosh(T(1)))) + end + + @testset "Arb - specific" begin + @test hypot(Arb(3), Arb(4)) == 5 + + @test Arblib.sqrtpos(Arb(4)) == 2 + @test Arblib.sqrtpos(Arb(-4)) == 0 + @test Arblib.sqrt1pm1(Arb(3)) == 1 + + @test atan(Arb(2), Arb(3)) ≈ atan(2, 3) + end +end diff --git a/test/minmax.jl b/test/minmax.jl new file mode 100644 index 00000000..9e250774 --- /dev/null +++ b/test/minmax.jl @@ -0,0 +1,108 @@ +@testset "MinMax" begin + @testset "$T" for T in (Mag, Arf, Mag) + @test min(T(1), T(2)) == T(1) + @test max(T(1), T(2)) == T(2) + @test minmax(T(1), T(2)) == minmax(T(2), T(1)) == (T(1), T(2)) + end + + @testset "Arb - specific" begin + @test Arblib.contains(min(Arb((0, 2)), Arb((-1, 3))), -1) + @test Arblib.contains(min(Arb((0, 2)), Arb((-1, 3))), 2) + @test !Arblib.contains(min(Arb((0, 2)), Arb((-1, 3))), 3) + @test Arblib.contains(max(Arb((0, 2)), Arb((-1, 3))), 0) + @test Arblib.contains(max(Arb((0, 2)), Arb((-1, 3))), 3) + @test !Arblib.contains(max(Arb((0, 2)), Arb((-1, 3))), -1) + @test all(Arblib.contains.(minmax(Arb((0, 2)), Arb((-1, 3))), (-1, 0))) + @test all(Arblib.contains.(minmax(Arb((0, 2)), Arb((-1, 3))), (2, 3))) + @test all(.!Arblib.contains.(minmax(Arb((0, 2)), Arb((-1, 3))), (3, -1))) + end + + @testset "minimum/maximum/extrema" begin + # See https://github.com/JuliaLang/julia/issues/45932 for + # discussions about the issues with the Base implementation + + # Currently there is no special implementation of extrema, the + # default implementation works well. But to help find future + # issues we test it here as well. + + @testset "$T" for T in (Mag, Arf, Mag) + A = T[10:20; 0:9] + @test minimum(A) == T(0) + @test maximum(A) == T(20) + @test extrema(A) == (T(0), T(20)) + end + + A = [Arb((i, i + 1)) for i = 0:20] + @test contains(minimum(A), Arb((0, 1))) + @test contains(minimum(reverse(A)), Arb((0, 1))) + @test contains(maximum(A), Arb((20, 21))) + @test contains(maximum(reverse(A)), Arb((20, 21))) + @test all(contains.(extrema(A), (Arb((0, 1)), Arb((20, 21))))) + @test all(contains.(extrema(reverse(A)), (Arb((0, 1)), Arb((20, 21))))) + + # Fails in Julia version < 1.8 with default implementation due + # to short circuiting in Base.mapreduce_impl + A = [setball(Arb, 2, 1); zeros(Arb, 257)] + @test iszero(minimum(A)) + @test iszero(maximum(-A)) + @test iszero(extrema(A)[1]) + @test iszero(extrema(-A)[2]) + # Before 1.8.0 these still fails and there is no real way to + # overload them + if VERSION < v"1.8.0-rc3" + @test_broken iszero(minimum(identity, A)) + @test_broken iszero(maximum(identity, -A)) + else + @test iszero(minimum(identity, A)) + @test iszero(maximum(identity, -A)) + end + # These work + @test iszero(extrema(identity, A)[1]) + @test iszero(extrema(identity, -A)[2]) + + # Fails with default implementation due to Base._fast + #A = [Arb(0); [setball(Arb, 0, i) for i in reverse(0:257)]] + A = [setball(Arb, 0, i) for i = 0:257] + @test contains(minimum(A), -257) + @test contains(maximum(A), 257) + @test contains(extrema(A)[1], -257) + @test contains(extrema(A)[2], 257) + @test contains(minimum(identity, A), -257) + @test contains(maximum(identity, A), 257) + @test contains(extrema(identity, A)[1], -257) + @test contains(extrema(identity, A)[2], 257) + + # Fails with default implementation due to both short circuit + # and Base._fast + A = [setball(Arb, 0, i) for i = 0:1000] + @test contains(minimum(A), -1000) + @test contains(maximum(A), 1000) + @test contains(extrema(A)[1], -1000) + @test contains(extrema(A)[2], 1000) + if VERSION < v"1.8.0-rc3" + @test_broken contains(minimum(identity, A), -1000) + @test_broken contains(maximum(identity, A), 1000) + else + @test contains(minimum(identity, A), -1000) + @test contains(maximum(identity, A), 1000) + end + @test contains(extrema(identity, A)[1], -1000) + @test contains(extrema(identity, A)[2], 1000) + + @test !Base.isbadzero(min, zero(Mag)) + @test !Base.isbadzero(min, zero(Arf)) + @test !Base.isbadzero(min, zero(Arb)) + + @test !Base.isbadzero(max, zero(Mag)) + @test !Base.isbadzero(max, zero(Arf)) + @test !Base.isbadzero(max, zero(Arb)) + + @test !Base.isgoodzero(min, zero(Mag)) + @test !Base.isgoodzero(min, zero(Arf)) + @test !Base.isgoodzero(min, zero(Arb)) + + @test !Base.isgoodzero(max, zero(Mag)) + @test !Base.isgoodzero(max, zero(Arf)) + @test !Base.isgoodzero(max, zero(Arb)) + end +end diff --git a/test/multi-argument.jl b/test/multi-argument.jl new file mode 100644 index 00000000..6429093b --- /dev/null +++ b/test/multi-argument.jl @@ -0,0 +1,49 @@ +@testset "Multi-argument" begin + @testset "Mag + and *" begin + @test Mag(6) <= +(Mag.((1, 2, 3))...) <= Mag(7) + @test Mag(10) <= +(Mag.((1, 2, 3, 4))...) <= Mag(11) + @test Mag(15) <= +(Mag.((1, 2, 3, 4, 5))...) <= Mag(16) + + @test Mag(6) <= *(Mag.((1, 2, 3))...) <= Mag(7) + @test Mag(24) <= *(Mag.((1, 2, 3, 4))...) <= Mag(25) + @test Mag(120) <= *(Mag.((1, 2, 3, 4, 5))...) <= Mag(121) + end + + @testset "$T + and *" for T in [Arf, Arb, Acb] + @test +(T.((1, 2, 3))...) == 6 + @test +(T.((1, 2, 3, 4))...) == 10 + @test +(T.((1, 2, 3, 4, 5))...) == 15 + + @test *(T.((1, 2, 3))...) == 6 + @test *(T.((1, 2, 3, 4))...) == 24 + @test *(T.((1, 2, 3, 4, 5))...) == 120 + end + + @testset "$T min and max" for T in [Mag, Arf, Arb] + @test min(T.((1, 2, 2))...) == T(1) + @test min(T.((1, 1, 2))...) == T(1) + @test min(T.((2, 2, 1))...) == T(1) + @test min(T.((1, 2, 2, 2))...) == T(1) + @test min(T.((1, 1, 2, 2))...) == T(1) + @test min(T.((2, 2, 1, 2))...) == T(1) + @test min(T.((2, 2, 2, 1))...) == T(1) + @test min(T.((1, 2, 2, 2, 2))...) == T(1) + @test min(T.((1, 1, 2, 2, 2))...) == T(1) + @test min(T.((2, 2, 1, 2, 2))...) == T(1) + @test min(T.((2, 2, 2, 1, 2))...) == T(1) + @test min(T.((2, 2, 2, 2, 1))...) == T(1) + + @test max(T.((2, 1, 1))...) == T(2) + @test max(T.((2, 2, 1))...) == T(2) + @test max(T.((1, 1, 2))...) == T(2) + @test max(T.((2, 1, 1, 1))...) == T(2) + @test max(T.((2, 2, 1, 1))...) == T(2) + @test max(T.((1, 1, 2, 1))...) == T(2) + @test max(T.((1, 1, 1, 2))...) == T(2) + @test max(T.((2, 1, 1, 1, 1))...) == T(2) + @test max(T.((2, 2, 1, 1, 1))...) == T(2) + @test max(T.((1, 1, 2, 1, 1))...) == T(2) + @test max(T.((1, 1, 1, 2, 1))...) == T(2) + @test max(T.((1, 1, 1, 1, 2))...) == T(2) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 9ddae11e..151d9273 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,9 +18,12 @@ using Arblib, Test, LinearAlgebra, Random, Serialization, SpecialFunctions include("promotion.jl") include("examples.jl") include("arithmetic.jl") + include("elementary.jl") + include("minmax.jl") include("rand.jl") include("float.jl") include("interval.jl") + include("multi-argument.jl") include("vector.jl") include("matrix.jl") include("eigen.jl")