## Logarithms (Part I?)

The traditional—but certainly not the best—way to compute the value of the logarithm of some number $x$ is to use a Taylor series, for example

$\displaystyle \ln x = (x-1)-\frac{1}{2}(x-1)^2+\frac{1}{3}(x-1)^3-\frac{1}{4}(x-1)^4+\cdots$

but that expansion is only valid for $0, or so, because it is the Taylor expansion of $\ln x$ "around 1", and the convergence radius of this particular expression isn't very large. Furthermore, it needs a great deal of terms before converging.

To accelerate convergence, Young and Gregory [1] propose to “compress” the series (by a short series of identities and manipulations):

$\displaystyle \ln x = 2\left( \frac{x-1}{x+1}+\frac{1}{3}\left(\frac{x-1}{x+1}\right)^3+\frac{1}{5}\left(\frac{x-1}{x+1}\right)^5+\cdots\right)$,

or

$\displaystyle \ln x = 2 \sum_{i=0}^{\infty} \frac{1}{2i+1}\left(\frac{x-1}{x+1}\right)^{2i+1}$.

This series, where the “around $a$” condition has been removed, should converge for any $x$. However, it does not so very fast. So the problem seems to be to have a good initial guess on $\ln x$ and use a Taylor series expansion around this good guess. Let’s revisit the first formula, valid for $x\approx 1$, but this time setting $x\approx a$. We get

$\displaystyle \ln x = \ln a + \frac{x-a}{a} -\frac{1}{2}\left(\frac{x-a}{a}\right)^2+\frac{1}{3}\left(\frac{x-a}{a}\right)^3-\cdots$,

or

$\displaystyle \ln x = \ln a+ \sum_{i=1}^{\infty} \frac{(-1)^{i+1}}{i}\left(\frac{x-a}{a}\right)^i$.

So, now we’re stuck with $\ln a$. The constant $a$ clearly depends on $x$ (since we have that $x\approx{a}$), and so can’t really be a constant. We must find a convenient $a$ from $x$. Finding $\ln z$ is easy when $z$ is of some special form, namely $z=e^w$. It would be especially easy to find if $z=2^k$ (because we could use some compiler built-in, say __builtin_clz, to find $k$), but we want $\ln x$, not $\log_2 x$. Well, no matter, we observe that since

$\ln x= \ln 2 \log_2 x$,

we can write

$\displaystyle \log_2 x = \log_2 a+ \frac{1}{\ln 2}\sum_{i=1}^{\infty} \frac{(-1)^{i+1}}{i}\left(\frac{x-a}{a}\right)^i$,

where $\ln 2=0.693147\ldots$, but since we choose $a=2^k$, we get

$\displaystyle \log_2 x = k + \frac{1}{\ln 2}\sum_{i=1}^{\infty} \frac{(-1)^{i+1}}{i}\left(\frac{x-2^k}{2^k}\right)^i$.

Lastly, since $x \approx 2^k$, $k$ can be approximated as the number of bits in $x$, something like $k=\lfloor log_2 x\rfloor$. This operation, of course, is to be carried not in floating point (because if we’re doing series, it’s likely because we do not have a log function in the first place) but with an efficient built-in, something like __builtin_clz.

*
* *

Let us now compare how fast each of the variants converge.

The series “around $a$” with a very good (and deusexmechanically chosen) guess for $\ln a$ converges very fast ($\ln 23$ was chosen). The second, in green, is the series with $k=[\log_2 x]$ (where the brackets stand for rounding). Lastly, Young’s and Gregory’s “compressed” series. It does converge, but it requires a lot of iterations.

Numerically, we also see rather rapid convergence: (with the columns in the same order: good initial guess, $2^k$ guess, and Young’s and Gregory’s):

*
* *

Lastly, we can ask ourselves how bad the initial guess can be given the number of iterations we are willing to perform. For now, I do not have a good characterization of it. But if we look at the following graph:

we find that, even as the number of terms grows, the convergence radius (which we should define more formally but that could be understood to be the region where the initial guess still leads to an approximation with a small error) doesn’t grow exponentially fast. More on this later, maybe.

*
* *

Choosing an initial guess using base 2 logarithms makes only sense if we have a very efficient way of finding integer $k$ in $x=2^k+b$ (where $b$ is not necessarily an integer. For this, you may use, as I said earlier, a compiler-specific extension, such as GCC’s __builtin_clz, that we used before, or something like it that maps to an efficient machine instruction. I also am not sure that Taylor series, and variants, are the best way of computing logarithms. However, except maybe for CORDIC algorithms, pretty much every textbook on numerical mathematics approaches the problem as we did here, except that they do not choose $a$ in “around $a$” adaptively.

[1] David M. Young, Robert Todd Gregory — A Survey of Numerical Mathematics, Vol 1 — Addison-Wesley (1972)