Author Archives: Sung Lee

About Sung Lee

I am a mathematician and a physics buff. I am interested in studying mathematics inspired by theoretical physics as well as theoretical physics itself. Besides, I am also interested in studying theoretical computer science.

Black Hole 1: Schwarzschild Black Hole

Leave a reply

In light of the recent exciting press release about the supermassive black hole at the center of our own Milky Way galaxy, I am going to discuss a brief introduction to black holes and other related concepts in a series of blog articles for those who wish to understand better these amazing monstrous heavenly creatures. For starter, this first article is about Schwarzschild black hole, the simplest black hole model and also the first exact solution of the Einstein’s Field Equation.

The Supermassive Black Hole (Sagittarius A*) at the Center of Milky Way Galaxy.

First, what is a black hole?

A black hole is a hypothetical (though its existence is no longer in doubt) region of spacetime where gravity is so immense that nothing, even light (the lightest matter in the universe) can escape from it.

How did we come to know about a black hole? Was it just an imagination of some mad scientists?

No imagination. We came to know about the theoretical notion of a black hole from Einstein’s general relativity. But long before general relativity came out, French mathematician Pierre-Simon Laplace had already thought about a celestial object whose gravity is so strong that even light can’t escape. It is discussed in his five-volume book Mécanique Céleste (Celestial Mechanics).

In general relativity, black holes are solutions of the Einstein’s Field Equations.

What are the Einstein’s Field Equations?

They are the fundamental equations of Einstein’s theory of gravitation. Such equations are expected to satisfy the following requirements:

  • The field equations should be independent of coordinate systems, i.e. they should be tensorial.
  • They should be partial differential equations of second order for the components $g_{ij}$ of the unknown metric tensor. If you want to know what a metric tensor is, go deeper here.
  • They are a relativistic generalization of the Poisson equation of Newtonian gravitational potential $$\nabla^2\Phi=4\pi G\rho$$ where $\rho$ is the mass density.
  • Since the energy-momentum tensor $T_{ij}$ is the special relativistic analogue of the mass density, it should be the source of the gravitational field.
  • If the space is flat, $T_{ij}$ should vanish.

The equations satisfying all these requirements are given by \begin{equation}\label{eq:EFT}G_{ij}=8\pi GT_{ij}\end{equation} where $G_{ij}=R_{ij}-\frac{1}{2}Rg_{ij}$. $G_{ij}$ is called the Einstein tensor. $R_{ij}$ and $R$ are, respectively, Ricci curvature tensor and scalar curvature. If you want to know what they are, go deeper here. Interestingly, the field equations \eqref{eq:EFT} were derived independently and almost simultaneously by Albert Einstein and David Hilbert in 1915.

So, how do we solve the Field Equations?

