Repository: JuliaLang/julia
File: base/intfuncs.jl (Julia)

Code before:

!!! note
    Bézout coefficients are *not* uniquely defined. `gcdx` returns the minimal
    Bézout coefficients that are computed by the extended Euclidean algorithm.
    (Ref: D. Knuth, TAoCP, 2/e, p. 325, Algorithm X.)
    For signed integers, these coefficients `u` and `v` are minimal in
    the sense that ``|u| < |b/d|`` and ``|v| < |a/d|``. Furthermore,
    the signs of `u` and `v` are chosen so that `d` is positive.
    For unsigned integers, the coefficients `u` and `v` might be near
    their `typemax`, and the identity then holds only via the unsigned
    integers' modulo arithmetic.
"""
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
    while b != 0
        q, r = divrem(a, b)
        a, b = b, r
        s0, s1 = s1, s0 - q*s1
        t0, t1 = t1, t0 - q*t1
    end
    a < 0 ? (checked_abs(a), -s0, -t0) : (a, s0, t0)
end
gcdx(a::Real, b::Real) = gcdx(promote(a,b)...)
gcdx(a::T, b::T) where T<:Real = throw(MethodError(gcdx, (a,b)))
gcdx(a::Real) = (gcd(a), signbit(a) ? -one(a) : one(a))
function gcdx(a::Real, b::Real, cs::Real...)
    # a solution to the 3-arg `gcdx(a,b,c)` problem, `u*a + v*b + w*c = gcd(a,b,c)`, can be
    # obtained from the 2-arg problem in three steps:
    #   1. `gcdx(a,b)`: solve `i*a + j*b = d′ = gcd(a,b)` for `(i,j)`
    #   2. `gcdx(d′,c)`: solve `x*gcd(a,b) + yc = gcd(gcd(a,b),c) = gcd(a,b,c)` for `(x,y)`
    #   3. return `d = gcd(a,b,c)`, `u = i*x`, `v = j*x`, and `w = y`
    # the N-arg solution proceeds similarly by recursion
    d, i, j = gcdx(a, b)
    d′, x, ys... = gcdx(d, cs...)
    return d′, i*x, j*x, ys...
end

gcdx(a::Integer, b::Integer) = gcdx(promote(a, b)...)
function gcdx(a::Signed, b::Unsigned)
    d, u, v = gcdx(_promote_mixed_signs(a, b)...)

Review comment (suggestion):
Wouldn't it be better / fit more with the rest of the code to use this:
```suggestion
    while !iszero(b)
```

Code after:

!!! note
    Bézout coefficients are *not* uniquely defined. `gcdx` returns the minimal
    Bézout coefficients that are computed by the extended Euclidean algorithm.
    (Ref: D. Knuth, TAoCP, 2/e, p. 325, Algorithm X.)
    For signed integers, these coefficients `u` and `v` are minimal in
    the sense that ``|u| < |b/d|`` and ``|v| < |a/d|``. Furthermore,
    the signs of `u` and `v` are chosen so that `d` is positive.
    For unsigned integers, the coefficients `u` and `v` might be near
    their `typemax`, and the identity then holds only via the unsigned
    integers' modulo arithmetic.
"""
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
    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
    # 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)))
gcdx(a::Real) = (gcd(a), signbit(a) ? -one(a) : one(a))
function gcdx(a::Real, b::Real, cs::Real...)
    # a solution to the 3-arg `gcdx(a,b,c)` problem, `u*a + v*b + w*c = gcd(a,b,c)`, can be
    # obtained from the 2-arg problem in three steps:
    #   1. `gcdx(a,b)`: solve `i*a + j*b = d′ = gcd(a,b)` for `(i,j)`
    #   2. `gcdx(d′,c)`: solve `x*gcd(a,b) + yc = gcd(gcd(a,b),c) = gcd(a,b,c)` for `(x,y)`
    #   3. return `d = gcd(a,b,c)`, `u = i*x`, `v = j*x`, and `w = y`
    # the N-arg solution proceeds similarly by recursion
    d, i, j = gcdx(a, b)
    d′, x, ys... = gcdx(d, cs...)
    return d′, i*x, j*x, ys...