INTERNAL/REVIEW: Calculating the greatest common divisor using IDL.

Topic:

This tech tip provides an example that will calculate the greatest common divisor of two numbers.Discussion:
The greatest common divisor of two integers a and b, written (a,b), is the largest divisor common to a and b. For example, (24,64) = 8. The GCD routine provided in this tech tip will return this information.

IDL> print, gcd(24,64)
8

Furthermore, given (a,b) = g, there exists integers x and y such that ax + by = g. Using the previous example, these values would be 3 and -1 (3*24 + (-1)*64 = 3). This information can optionally be returned via the XY keyword to the GCD function.

IDL> print, gcd(24,64, xy = xy)
8
IDL> print, xy
3 -1
IDL> print, 24*xy[0] + 64*xy[1]
8
Solution:

function gcd, a, b, xy = xy
; The GCD routine returns the greatest
; common divisor of A and B.
;
; The XY keyword is to be set to a named
; variable which will contain the integer
; array XY such that:
;
;    A * XY[0] + B * XY[1] = GCD(A,B)
 
  if n_params() ne 2 then return, -1
  aa = max([a,b], min = bb)
  aa = abs(aa) & bb = abs(bb)
 
  xy = [0ll,0] & xy_1 = [0ll,1] & xy_2 = [1ll,0]
 
  while bb ne 0 do begin
    q = aa/bb
    r = aa mod bb
    xy = xy_2 - q * xy_1
    aa = bb & bb = r
    xy_2 = xy_1 & xy_1 = xy
  endwhile
 
  if a gt b then $
    xy = xy_2 $
  else $
    xy = reverse(xy_2, /over)
 
  if (a lt 0) then xy = xy*[-1,1]
  if (b lt 0) then xy = xy*[1,-1]
 
  return, aa
end