The field equations are nonlinear partial differential equations. Partial differential equations are, in general,  difficult to solve with equations alone. We impose additional (physical) assumptions and conditions (such as boundary conditions and initial conditions) so the equations can becomes simpler enough for us to be able to solve. First, we have no information on the source of gravity, i.e. no information on the energy-momentum tensor $T_{ij}$. Instead, we consider the field equations outside the field-producing masses. Thus, we obtain the vacuum field equations $$R_{ij}-\frac{1}{2}Rg_{ij}=0$$ Without loss of generality, we may assume that the metric tensor $g_{ij}$ is diagonal, i.e. $g_{ij}=0$ if $i\ne j$. This results in $R_{ij}=0$ if $i\ne j$ from the vacuum field equations. For $i=j$, since the scalar curvature $R$ is given by $R=\sum_ig^{ii}R_{ii}$, from the vacuum field equations, we have $(n-2)R=0$. For $n=4$, $R=0$ and therefore the vacuum field equations turn into \begin{equation}\label{eq:ricciflat}R_{ij}=0\end{equation} i.e. vanishing Ricci curvature tensor. \eqref{eq:ricciflat} is simpler than the original vacuum field equations but by no means it is easy to solve. The simplest solution of \eqref{eq:ricciflat} was obtained in 1915 by a German physicist Karl Schwarzschild while serving in the army during World War I. This is also the first exact solution to the Einstein’s field equations. Schwarzschild assumed a static spherical symmetric metric as a solution of the vacuum equations \eqref{eq:ricciflat}. The ansatz for a static isotropic metric is given by \begin{equation}\label{eq:ansatz}ds^2=-A(r)dt^2+B(r)dr^2+r^2(d\theta^2+\sin^2\theta d\phi^2)\end{equation} In addition, we want the solution to become the flat spacetime far away from the masses i.e. the solution is asymptotically flat, so we require $$\lim_{r\to\infty}A(r)=\lim_{r\to\infty}B(r)=1$$ Now, we are ready to solve the vacuum field equations \eqref{eq:ricciflat}. From the ansatz \eqref{eq:ansatz}, we have $$g_{tt}=-A(r),\ g_{rr}=B(r),\ g_{\theta\theta}=r^2,\ g_{\phi\phi}=r^2\sin^2\theta$$ The nonzero Christoffel symbols are computed to be \begin{align*}\Gamma_{rr}^r&=\frac{B’}{2B},\ \Gamma_{tt}^r=\frac{A’}{2B},\ \Gamma_{\theta\theta}^r=-\frac{r}{B}\\\Gamma_{\phi\phi}^r&=-\frac{r\sin^2\theta}{B},\ \Gamma_{\theta r}^\theta=\Gamma_{\phi r}^\phi=\frac{1}{r},\\\Gamma_{tr}^t&=\frac{A’}{2A},\ \Gamma_{\phi\phi}^\theta=-\sin\theta\cos\theta,\ \Gamma_{\phi\theta}^\phi=\cot\theta\end{align*} Here, ‘ stands for $\frac{d}{dr}$. Using these Christoffel symbols, the Ricci tensors are computed to be \begin{align*}R_{tt}&=\frac{A^{\prime\prime}}{2B}-\frac{A’}{4B}\left(\frac{A’}{A}+\frac{B’}{B}\right)+\frac{A’}{rB}\\R_{rr}&=-\frac{A^{\prime\prime}}{2A}+\frac{A’}{4A}\left(\frac{A’}{A}+\frac{B’}{B}\right)+\frac{B’}{rB}\\R_{\theta\theta}&=1-\frac{1}{B}-\frac{r}{2B}\left(\frac{A’}{A}-\frac{B’}{B}\right)\\R_{\phi\phi}&=\sin^2\theta R_{\theta\theta}\end{align*} From the first two equations of the Ricci tensors above and by requirung them to be 0, $$BR_{tt}+AR_{rr}=\frac{1}{rB}(A’B+AB’)=\frac{1}{rB}(AB)’=0$$ So, we have $A(r)B(r)$=constant. Asymptotic flatness implies that $A(r)B(r)=1$. $R_{\theta\theta}=0$ along with $B(r)=\frac{1}{A(r)}$ and $A’B+AB’=0$ results in $A+rA’=(rA)’=1$. Hence, $A(r)=1+\frac{C}{r}$, where $C$ is a constant. By weak field approximation (Newtonian limit) , we obtain \begin{align*}A(r)=-g_{00}&=1+\frac{2\Phi}{c^2}\\&=1-\frac{2MG}{c^2r}\end{align*} where $\Phi=-\frac{GM}{r}$. This determines the constant $C=-\frac{2MG}{c^2r}$. Finally, the Schwarzschild solution is given by \begin{equation}\label{eq:bh}ds^2=-\left(1-\frac{2MG}{c^2r}\right)c^2dt^2+\left(1-\frac{2MG}{c^2r}\right)^{-1}dr^2+r^2(d\theta^2+\sin^2\theta d\phi^2)\end{equation} As I mentioned in the beginning, this is the simplest black hole model and also is the first exact solution of the Einstein’s Field Equations. It turns out that Schwarzschild solution is the only spherically symmetric solution to the vacuum Einstein equation \eqref{eq:ricciflat} due to

Birkhoff’s Theorem

In 1923, George David Birkhoff proved:

Theorem. Any spherically symmetric solution of the vacuum Einstein’s field equations \eqref{eq:ricciflat} must be static and asymptotically flat.

The Shape of a Black Hole

Is the converse of Birkhoff’s theorem true? Namely, is any static and asymptotically flat solution of the vacuum field equations necessarily spherically symmetric? I don’t know the answer, though I assume the answer is known. While the notion of a black hole originated from general relativity as a solution of the field equations, general relativity does not tell us how a black hole can be actually formed. A mechanism of black hole formation was provided by nuclear physics as a consequence of the gravitational collapse of a dead star. From our observations of the shape of massive objects , one would think the shape of a black hole is a sphere also. Since we don’t know (I mean I don’t know) for sure, it would be interesting to ponder whether the shape of a black hole different from a sphere can exist, for example a toroidal (doughnut shaped) black hole. If not for Schwarzschild black hole, how about for a rotational black hole?

