## digitalmars.D - greatest common divisor implementation

- Matthias Walter <xammy xammy.homelinux.net> Feb 07 2011
- bearophile <bearophileHUGS lycos.com> Feb 07 2011
- Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> Feb 07 2011
- bearophile <bearophileHUGS lycos.com> Feb 14 2011
- Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> Feb 15 2011
- Don <nospam nospam.com> Feb 19 2011
- Matthias Walter <xammy xammy.homelinux.net> Feb 14 2011
- Matthias Walter <xammy xammy.homelinux.net> Feb 14 2011
- Matthias Walter <xammy xammy.homelinux.net> Feb 13 2011

Hi everyone, as I'm currently working on a C++ project which involves gcd computations I had a quick look at phobos' implementation. 1. First thing I saw is that gcd(-3,6) raises an exception, although mathematically it is just 2. Of course checking the sign of the arguments takes computation time, but for my stuff I'd definitely need it. If it's really too expensive then there should be at least another routine. 2. I'd be happy to implement the Extended Euclidean algorithm for phobos which also gives multipliers s and t such that gcd(x,y) = s*x + t*y. Any comments before I start working on it? Matthias

Feb 07 2011

Matthias Walter:1. First thing I saw is that gcd(-3,6) raises an exception, although mathematically it is just 2. Of course checking the sign of the arguments takes computation time, but for my stuff I'd definitely need it. If it's really too expensive then there should be at least another routine.

I agree. I my code I have not used Phobos gcd() because of that.Any comments before I start working on it?

See also: http://d.puremagic.com/issues/show_bug.cgi?id=4125 Bye, bearophile

Feb 07 2011

On 2/7/11 3:18 PM, Matthias Walter wrote:Hi everyone, as I'm currently working on a C++ project which involves gcd computations I had a quick look at phobos' implementation. 1. First thing I saw is that gcd(-3,6) raises an exception, although mathematically it is just 2. Of course checking the sign of the arguments takes computation time, but for my stuff I'd definitely need it. If it's really too expensive then there should be at least another routine. 2. I'd be happy to implement the Extended Euclidean algorithm for phobos which also gives multipliers s and t such that gcd(x,y) = s*x + t*y. Any comments before I start working on it? Matthias

Please do. I implemented gcd in a hurry because I needed it for unsigned numbers, so it's not the greatest it could be. A github pull request would be much appreciated! Andrei

Feb 07 2011

Matthias Walter: You have done lot of nice work :-)1. Are there any further suggestions on the implementations / Did I forget something?

Are benchmarks done with "BigInt" and "long" too? (If you test bigints you need bigger numbers too, and to test that the results are correct).2. What do you think, we should finally use?

It seems more complex algorithms don't pay much for int values. Bye, bearophile

Feb 14 2011

On 2/14/11 10:38 PM, Matthias Walter wrote:1. Are there any further suggestions on the implementations / Did I forget something?

