`isapprox` can give absurd result in case of integer overflow

[nalimilan]

I was surprised to realize this:

julia> Int8(-68) ≈ Int8(60)
true

This is because of the way the isapprox fallback is defined, and more precisely to:

julia> abs(Int8(-68) - Int8(60))
-128

Shouldn't we define a custom method for integers that is robust to overflow?

[oscardssmith]

yes

[nalimilan]

OK. I'm not sure how to handle overflow though. Using widen is the easy way, but the performance impact will be large for Int64/UInt64 and even more Int128/UInt128. The best I could find is a hybrid approach which simply calls == in the most common (default) case and relies on widen for other situations:

function isapprox(x::Integer, y::Integer;
                  atol::Real=0, rtol::Real=rtoldefault(x,y,atol),
                  nans::Bool=false, norm::Function=abs)
    if norm === abs && atol < 1 && rtol == 0
        return x == y
    else
        return norm(widen(x) - widen(y)) <= max(atol, rtol*max(norm(widen(x)), norm(widen(y))))
    end
end

Unfortunately not all integer types define widen (notably not Bool) so that method would have to be used only for types defined in Base.

This method will also have to be used/adapted for cases mixing integers and floats. That case is more complex as rtol will be nonzero by default. For performance, an intermediate branch could be added to check whether x - y will overflow and whether x == typemin(x) or y == typemin(y) (which makes abs return a negative value). If that's not the case (99% of cases) we can call the formula without widen. (An alternative solution would be to call float on the integer value, but for Int64/UInt64 and Int128/UInt128 this would give incorrect results for large numbers above maxintfloat(Float64), which is worse than the current overflow risk.)

[oscardssmith]

can't you just check use a subtract with overflow checking?

[nalimilan]

Well that's no enough. First, norm(x) can also pose problems. For the default norm=abs we can check whether x == typemin(x). But that doesn't work for arbitrary user functions. Maybe we can say that it's the user's responsibility to pass a function that checks for overflow?

The second issue is what to do if overflow happens. In general that means that x and y are not approximately equal, but nothing prevents the user from passing atol=typemax(Int) or rtol=100, so for full correctness we have to compute the full formula in a wider type AFAICT, at least in a fallback branch which is almost never taken.

[vtjnash]

I guess part of the confusion here is that (except for exactly typemin as shown in the first answer) most of the answers actually are close, in the mod Int8 arithmetic base that is being expressed here.

julia> isapprox(Int8(127), Int8(-118), atol=10)
false

julia> isapprox(Int8(127), Int8(-119), atol=10)
true

[nalimilan]

Yes, but do we really want isapprox to use modular arithmetic? That sounds confusing, and for the most common case we can implement a solution robust to overflow with a very small performance cost.

[LilithHafner]

Here's another overflow:

julia> isapprox(UInt8(60), UInt8(61), atol=1)
false

julia> isapprox(UInt8(60), UInt8(59), atol=1)
true

Technically, the overflow behavior described in this thread is correct because the docstring defines the semantics in terms of the implementation:

isapprox returns true if norm(x-y) <= max(atol, rtol*max(norm(x), norm(y)))

But I still think some of this behavior is bad enough to warrant changing (and changing the documentation to reflect that)

[nalimilan]

#50730 does fix the case described in the last comment as it only covers atol < 1. I can try proposing more code to cover this kind of overflow, but it seems almost impossible to cover all situations, notably when mixing signed and unsigned.