Update: Now I know the answers to questions about the converse of Birkhoff’s theorem and the shape of a black hole. I wrote about them here. At the time of writing this article I knew even lesser than what I know now about GR, which is still pretty much close to nothing.

The Event Horizon

I conclude this article by mentioning that the Schwarzschild solution in \eqref{eq:bh} only describes the outside of the black hole. What determines the inside and the outside of a black hole? It is the event horizon. The event horizon is where $r=r_g=\frac{2MG}{c^2}$. $r_g$ is called the Schwarzschild radius.The Schwarzschild solution requires that $r>r_g$. The Schwarzschild solution unfortunately does not tell anything about whatever happens inside the black hole. We will never know for sure what happens once an object falls into a black hole past its event horizon unless we have information on the source of its gravity, i.e. the stress-energy tensor $T_{ij}$.

The Diagonalization of a Riemannian Metric

Leave a reply

Let $(M^n, ds^2)$ denote an $n$-dimensional Riemannian manifold. The Riemannian metric $ds^2$, in general, is given, in terms of local coordinates $x’_1,\cdots,x’_n$, by $$ds^2=\sum_{i,j}h_{ij}(x’_1,\cdots,x’_n)dx’_idx’_j$$ Since the metric tensor $(h_{ij})$ is symmetric, by the finite dimensional spectral theorem, there exists an orthogonal matrix $O$ such that $O^t(h_{ij})O=(g_{ii})$. Let $$\begin{pmatrix}dx_1\\\vdots\\dx_n\end{pmatrix}=O^t\begin{pmatrix}dx’_1\\\vdots\\dx’_n\end{pmatrix}$$ Then \begin{align*}ds^2&=\sum_{i,j}h_{ij}(x’_1,\cdots,x’_n)dx’_idx’_j\\&=\begin{pmatrix}dx’_1\\\vdots\\dx’_n\end{pmatrix}^t(h_{ij})\begin{pmatrix}dx’_1\\\vdots\\dx’_n\end{pmatrix}\\&=\left[O\begin{pmatrix}dx_1\\\vdots\\dx_n\end{pmatrix}\right]^tO(g_{ii})O^t\left[O\begin{pmatrix}dx_1\\\vdots\\dx_n\end{pmatrix}\right]\\&=\begin{pmatrix}dx_1\\\vdots\\dx_n\end{pmatrix}^t(g_{ii})\begin{pmatrix}dx_1\\\vdots\\dx_n\end{pmatrix}\\&=\sum_{i}g_{ii}dx_i^2\end{align*} Therefore, without loss of generality, we may assume that a Riemannian metric is of the diagonal form $$ds^2=\sum_{i}g_{ii}(x_1,\cdots,x_n)dx_i^2$$

An Easier Way to Solve Quadratic Equations

Leave a reply

The Greek mathematician Heron (100? B.C.) is credited by some as the first (historically known) person who studied quadratic equations. He considered questions like: “There is a square. The sum of its area and circumference is 896. What is its length (or width)?” But some simple quadratic equations were already known to Mesopotamians 4,000 years ago. Along with linear equations, quadratic equations are among the oldest subjects of mathematics that have been studied. So, it is amazing and even somewhat amusing that there is still something new to find out about the solution of quadratic equations. I recently came across an article on the arXiv titled “A Simple Proof of the Quadratic Formula” by Po-Shen Loh at Carnegie Mellon. His proof also provides an alternative way to solve quadratic equations without dealing with factoring quadratic polynomials or having to remember the quadratic formula. I believe students will find it much easier than the conventional methods (factoring, completing the square, or the quadratic formula). First, here is the proof. Any quadratic equation can be written as a monic equation $$x^2+bc+c=0$$ Suppose that its solutions are $\alpha$ and $\beta$. Then $$x^2+bx+c=(x-\alpha)(x-\beta)=x^2-(\alpha+\beta)x+\alpha\beta=0$$ By comparing the coefficients, we have $\alpha+\beta=-b$ and $\alpha\beta=c$. Thus, $\alpha$ and $\beta$ can be written as $u-\frac{b}{2}$ and $-u-\frac{b}{2}$, respectively, where $u$ is an unknown quantity that is to be determined. Since $\alpha\beta=c$, we have $u^2-\frac{b^2}{4}=-c$. Solving this equation for $u$, we obtain $u=\pm\frac{\sqrt{b^2-4ac}}{2}$ and hence $\alpha$ and $\beta$ are given by $-\frac{b}{2}\pm\frac{\sqrt{b^2-4ac}}{2}$.

