Skip to content

Commit

Permalink
Fix implementation of minimum and maximum
Browse files Browse the repository at this point in the history 
This addresses some of the issues discussed in
JuliaLang/julia#45932.

For Julia 1.8.0 it completely fixes minimum and maximum, for earlier
versions of Julia it only partially fixes some of the issues.

The change completely removes the specialised implementation of
extrema. It also completely removes the specialised implementation of
minimum and maximum for Mag and Arf. The default implementations work
well, though are slower due to more allocations. In the future it
could maybe be an option to have a faster version using
MutableArithmetics.jl.

The Base implementation of minimum and maximum is fixed for Arb by
overloading Base._fast to make it correct for Arb. Before 1.8.0 this
is not enough to fix all problems and it thus keeps the specialised
implementation for Julia before 1.8.0-rc3 (latest release candidate as
of writing), for later versions the specialised implementation can be
removed.
  • Loading branch information
@Joel-Dahne
Joel-Dahne committed Jul 29, 2022
1 parent f6a07c3 commit 5e81504
Show file tree
Hide file tree
Showing 3 changed files with 145 additions and 70 deletions.
32 changes: 9 additions & 23 deletions README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -173,29 +173,15 @@ TODO: Come up with more examples
### Implementation details
In some cases the implementation in Julia implicitly makes certain
assumptions to improve performance and this can lead to issues. For
example the `maximum` method in Julia checks for `NaN` results (on
which is short fuses) using `x == x`, which works for most numerical
types but not for `Arb` (`x == x` is only true if the radius is zero).
See <https://github.com/JuliaLang/julia/issues/36287> for some more
details. Arblib implements its own `maximum` method which gives
rigorous results, but it only covers the case
`maximum(AbstractFloat{Arb})`.

``` julia
julia> f = i -> Arb((i, i + 1));

julia> A = f.(0:1000);

julia> maximum(A)
[1.00e+3 +/- 1.01]

julia> maximum(A, dims = 1)
1-element Array{Arb,1}:
[+/- 1.01]

julia> maximum(f, 0:1000)
[+/- 1.01]
```
example, prior to Julia version 1.8 the `minimum` and `maximum`
methods in Julia checked for `NaN` results (on which is short fuses)
using `x == x`, which works for most numerical types but not for `Arb`
(`x == x` is only true if the radius is zero). See
<https://github.com/JuliaLang/julia/issues/36287> and in particular
<https://github.com/JuliaLang/julia/issues/45932> for more details.
Since Julia version 1.8 the `minimum` and `maximum` methods work
correctly for `Arb`, for earlier versions of Julia it only works
correctly in some cases.

These types of problems are the hardest to find since they are not
clear from the documentation but you have to read the implementation,
Expand Down
75 changes: 49 additions & 26 deletions src/arithmetic.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,10 @@ 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
Expand Down Expand Up @@ -74,6 +78,10 @@ function root(x::ArfOrRef, k::Integer)
return y
end

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))

### Arb and Acb
for (jf, af) in [(:+, :add!), (:-, :sub!), (:*, :mul!), (:/, :div!)]
@eval Base.$jf(x::ArbOrRef, y::Union{ArbOrRef,ArfOrRef,_BitInteger}) =
Expand Down Expand Up @@ -174,6 +182,10 @@ function sinhcosh(x::Union{ArbOrRef,AcbOrRef})
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)
Expand All @@ -194,38 +206,49 @@ 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)

### min and max
for T in [MagOrRef, ArfOrRef, ArbOrRef]
for op in [:min, :max]
# min/max
@eval Base.$op(x::$T, y::$T) = $(Symbol(op, :!))(zero(x), x, y)
# minimum/maximum
# The default implementation of minimum/maximum in Julia is
# not correct for Arb, we implemented a correct version when
# no specific dimension is given. The default give correct
# results for Mag and Arf but wrong results for Arb in some
# cases, be careful!
@eval function Base.$(Symbol(op, :imum))(A::AbstractArray{<:$T})
isempty(A) &&
throw(ArgumentError("reducing over an empty collection is not allowed"))
res = copy(first(A))
for x in Iterators.drop(A, 1)
$(Symbol(op, :!))(res, res, x)
end
return res
### 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
@eval Base.minmax(x::$T, y::$T) = (min(x, y), max(x, y))

@eval function Base.extrema(A::AbstractArray{<:$T})
function Base.maximum(A::AbstractArray{<:ArbOrRef})
isempty(A) &&
throw(ArgumentError("reducing over an empty collection is not allowed"))
l = copy(first(A))
u = copy(first(A))
res = copy(first(A))
for x in Iterators.drop(A, 1)
min!(l, l, x)
max!(u, u, x)
Arblib.max!(res, res, x)
end
return (l, u)
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
108 changes: 87 additions & 21 deletions test/arithmetic.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,9 +29,6 @@
@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))
@test minimum(Mag[10:20; 0:9]) == Mag(0)
@test maximum(Mag[10:20; 0:9]) == Mag(20)
@test extrema(Mag[10:20; 0:9]) == (Mag(0), Mag(20))
end

@testset "Arf" begin
Expand Down Expand Up @@ -74,9 +71,6 @@
@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 minimum(Arf[10:20; 0:9]) == Arf(0)
@test maximum(Arf[10:20; 0:9]) == Arf(20)
@test extrema(Arf[10:20; 0:9]) == (Arf(0), Arf(20))
end

@testset "$T" for T in [Arb, Acb]
Expand Down Expand Up @@ -188,21 +182,6 @@
@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)))
@test minimum(Arb[10:20; 0:9]) == Arb(0)
@test maximum(Arb[10:20; 0:9]) == Arb(20)
@test extrema(Arb[10:20; 0:9]) == (Arb(0), Arb(20))
A = [Arb((i, i + 1)) for i = 0:10]
@test Arblib.contains(minimum(A), Arb((0, 1)))
@test Arblib.contains(minimum(reverse(A)), Arb((0, 1)))
@test Arblib.contains(maximum(A), Arb((10, 11)))
@test Arblib.contains(maximum(reverse(A)), Arb((10, 11)))
@test all(Arblib.contains.(extrema(A), (Arb((0, 1)), Arb((10, 11)))))
# These fails with the default implementation of minimum and maximum
@test Arblib.contains(
minimum([Arb((-i, -i + 1)) for i = 0:1000]),
Arb((-1000, -999)),
)
@test Arblib.contains(maximum([Arb((i, i + 1)) for i = 0:1000]), Arb((1000, 1001)))
end

@testset "Acb - specific" begin
Expand Down Expand Up @@ -243,4 +222,91 @@
@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])

# 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

0 comments on commit 5e81504

Please sign in to comment.