gcd improved & optimized + questions

[Jaroslav Hajek]

hi all,

with the changeset
http://hg.savannah.gnu.org/hgweb/octave/rev/0ee74d581c00

I have rewritten gcd to become more efficient and robust.

1. It now works reliably for large 64 bit integers:

With a recent tip:
address@hidden:~/devel/octave$ octave -q
octave:1> gcd (int64(2^54) + 5, int64(2^54) + 10427872)
ans =  68
---> completely wrong

with the new patch:
address@hidden:~/devel/octave$ ./run-octave -q
octave:1> gcd (int64(2^54) + 5, int64(2^54) + 10427872)
ans = 10427867
---> correct answer (as given by PARI)

2. performance is significantly better. here's a benchmark:

n = 1e6;
a = ceil (5e4*rand (n, 1));
b = ceil (2e4*rand (n, 1));

tic; gcd (a, b); toc
tic; [g, c, d] = gcd (a, b); toc

a = int16 (a);
b = int16 (b);

tic; gcd (a, b); toc
tic; [g, c, d] = gcd (a, b); toc

with a recent tip (Core 2 Duo @ 2.83GHz, g++ -O3 -march=native), I get

Elapsed time is 5.18626 seconds.
Elapsed time is 4.74562 seconds.
Elapsed time is 4.85944 seconds.
Elapsed time is 4.71784 seconds.

with the new patch, I get:

Elapsed time is 0.205938 seconds.
Elapsed time is 0.452337 seconds.
Elapsed time is 0.132734 seconds.
Elapsed time is 0.137132 seconds.

I retained the useful ability to specify more than two args to gcd,
even though Matlab apparently doesn't allow it. However, the old
implementation also featured a special case when a single input was
given - gcd of all elements was computed. This was marked as being for
backward compatibility. I didn't implement it because I thought that
now may be a good time to drop it. I think it was confusing because if
gcd is a mapper function, then the natural definition is gcd(X) = X.
Chances are good that Matlab will do the same if they ever extend gcd
to allow variable argument list. Octave's gcd was behaving a bit like
a reduction function in this special case (although it couldn't reduce
along a dimension).
So, to avoid confusion (and save myself some work), I simply forbade
the case.  It didn't seem to be used anywhere else in Octave or
OctaveForge anyway. Opinions? Is this OK? Should it be mentioned in
NEWS?

enjoy

-- 
RNDr. Jaroslav Hajek, PhD
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz