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

1. 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$.
2. If b == 0 then gcd(a,b) returns 0 because every nonzero integer divides 0.
3. 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.
4. 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$.

1. Find the GCD of $a$ and $b$ by running the Python program gcd(a,b).
2. 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