• Matthias Walter (12/12) Feb 07 2011 Hi everyone,
• bearophile (6/12) Feb 07 2011 See also:
• Andrei Alexandrescu (5/17) Feb 07 2011 Please do. I implemented gcd in a hurry because I needed it for unsigned...
• Matthias Walter (38/59) Feb 13 2011 I re-implemented [1] several things like recursive / iterative or modulo
• bearophile (6/9) Feb 14 2011 Are benchmarks done with "BigInt" and "long" too? (If you test bigints y...
• Matthias Walter (48/53) Feb 14 2011 Yeah, so I did tests for long now, too. Unfortunately with gdc I was
• Andrei Alexandrescu (3/56) Feb 15 2011 Sounds great. Thanks!
• Matthias Walter (13/16) Feb 14 2011 I'd like to work on the BigInt things but there are a couple of blockers...
• Don (2/22) Feb 19 2011 There's a (very complicated) O(n log n) algorithm for GCD for large numb...

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

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

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

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!

I re-implemented [1] several things like recursive / iterative or modulo
/ 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

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

Matthias Walter <xammy xammy.homelinux.net> writes:

 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).

Yeah, so I did tests for long now, too. Unfortunately with gdc I was
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?

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

I agree on that and suggest to take Boost's implementation which is very
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

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?

 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).

 Yeah, so I did tests for long now, too. Unfortunately with gdc I was
 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?

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

 I agree on that and suggest to take Boost's implementation which is very
 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

Matthias Walter <xammy xammy.homelinux.net> writes:

 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).

I'd like to work on the BigInt things but there are a couple of blockers:

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

Don <nospam nospam.com> writes:

Matthias Walter wrote:
 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).

 I'd like to work on the BigInt things but there are a couple of blockers:
 
 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.