Uncaught overflow for gcd(typemin(T), zero(T)) where T<:Signed

[ivirshup]

Hi! Here's a quick example of the issue:

julia> gcd(typemin(Int8), zero(Int8))
-128

Since gcd(typemin(T), typemin(T)) where T<:Signed (per #15228) is an overflow error, and results should be positive I think this is a bug. I'd like to open a PR to fix this (I think an abs to be swapped out with a checked_abs). Does this sound reasonable?

[tvrtadx]

if we return 128.it not a type of Int8 .
So gcd is not a "type-stable" functions

[ivirshup]

For sure. It shouldn't return -128 either, since the return value is expected to be positive. I think gcd(zero(T), typemin(T)) should do the same thing as gcd(typemin(T), typemin(T)), which is currently throwing an OverflowError:

julia> gcd(typemin(Int8), typemin(Int8))
ERROR: OverflowError: gcd(-128, -128) overflows

[tvrtadx]

In general,we don’t check for overflow in integer arithmetic for performance.Such as:

julia> x = Int8(-128)
julia> abs(x) == -x == x
julia> true
julia> gcd(Int8(-128),Int8(0)) == abs(Int8(-128))
julia> true

This is a compromise that has to be done for efficiency.

[ivirshup]

I thought about this, but I think there are two conflicting points. First, consistency. The documentation says the result is positive and the function already checks for overflow in a very similar case:

julia> gcd(typemin(Int8), typemin(Int8))
ERROR: OverflowError: gcd(-128, -128) overflows
Stacktrace:
 [1] __throw_gcd_overflow(::Int8, ::Int8) at ./intfuncs.jl:50
 [2] gcd(::Int8, ::Int8) at ./intfuncs.jl:47
 [3] top-level scope at REPL[4]:1

Second, I don't think it's actually a performance hit, or at least I can't consistently get a measure of one. Here's the definition I was using:

checked_gcd definition

using Base: checked_abs

function checked_gcd(a::T, b::T) where T<:Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128}
    if a == 0
        return checked_abs(b)
    end
    if b == 0
        return checked_abs(a)
    end
    # a == 0 && return abs(b)
    # b == 0 && return abs(a)
    za = trailing_zeros(a)
    zb = trailing_zeros(b)
    k = min(za, zb)
    u = unsigned(abs(a >> za))
    v = unsigned(abs(b >> zb))
    while u != v
        if u > v
            u, v = v, u
        end
        v -= u
        v >>= trailing_zeros(v)
    end
    r = u << k
    # T(r) would throw InexactError; we want OverflowError instead
    r > typemax(T) && __throw_gcd_overflow(a, b)
    r % T
end
@noinline __throw_gcd_overflow(a, b) = throw(OverflowError("gcd($a, $b) overflows"))

Interestingly, elseif instead of two separate if statements seemed to cause a slowdown. There was sometimes a speed up for using an if statement instead of && when gcd was broadcast over arrays.

All variations in minimum time (using @benchmark) were within 2% of each other – except for differences switching between && and if during a broadcasted call.

No code in the normal path changes, so my assumption performance issues would come more from branch prediction. Also if performance was harmed by a check, shouldn't the overflow error from above be removed?

[tvrtadx]

So, you can only add another checked_gcd method.Otherwise,gcd will not get correct result for arrays.

julia> a = Int8[-128,126,12]
julia> gcd(a)
julia> 2

Using your checked_gcd methods in stead of original gcd:

julia>  a = Int8[-128,126,12]
julia> gcd(a)
julia> ERROR: LoadError: OverflowError: checked arithmetic: cannot compute |x| for x = -128::Int8

[ivirshup]

Can't that error already occur from gcd(Int8[-128, -128, 1]) or any leading combination of two or more typemins and zeros?

I'd also agree with just not throwing overflow errors from gcd. But it currently throws one. I just think the behavior should be consistent about whether typemin is a valid result from gcd or not.

[ivirshup]

