gcl-devel
[Top][All Lists]

## [Gcl-devel] gcd proposal

 From: Andrei Zorine Subject: [Gcl-devel] gcd proposal Date: Sat, 10 Jan 2004 05:55:54 +0300 User-agent: Mozilla/5.0 (Windows; U; Win98; en-US; rv:1.0rc1) Gecko/20020417

```hello,
```
i have a small proposal about gcd computation. I propose to use a binary gcd algorithm in the get-gcd function (num_arithm.c). It gives speed improvement over the existing algorithm. Maxima source tree contains a commac.lisp file with gcd redefined and commented out. This way it works.
```--
Andrei Zorine
```
```object
get_gcd(object x, object y)
{
int     i, j, k;
long k1,t;
object  q, r;

if (number_minusp(x))
x = number_negate(x);
if (number_minusp(y))
y = number_negate(y);

L:
if (type_of(x) == t_fixnum && type_of(y) == t_fixnum) {
i = fix(x);
j = fix(y);
/* binary gcd starts here */
k1=0;
while(i%2==0 && j%2==0)
{
k1++;
i=i>>1;
j=j>>1;
};
if(i%2==1) t=-j;
else {
t=i;
t=t>>1;
};
B3:
while(t%2==0) t=t>>1;
if(t>0) i=t; else j=-t;
t=i-j;
if(t!=0) goto B3;
return(make_fixnum(i<<k1));
} /*binary gcd ends here */

if (number_compare(x, y) < 0) {
r = x;
x = y;
y = r;
}
if (type_of(y) == t_fixnum && fix(y) == 0) {
return(x);
}
integer_quotient_remainder_1(x, y, &q, &r);
x = y;
y = r;
goto L;
}
```
```#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <sys/times.h>

long gcl_gcd(long i, long j)
{
long k;
LL:
if (i < j) {
k = i;
i = j;
j = k;
}
if (j == 0) {
return(i);
}
k = i % j;
i = j;
j = k;
goto LL;
}

long binary_gcd(long u, long v)
{
long k,t;
k=0;
while(u%2==0 && v%2==0)
{
k++;
u=u/2;
v=v/2;
};
if(u%2==1) t=-v;
else {
t=u;
t=t/2;
};
B3:
while(t%2==0) t=t/2;
if(t>0) u=t; else v=-t;
t=u-v;
if(t!=0) goto B3;
return u<<k;
};

int main(int argc, char** argv)
{
long i,j,c,cmax=45000,d,modulus=0;
clock_t start, end;
struct tms beg, fin;
double used1, used2;

if(argc>1) modulus = atoi(argv[1]);

srand( time(0) );
for(;;)
{
i = rand();
j = rand();
if(modulus>0)
{
i=i%modulus;
j=j%modulus;
};
printf("gcd(%d, %d) = 1) %d, 2) %d\t",i,j,gcl_gcd(i,j), binary_gcd(i,j));
start = times(&beg);
for(c=0;c<cmax;c++)
d=gcl_gcd(i,j);
end = times(&fin);
used1 = ((double) fin.tms_utime - (double) beg.tms_utime);
start = times(&beg);
for(c=0;c<cmax;c++)
d=binary_gcd(i,j);
end = times(&fin);
used2 = ((double) fin.tms_utime - (double) beg.tms_utime);
printf("gcl: %10.8f, binary: %10.8f, gcl-binary = %10.8f\n",
used1/(double)cmax, used2/(double)cmax, (used1-used2)/(double)cmax);
}
return 0;
}
```