Re: [Gcl-devel] gcd proposal

[Camm Maguire]

Greetings, and thanks for your suggestion!  Tested and installed.
About a 10% speedup on a pIII, ~16% on a P4/xeon.

Take care,

Andrei Zorine <address@hidden> writes:

> 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;
> }
> _______________________________________________
> Gcl-devel mailing list
> address@hidden
> http://mail.gnu.org/mailman/listinfo/gcl-devel

-- 
Camm Maguire                                            address@hidden
==========================================================================
"The earth is but one country, and mankind its citizens."  --  Baha'u'llah