unable to generate uniform long numbers ("UniformIntGenerator.popFront not implemented for large ranges"). Interestingly, with DMD it works... The dmd-only results of long as type, with random long values: SimpleRecursive!(long): 0.134 SimpleIterative!(long): 0.127 BinaryRecursive!(long): 1.727 BinaryIterative!(long): 0.345 BoostIterative!(long): 0.136 BoostBinary!(long): 0.575 gcdSteinTable!(long, 1): 0.594 gcdSteinTable!(long, 2): 0.527 gcdSteinTable!(long, 3): 0.489 gcdSteinTable!(long, 4): 0.536 gcdSteinTable!(long, 5): 0.489 gcdSteinTable!(long, 6): 0.531 gcdSteinTable!(long, 7): 0.498 gcdSteinTable!(long, 8): 0.542 gcdSteinTable!(long, 9): 0.491 gcdSteinTable!(long,10): 0.539 gcdSteinTable!(long,11): 0.494 gcdSteinTable!(long,12): 0.538 The results for dmd and gdc with long as type, but random int values: SimpleRecursive!(long): 0.110 0.143 SimpleIterative!(long): 0.100 0.136 BinaryRecursive!(long): 0.661 0.208 BinaryIterative!(long): 0.166 0.133 BoostIterative!(long): 0.104 0.142 BoostBinary!(long): 0.277 0.174 gcdSteinTable!(long, 1): 0.275 0.134 gcdSteinTable!(long, 2): 0.259 0.130 gcdSteinTable!(long, 3): 0.249 0.108 gcdSteinTable!(long, 4): 0.253 0.095 gcdSteinTable!(long, 5): 0.243 0.089 gcdSteinTable!(long, 6): 0.243 0.086 gcdSteinTable!(long, 7): 0.245 0.084 gcdSteinTable!(long, 8): 0.245 0.082 gcdSteinTable!(long, 9): 0.246 0.081 gcdSteinTable!(long,10): 0.245 0.081 gcdSteinTable!(long,11): 0.246 0.080 gcdSteinTable!(long,12): 0.250 0.0812. What do you think, we should finally use?

readable, too. I will probably do some tests with BigInt, too. But I suspect that the result looks totally different which means that we'll use different specializations. The same holds true for the extended algorithm. Andrei: As soon as I'm done, I'll submit a pull request to github. Matthias

Sounds great. Thanks! Andrei

Feb 15 2011

Matthias Walter wrote:1. Are there any further suggestions on the implementations / Did I forget something?

1. abs() does not work, because opCmp for BigInts is not pure, yet. IIRC Don waits / waited for a pure/nothrow bug to be fixed before he'd continue working on BigInt purity, etc. 2. I don't think the current random number implementation has code to generate BigInts easily via uniform() method. There was a fixed number of bits that are used per call. So I suggest to work on the Extended Euclidean algorithm implementations and get results for int/long for those. If this is done I'll submit a pull request with nice code for int+long and code for BigInt that just works but is not benchmarked. Then we should change the gcd bug to BigInt-only (or open another one...) Matthias

There's a (very complicated) O(n log n) algorithm for GCD for large numbers.

Feb 19 2011

1. Are there any further suggestions on the implementations / Did I forget something?

unable to generate uniform long numbers ("UniformIntGenerator.popFront not implemented for large ranges"). Interestingly, with DMD it works... The dmd-only results of long as type, with random long values: SimpleRecursive!(long): 0.134 SimpleIterative!(long): 0.127 BinaryRecursive!(long): 1.727 BinaryIterative!(long): 0.345 BoostIterative!(long): 0.136 BoostBinary!(long): 0.575 gcdSteinTable!(long, 1): 0.594 gcdSteinTable!(long, 2): 0.527 gcdSteinTable!(long, 3): 0.489 gcdSteinTable!(long, 4): 0.536 gcdSteinTable!(long, 5): 0.489 gcdSteinTable!(long, 6): 0.531 gcdSteinTable!(long, 7): 0.498 gcdSteinTable!(long, 8): 0.542 gcdSteinTable!(long, 9): 0.491 gcdSteinTable!(long,10): 0.539 gcdSteinTable!(long,11): 0.494 gcdSteinTable!(long,12): 0.538 The results for dmd and gdc with long as type, but random int values: SimpleRecursive!(long): 0.110 0.143 SimpleIterative!(long): 0.100 0.136 BinaryRecursive!(long): 0.661 0.208 BinaryIterative!(long): 0.166 0.133 BoostIterative!(long): 0.104 0.142 BoostBinary!(long): 0.277 0.174 gcdSteinTable!(long, 1): 0.275 0.134 gcdSteinTable!(long, 2): 0.259 0.130 gcdSteinTable!(long, 3): 0.249 0.108 gcdSteinTable!(long, 4): 0.253 0.095 gcdSteinTable!(long, 5): 0.243 0.089 gcdSteinTable!(long, 6): 0.243 0.086 gcdSteinTable!(long, 7): 0.245 0.084 gcdSteinTable!(long, 8): 0.245 0.082 gcdSteinTable!(long, 9): 0.246 0.081 gcdSteinTable!(long,10): 0.245 0.081 gcdSteinTable!(long,11): 0.246 0.080 gcdSteinTable!(long,12): 0.250 0.0812. What do you think, we should finally use?

