fix powermod problems, add test coverage

base/gmp.jl

@@ -307,7 +307,7 @@ function powermod(x::BigInt, p::BigInt, m::BigInt)
     ccall((:__gmpz_powm, :libgmp), Void,
           (Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}, Ptr{BigInt}),
           &r, &x, &p, &m)
-    return r
+    return m < 0 && r > 0 ? r + m : r # choose sign conistent with mod(x^p, m)
 end
 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))

base/intfuncs.jl

@@ -115,7 +115,9 @@ end
 # b^p mod m
 function powermod{T}(b::Integer, p::Integer, m::T)
     p < 0 && throw(DomainError())
-    p == 0 && return one(b)
+    p == 0 && return m == 0 ? throw(DivideError()) : m < 0 ? 
+    (m == -1 ? zero(promote_type(typeof(b),typeof(m))) : one(b)+m) :
+    (m == 1 ? zero(promote_type(typeof(b),typeof(m))) : one(promote_type(typeof(b),typeof(m))))
     b = oftype(m,mod(b,m))
     t = prevpow2(p)
     local r::T

test/numbers.jl

@@ -1524,3 +1524,11 @@ for i = -20:20, j = -20:20
     @test gcdx(i,j) == gcdx(ib,jb)
     @test d != 1 || invmod(i,j) == invmod(ib,jb)
 end
+
+# check powermod function against GMP
+for i = -10:10, p = 0:5, m = -10:10
+    if m != 0
+        @test powermod(i,p,m) == powermod(i,p,big(m)) == powermod(big(i),big(p),big(m))
+        @test mod(i^p,m) == powermod(i,p,m) == mod(big(i)^p,big(m))
+    end
+end