05.30.06

Computing gcd’s was never this hard!

Posted in Mathematics, ruby/rails at 1:11 pm by Haris

This post has gotten a bit longer than expected, so a quick summary is in order here: If you run Ruby 1.8.2 or earlier, be very careful about using rational.rb and mathn.rb. The algorithms in there for computing the gcd of numbers, necessary in almost every operation, are from kind of bad to utterly horrible, and in particular in one of my projects 30% of the program time was spent in computing gcds.

Ok, and now for our feature presentation:

We all know how to compute the greatest common divisor (gcd) of two numbers, right? It is very easy and fast to do. For instance in Ruby it would be done via something like this:

def mygcd2(one,other)
min = one.abs
max = other.abs
while min > 0
min, max = max % min, min
end
max
end

For those not familiar with the idioms of Ruby, every method in Ruby returns the last expression evaluated, if there is no explicit return, hence the max at the last line. Also, you can do assignments in pairs as in min, max = max % min, min. This is almost equivalent to:

def mygcd(one,other)
min = one.abs
max = other.abs
while min > 0
tmp = min
min = max % min
max = tmp
end
max
end

There is a subtle difference, in that for the double assignment that takes place in mygcd2 a temporary array is created each time, which as we will see in a bit offers some slowdown.

Ok, so, what is the big deal, we all know this, right?

I was profiling a Series class I have constructed to do manipulations with power series, and found the following:

%   cumulative   self              self     total
time   seconds   seconds    calls  ms/call  ms/call  name
28.39   410.17    410.17    39850    10.29    15.95  Integer#gcd

For those of you not familiar with the output of the Ruby profiler, my application was spending 30% of its time running the method that computes the gcd of two integers!!!! And each call to the method was taking on average 10ms!!!!

WHAT!!! HOW CAN THIS BE????

I mean it has to do a lot of operations with rational numbers, so gcd is probably called a number of times (in fact every time a new rational number is defined, hence the 40000 calls), but still….

Now before I continue I should point out that I am using Ruby 1.8.2, and I understand that these things are much more improved in future versions. So I will criticize this particular version of Ruby, along with version 0.5 of mathn, which is what ships with MacOSX 10.4.

Ok, so how does Ruby compute the gcd of two numbers? Well, the beauty of it is that you can just look at the source, which in my case is located at file://localhost/usr/lib/ruby/1.8/rational.rb What you see there is the following:

def gcd(n)
m = self.abs
n = n.abs

return n if m == 0
return m if n == 0

b = 0
while n[0] == 0 && m[0] == 0
b += 1; n >>= 1; m >>= 1
end
m >>= 1 while m[0] == 0
n >>= 1 while n[0] == 0
while m != n
m, n = n, m if n > m
m -= n; m >>= 1 while m[0] == 0
end
m << b
end

def gcd2(int)
a = self.abs
b = int.abs

a, b = b, a if a < b

while b != 0
void, a = a.divmod(b)
a, b = b, a
end
return a
end

Hm, this looks a bit complicated, but maybe it performs well, we’ll see. If I understand what it does, it effectively first eliminates any power of 2 that the two numbers might have in common, by doing bit-shifting, which one might think is fast. But the real question, which I still have and do not understand is: WHAT IS gcd2? Where is it used? It looks close to my algorithms above, though still more inefficient. But what is it for? Why do we have two gcd’s???

Please, if anyone knows the answer to that, tell me because I desperately need to know.

Anyway, let’s see how these perform. This was done using the Benchmark module in a standard way, and you can find the code here. Here are the results:

Computing gcd of 100x100 pairs starting at 0
user     system      total        real
gcd     0.580000   0.010000   0.590000 (  0.654295)
gcd2    0.310000   0.000000   0.310000 (  0.366262)
mygcd   0.180000   0.010000   0.190000 (  0.225134)
mygcd2  0.240000   0.000000   0.240000 (  0.272406)

Computing gcd of 1000x1000 pairs starting at 0
user     system      total        real
gcd    86.790000   0.580000  87.370000 ( 98.705976)
gcd2   43.550000   0.440000  43.990000 ( 53.308317)
mygcd  24.120000   0.150000  24.270000 ( 26.866034)
mygcd2 31.770000   0.310000  32.080000 ( 38.794325)

