Here, we study some computing problems that are related to Divisibility and the Division Algorithm.

#### Maxima

**GCD:** Maxima has a function to compute the GCD (Greatest Common Divisor) `gcd(a,b)`

.

*Example*. Find the GCD of 158 and 188.

`(%i1) gcd(158,188);`

`(%o1) 2`

**LCM:** Maxima has a function to compute the LCM (Least Common Multiple) `lcm(express_1,...express_n)`

. This function needs to be loaded by running the command `load("functs")`

.

*Example*. Find the LCD of 6, 8, 15, 32.

`(%i1) load("functs")$`

`(%i2) lcm(6,8,15,32);`

`(%o2) 480`

Bézout’s Lemma (see here) says that given two integers $a,b$, there is an integer solution $(x,y)$ to the equation $ax+by=(a,b)$. The Maxima command `gcdex`

implements the Euclidean algorithm. `gcdex(a,b)`

returns `[x, y, u]`

where $u=(a,b)$ and $ax+by=u$. So, one can find a solution of the equation $ax+by=(a,b)$ using `gcdex`

.

*Example*. Solve $2=158x+188y$.

We saw earlier that $(158,188)=2$.

`(%i2) gcdex(158,188);`

`(%o2) [25, - 21, 2]`

Hence, $x=25$ and $y=-21$.

One can also use the command `igcdex`

. `gcdex`

and `igcdex`

both implement the Euclidean algorithm. The only difference between them is that the argument of `igcdex`

must be integers while the arguments of `gcdex`

can be also polynomials. For instance,

(%i7) gcdex(x^2 + 1,x^3 + 4);

(%o7)/R/ [$\displaystyle\frac{x^2+4x-1}{17}$, $\displaystyle\frac{x+4}{17}$, 1]

#### Python

**GCD:** `Algorithm`

- The name of the Python code is
`gcd(a,b)`

and it implements the Euclidean algorithm to find the GCD of two integers $a$ and $b$. - If
`b == 0`

then `gcd(a,b)`

returns `0`

because every nonzero integer divides 0. - Otherwise, it returns
` gcd(b,a % b)`

. % is the *modulo* or *remainder operator*. `a % b`

returns the remainder when $a$ is divided by $b$ if $a\geq b$ or $a$ if $a<b$ because $a=b\cdot 0+a$. For example, `12345 % 2`

returns 1. - The loop continues until
`b==0`

.

Use your favorite editor to write the following python code (it is called a Python module) and save it as, for example, gcd.py. The Python interface IDLE has own editor. I also recommend Emacs and Geany for Python editors if you are running Python in command shell with iPython.

# Computing the GCD of two integers a and b recursively using the Euclidean algorithm
def gcd(a,b):
if (b == 0):
return a
return gcd(b, a%b)

On IDLE, the Python module can be run from its pull-down (or drop-down) menu under Run. Its shortcut key is F5. To run it on iPython, first change your current folder location to the one that contains the Python module gcd.py using `cd`

(change directory) command in command shell and then start iPython by typing `ipython`

(you can first start iPython and then change your current folder location as well using `cd`

within iPython) and pressing the enter key. Once iPython is loaded, do

$ ipython
Python 3.7.9 (default, Feb 4 2021, 01:17:56)
Type 'copyright', 'credits' or 'license' for more information
IPython 7.19.0 -- An enhanced Interactive Python. Type '?' for help.
In [1]: run gcd.py
In [2]: gcd(98,56)
Out[2]: 14

**The Euclidean Algorithm:** The Euclidean algorithm itself is a perfect computer algorithm. A Python code of the Euclidean algorithm is

# This Python module runs the Euclidean Algorithm to find the GCD of two integers a and b
def euclid(a,b):
r = a % b
while r > 0:
print (a, b, r)
a, b, r = b, r, b % r
return b

Running this Python module (named euclid.py) on iPython shows

$ ipython
Python 3.7.9 (default, Feb 4 2021, 01:17:56)
Type 'copyright', 'credits' or 'license' for more information
IPython 7.19.0 -- An enhanced Interactive Python. Type '?' for help.
In [1]: run euclid.py
In [2]: euclid(98,56)
98 56 42
56 42 14
Out[2]: 14

**LCM:** The LCM $[a,b]$ can be obtained by $\frac{ab}{(a,b)}$ (see here).

`Algorithm`

: The name of the Python code is `lcm(a,b)`

. It computes the LCM of two integers $a$ and $b$.

- Find the GCD of $a$ and $b$ by running the Python program
`gcd(a,b)`

. - Compute the LCM of $a$ and $b$ using the formula
`a*b/gcd(a,b)`

.

Its Python code is

# Compute the LCM of two integers a and b
Recursive function to return the GCD of a and b
def gcd(a,b):
if (b == 0):
return a
return gcd(b, a%b)
# Function to return the LCM of a and b
def lcm(a,b):
return a*b/gcd(a,b)

Running this Python code (named lcm.py) on iPython shows

In [3]: run lcm.py
In [4]: lcm(15,20)
Out[4]: 60.0

**Linear Diophantine Equations:** As seen here, a linear diophantine equation

\begin{equation}

\label{eq:lindiophant}

ax+by=c

\end{equation}

is solvable if and only if $(a,b)|c$. If $b|a$ then $(a,b)=b$ and

$$x=0,\ y=\frac{c}{b}$$

is a solution. If $b\not|a$, then write $a=bq+r$ and substitute into \eqref{eq:lindiophant}:

$$(bq+r)x+by=c$$

This can be rearranged as

$$b(qx+y)+rx=c$$

Setting $u=qx+y$ and $v=x$, we obtain the equation

\begin{equation}

\label{eq:lindiophant2}

bu+rv=c

\end{equation}

with smaller coefficients. If we can solve \eqref{eq:lindiophant2}, we recover a solution of the original equation as

$$x=v,\ y=u-qv$$

The process terminates (according to the Euclidean algorithm) with an equation where the second coefficient divides the first. This scheme is embodied in the following Python code

# Solve a linear diophantine equation using the Euclidean algorithm
# The function divmod returns a pair of numbers, the quotient and remainder, which we store in q, r.
def diosolve(a,b,c):
q, r = divmod(a,b)
if r == 0:
return( [0,c/b] )
else:
sol = diosolve( b, r, c )
u = sol[0]
v = sol[1]
return( [ v, u-q*v ] )

Running this Python code (named diosolve.py) on iPython to solve the linear diophantine equation $12345x+54321y=3$:

In [9]: run diosolve.py
In [10]: diosolve(12345, 54321, 3)
Out[10]: [3617.0, -822.0]

Note that $(12345,54321)=3$, so of course the diophantine equation is solvable.

*Refernces*.

[1] Jim Carlson, A Short Course in Python for Number Theory, Draft of May 21, 2004

[2] Maximal 5.29.0 Manual, 29. Number Theory