A matrix formulation of extended Euclid algorithm.

I have written about Extended Euclid Algorithm in the past. However more recently while solving a problem I had to express it in a matrix form which I thought is something worthwhile.

Recall the problem statement: Given two non-zero integers $a,b$ (assume $a\geq b$), the problem is to determine integers $x,y$ such that $\mbox{gcd}(a,b) = a\times x + b \times y$. Notice that the original Euclid's algorithm can compute the gcd but not the coefficients $x,y$ whose linear combination yields the gcd. The idea here is to use the quotient instead of the remainder as in the original algorithm. We express the operands in each iteration a linear combination of the original inputs and use the following to update as follows. The quotient in $i^{th}$ iteration is denoted by $q_i$ leading to the following formulation.

$$\begin{array}{cc} \mbox{gcd}(a,b) = \mbox{ExGcd}(\left[\begin{array}{c} a \\ b \end{array}\right]) &\\ \left [ \begin{array}{cc} x_{0,0} & y_{0,0}\\ x_{1,0} & y_{1,0} \end{array} \right ] = \left [ \begin{array}{cc}1 & 0\\ 0 & 1\end{array} \right ] & \mbox{initialization } i=0\\ \mbox{ExGcd}(\left [\begin{array}{cc} a & b \end{array}\right ] \left [ \begin{array}{cc} x_{0,i} & y_{0,i} \\ x_{1.i} & y_{1,i} \end{array} \right]) = \mbox{ExGcd}(\left[\begin{array}{cc} a & b \end{array}\right] \left [ \begin{array}{cc} x_{0,i-1} & y_{0,i-1} \\ x_{1,i-1} & y_{1,i-1} \end{array} \right ] \left [ \begin{array}{cc} 0 & 1 \\ 1 & -q_{i-1} \end{array} \right ]) & i > 0 \\ q_i = \mbox{Quotient}(x_{0,i}a+x_{1,i}b, y_{0,i}a+y_{1,i}b) & \\ \end{array}$$


Following is an implementation using std::div:
// Input: a,b
// Output: x, y (y0,y1)

void ExtendedEuclid(long a, long b, long& y0, long& y1)
{
 long x0, x1, t0, t1;

 if (a < b) { std::swap(a, b); }
 // Initialize:
 x0 = 1; x1 = 0; // a = 1 * a + 0 * b
 y0 = 0; y1 = 1; // b = 0 * a + 1 * b
 
 std::ldiv_t div_res = std::div(a, b);
 while (div_res.rem) {
  t0 = y0; t1 = y1;

  y0 = x0 - (div_res.quot * y0);
  y1 = x1 - (div_res.quot * y1);
  x0 = t0; x1 = t1;

  // Invariant: (x0*a+x1*b >= y0*a+y1*b)
  div_res = std::div(x0*a+x1*b, y0*a+y1*b);
 }
}

Following are few examples obtained by calling ExtendedEuclid(a,b,x,y):

a=25375,b=667 x=1, y=-38, gcd=29 
29 = 25375 * 1 + 667*38

a=17765,b=13481 x=-343, y=452, gcd=17
17 = 17765 * -343 + 13481 * 452

a=19109,b=18321 x=-23, y=24, gcd=197
197 = 19109 * -23 + 18321 * 24