make gcdx type-stable & 20x faster, fix invmod to work with negative …

…arguments, use GMP for BigInt invmod

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)

base/intfuncs.jl

@@ -52,21 +52,24 @@ 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, 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 ? m + x : x) : 
+    g == -1 ? (x > 0 ? abs(m) - x : -x) : error("no inverse exists")
 end
 
 # ^ for any x supporting *