Merge pull request JuliaLang#17615 from eschnett/eschnett/invmod

Correct `invmod` for unsigned integers and for negative moduli

base/docs/helpdb/Base.jl

@@ -6960,7 +6960,9 @@ qr
 """
  invmod(x,m)
 
-Take the inverse of `x` modulo `m`: `y` such that ``xy = 1 \\pmod m``.
+Take the inverse of `x` modulo `m`: `y` such that ``x y = 1 \\pmod m``,
+with ``div(x,y) = 0``. This is undefined for ``m = 0``, or if
+``gcd(x,m) \\neq 1``.
 """
 invmod
 

base/gmp.jl

@@ -260,14 +260,21 @@ for (fJ, fC) in ((:+, :add), (:-,:sub), (:*, :mul),
 end
 
 function invmod(x::BigInt, y::BigInt)
- z = BigInt()
- y = abs(y)
- if y == 1
- return big(0)
+ z = zero(BigInt)
+ ya = abs(y)
+ if ya == 1
+ return z
+ end
+ if (y==0 || ccall((:__gmpz_invert, :libgmp), Cint, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &x, &ya) == 0)
+ throw(DomainError())
  end
- if (y==0 || ccall((:__gmpz_invert, :libgmp), Cint, (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}), &z, &x, &y) == 0)
- error("no inverse exists")
+ # GMP always returns a positive inverse; we instead want to
+ # normalize such that div(z, y) == 0, i.e. we want a negative z
+ # when y is negative.
+ if y < 0
+ z = z + y
  end
+ # The postcondition is: mod(z * x, y) == mod(big(1), m) && div(z, y) == 0
  return z
 end
 

base/intfuncs.jl

@@ -53,7 +53,7 @@ lcm{T<:Integer}(abc::AbstractArray{T}) = reduce(lcm,abc)
 
 Computes the greatest common (positive) divisor of `x` and `y` and their Bézout
 coefficients, i.e. the integer coefficients `u` and `v` that satisfy
-``ux+vy = d = gcd(x,y)``.
+``ux+vy = d = gcd(x,y)``. ``gcdx(x,y)`` returns ``(d,u,v)``.
 
 ```jldoctest
 julia> gcdx(12, 42)
@@ -67,15 +67,20 @@ julia> gcdx(240, 46)
 
 !!! note
  Bézout coefficients are *not* uniquely defined. `gcdx` returns the minimal
- Bézout coefficients that are computed by the extended Euclid algorithm.
- (Ref: D. Knuth, TAoCP, 2/e, p. 325, Algorithm X.) These coefficients `u`
- and `v` are minimal in the sense that ``|u| < |\\frac y d|`` and ``|v| <
- |\\frac x d|``. Furthermore, the signs of `u` and `v` are chosen so that `d`
- is positive.
+ 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| < |y/d|`` and ``|v| < |x/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.
 """
 function gcdx{T<:Integer}(a::T, b::T)
+ # a0, b0 = a, b
  s0, s1 = one(T), zero(T)
  t0, t1 = s1, s0
+ # The loop invariant is: s0*a0 + t0*b0 == a
  while b != 0
  q = div(a, b)
  a, b = b, rem(a, b)
@@ -87,13 +92,16 @@ end
 gcdx(a::Integer, b::Integer) = gcdx(promote(a,b)...)
 
 # multiplicative inverse of n mod m, error if none
-function invmod(n, m)
+function invmod{T<:Integer}(n::T, m::T)
  g, x, y = gcdx(n, m)
- if g != 1 || m == 0
- error("no inverse exists")
- end
- x < 0 ? abs(m) + x : x
+ (g != 1 || m == 0) && throw(DomainError())
+ # Note that m might be negative here.
+ # For unsigned T, x might be close to typemax; add m to force a wrap-around.
+ r = mod(x + m, m)
+ # The postcondition is: mod(r * n, m) == mod(T(1), m) && div(r, m) == 0
+ r
 end
+invmod(n::Integer, m::Integer) = invmod(promote(n,m)...)
 
 # ^ for any x supporting *
 to_power_type(x::Number) = oftype(x*x, x)

doc/stdlib/math.rst

@@ -1266,7 +1266,7 @@ Mathematical Functions
 
  .. Docstring generated from Julia source
 
- Computes the greatest common (positive) divisor of ``x`` and ``y`` and their Bézout coefficients, i.e. the integer coefficients ``u`` and ``v`` that satisfy :math:`ux+vy = d = gcd(x,y)`\ .
+ Computes the greatest common (positive) divisor of ``x`` and ``y`` and their Bézout coefficients, i.e. the integer coefficients ``u`` and ``v`` that satisfy :math:`ux+vy = d = gcd(x,y)`\ . :math:`gcdx(x,y)` returns :math:`(d,u,v)`\ .
 
  .. doctest::
 
@@ -1279,7 +1279,7 @@ Mathematical Functions
  (2,-9,47)
 
  .. note::
- Bézout coefficients are *not* uniquely defined. ``gcdx`` returns the minimal Bézout coefficients that are computed by the extended Euclid algorithm. (Ref: D. Knuth, TAoCP, 2/e, p. 325, Algorithm X.) These coefficients ``u`` and ``v`` are minimal in the sense that :math:`|u| < |\frac y d|` and :math:`|v| < |\frac x d|`\ . Furthermore, the signs of ``u`` and ``v`` are chosen so that ``d`` is positive.
+ 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 :math:`|u| < |y/d|` and :math:`|v| < |x/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.
 
 
 .. function:: ispow2(n) -> Bool
