Content deleted Content added
→Implementations: added pseudocode |
No edit summary Tag: repeating characters |
||
Line 145:
I didn't understand why the case where u and v are powers of two needed fixing, could you please elaborate? -- [[User:Regregex|Regregex]] ([[User talk:Regregex|talk]]) 17:02, 10 October 2008 (UTC)
:I don't remember, probably I just screwed up. :-) [[User:Balabiot|Balabiot]] ([[User talk:Balabiot|talk]]) 17:03, 23 October 2008 (UTC)
I ran 10 000 000 random tests on Mac OS X.6.4 with C Code and got the following results:
*Euclide : 6180 ms
*Stein : 9038 ms
<source lang="text">
#include <stdio.h>
#include <sys/time.h>
#include <stdlib.h>
#define N 10000000
#define MAX_UINTplus1 0x100000000LLU
#define uint64 unsigned long long
#define INIT_RANDOM srandom(3)
uint64 pgcdEuclide(uint64,uint64);
uint64 pgcdStein(uint64,uint64);
int main(void) {
int i;
uint64 a,b,d,t0,t1;
struct timeval tv;
// Version Euclide
INIT_RANDOM;
gettimeofday(&tv,0);
t0 = tv.tv_sec* 1000 *1000+tv.tv_usec;
for(i=0;i<N;i++){
a=random()*MAX_UINTplus1+random();
b=random()*MAX_UINTplus1+random();
d=pgcdEuclide(a,b);
}
gettimeofday(&tv,0);
t1 = tv.tv_sec*1000*1000+tv.tv_usec;
printf("Euclide : %llu ms\n",(t1-t0)/1000);
// Version Stein
INIT_RANDOM;
gettimeofday(&tv,0);
t0 = tv.tv_sec* 1000 *1000+tv.tv_usec;
for(i=0;i<N;i++){
a=(uint64)random()*MAX_UINTplus1+random();
b=(uint64)random()*MAX_UINTplus1+random();
d=pgcdStein(a,b);
}
gettimeofday(&tv,0);
t1 = tv.tv_sec*1000*1000+tv.tv_usec;
printf("Stein : %llu ms\n",(t1-t0)/1000);
return 0;
}
uint64 pgcdEuclide(uint64 a,uint64 b){
if(b==0) return a;
return pgcdEuclide(b,a%b);
}
uint64 pgcdStein(uint64 u,uint64 v){
int shift;
/* GCD(0,x) := x */
if (u == 0 || v == 0)
return u | v;
/* Let shift := lg K, where K is the greatest power of 2
dividing both u and v. */
for (shift = 0; ((u | v) & 1) == 0; ++shift) {
u >>= 1;
v >>= 1;
}
while ((u & 1) == 0)
u >>= 1;
/* From here on, u is always odd. */
do {
while ((v & 1) == 0) /* Loop X */
v >>= 1;
/* Now u and v are both odd, so diff(u, v) is even.
Let u = min(u, v), v = diff(u, v)/2. */
if (u < v) {
v -= u;
} else {
uint64 diff = u - v;
u = v;
v = diff;
}
v >>= 1;
} while (v != 0);
return u << shift;
}
</source>
| |||