Thought about this a bit more, and drafted an implementation of gcd(a::AbstractArray{<:Integer}) which should be more stable. This won't throw errors unless the result would actually overflow (i.e. gcd(Int8[-128, -128, 1]) wouldn't error.

function checked_gcd(abc::AbstractArray{T}) where T<:Integer
    idx = 1
    a = zero(T)
    n = length(abc)
    while idx <= n && (abc[idx] == typemin(T) || abc[idx] == zero(T))
        a = min(a, abc[idx])
        idx += 1
    end
    while idx <= n
        a = checked_gcd(a, abc[idx])
        if a == one(T)
            return a
        end
        idx += 1
    end
    return checked_abs(a)
end

[simonbyrne]

To add to this issue: currently gcd(x::Int8, m::Int8) only overflows for (x,m) = (Int8(-128),Int8(-128)), however gcdx(x::Int8, m::Int8) overflows only for (x,m) = (Int8(-128),Int8(-1)) and (x,m) = (Int8(-1),Int8(-128)).

My vague preference is that neither should throw an error, but I don't feel particularly strongly other than that these should be consistent.

[KlausC]

IMO, r = gcd(a::T,b::T) should have following properties:

1. r isa T
2. r >= 0
3. There are integers c and d with r * c == a and r * d == b
4. If a == b == 0 r == 0, otherwise r is the greatest integer with properties EDIT: 2.-3. (not 1.)
   If there is no such r for certain a, b, the call gcd(a,b) should throw an OverflowError.

That corresponds to the documentation of gcd.

I think, we should prefer appropriate checks over performance, as we do in a similar situation comment.

Concerning gcdx, that should consistently throw the same exceptions for the same arguments as gcd.
gcdx(0, 0) == (0, 1, 0).
The julia docu of gcdx describes the expected output correctly.

[KlausC]

The currently implemented algorithms deviate from the previous claims only in these cases:

1. a == b == typemin(T)
2. a == T(-1) and b == typemin(T)
3. a = typemin(T) and b == T(-1)
   for all signed integer bit types. For unsigned and big integers everything looks fine.

[KlausC]

Interestingly, as the domain of gcd and gcdx has been extended to rationals, both functions throw exceptions for a == b == typemin(T) // T(1)

wrap(f, x, y) = try f(x, y); catch ex; return ex; end
function foo(T, rat=false)
         for x in typemin(T):typemax(T)
           for y in typemin(T):typemax(T)
               if rat
                   x = x // T(1)
                   y = y // T(1)
               end
               g = wrap(gcd, x, y)
               gx = wrap(gcdx, x, y)
               if g isa Exception || gx isa Exception
                   println("exception gcd($x,$y) == $g != gcdx() == $gx")
               elseif g >= 0
                   h, u, v = gx
                   h != g && println("failure: gcd($x,$y) == $g != gcdx() == $gx")
               else
                   println("failure: gcd($x,$y) == $g < 0")
               end
           end
         end
       end
foo (generic function with 2 methods)

julia> foo(Int8)
exception gcd(-128,-128) == OverflowError("gcd(-128, -128) overflows") != gcdx() == (-128, 0, -1)
exception gcd(-128,-1) == 1 != gcdx() == DivideError()
failure: gcd(-128,0) == -128 < 0
exception gcd(-1,-128) == 1 != gcdx() == DivideError()
failure: gcd(0,-128) == -128 < 0

julia> foo(Int8, true)
exception gcd(-128//1,-128//1) == OverflowError("gcd(-128, -128) overflows") != gcdx() == OverflowError("gcd(-128, -128) overflows")
exception gcd(-128//1,-1//1) == 1//1 != gcdx() == DivideError()
failure: gcd(-128//1,0//1) == -128//1 < 0
exception gcd(-1//1,-128//1) == 1//1 != gcdx() == DivideError()
failure: gcd(0//1,-128//1) == -128//1 < 0

[nhz2]

This issue has been fixed, I think in #35451