Thursday, 11 February 2016

Finding gcd by Euclidean algorithm and extended Euclidean algorithm

Task: find the greatest common divisor of two integers a and b using python

The first program relies on the fact that if $d\mid a$ and $d\mid b$, then $d|(a-b)$. In the third program, we need to ensure that b is not zero because division by zero is not defined. Line 21 follows from Euclidean algorithm. (In fact, it does not matter whether a or b is the larger number. Here is an example for illustration. When a = 210 and b = 7, 210 is assigned to 7 and 7 to 0. The program then stops and returns 7. When a = 7 and b = 210, 7 is assigned to 210 and 210 is assigned to 7 % 210. 7 % 210 is 0 with remainder 7. Now, 210 is assigned to 7 and 7 to 0. The program terminates and returns 7.)

The fourth program is known as the extended Euclidean algorithm and it takes some time to understand. Note that for integers $a,b$ where $a>b$, $a=qb+r=\lfloor a/b \rfloor b+a-\lfloor a/b \rfloor b$, or in python $a=a//b*b+a\%b$. Recall that each iteration in the Euclidean algorithm replaces $(a,b)$ by $(b,a\;\text{mod}\; b)$. We can formulate this as a matrix multiplication: $$\begin{pmatrix}b & a\;\text{mod}\;b\end{pmatrix}=\begin{pmatrix}a&b\end{pmatrix}\begin{pmatrix}0 & 1\\ 1 & -q\end{pmatrix},\qquad q=\lfloor a/b \rfloor.$$Suppose that the algorithm terminates after $k$ iterations, then the operations performed on $\text{gcd}(a,b)$ can be summarised in matrix notation by $$\begin{pmatrix}\text{gcd}(a,b)&0\end{pmatrix}=\begin{pmatrix}a&b\end{pmatrix}\begin{pmatrix}0 & 1\\ 1 & -q_1\end{pmatrix}\begin{pmatrix}0 & 1\\ 1 & -q_2\end{pmatrix}\cdots\begin{pmatrix}0 & 1\\ 1 & -q_k\end{pmatrix}$$ or $$\begin{pmatrix}\text{gcd}(a,b)&0\end{pmatrix}=\begin{pmatrix}a&b\end{pmatrix}\begin{pmatrix}x & u\\ y & v\end{pmatrix}.$$Since we want to find the triple $(d,x,y)$, namely, the gcd and the first column of the resulting matrix, we want to keep track of the product of matrices along with the remainder computation. Therefore, we start with the matrix $$\begin{pmatrix}x & u\\ y & v\\ a & b\end{pmatrix}=\begin{pmatrix}1 & 0\\ 0 & 1\\ a & b\end{pmatrix},$$ or set $x,y,u,v=1,0,0,1$. Note that $$\begin{pmatrix}x & u\\ y & v\\ a & b\end{pmatrix}\begin{pmatrix}0 & 1\\ 1 & -q\end{pmatrix}=\begin{pmatrix}u & x-uq\\ v & y-vq\\ b & a-bq\end{pmatrix}$$ It is then natural to assign $x$ to $u$, $y$ to $v$; $u$ to $x-uq$, $v$ to $y-vq$. The rest is similar to program $3$. Here's a second interpretation. The algorithm aims to compute the tuple $(d,x,y)$ such that $ax+by=d$. We know that $$a\cdot 1+b\cdot 0=a\\ a\cdot 0+b\cdot 1=b.$$ We then set $$a\cdot x+b\cdot y=a\quad(1)\\ a\cdot u+b\cdot v=b\quad(2)$$ where $x,y,u,v$ are $1,0,0,1$. Note that $\text{gcd}(a,b)=\text{gcd}(b,r=a\;\text{mod}\; b)$. We thus have a new pair of equations: $$(1)-(2)\cdot q,\quad a\cdot(x-uq)+b\cdot(y-vq)=r,\quad q=\lfloor a/b \rfloor. \\ a\cdot x+b\cdot y=b$$ Assigning $a$ to $b$, $b$ to $r$ and changing the values of $x,y,u,v$ accordingly, the program will eventually terminate when $b$ becomes $0$. It follows that $$a\cdot u+b\cdot v=0\\ a\cdot x+b\cdot y=\text{gcd}(a,b).$$Then we take the second equation and we have $(d,x,y)$.

In the sixth program, we are still relying on the fact that $\text{gcd}(a,b)=\text{gcd}(b,r=a\;\text{mod}\; b)$. We have $$ax+by=\text{gcd}(a,b)\\ bx'+ry'=\text{gcd}(a,b).$$ Put $x=y',y=x'$. Then $$\begin{align}ax+by&=ay'+bx'\\ &=(qb+r)y'+bx'\\ &=qby'+(ry'+bx')\\ ay'+bx'&=qby'+\text{gcd}(a,b)\\ \Rightarrow ay'+b(x'-qy')&=\text{gcd}(a,b).\end{align}$$ Finally, we put $x=y',y=x'-qy'$, where $x,y$ are the values we want.

Lessons learnt:
  • = vs ==: The first one assigns a value to the variable on the left; the second one means that the value of the variable is that on the right.
  • % // operators
  • while loop
  • iterative (loop) vs recursive (a function calling itself repeatedly)
  • multiple assignment: In python, we don't need a swap function because a, b = b, a works sufficiently without defining a temp variable, c. This is useful when we want a to be larger than b:
if a < b:
    a, b = b, a
else:
    pass

More to know:
subtractive algorithm 123
binary gcd algorithm 1
extended Euclidean algorithm

Reference:
program 4 123
program 6 123

No comments:

Post a Comment