clear up several more inconsistencies among gcd , gcdx and invmod (#59678)

test suite of invariances satisfied among `gcd`, `gcdx`, `invmod`:

```
function can_compute_gcd(::Type{T}, a, b) where {T}
    (a > b) && return can_compute_gcd(T, b, a)
    if T <: Signed
        return a != typemin(T) || (b != typemin(T) && b != 0)
    else
        return (a > typemin(a)) || iszero(a)
    end
end

check_bezout(T::Type{U}, a, b, d, u, v) where U = u * a + v * b == d
function check_bezout(T::Type{S}, a, b, d, u, v) where {S<:Signed}
    W = widen(T)
    return W(u) * W(a) + W(v) * W(b) == W(d)
end

bad_values = Tuple{Integer, Integer}[]

let
    Ts_signed   = (Int8, Int16, Int32, Int64)
    Ts_unsigned = (UInt8, UInt16, UInt32, UInt64)
    Ts = (Ts_signed..., Ts_unsigned...)

    T = UInt128
    test_values = T[]
    nbits = (sizeof(T) << 3)
    append!(test_values, [-one(T), zero(T), one(T), typemin(T), typemax(T)])
    append!(test_values, reduce(vcat, [
        [T(1) << i, (T(1) << i) - 1, (T(1) << i) + 1]
        for i in 0:nbits
    ]))
    append!(test_values, reduce(vcat, [
        [T(3)^div(i,2), T(3)^div(i,2) - 1, T(3)^div(i,2) + 1]
        for i in 0:2:nbits
    ]))

    function do_tests(a, b)
        T = Base.promote_typeof(a, b)
        if !can_compute_gcd(T, a, b)
            return
        else
            try
                gcd(a, b)
                gcdx(a, b)
            catch
                push!(bad_values, (a, b))
                return
            end
        end

        g = gcd(a, b)
        d, u, v = gcdx(a, b)

        properties = (
            ispositive(g) &&
            g === d &&
            (typeof(g) == typeof(d)) &&
            check_bezout(T, a, b, d, u, v) &&
            mod(a, g) == mod(b, g) == 0
        )

        if isone(g) && !iszero(b)
            ai = invmod(a, b)
            properties = properties && (
                (mod(one(b), b) == mod(widemul(a, ai), b)) &&
                ((sign(ai) == sign(b)) || (sign(ai) == 0)) &&
                (ai === mod(ai, b)) &&
                (typeof(ai) == typeof(mod(a, b))) &&
                iszero(div(ai, b))
            )
        end
        if !properties
            println("$(typeof(a))($a), $(typeof(b))($b)")
            push!(bad_values, (a, b))
        end
    end

    for T1 in Ts
        As = test_values .% T1
        for T2 in Ts
            Bs = test_values .% T2

            for a in As
                for b in Bs
                    iszero(a) && iszero(b) && continue
                    do_tests(a, b)
                end
            end
        end
    end
end
```
---------

Co-authored-by: Max Horn <max@quendi.de>
Co-authored-by: Sukera <11753998+Seelengrab@users.noreply.github.com>

Author: 3 people

base/intfuncs.jl

@@ -143,13 +143,19 @@ function lcm(a::T, b::T) where T<:Integer
     end
 end
 
+function _promote_mixed_signs(a::Signed, b::Unsigned)
+    # handle the case a == typemin(typeof(a)) if R != typeof(a)
+    R = promote_typeof(a, b)
+    promote(abs(a % signed(R)), b)
+end
+
 gcd(a::Integer) = checked_abs(a)
 gcd(a::Rational) = checked_abs(a.num) // a.den
 lcm(a::Union{Integer,Rational}) = gcd(a)
-gcd(a::Unsigned, b::Signed) = gcd(promote(a, abs(b))...)
-gcd(a::Signed, b::Unsigned) = gcd(promote(abs(a), b)...)
+gcd(a::Unsigned, b::Signed) = gcd(b, a)
+gcd(a::Signed, b::Unsigned) = gcd(_promote_mixed_signs(a, b)...)
 lcm(a::Unsigned, b::Signed) = lcm(promote(a, abs(b))...)
-lcm(a::Signed, b::Unsigned) = lcm(promote(abs(a), b)...)
+lcm(a::Signed, b::Unsigned) = lcm(_promote_mixed_signs(a, b)...)
 gcd(a::Real, b::Real) = gcd(promote(a,b)...)
 lcm(a::Real, b::Real) = lcm(promote(a,b)...)
 gcd(a::Real, b::Real, c::Real...) = gcd(a, gcd(b, c...))
@@ -223,22 +229,31 @@ julia> gcdx(15, 12, 20)
     their `typemax`, and the identity then holds only via the unsigned
     integers' modulo arithmetic.
 """
-Base.@assume_effects :terminates_locally function gcdx(a::Integer, b::Integer)
-    T = promote_type(typeof(a), typeof(b))
-    a == b == 0 && return (zero(T), zero(T), zero(T))
+Base.@assume_effects :terminates_locally function gcdx(a::T, b::T) where {T<:Integer}
+    if iszero(a) && iszero(b)
+        return (zero(T), zero(T), zero(T))
+    elseif isone(abs(b))
+        # handles (typemin(::Signed), -1)
+        return (one(T), zero(T), b)
+    elseif isone(abs(a))
+        return (one(T), a, zero(T))
+    end
     # a0, b0 = a, b
     s0, s1 = oneunit(T), zero(T)
     t0, t1 = s1, s0
     # The loop invariant is: s0*a0 + t0*b0 == a && s1*a0 + t1*b0 == b
-    x = a % T
-    y = b % T
-    while y != 0
-        q, r = divrem(x, y)
-        x, y = y, r
+    while !iszero(b)
+        q, r = divrem(a, b)
+        a, b = b, r
         s0, s1 = s1, s0 - q*s1
         t0, t1 = t1, t0 - q*t1
     end
-    x < 0 ? (-x, -s0, -t0) : (x, s0, t0)
+    # for cases like abs(Int8(-128))
+    if isnegative(a) && isnegative(abs(a))
+        throw(DomainError((a, b), LazyString("gcd not representable in ", T)))
+    else
+        return isnegative(a) ? (abs(a), -s0, -t0) : (a, s0, t0)
+    end
 end
 gcdx(a::Real, b::Real) = gcdx(promote(a,b)...)
 gcdx(a::T, b::T) where T<:Real = throw(MethodError(gcdx, (a,b)))
@@ -254,11 +269,12 @@ function gcdx(a::Real, b::Real, cs::Real...)
     d′, x, ys... = gcdx(d, cs...)
     return d′, i*x, j*x, ys...
 end
+
 function gcdx(a::Signed, b::Unsigned)
-    R = promote_type(typeof(a), typeof(b))
-    _a = a % signed(R) # handle the case a == typemin(typeof(a)) if R != typeof(a)
-    d, u, v = gcdx(promote(abs(_a), b)...)
-    d, flipsign(u, a), v
+    R = promote_typeof(a, b)
+    d, u, v = gcdx(promote(abs(a % signed(R)), b)...)
+    flip_typemin = isnegative(a) & (R <: Signed)
+    d, flipsign(u, a - flip_typemin), v
 end
 function gcdx(a::Unsigned, b::Signed)
     d, v, u = gcdx(b, a)
@@ -287,24 +303,38 @@ julia> invmod(5, 6)
 ```
 """
 function invmod(n::Integer, m::Integer)
+    # The postcondition is: mod(widemul(result, n), m) == mod(one(T), m) && iszero(div(result, m))
     iszero(m) && throw(DomainError(m, "`m` must not be 0."))
-    if n isa Signed && hastypemax(typeof(n))
-        # work around inconsistencies in gcdx
-        # https://github.com/JuliaLang/julia/issues/33781
-        T = promote_type(typeof(n), typeof(m))
-        n == typemin(typeof(n)) && m == typeof(n)(-1) && return T(0)
-        n == typeof(n)(-1) && m == typemin(typeof(n)) && return T(-1)
+    R = promote_typeof(n, m)
+    if R <: Signed
+        x = _bezout_coef(n, m)
+        return mod(x, m)
+    else
+        S = signed(R)
+        if !hastypemax(S) || (n <= typemax(S)) && (m <= typemax(S))
+            x = _bezout_coef(n % S, m % S)
+
+            # this branch is only hit if R <: Unsigned, so we don't have
+            # to worry about abs(typemin(::Signed)) overflow. If `m` is
+            # signed then `x` must be unsigned, and thus never negative
+            isnegative(x) && (x += abs(m))
+            return mod(x % R, m)
+        else
+            # since gcdx only promises bezout w.r.t overflow for unsigned ints,
+            # we have to widen to a signed type
+            W = widen(S)
+            x = _bezout_coef(n % W, m % W)
+            t = mod(x, m % W)
+            isnegative(m) && (t -= m)
+            return mod(t % R, m)
+        end
     end
-    g, x, y = gcdx(n, m)
+end
+
+function _bezout_coef(n, m)
+    g, x, _ = gcdx(n, m)
     g != 1 && throw(DomainError((n, m), LazyString("Greatest common divisor is ", g, ".")))
-    # Note that m might be negative here.
-    if x isa Unsigned && hastypemax(typeof(x)) && x > typemax(x)>>1
-        # x might have wrapped if it would have been negative
-        # adding back m forces a correction
-        x += m
-    end
-    # The postcondition is: mod(result * n, m) == mod(T(1), m) && div(result, m) == 0
-    return mod(x, m)
+    return x
 end
 
 """

test/intfuncs.jl

@@ -217,6 +217,19 @@ end
     d, u, v = gcdx(x, y)
     @test x*u + y*v == d
 
+    for T in (Int8, Int16, Int32, Int64, Int128)
+        @test_throws DomainError gcdx(typemin(T), typemin(T))
+        @test_throws DomainError gcdx(typemin(T), T(0))
+        @test_throws DomainError gcdx(T(0), typemin(T))
+        d, u, v = gcdx(typemin(T), T(-1))
+        @test d == T(1)
+        @test typemin(T) * u + T(-1) * v == T(1)
+        @test gcdx(T(-1), typemin(T)) == (d, v, u)
+        d, u, v = gcdx(typemin(T), T(1))
+        @test d == T(1)
+        @test typemin(T) * u + T(1) * v == T(1)
+        @test gcdx(T(1), typemin(T)) == (d, v, u)
+    end
 end
 
 # issue #58025
@@ -244,7 +257,7 @@ end
 
     @test gcdx(Int16(-32768), Int8(-128)) === (Int16(128), Int16(0), Int16(-1))
     @test gcdx(Int8(-128), UInt16(256)) === (0x0080, 0xffff, 0x0000)
-    @test_broken gcd(Int8(-128), UInt16(256)) === 0x0080
+    @test gcd(Int8(-128), UInt16(256)) === 0x0080
 end
 
 @testset "gcd/lcm/gcdx for custom types" begin
@@ -294,10 +307,20 @@ end
     # Verify issue described in PR 58010 is fixed
     @test invmod(UInt8(3), UInt16(50000)) === 0x411b
 
+    @test invmod(0x00000001, Int8(-128)) === Int32(-127)
+    @test invmod(0xffffffff, Int8(-38)) === Int32(-15)
+    @test invmod(Int8(-1), 0xffffffff) === 0xfffffffe
+    @test invmod(Int32(-1), typemin(Int64)) === Int64(-1)
+    @test invmod(0x3e81, Int16(-5716)) === Int16(-2407)
+
     for T in (Int8, UInt8)
         for x in typemin(T):typemax(T)
             for m in typemin(T):typemax(T)
-                if m != 0 && try gcdx(x, m)[1] == 1 catch _ true end
+                if !(
+                    iszero(m) ||
+                    iszero(mod(x, m)) && !isone(abs(m)) ||
+                    !isone(gcd(x, m))
+                )
                     y = invmod(x, m)
                     @test mod(widemul(y, x), m) == mod(1, m)
                     @test div(y, m) == 0