How to calculate logrithm
Simple And Interesting Numerical Method in TextBook.
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.
Another Algorithm in Excersise.
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\).
- [Initialize.] Set \(y \gets 0, z \gets x\ shifted\ right\ 1, k \gets 1\).
- [Test for end.] if \(x=1\), stop.
- [Compare.] if \(x - z < 1\), set \(z \gets z\ shifted\ right\ 1, k \gets k + 1\), and repeat this step.
- [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\).