Example. The above question by Heron is equivalent to solving the quadratic equation $x^2+4x-896=0$. $\frac{b}{2}=2$ so we have $u^2-4=896$ i.e. $u^2=900$. Therefore, $u=\pm 30$ and the two solutions are $x=-2\pm 30$. Since $x>0$, the answer is $x=28$.

The method still applies to quadratic equations that have complex solutions. For example,

Example. Consider $x^2+x+1=0$. $\frac{b}{2}=\frac{1}{2}$, so we have the equation $u^2-\frac{1}{4}=-1$ i.e. $u^2=-\frac{3}{4}$. Therefore, $u=\pm\frac{\sqrt{3}i}{2}$ and the solutions are $x=-\frac{1}{2}\pm\frac{\sqrt{3}i}{2}$.

Shannon Entropy

Leave a reply

Let $X$ be a discrete random variable and suppose it takes values $x_1,\cdots,x_n$. Let $p$ be the probability distribution function (pdf in short) of the random variable $X$. For each $i=1,\cdots,n$, denote by $p_i$ $p(x_i)$, i.e.
$$p_i=P\{X=x_i\}$$
Here is a question. Is there a unique way to measure the amount or degree of uncertainty represented by the probability distribution? The answer is affirmative and Claude E. Shannon, in his landmark paper [2], was able to obtain such measurement of uncertainty. It is denoted by $H(p_1,\cdots,p_n)$ and called the entropy. Regardless of the name, there is no reason to expect that this notion of entropy is related to thermodynamic entropy because it is obtained independently from any physical hypotheses. Remarkably however, it coincides with thermodynamic entropy except for the presence of Boltzmann constant. To distinguish it from thermodynamic entropy, it is often called the information entropy or Shannon entropy. Shannon imposed the following three conditions on $H$:

  1. $H$ is a continuous function of the $p_i$.
  2. If all $p_i$ are equal, the quantity $$A(n)=H\left(\frac{1}{n},\cdots,\frac{1}{n}\right)$$ is a monotone increasing function of $n$. An increase in $n$ results in a decrease in the probability, therefore the uncertainty (entropy) increases.
  3. Instead of considering of the events $x_1,\cdots,x_n$ with the respective probabilities $p_1,\cdots,p_n$, we may group the first $k$ of them $x_,\cdots,x_k$ as a single composite event and is given the probability $w_1=p_1+\cdots+p_k$. Then next $m$ events $x_{k+1},\cdots,x_{k+m}$ are gouped together as a single composite event and is given the probability $w_2=p_{k+1}+\cdots+p_{k+m}$, and so on so forth. The entropy regarding these composite events is $H(w_1,\cdots,w_r)$. In addition, since the composite events have occurred, the events $x_1,\cdots,x_k$ should be given the conditional probabilities $\frac{p_1}{w_1},\cdots,\frac{p_k}{w_1}$, the events $x_{k+1},\cdots,x_{k+m}$ should be given the conditional probabilities $\frac{p_{k+1}}{w_2},\cdots,\frac{p_{k+m}}{w_2}$, and so on so forth. In order for our information measurement is to be consistent, the ultimate state of knowledge we obtained this way should be the same as one obtained by considering the event $x_1,\cdots,x_n$ with probabilities $p_1,\cdots,p_n$. The argument is summed up as \begin{equation}\begin{aligned}H(p_1,\cdots,p_n)=&H(w_1,\cdots,w_r)+w_1H\left(\frac{p_1}{w_1},\cdots,\frac{p_k}{w_1}\right)\\&+w_2H\left(\frac{p_{k+1}}{w_2},\cdots,\frac{p_{k+m}}{w_2}\right)+\cdots\end{aligned}\label{eq:composition}\end{equation}The appearance of the weighting factors $w_i$ is of course due to the fact that, for example, we run into the additional uncertainty $H\left(\frac{p_1}{w_1},\cdots,\frac{p_k}{w_1}\right)$ only with the probability $w_1$. \eqref{eq:composition} is called the composition law.

Example. $H\left(\frac{1}{2},\frac{1}{3},\frac{1}{6}\right)=H\left(\frac{1}{2},\frac{1}{2}\right)+\frac{1}{2}H\left(\frac{2}{3},\frac{1}{3}\right)$ by grouping the events with probabilities $\frac{1}{3}$ and $\frac{1}{6}$ as a composite event.

