"Ben Pfaff" <bl*@cs.stanfor d.eduwrote in message
news:87******** ****@benpfaff.o rg...
"Dann Corbit" <dc*****@connx. comwrites:
>/*
Now that I actually understand how it works, I like the subtraction trick
a
lot better. In this instance, (lots of small repeated factors) it takes
the
same number of iterations as the modulo version
*/
Are we talking about the same binary GCD algorithm that's in
Knuth? There's a wikipedia page about it:
http://en.wikipedia.org/wiki/Binary_GCD_algorithm
(Every algorithm is in Knuth if you look hard enough.)
Tom's version seems somewhat faster than the web article version.
The timings below are after profile guided optimization.
/*
Now that I actually understand how it works, I like the subtraction trick a
lot better. In this instance, (lots of small repeated factors) it takes the
same number of iterations as the modulo version
*/
unsigned long long modulos = 0;
unsigned long long subtractions = 0;
unsigned long long gcdm(unsigned long long a, unsigned long long b)
{
unsigned long long k, t;
k = 0;
while (!(a & 1 || b & 1)) {
++k; a >>= 1; b >>= 1;
}
while (!(a & 1)) { a >>= 1; }
while (!(b & 1)) { b >>= 1; }
while (b) {
if (a b) { t = a; a = b; b = t; }
b = b % a;
while (b && !(b & 1)) { b >>= 1; }
}
return a << k;
}
unsigned long long gcds(unsigned long long a, unsigned long long b)
{
unsigned long long k, t;
k = 0;
while (!(a & 1 || b & 1)) {
++k; a >>= 1; b >>= 1;
}
while (!(a & 1)) { a >>= 1; }
while (!(b & 1)) { b >>= 1; }
while (b) {
if (a b) { t = a; a = b; b = t; }
b = b - a;
while (b && !(b & 1)) { b >>= 1; }
}
return a << k;
}
unsigned long long gcd(unsigned long long u, unsigned long long v) {
unsigned long long k = 0;
if (u == 0)
return v;
if (v == 0)
return u;
while ((u & 1) == 0 && (v & 1) == 0) { /* while both u and v are even
*/
u >>= 1; /* shift u right, dividing it by 2 */
v >>= 1; /* shift v right, dividing it by 2 */
k++; /* add a power of 2 to the final result */
}
/* At this point either u or v (or both) is odd */
do {
if ((u & 1) == 0) /* if u is even */
u >>= 1; /* divide u by 2 */
else if ((v & 1) == 0) /* else if v is even */
v >>= 1; /* divide v by 2 */
else if (u >= v) /* u and v are both odd */
u = (u-v) >1;
else /* u and v both odd, v u */
v = (v-u) >1;
} while (u 0);
return v << k; /* returns v * 2^k */
}
unsigned long long big=12974633789 0625;
unsigned long long med=86497558593 75*7;
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
unsigned long long randvals[1000000];
int main(void)
{
clock_t start;
clock_t end;
unsigned long long answer=0;
size_t index;
for (index = 0; index < 1000000; index++)
{
randvals[index] = rand();
randvals[index] <<= 32;
randvals[index] += rand();
}
start = clock();
for (index = 0; index < 1000000-1; index++)
answer += gcdm(randvals[index],randvals[index+1]);
end = clock();
printf("clocks to do one million modular gcd's = %u, summed GCD =
%llu\n", (unsigned) end-start, answer) ;
answer = 0;
start = clock();
for (index = 0; index < 1000000-1; index++)
answer += gcds(randvals[index],randvals[index+1]);
end = clock();
printf("clocks to do one million subtraction gcd's = %u, summed GCD =
%llu\n", (unsigned) end-start, answer) ;
answer = 0;
start = clock();
for (index = 0; index < 1000000-1; index++)
answer += gcd(randvals[index],randvals[index+1]);
end = clock();
printf("clocks to do one million subtraction gcd's {web article version} =
%u, summed GCD = %llu\n", (unsigned) end-start, answer) ;
return 0;
}
/*
clocks to do one million modular gcd's = 1031, summed GCD = 9203554
clocks to do one million subtraction gcd's = 766, summed GCD = 9203554
clocks to do one million subtraction gcd's {web article version} = 1109,
summed GCD = 9203554
Press any key to continue . . .
*/
--
int main(void){char
p[]="ABCDEFGHIJKLM NOPQRSTUVWXYZab cdefghijklmnopq rstuvwxyz.\
\n",*q="kl BIcNBFr.NKEzjwC IxNJC";int i=sizeof p/2;char *strchr();int
putchar(\
);while(*q){i+= strchr(p,*q++)-p;if(i>=(int)si zeof p)i-=sizeof
p-1;putchar(p[i]\
);}return 0;}
