www.digitalmars.com         C & C++   DMDScript  

digitalmars.D - greatest common divisor implementation

reply Matthias Walter <xammy xammy.homelinux.net> writes:
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
next sibling parent bearophile <bearophileHUGS lycos.com> writes:
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
prev sibling next sibling parent reply Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
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
next sibling parent reply bearophile <bearophileHUGS lycos.com> writes:
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
next sibling parent Andrei Alexandrescu <SeeWebsiteForEmail erdani.org> writes:
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.081
 2. 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
prev sibling parent Don <nospam nospam.com> writes:
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
prev sibling next sibling parent Matthias Walter <xammy xammy.homelinux.net> writes:
 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.081
 2. 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
prev sibling parent Matthias Walter <xammy xammy.homelinux.net> writes:
 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
prev sibling parent Matthias Walter <xammy xammy.homelinux.net> writes:
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