Merge pull request #4811 from stevengj/invmod

fix gcd & lcm sign and optimize gcdx and invmod

Author: JeffBezanson

base/gmp.jl

@@ -6,7 +6,7 @@ import Base: *, +, -, /, <, <<, >>, >>>, <=, ==, >, >=, ^, (~), (&), (|), ($),
              binomial, cmp, convert, div, divrem, factorial, fld, gcd, gcdx, lcm, mod,
              ndigits, promote_rule, rem, show, isqrt, string, isprime, powermod,
              widemul, sum, trailing_zeros, trailing_ones, count_ones, base, parseint,
-             serialize, deserialize, bin, oct, dec, hex, isequal
+             serialize, deserialize, bin, oct, dec, hex, isequal, invmod
 
 type BigInt <: Integer
     alloc::Cint
@@ -152,7 +152,7 @@ deserialize(s, ::Type{BigInt}) = parseint(BigInt, deserialize(s), 62)
 # Binary ops
 for (fJ, fC) in ((:+, :add), (:-,:sub), (:*, :mul),
                  (:fld, :fdiv_q), (:div, :tdiv_q), (:mod, :fdiv_r), (:rem, :tdiv_r),
-                 (:gcd, :gcd), (:lcm, :lcm),
+                 (:gcd, :gcd), (:lcm, :lcm), (:invmod, :invert),
                  (:&, :and), (:|, :ior), (:$, :xor))
     @eval begin
         function ($fJ)(x::BigInt, y::BigInt)
@@ -314,6 +314,9 @@ powermod(x::BigInt, p::Integer, m::BigInt) = powermod(x, BigInt(p), m)
 powermod(x::BigInt, p::Integer, m::Integer) = powermod(x, BigInt(p), BigInt(m))
 
 function gcdx(a::BigInt, b::BigInt)
+    if b == 0 # shortcut this to ensure consistent results with gcdx(a,b)
+        return a < 0 ? (-a,-one(BigInt),zero(BigInt)) : (a,one(BigInt),zero(BigInt))
+    end
     g = BigInt()
     s = BigInt()
     t = BigInt()

base/intfuncs.jl

@@ -33,16 +33,16 @@ abs(x::Signed) = flipsign(x,x)
 ## number-theoretic functions ##
 
 function gcd{T<:Integer}(a::T, b::T)
-    neg = a < 0
     while b != 0
         t = b
         b = rem(a, b)
         a = t
     end
-    g = abs(a)
-    neg ? -g : g
+    abs(a)
 end
-lcm{T<:Integer}(a::T, b::T) = a * div(b, gcd(b,a))
+
+# explicit a==0 test is to handle case of lcm(0,0) correctly
+lcm{T<:Integer}(a::T, b::T) = a == 0 ? a : abs(a * div(b, gcd(b,a)))
 
 gcd(a::Integer) = a
 lcm(a::Integer) = a
@@ -52,21 +52,23 @@ gcd(a::Integer, b::Integer...) = gcd(a, gcd(b...))
 lcm(a::Integer, b::Integer...) = lcm(a, lcm(b...))
 
 # return (gcd(a,b),x,y) such that ax+by == gcd(a,b)
-function gcdx(a, b)
-    if b == 0
-        (a, 1, 0)
-    else
-        m = rem(a, b)
-        k = div((a-m), b)
-        (g, x, y) = gcdx(b, m)
-        (g, y, x-k*y)
+function gcdx{T<:Integer}(a::T, b::T)
+    s0, s1 = one(T), zero(T)
+    t0, t1 = s1, s0
+    while b != 0
+        q = div(a, b)
+        a, b = b, rem(a, b)
+        s0, s1 = s1, s0 - q*s1
+        t0, t1 = t1, t0 - q*t1
     end
+    a < 0 ? (-a, -s0, -t0) : (a, s0, t0)
 end
+gcdx(a::Integer, b::Integer) = gcdx(promote(a,b)...)
 
-# multiplicative inverse of x mod m, error if none
+# multiplicative inverse of n mod m, error if none
 function invmod(n, m)
     g, x, y = gcdx(n, m)
-    g != 1 ? error("no inverse exists") : (x < 0 ? m + x : x)
+    g == 1 ? (x < 0 ? abs(m) + x : x) : error("no inverse exists")
 end
 
 # ^ for any x supporting *

base/rational.jl

@@ -6,7 +6,7 @@ immutable Rational{T<:Integer} <: Real
         if num == 0 && den == 0
             error("invalid rational: 0//0")
         end
-        g = gcd(den, num)
+        g = den < 0 ? -gcd(den,num) : gcd(den, num)
         new(div(num, g), div(den, g))
     end
 end

doc/helpdb.jl

@@ -3885,19 +3885,19 @@ popdisplay(d::Display)
 
 ("Mathematical Functions","Base","gcd","gcd(x, y)
 
-   Greatest common divisor
+   Greatest common (positive) divisor (or zero if x and y are both zero).
 
 "),
 
 ("Mathematical Functions","Base","lcm","lcm(x, y)
 
-   Least common multiple
+   Least common (non-negative) multiple.
 
 "),
 
 ("Mathematical Functions","Base","gcdx","gcdx(x, y)
 
-   Greatest common divisor, also returning integer coefficients \"u\"
+   Greatest common (positive) divisor, also returning integer coefficients \"u\"
    and \"v\" that solve \"ux+vy == gcd(x,y)\"
 
 "),

doc/stdlib/base.rst

@@ -2632,15 +2632,15 @@ Mathematical Functions
 
 .. function:: gcd(x,y)
 
-   Greatest common divisor
+   Greatest common (positive) divisor (or zero if x and y are both zero).
 
 .. function:: lcm(x,y)
 
-   Least common multiple
+   Least common (non-negative) multiple.
 
 .. function:: gcdx(x,y)
 
-   Greatest common divisor, also returning integer coefficients ``u`` and ``v`` that solve ``ux+vy == gcd(x,y)``
+   Greatest common (positive) divisor, also returning integer coefficients ``u`` and ``v`` that solve ``ux+vy == gcd(x,y)``
 
 .. function:: ispow2(n)
 

test/numbers.jl

@@ -1511,3 +1511,16 @@ end
 # overflow in rational comparison
 @test 3//2 < typemax(Int)
 @test 3//2 <= typemax(Int)
+
+# check gcd and related functions against GMP
+for i = -20:20, j = -20:20
+    local d = gcd(i,j)
+    @test d >= 0
+    @test lcm(i,j) >= 0
+    local ib = big(i)
+    local jb = big(j)
+    @test d == gcd(ib,jb)
+    @test lcm(i,j) == lcm(ib,jb)
+    @test gcdx(i,j) == gcdx(ib,jb)
+    @test d != 1 || invmod(i,j) == invmod(ib,jb)
+end