Computing gcd of all pairs up to 100
user     system      total        real
gcd      0.560000   0.010000   0.570000 (  0.639058)
gcd2     0.280000   0.000000   0.280000 (  0.348463)
mygcd    0.170000   0.010000   0.180000 (  0.209354)
mygcd2   0.220000   0.000000   0.220000 (  0.243531)

Hm, both the gcd and gcd2 methods perform pretty slowly compared to my more naive versions. Could they perhaps be better for large numbers? Let’s see:

Computing gcd of 100x100 pairs starting at 100000
user     system      total        real
gcd      1.460000   0.020000   1.480000 (  1.649013)
gcd2     0.390000   0.000000   0.390000 (  0.437032)
mygcd    0.220000   0.000000   0.220000 (  0.245309)
mygcd2   0.290000   0.010000   0.300000 (  0.345750)

Computing gcd of 100x100 pairs starting at 100000000
user     system      total        real
gcd      2.210000   0.020000   2.230000 (  2.456261)
gcd2     0.390000   0.000000   0.390000 (  0.434390)
mygcd    0.230000   0.000000   0.230000 (  0.243268)
mygcd2   0.290000   0.010000   0.300000 (  0.330893)

Computing gcd of 100x100 pairs starting at 100000000000
user     system      total        real
gcd     3.410000   0.050000   3.460000 (  3.941590)
gcd2    0.580000   0.010000   0.590000 (  0.666822)
mygcd   0.420000   0.000000   0.420000 (  0.491752)
mygcd2  0.500000   0.010000   0.510000 (  0.560805)

Nope, in fact if anything they perform worse, especially gcd.

Ok, so I guess at this point I have three questions:

1. Why not implement the standard method, which is the first program one learns to write in Computer Science, in fact it is exactly how one learns to compute the gcd at school?
2. Why not benchmark it against the standard algorithm?
3. Is there something I am missing, that makes the implemented gcd method preferable?
4. What on earth does gcd2 do?

Hm, but I mentioned mathn.rb at some point and haven’t talked about it yet. First of all, what mathn does is to make the various number classes work better together. For instance, 1/2 actually becomes the rational number 1/2, instead of 0. So as you can imagine it is absolutely critical for what I wanted to do. Let’s look at some of the code implemented there, notably a reimplementation of gcd2, which is most certainly a criminal offense:

def gcd2(int)
a = self.abs
b = int.abs
a, b = b, a if a < b

pd_a = a.prime_division
pd_b = b.prime_division

gcd = 1
for pair in pd_a
as = pd_b.assoc(pair[0])
if as
gcd *= as[0] ** [as[1], pair[1]].min
end
end
return gcd
end

Yes, you guessed it right, this method factors a and b into their prime components, and then assembles together the common factors.

Why is this a criminal offense you ask? Well…

1. The algorithm in mygcd is O(logn), i.e. its running time is linear in the number of digits of the two numbers. It just divides each number into the other and so on, which is a relatively quick operation.
2. This new gcd2 effectively computes all divisors of the two numbers, and then from there obtains the greatest common divisor. To get an idea of why this is bad, note that until 3 years ago we did not even know how to determine in polynomial time whether a number is prime, and even now the running time follows the 6th power of logn at best, far from linear. What is more, it is believed that one cannot factor numbers in a reasonable amount of time. In fact, a number of encryption algorithms rely on the fact that we cannot effectively factor the product of two reasonably large prime numbers, let alone a general number with many more factors! This problem is not even solvable in polynomial time, let alone not being linear!!!
3. Did I mention that a new Prime number generator is used each time this gcd2 runs? Let me emphasize this: THE LIST OF PRIMES IS GENERATED ANEW EACH TIME THE GCD OF TWO NUMBERS IS TO BE COMPUTED. IT IS NOT CACHED

I don’t even need to show you any benchmarks here, they are absolutely laughable the moment the numbers have more than 4 digits. I guess my main questions here are:

1. How is it possible that the person implementing the code that makes the various math classes in Ruby work nicely together implements a gcd algorithm (used each time a rational number is created) that is not just slow, but not even believed to be running in a reasonable amount of time?
2. Why oh why did they feel the need to reimplement gcd2, which is already implemented much better in rational.rb, which is included?
3. Didn’t they test their algorithm for relatively large numbers?
4. How is it possible that this was not spotted earlier?

Luckily, as I understand it, this has been fixed in more recent versions on Ruby, so that’s great. Ok, enough ranting on my part.

Later

Permalink