/*
  Author:  Pate Williams (c) 1997

  "Algorithm 1.3.5 (Binary GCD). Given two
  non-negative integers a and b, this algorithm
  finds their GCD." -Henri Cohen- See "A Course
  in Computational Algebraic Number Theory" by
  Henri Cohen pages 14-15.
*/

#include <stdio.h>

long binary_gcd(long a, long b)
{
  long k = 0, r, t;

  if (a < b) t = a, a = b, b = t;
  if (b == 0) return a;
  r = a % b, a = b, b = r;
  printf("%ld %ld %ld\n", a, b, r);
  if (b == 0) return a;
  while (!(a & 1) && !(b & 1))
    k++, a >>= 1, b >>= 1;
  while (!(a & 1)) a >>= 1;
  while (!(b & 1)) b >>= 1;
  printf("%ld %ld\n", a, b);
  L4:
  t = (a - b) >> 1;
  printf("%ld %ld %ld\n", a, b, t);
  if (t == 0) return (1 << k) * a;
  while (!(t & 1)) t >>= 1;
  if (t > 0) a = t; else b = - t;
  goto L4;
}

int main(void)
{
  long a = 1764, b = 868, ap = 33, bp = 77;

  printf("a = %ld b = %ld\n", a, b);
  printf("gcd(a, b) = %ld\n", binary_gcd(a, b));
  printf("a = %ld b = %ld\n", ap, bp);
  printf("gcd(a, b) = %ld\n", binary_gcd(ap, bp));
  return 0;
}