Binary GCD (Stein’s Algorithm) in C

1. Algorithm Description

Binary GCD also known as Stein’s Algorithm is an algorithm that computes the greatest common divisor of two (positive) numbers . Discovered in 1967 by the Israeli programmer Josef Stein, it’s an alternative to the classical Euclid’s Algorithm, and is considered to be more efficient than this as it’s replacing divisions and multiplications with bitwise operations . The algorithm is recursive by nature, but loops can be used instead of recursion .

The algorithm can be described by the following rules (here and here) . Note that by B_GCD we will refer to a function that returns the greatest common divisor of two positive numbers .

1. B_GCD(0, 0) is not defined, but for convenience we will consider it 0 (zero), so B_GCD(0,0) = 0 ;
2. B_GCD(num1, 0) = num1 and B_GCD(0, num2) = num2 . The reason for this is beacuse zero is divided by everything ;
3. If num1 and num2 are even, B_GCD(num1, num2) = 2 * B_GCD(num1/2, num2/2), as 2 is a common divisor .
4. If num1 is even and num2 is odd, B_GCD(num1, num2) = B_GCD(num1 /2, num2), as 2 is not a common divisor . The steps are the same if num1 is odd and num2 is even : B_GCD(num1, num2) = B_GCD(num1, num2/2) .
5. If both num1 and num2 are odd, then if num1 >= num2 the B_GCD(num1, num2) = B_GCD((num1-num2)/2, num2), else B_GCD(num1, num2) = B_GCD((num2-num1)/2, num1) .
6. Step 4 and 5 are repeated until num1 = num2, or num1 = 0 .

We can also use pseudo code to describe the above algorithm .

Recursive Version of Binary GCD (Stein Algorithm)

FUNCTION  B_GCD(num1, num2)
  IF num1 = num2 THEN
    RETURN num1
  IF num1 = 0 AND num2 = 0 THEN
    RETURN 0
  IF num1 = 0 THEN
    RETURN num2
  IF num2 = 0 THEN
    RETURN num1
  IF num1 IS EVEN AND num2 IS EVEN THEN
    RETURN (B_GCD(num1/2, num2/2) * 2)
  IF num1 IS EVEN AND num2 IS ODD THEN
    RETURN B_GCD(num1/2, num2)
  IF num1 IS ODD AND num2 IS EVEN THEN
    RETURN B_GCD(num1, num2/2)
  IF num1 IS ODD AND num2 IS ODD THEN
    IF num1 >= num2 THEN
      RETURN B_GCD((num1-num2)/2, num2)
    ELSE
      RETURN B_GCD((num2-num1)/2, num1)

The loop-version of the Binary GCD Algorithm

FUNCTION B_GCD(num1, num2)
  power_of_two := 0
  IF (num1 = 0 OR num2 = 0) THEN
    RETURN num1 | num2
  WHILE ((num1 IS EVEN) AND (num2 IS EVEN))
    num1 := num1 / 2
    num2 := num2 / 2
    power_of_two := power_of_two + 1
  DO
    WHILE(num1 IS EVEN)
      num1 := num1 / 2
    WHILE(num2 IS EVEN)
      num2 := num2 / 2
    IF (num1 >= num2) THEN
      num1 := (num1 - num2) / 2
    ELSE
      tmp  := num1
      num1 := (num2 - num1) / 2
      num2 := tmp
  WHILE NOT ((num1 = num2) OR (num1 = 0))
  RETURN num2 * power_of_two

Euclid’s Algorithm: Reducing fraction to lowest terms

In the last article I’ve described Euclid’s algorithm for finding the greatest common divisor of two given numbers .A simple programming exercise I found in the book called “Algorithms in C” (Sedgewick) asks us to reduce a given fraction to lowest terms, using the Euclid’s Algorithm .

Eg.
Fraction = 36 / 120
Fraction = Reduce(Fraction) = 3 / 12

1. C Implementation

Step 1:

We will start by defining the data-structure we will use for encapsulating the fraction information (numerator and denominator):

/**
* Fraction data-structure:
* f = numerator / denominator
*/
typedef struct s_fraction {
  int numerator;
  int denominator;
} Fraction;

Step 2:

Usually when we play with C structures we should have two functions, for allocating and de-allocating memory associated with a structure (if you have an object-oriented background, this step will be equivalent with writing a constructor / destructor pair ) .

/**
* Allocates memory for a fraction structure .
* Initialize fraction structure .
* @param numerator
* @param denominator
* @return A pointer to the fraction structure on heap . (memmory should
be de-allocated manually)
*/
Fraction *fraction_new(int numerator, int denominator)
{
  Fraction *f = NULL;
  if (!denominator) {
    fprintf(stderr, "Invalid fraction. Denominator cannot be 0 (zero).n");
    return (f);
  }
  f = malloc(sizeof(*f));
  if (!f) {
    MALLOC_FAIL;
  }
  f->numerator = numerator;
  f->denominator = denominator;
  return (f);
}

/**
* De-allocates the memory associated with a given fraction .
* @param f The fraction to be free'd .
*/
void fraction_delete(Fraction *f)
{
  if (f) {
    free(f);
  }
}

An important aspect is that the fraction_new() function fails if a 0 (zero) denominator is given .

Step 3:

On this step we will write a function that reduces the fraction to lowest terms . For this we will also need to determine the greatest common divisor of the fraction’s denominator and numerator in order to divide the fraction

Euclid’s Algorithm

Recently I’ve started to implement (or reimplement) the most common algorithms a software developer should know . One of the nicest books I found on this topic is Algorithms in C (Robert Sedgewick) . Of course, there is this and this, but lately I am more interested on the “implementation” side of things than on maths and theory .

1. Algorithm Description

Euclid’s Algorithm is an efficient method for calculating the greatest common divisor of two numbers (aka GCD) . The  GCD of two numbers is considered the largest number that divides both of them (without leaving a reminder) .

Euclid’s algorithm is based on the principle that given two numbers a and b, if a is greater than b than the greatest common divisor of a and b is the same with the common divisor of b and (b – a) . If we transform the above property into “code” (actually pseudo-code) the algorithm looks like this:

FUNCTION GCD(num1, num2)
  WHILE num1 > 0
    IF num1 < num2
      SWAP (num1, num2)
    num1 := num1 - num2
  RETURN num2

The above pseudo-code is called the subtraction-variant . We can of course replace the repetitive subtractions with one division . The division-based version looks like this (in pseudo-code):

FUNCTION GCD(num1, num2)
  WHILE num1 > 0
    tmp  := num1
    num1 := num2 MOD num1
    num2 := tmp
  RETURN num2

There is also a recursive version of the algorithm:

FUNCTION GCD(num1, num2)
  IF num1 <> 0 THEN
    RETURN GCD(num2 MOD num1, num1)
  RETURN num2