fix gcd & lcm sign and optimize gcdx and invmod

[stevengj]

The attached patch fixes several bugs (including #4810):

• gcd, lcm, and gcdx, returned inconsistent signs between BigInt and other integer types. Now they always return a positive gcd and lcm.
• invmod was broken for negative arguments, since it was incorrectly interpreting a gcd of -1 as meaning that the inverse did not exist.
• Rational construction was broken for BigInt because of the gcd sign inconsistency.
• gcdx was not type-stable (it promoted a result to Int in one case).

I also rewrote gcdx non-recursively to make it about 20x faster on my machine, and use the GMP invmod function directly for BigInt.

[StefanKarpinski]

So my take on this is that we should cherry-pick the rational constructor fix for 0.2 and wait on the rest.

[ivarne]

You should also document the behaviour of gcd to always return a positive number.

[JeffBezanson]

I'd prefer to merge this whole PR. In any case, gcd was not documented in enough detail before.

[stevengj]

Another omission is that we don't seem to have any test coverage for gcdx etc. Should there be a new test/intfuncs.jl file or would this go into one of the existing test files?

[StefanKarpinski]

I guess we could merge this, although I worry that someone else may already be depending on the old sign behavior of gcd and lcm.

[stevengj]

I'm adding a test case to test/numbers.jl, and found another bug: lcm(0,0) throws an exception, as opposed to returning the correct result (0). Patch forthcoming.

[stevengj]

@StefanKarpinski, since the old sign behavior was undocumented, it's not unreasonable to treat it as a bug and to say that people who relied on it are on their own. (In practice, most people doing gcd-like things are probably working with nonnegative integers, though.)

[stevengj]

Hmm, another inconsistency: powermod(3,7,-5) gives -3, but with BigInt it gives 2.

Looks like this can be fixed simply by using rem instead of mod in the powermod implementation.

[stevengj]

More problems with powermod: powermod(1,0,1) gives 1, not zero, and powermod(1,0,0) also gives 1 rather than throwing an exception. The lack of systematic test coverage really shows.

[JeffBezanson]

The issue is that GMP defines mod to be non-negative, while our mod has the same sign as the divisor.

[stevengj]

@JeffBezanson, yes I noticed. The right thing is to just change the sign in the BigInt implementation, to enforce powermod(a,b,m) == mod(a^b, m).

Another bug: the powermod implementation is not currently type-stable: e.g. returning one(b) for an exponent of 0 is wrong, since it should be one(promote_type(typeof(b),typeof(m))).

[JeffBezanson]

I'm going to merge this and then fix these powermod problems.

[stevengj]

I have a patch to fix the powermod problems, I'll push it.

[staticfloat]

Merging this has broken Travis. Here's the test failure log, looks like a test of gcdx hasn't been updated/was broken:

test failed: :((gcdx(i,j)==gcdx(ib,jb)))
 in error at error.jl:21
 in default_handler at test.jl:19
 in do_test at test.jl:39
 in anonymous at no file:1524
 in runtests at /tmp/julia/share/julia/test/testdefs.jl:5
 in anonymous at multi.jl:834
 in run_work_thunk at multi.jl:575
 in anonymous at task.jl:834
at numbers.jl:1526

[JeffBezanson]

This must be due to a GMP version difference. For gcdx(-20,-20), my GMP gives (20,0,-1) but the one on Travis gives (20,-1,0). Maybe we should extend the b==0 special case to also get called for a==b?

[JeffBezanson]

A further problem is that since our mod gives negative results when the second argument is negative, no inverses exist mod negative numbers.

[stevengj]

Inverses still exist. e.g. mod(7 * 16, -37) gives -36, and -36 is indeed a multiplicative identity element modulo -37.

[JeffBezanson]

In that case we are just up against GMP's incompatible definition of mod.

[staticfloat]

Looks like this was fixed by c92e98d. Thanks Jeff!