/*---------------------------------------------------------------------- File : gcd.c Contents: compute greatest common divisor with school algorithm (i.e. using a common prime factors approach) Author : Christian Borgelt History : 09.12.1997 file created 18.12.1997 lecture algorithm and its improvement added 05.01.1998 second improvement of lecture algorithm added ----------------------------------------------------------------------*/ #include #include #ifndef VERSION #define VERSION 5 /* default: compile fifth version */ #endif /*---------------------------------------------------------------------- Functions ----------------------------------------------------------------------*/ #if (VERSION == 1) /* if to compile first version */ int power (int p, int q) { /* --- compute p to the power of q */ int t; /* temporary buffer */ if (q <= 1) /* p^0 and p^1 can be */ return (q <= 0) ? 1 : p; /* computed directly */ t = power(p, q /2); /* compute p^|_q/2_| */ if (q %2 == 0) return t*t; /* even q: p^q = (p^(q/2))^2 */ else return t*t*p; /* odd q: p^q = (p^((q-1)/2))^2 *p */ } /* power() */ /*---------------------------------------------------------------------- Note that the function above is already optimized (a very similar function to compute p^q has been used in the program `euler.c'). A naive approach to compute p^q is: int power (int p, int q) { int t = 1; while (--q >= 0) t *= p; return t; } ----------------------------------------------------------------------*/ int gcd (int m, int n) { /* --- lecture algorithm */ int i; /* to traverse test divisors */ int si, ti; /* exponents of i in m and n */ int ui; /* minimum of the exponents ti and si */ int d = 1; /* greatest common divisor */ /* 1: i = 2, d = 1. (see also steps 5 and 6) */ for (i = 2; i <= n; i++) { /* traverse test divisors */ /* 2: find maximal si, such that i^si divides m. */ for (si = 0; m % power(i, si+1) == 0; si++); /* 3: find maximal ti, such that i^ti divides n. */ for (ti = 0; n % power(i, ti+1) == 0; ti++); /* 4: ui = min(si,ti), d = d*i^ui. */ ui = (si < ti) ? si : ti; /* set ui = min(si,ti) */ d *= power(i, ui); /* multiply common factor i^ui */ /* into the gcd (stored in d) */ /* 5: m = m/i^si, n = n/i^ti, i = i+1. */ m /= power(i, si); /* divide out the factors i^si */ n /= power(i, ti); /* and i^ti from m and n, resp. */ } /* 6: if i <= n, then goto 2, otherwise return d. */ return d; /* return the computed gcd */ } /* gcd() */ #endif /*---------------------------------------------------------------------- The above functions implement directly the algorithm discussed in the lecture (the steps of this algorithm are given in the comments; note that step 1, the `i = i+1' of step 5, and the conditional `goto' in step 6 have all be combined into a for-statement). This makes the function `gcd' easy to understand, but not very efficient, since in this function the expensive function `power' is called several times. ----------------------------------------------------------------------*/ #if (VERSION == 2) /* if to compile second version */ int gcd (int a, int b) { /* --- improved lecture algorithm 1 */ int d; /* test divisor */ int e; /* exponent of test divisor */ int r = 1; /* result (greatest common divisor) */ for (d = 2; d <= a; d++) { /* traverse test divisors */ for (e = 0; a %d == 0; e++) /* determine the exponent of d in a */ a /= d; /* and divide out the factor d from a */ while (b %d == 0) { /* while b is (still) divisible by d, */ b /= d; /* divide out the factor d from b */ if (--e >= 0) r *= d; /* if the exponent of d in m is not */ } /* already used up, multiply d into r */ } /* (i.e. multiply r by d^min(ea,eb)) */ return r; /* return the greatest common divisor */ } /* gcd() */ #endif /*---------------------------------------------------------------------- This is simply the lecture algorithm without any calls to the function `power'. This optimization can be achieved by noting that the exponent of the test divisor d in a (or b) can be determined by counting how often a (or b) can be divided by d without a remainder. This is done for a in the inner for-loop, where e counts the number of times a can be divided (i.e. after the loop e contains the exponent of d in the original a). The same is done for b in the while loop. The only difference is that we do not count the number of times we can divide b by d, but simply divide out the factor. Each time we do so, we multiply d into the result r, if the remaining exponent of d in a still allows us to. Note that by dividing out the factors, the outer for-loop is executed less often, provided a is not prime. Note also that in the while loop, it is possible to shift the condition of the `if' into the condition of the `while', i.e. it is possible to write: while ((--e >= 0) && (b %d == 0)) { b /= d; r *= d; } In this case not all factors d are divided out from b, but this does not invalidate the result (see also version 4, especially the exercise). ----------------------------------------------------------------------*/ #if (VERSION == 3) /* if to compile third version */ int gcd (int a, int b) { /* --- improved lecture algorithm 2 */ int d; /* test divisor */ int fa, fb; /* factors of a and b */ int r = 1; /* result (greatest common divisor) */ for (d = 2; d <= a; d++) { /* traverse test divisors */ for (fa = 1; a %d == 0; fa *= d) /* compute largest factor of a */ a /= d; /* that is a power of d */ for (fb = 1; b %d == 0; fb *= d) /* compute largest factor of b */ b /= d; /* that is a power of d */ r *= (fa < fb) ? fa : fb; /* multiply smaller factor into */ } /* the greatest common divisor */ return r; /* return computed gcd */ } /* gcd() */ #endif /*---------------------------------------------------------------------- This is also the lecture algorithm without any calls to the function `power'. It computes the factors of a and b that are powers of a given test divisor d in the variables fa and fb, i.e. it determines directly the values of fa = d^ea and fb = d^eb, where ea and eb are the maximal values, such that d^ea divides a and d^eb divides b. In the same loop these factors are divided out from the numbers a and b. In this way the algorithm combines step 2, step 3, the second part of step 4, and step 5 into two for-loops. The smaller of the factors fa and fb is then multiplied into the gcd (first part of step 4). ----------------------------------------------------------------------*/ #if (VERSION == 4) /* if to compile fourth version */ int gcd (int a, int b) { /* --- O(min(a,b)) */ int d; /* test divisor, temporary buffer */ int r = 1; /* result (greatest common divisor) */ if (a > b) { /* if a is larger than b, exchange */ d = a; a = b; b = d; } /* (to make a the smaller of the two) */ for (d = 2; d <= a; d++) { /* traverse test divisors */ while ((a %d == 0) /* while both a and b are */ && (b %d == 0)) { /* divisible by d, divide out */ a /= d; b /= d; r *= d; } /* this common factor and */ } /* multiply it into the gcd */ return r; /* return greatest common divisor */ } /* gcd() */ #endif /*---------------------------------------------------------------------- This function tests divisors up to the smaller of the two numbers. Hence, for numbers that are prime relative to each other, the time complexity is O(min(a,b)). Note that exchanging the two arguments in case a > b is done only to speed up the computation. The function would also work fine without this exchange, but would take longer to finish. Exercise: Note that in the loop we do no longer divide out all factors from a and b. Nevertheless we get a correct result. Why? ----------------------------------------------------------------------*/ #if (VERSION == 5) /* if to compile fifth version */ int gcd (int a, int b) { /* --- O(sqrt(min(a,b))) */ int d; /* test divisor, temporary buffer */ int i; /* increment of test divisor */ int e; /* exponent of prime factor in a */ int r = 1; /* result (greatest common divisor) */ if (a > b) { /* if a is larger than b, exchange */ d = a; a = b; b = d; } /* (to make a the smaller of the two) */ for (e = 0; a %2 == 0; e++) /* determine exponent of 2 in a */ a /= 2; /* and divide out prime factor 2 */ while ((--e >= 0) && (b %2 == 0)) { b /= 2; r *= 2; } /* compute factor 2 part of gcd */ for (e = 0; a %3 == 0; e++) /* determine exponent of 3 in a */ a /= 3; /* and divide out prime factor 3 */ while ((--e >= 0) && (b %3 == 0)) { b /= 3; r *= 3; } /* compute factor 3 part of gcd */ for (d = 5, i = 2; d*d <= a; d += i, i ^= 6) { for (e = 0; a %d == 0; e++) /* determine exponent of d in a */ a /= d; /* and divide out prime factor d */ while ((--e >= 0) && (b %d == 0)) { b /= d; r *= d; } /* compute factor d part of gcd */ } if (b %a == 0) /* test possible last prime factor */ r *= a; /* and multiply it into the gcd */ return r; /* return greatest common divisor */ } /* gcd() */ #endif /*---------------------------------------------------------------------- This function is an improvement over version 2 derived directly from the improvements of the prime number test. It tests for divisors only up to the square root of the smaller of the two numbers and skips multiples of two and three. The former brings the time complexity down to O(sqrt(min(a,b))), whereas the second only reduces the constant factor. ----------------------------------------------------------------------*/ int main (int argc, char *argv[]) { /* --- main function */ int a, b; /* numbers to compute the gcd of */ if (argc != 3) { /* if wrong number of arguments given */ printf("usage: %s a b\n", argv[0]); printf("compute the greatest common divisor of a and b\n"); return 0; /* print a usage message */ } /* and abort the function */ a = atoi(argv[1]); /* convert both arguments */ b = atoi(argv[2]); /* to integer numbers */ if ((a <= 0) || (b <= 0)) { /* test for correct arguments */ printf("%s: illegal arguments\n", argv[0]); return 1; } printf("%i\n", gcd(a, b)); /* compute and print gcd */ return 0; /* return 'ok' */ } /* main() */