Skip to content

Commit

Permalink
add vararg gcdx(a, b, c...)
Browse files Browse the repository at this point in the history
  • Loading branch information
thchr committed Sep 30, 2024
1 parent 17445fe commit a6f0962
Show file tree
Hide file tree
Showing 2 changed files with 40 additions and 2 deletions.
26 changes: 24 additions & 2 deletions base/intfuncs.jl
Original file line number Diff line number Diff line change
Expand Up @@ -165,24 +165,34 @@ end

# return (gcd(a, b), x, y) such that ax+by == gcd(a, b)
"""
gcdx(a, b)
gcdx(a, b...)
Computes the greatest common (positive) divisor of `a` and `b` and their Bézout
coefficients, i.e. the integer coefficients `u` and `v` that satisfy
``ua+vb = d = gcd(a, b)``. ``gcdx(a, b)`` returns ``(d, u, v)``.
``u*a + v*b = d = gcd(a, b)``. ``gcdx(a, b)`` returns ``(d, u, v)``.
For more arguments than two, i.e., `gcdx(a, b, c, ...)` the Bézout coefficients are computed
recursively, returning a solution `(d, u, v, w, ...)` to
``u*a + v*b + w*c + ... = d = gcd(a, b, c, ...)``.
The arguments may be integer and rational numbers.
!!! compat "Julia 1.4"
Rational arguments require Julia 1.4 or later.
!!! compat "Julia 1.11"
More or fewer arguments than two require Julia 1.12 or later.
# Examples
```jldoctest
julia> gcdx(12, 42)
(6, -3, 1)
julia> gcdx(240, 46)
(2, -9, 47)
julia> gcdx(15, 12, 20)
(1, 7, -7, -1)
```
!!! note
Expand Down Expand Up @@ -215,6 +225,18 @@ Base.@assume_effects :terminates_locally function gcdx(a::Integer, b::Integer)
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

# multiplicative inverse of n mod m, error if none

Expand Down
16 changes: 16 additions & 0 deletions test/rational.jl
Original file line number Diff line number Diff line change
Expand Up @@ -702,6 +702,22 @@ end
end
end

@testset "gcdx for 1 and 3+ arguments" begin
# one-argument
@test gcdx(7) == (7, 1)
@test gcdx(-7) == (7, -1)
@test gcdx(1//4) == (1//4, 1)

# 3+ arguments
@test gcdx(2//3) == gcdx(2//3) == (2//3, 1)
@test gcdx(15, 12, 20) == (1, 7, -7, -1)
@test gcdx(60//4, 60//5, 60//3) == (1//1, 7, -7, -1)
abcd = (105, 1638, 2145, 3185)
d, uvwp... = gcdx(abcd...)
@test d == sum(abcd .* uvwp) # u*a + v*b + w*c + p*d == gcd(a, b, c, d)
@test (@inferred gcdx(1, 2, 3, 4, 5, 6, 7, 8, 9, 10)) isa NTuple{11, Int}
end

@testset "Binary operations with Integer" begin
@test 1//2 - 1 == -1//2
@test -1//2 + 1 == 1//2
Expand Down

0 comments on commit a6f0962

Please sign in to comment.