readable, too. I will probably do some tests with BigInt, too. But I suspect that the result looks totally different which means that we'll use different specializations. The same holds true for the extended algorithm. Andrei: As soon as I'm done, I'll submit a pull request to github. Matthias

Feb 14 2011

1. Are there any further suggestions on the implementations / Did I forget something?

1. abs() does not work, because opCmp for BigInts is not pure, yet. IIRC Don waits / waited for a pure/nothrow bug to be fixed before he'd continue working on BigInt purity, etc. 2. I don't think the current random number implementation has code to generate BigInts easily via uniform() method. There was a fixed number of bits that are used per call. So I suggest to work on the Extended Euclidean algorithm implementations and get results for int/long for those. If this is done I'll submit a pull request with nice code for int+long and code for BigInt that just works but is not benchmarked. Then we should change the gcd bug to BigInt-only (or open another one...) Matthias

Feb 14 2011

On 02/07/2011 05:13 PM, Andrei Alexandrescu wrote:On 2/7/11 3:18 PM, Matthias Walter wrote:Hi everyone, as I'm currently working on a C++ project which involves gcd computations I had a quick look at phobos' implementation. 1. First thing I saw is that gcd(-3,6) raises an exception, although mathematically it is just 2. Of course checking the sign of the arguments takes computation time, but for my stuff I'd definitely need it. If it's really too expensive then there should be at least another routine. 2. I'd be happy to implement the Extended Euclidean algorithm for phobos which also gives multipliers s and t such that gcd(x,y) = s*x + t*y. Any comments before I start working on it? Matthias

Please do. I implemented gcd in a hurry because I needed it for unsigned numbers, so it's not the greatest it could be. A github pull request would be much appreciated!

/ binary implementations. For the latter I found one with a lookup-table that pre-computes partial results. The testbed generated a sequence of 10000 int pairs, computing the gcd of a pair of ints and adding the sum of all (as checksum and to avoid removal by the compiler). This was repeated 100 times using the same sequence. Time is per gcd call and given in nanoseconds. First column: dmd -O -release -inline Second column: gdc -O3 SimpleRecursive: 0.096 0.086 SimpleIterative: 0.090 0.083 BinaryRecursive: 0.240 0.213 BinaryIterative: 0.141 0.133 BoostIterative: 0.087 0.085 BoostBinary: 0.219 0.177 gcdSteinTable!(int, 1): 0.161 0.136 gcdSteinTable!(int, 2): 0.133 0.130 gcdSteinTable!(int, 3): 0.107 0.108 gcdSteinTable!(int, 4): 0.093 0.094 gcdSteinTable!(int, 5): 0.089 0.088 gcdSteinTable!(int, 6): 0.087 0.083 gcdSteinTable!(int, 7): 0.086 0.081 gcdSteinTable!(int, 8): 0.085 0.081 gcdSteinTable!(int, 9): 0.085 0.082 gcdSteinTable!(int,10): 0.085 0.080 gcdSteinTable!(int,11): 0.086 0.082 gcdSteinTable!(int,12): 0.086 0.080 Remark: The 2nd template parameter of the gcdSteinTable implementation means the number of bits for the lookup-table. It therefore uses 2^bits bytes of memory! 1. Are there any further suggestions on the implementations / Did I forget something? 2. What do you think, we should finally use? 3. I will also implement the extended gcd (i.e. having gcd(a, b) = z = x * a + y * b we are interested in x,y,z) Matthias [1] http://pastebin.com/sXpsdK5S

Feb 13 2011