How to calculate logrithm

I haven't talked about math for a long long time.

But this week, I came across some materials in TAOCP. which it provides a numerical way to calculate logrithm.

Long story short, We can have a numerical way to calculate \(log_{10}x\) for any real number \(x\).

If we regard the decimal part of the result in the binary system, so \(b_k\) is in \({0, 1}\), \(k\) is any positive integer. We have:

\begin{equation} % arbitrary environments, \label{orga9ce4a1} log_{10}x = n + \frac{b_1}{2} + \frac{b_2}{2^2} + ... + \frac{b_k}{2^k} + ... \end{equation}

Where \(n\) is the integer part of the result.

Consider all number \(t\) larger than 10, we have:

\begin{equation} \label{org9631e53} log_{10}\frac{t}{10} = log_{10}t - 1 \end{equation}

Such every integer number \(t\) greater than 10 can be inducted to number between \([1, 10)\).

So, we can only consider \(x\) in \([1, 10)\), back to \eqref{orga9ce4a1}. We can conclude that \(log_{10}x < 1\), so \(n\) can be ignored.

\begin{equation} \label{orgb380669} log_{10}x = \frac{b_1}{2} + \frac{b_2}{2^2} + ... + \frac{b_k}{2^k} + ... \end{equation}

So let's consider \(log_{10}x^2\).

\begin{equation} \label{org37ba08f} log_{10}x^2 = 2log_{10}x = b_1 + \frac{b_2}{2} + ... + \frac{b_k}{2^{k-1}} + ... \end{equation}

There are two possibilities:

• \(x^2 < 10\)
• \(x^2 >= 10\)

If \(1 \leq x^2 < 10\), then \(0 \leq log_{10}x^2 < 1\), so \(b_1\) must be zero.

If \(x^2 \geqslant 10\), \(\frac{x^2}{10} < 10\) must be true. We can caculate \(log_{10}\frac{x^2}{10}\), combining \eqref{org37ba08f} and \eqref{org9631e53} we got:

\begin{equation} \label{org750681c} log_{10}\frac{x^2}{10} = 2log_{10}x - 1 = b_1 + \frac{b_2}{2} + ... + \frac{b_kac}{2^{k-1}} + ... - 1 \end{equation}

Because

\begin{equation} \label{org9c5c228} \frac{b_2}{2} + ... + \frac{b_k}{2^{k-1}}+... \leqslant \frac{1}{2} + ... \frac{1}{2^k} + ... = 1-(\frac{1}{2})^k = 1 \end{equation}

Only \(b_1 = 1\) could \eqref{org750681c} be in \([0, 1)\).

Now, we have this:

for all integer \(x \in [0, 10)\)

\begin{equation} \label{org604c6f4} \begin{split} b_1 &= 1\ if\ x^2 \geqslant 10 \\ &= 0\ if\ x^2 \le 10 \end{split} \end{equation}

\eqref{org9c5c228} actually provides a way to calculate \(b_1\) for any \(x \in [0, 10)\).

Whenever we have \(b_1\), we can iteratively square \(x\) again and again to get \(b_k\) in \eqref{org37ba08f}

\begin{equation} \label{orge002483} log_{10}\frac{x^{2^{k-1}}}{10^{2^k(\frac{b_1}{2}+....\frac{b_k}{2^k}}} = b_k + \frac{b_{k+1}}{2}... - 1 \end{equation}

Because \(1 \leq \frac{x^{2^{k-1}}}{10^{2^k(\frac{b_1}{2}+....\frac{b_k}{2^k}}} < 10\)

If and only if \(\frac{x^{2^{k-1}}}{10^{2^k(\frac{b_1}{2}+....\frac{b_k}{2^k}}} \geqslant 10\) then \(b_k = 1\).

We have \eqref{org9c5c228} And we also have \eqref{org9631e53}. So for any \(1 \leqslant \frac{x}{10^n} < 10\)

\begin{equation} \label{org93ecbac} log_{10}{x} = log_{10}\frac{x}{10^n} + n \end{equation}

And finally for \(x\) is real number, we can calculate \(b_k\) reduce to \(x \in [1, 10)\) for \eqref{orgb380669} and finally apply to \eqref{orga9ce4a1} as well.

if we'd like to proceed to precision that after digits to \(k\). we can just calculate it into \(k\) and drop the later.

\begin{equation} \label{org6be8af9} log_{10}x_k = \frac{b_1}{2} + \frac{b_2}{2^2} + ... + \frac{b_k}{2^k} \end{equation}

That's why the simple numerical logrithm caculation algrithm works.

In exercise 25, Knuth present a differnt approach dating back in essence to Henry Briggs.

Suppose a real number \(x\), \(1 \leqslant x < 2\).

To calculate \(y=log_bx\).

1. [Initialize.] Set \(y \gets 0, z \gets x\ shifted\ right\ 1, k \gets 1\).
2. [Test for end.] if \(x=1\), stop.
3. [Compare.] if \(x - z < 1\), set \(z \gets z\ shifted\ right\ 1, k \gets k + 1\), and repeat this step.
4. [Reduce values.] Set \(x \gets x - z, z \gets x\ shifted\ right\ k, y \gets y + log_b(\frac{2^k}{2^k-1})\), and go to 2.

How this could be work?

Luckily, I found a good explanation on the internet.

If we express \(x\) as product \(x = \Pi a_i\), where \(1 < a_{i+1} \leq a_i\), then

\begin{equation} log_{b}x = \sum log_b a_i \end{equation}

if \(a_i\) decrease fast enough, then we can get a good approximation by taking partial sums, assumming that we have precomputed \(log_b a_i\)

Knuth's algorithm choose this

\begin{equation} \label{orgff224d8} log_{b}x = \sum log_b \frac{2^{k_i}}{2^{k_i}-1} \end{equation}

In 3, this algorithm finds the smallest \(k\) where \(k \geqslant 1\) such that \(x - z = x(1-2^{-k}) > 1\), which is equivalent to \(x > (1-2^{-k})^{-1} = \frac{2^k}{2^k-1} = a_i\).

Because \(\frac{2^k}{2^k-1}\) decrease when \(k\) increase when \(k \geqslant 1\). So every time in 4, we found the largest \(\frac{2^k}{2^k-1}\) that below \(x\). This way, \(a_i\) can decrease fast.

We then replace \(x\) with \(x/a_i\) and continue recursively, knowing that \(log_b x = log_b a_i + log_b x / a_i\).

As a result, \(a_i\) decrease fast to the precision untouchable, finally, \(x\) becomes \(1\).