Add extended Euclid Algorithm

learning-lua/algo/Euclid-tests.lua

@@ -1,21 +1,52 @@
 local gcd = require "algo.Euclid".gcd
+local xgcd = require "algo.Euclid".xgcd
 
-local function negative_test(a, b)
+local DEBUG = require "algo.Debug":new(false)
+local debugf = DEBUG.debugf
+
+
+local function negative_test(f, a, b)
   assert(a < 0 or b < 0)
-  local ok, err = pcall(gcd, a, b)
+  local ok, err = pcall(f, a, b)
   assert(not ok)
   assert(string.find(err, "input numbers must not be negative"))
 end
 
--- edge cases
+local function verify_xgcd(s, t, d, m, n)
+  debugf("%d * %d + %d * %d = %d", m, s, n, t, d)
+  assert(m * s + n * t == d)
+end
+
+-- gcd edge cases
 assert(gcd(0, 0) == 0)
 assert(gcd(1, 0) == 1)
 assert(gcd(0, 1) == 1)
 
--- negative is not allowed
-negative_test(0, -1)
-negative_test(-1, 0)
+-- negative is not allowed for gcd
+negative_test(gcd, 0, -1)
+negative_test(gcd, -1, 0)
 
--- positive tests
+-- gcd positive tests
 assert(gcd(123, 54) == 3)
 assert(gcd(54, 123) == 3)
+
+-- negative is not allowed for xgcd
+negative_test(xgcd, 0, -1)
+negative_test(xgcd, -1, 0)
+
+-- xgcd edge cases
+local d, m, n = xgcd(0, 0)
+verify_xgcd(0, 0, d, m, n)
+assert(d == 0 and m == 0 and n == 0)
+
+verify_xgcd(1, 0, xgcd(1, 0))
+verify_xgcd(0, 1, xgcd(0, 1))
+
+-- xgcd positive tests
+verify_xgcd(240, 46, xgcd(240, 46))
+verify_xgcd(46, 240, xgcd(46, 240))
+verify_xgcd(123, 54, xgcd(123, 54))
+verify_xgcd(54, 123, xgcd(54, 123))
+
+verify_xgcd(84, 33, xgcd(84, 33))
+verify_xgcd(270, 192, xgcd(270, 192))

learning-lua/algo/Euclid.lua

@@ -10,9 +10,50 @@ end
 function Euclid.gcd(a, b)
   assert(a >= 0 and b >= 0, "input numbers must not be negative")
   if a < b then
-    return _gcd(b, a)
+    return _gcd(b, a) -- tail recursive
   else
-    return _gcd(a, b)
+    return _gcd(a, b) -- tail recursive
+  end
+end
+
+---(See Euclid.pdf on the idea behind this specific implementation.)
+---@param n number
+---@param d number
+---@param npair table the pair of numbers for `n`.
+---@param dpair table the pair of numbers for `d`.
+---@return number d greatest common divisor of the `a` and `b` parameters passed to `Euclid.xgcd(a, b)`.
+---@return number s integer such that `s*a + t*b = d`.
+---@return number t integer such that `s*a + t*b = d`.
+local function _xgcd(n, d, npair, dpair)
+  if d == 0 then
+    return n, npair[1], npair[2]
+  end
+  local q = math.floor(n / d)
+  local r = n % d
+  local dpair_1, dpair_2 = table.unpack(dpair)
+  dpair[1], dpair[2] = npair[1] - q * dpair_1, npair[2] - q * dpair_2
+  npair[1], npair[2] = dpair_1, dpair_2
+  return _xgcd(d, r, npair, dpair) -- tail recursive
+end
+
+---Computes the greatest common divisor of `a` and `b` using Extended Euclid algorithm.
+---Returns the greatest common divisor of `a` and `b`,
+---and integers `s` and `t` such that `s*a + t*b = d`.
+---@param a number non-negative integer.
+---@param b number non-negative integer.
+---@return number d greatest common divisor of `a` and `b`.
+---@return number s integer such that `s*a + t*b = d`.
+---@return number t integer such that `s*a + t*b = d`.
+function Euclid.xgcd(a, b)
+  assert(a >= 0 and b >= 0, "input numbers must not be negative")
+  if a == 0 and b == 0 then
+    return 0, 0, 0 -- edge case.
+  end
+  if a < b then
+    local d, m, n = _xgcd(b, a, {1, 0}, {0, 1})
+    return d, n, m
+  else
+    return _xgcd(a, b, {1, 0}, {0, 1}) -- tail recursive
   end
 end