Example. Let us consider three events $x_1,x_2,x_3$ with respective probabilities $p_1=\frac{3}{9}$, $p_2=\frac{4}{9}$, $p_3=\frac{2}{9}$. We may view these three events as three composite events resulted from nine events with equal probability $\frac{1}{9}$ by grouping the first 3 events, then the next 4 events, and finally the remaining 2 events together. Using \eqref{eq:composition}, we obtain
\begin{equation}
\label{eq:entropy}
H\left(\frac{3}{9},\frac{4}{9},\frac{2}{9}\right)+\frac{3}{9}A(3)+\frac{4}{9}A(4)+\frac{2}{9}A(2)=A(9)
\end{equation}

In exactly the same manner, if we write $p_i=\frac{n_i}{\sum_in_i}$, \eqref{eq:entropy} is readily generalized to
\begin{equation}
\label{eq:entropy2}
H(p_1,\cdots,p_n)+\sum_ip_iA(n_i)=A\left(\sum_in_i\right)
\end{equation}
As a special case, if we choose $n_i=m$, then \eqref{eq:entropy2} becomes
\begin{equation}
\label{eq:entropy3}
A(m)+A(n)=A(mn)
\end{equation}
Obviously,
\begin{equation}
\label{eq:entropy4}
A(n)=K\ln n
\end{equation}
(where $K>0$ due to condition 2) is a solution. But is that a unique solution? To see that, let us consider a continuum extension $A(x)$ of $A(n)$ which is at least twice differentiable. Then \eqref{eq:entropy3} can be written as the functional equation
\begin{equation}
\label{eq:entropy5}
A(x)+A(y)=A(xy)
\end{equation}
Differentiating \eqref{eq:entropy5} with respect to $x$, we obtain
\begin{equation}
\label{eq:entropy6}
A'(x)=yA'(xy)
\end{equation}
Differentiating \eqref{eq:entropy6} with respect to $y$, we obtain
\begin{equation}
\label{eq:entropy7}
0=A'(xy)+xyA^{\prime\prime}(xy)
\end{equation}
By introducing a new variable $z=xy$, \eqref{eq:entropy7} can be written as the second-order linear differential equation
$$zA^{\prime\prime}(z)+A'(z)=0$$
whose solution is
$$A(z)=K\ln(Cz)$$
That is, we have $A(n)=K\ln(Cn)$ where $K>0$ and $C>0$ are constants. Hence, $A(n)$ is unique up to constant multiples. Since entropy is a qualitative measurement of uncertainty rather than a quantitative one, we may choose $C=1$. For the same reason, we may choose $K=1$ as well but for now we will keep the constant multiple $K$. Substituting \eqref{eq:entropy4} into \eqref{eq:entropy2}, we finally obtain Shannon entropy
\begin{equation}
\label{eq:shannon}
\begin{aligned}
H(p_1,\cdots,p_n)&=K\ln\left(\sum_in_i\right)-K\sum_ip_i\ln n_i\\
&=-K\sum_ip_i\ln\left(\frac{n_i}{\sum_in_i}\right)\\
&=-K\sum_ip_i\ln p_i
\end{aligned}
\end{equation}

Definition. There is a generalization of Shannon entropy called the Rényi entropy
\begin{equation}
\label{eq:renyi}
H_r(p_1,\cdots,p_n)=\frac{1}{1-r}\ln\left(\sum_ip_i^r\right)
\end{equation}
for $0<r<\infty$ and $r\ne 1$.

Proposition 1. $$\lim_{r\to 1}H_r(p_1,\cdots,p_n)=-\sum_ip_i\ln p_i$$
i.e., Shannon entropy (with $K=1$) is the limit of the Rényi entropy as $r\to 1$.

Proof. Let $f(r)=\ln\left(\sum_ip_i^r\right)$. Then it is readily seen that the limit on the left is simply $f'(1)$.

Definition. Tsallis entropy is another generalization of Shannon entropy. It is defined by
\begin{equation}
\label{eq:tsallis}
S_r(p_1,\cdots,p_n)=\frac{1}{r-1}\left(1-\sum_ip_i^r\right)
\end{equation}

One can easily show that the limit of Tsallis entropy as $r\to 1$ is Shannon entropy (with $K=1$) similarly to the proof of proposition 1.

References:

  1. E. T. Jaynes, Information Theory and Statistical Mechanics, The Physical Review, Vol. 106, No. 4, 620-630, May 15, 1957
  2. Claude E. Shannon, A mathematical theory of communication, Bell Sys. Tech. Journal, 27: 379-423, 623-656, 1948

Computing for Number Theory 1: GCD, LCD, Linear Diophantine Equations

Leave a reply

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