*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