@@ -1322,7 +1322,7 @@ Mathematical Functions
 
  .. Docstring generated from Julia source
 
- Take the inverse of ``x`` modulo ``m``\ : ``y`` such that :math:`xy = 1 \pmod m`\ .
+ Take the inverse of ``x`` modulo ``m``\ : ``y`` such that :math:`x y = 1 \pmod m`\ , with :math:`div(x,y) = 0`\ . This is undefined for :math:`m = 0`\ , or if :math:`gcd(x,m) \neq 1`\ .
 
 .. function:: powermod(x, p, m)
 

test/examples.jl

@@ -18,7 +18,7 @@ x = ModInts.ModInt{256}(13)
 y = inv(x)
 @test y == ModInts.ModInt{256}(197)
 @test x*y == ModInts.ModInt{256}(1)
-@test_throws ErrorException inv(ModInts.ModInt{8}(4))
+@test_throws DomainError inv(ModInts.ModInt{8}(4))
 
 include(joinpath(dir, "ndgrid.jl"))
 r = repmat(1:10,1,10)

test/intfuncs.jl

@@ -35,10 +35,13 @@ end
 @test gcdx(5, -12) == (1, 5, 2)
 @test gcdx(-25, -4) == (1, -1, 6)
 
-@test invmod(6, 31) == 26
-@test invmod(-1, 3) == 2
-@test invmod(-1, -3) == 2
-@test_throws ErrorException invmod(0, 3)
+@test invmod(6, 31) === 26
+@test invmod(-1, 3) === 2
+@test invmod(1, -3) === -2
+@test invmod(-1, -3) === -1
+@test invmod(0x2, 0x3) === 0x2
+@test invmod(2, 0x3) === 2
+@test_throws DomainError invmod(0, 3)
 
 @test powermod(2, 3, 5) == 3
 @test powermod(2, 3, -5) == -2

test/numbers.jl

@@ -2167,11 +2167,12 @@ for T in (Int32,Int64), ii = -20:20, jj = -20:20
  @test d == gcd(ib,jb)
  @test lcm(i,j) == lcm(ib,jb)
  @test gcdx(i,j) == gcdx(ib,jb)
- if j == 0
- @test_throws ErrorException invmod(i,j)
- @test_throws ErrorException invmod(ib,jb)
- elseif d == 1
+ if j == 0 || d != 1
+ @test_throws DomainError invmod(i,j)
+ @test_throws DomainError invmod(ib,jb)
+ else
  n = invmod(i,j)
+ @test div(n, j) == 0
  @test n == invmod(ib,jb)
  @test mod(n*i,j) == mod(1,j)
  end