gcl-devel
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[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;
}

reply via email to

[Prev in Thread] Current Thread [Next in Thread]