NOTE & BLOG

Brief Introduction to the Implementation Approach in glibc:

For the logarithm of any base, we can utilize the change of base formula:

The term $\ln k$ in equation (1) is a constant, so for logarithms of any base, we only need to calculate $\ln x$.

Computation of $\ln x$

In advanced calculus, we learned the following Taylor series:

Suppose we want to find $\ln x$; naturally, we can set $t=x-1$, then:

However, expanding the series to infinity in practical calculations is often not feasible. If we expand it to the $n$th term, the Taylor series with the Peano remainder term can be expressed as:

Note that the Peano remainder term here is $o(t^n)$, representing the error term. If $t>1$ (or $x > 2$), the error in the above formula may become significantly large with the increase of $t$ or $x$.

However, another formula is available:

Here's an informal proof for equation $(4)$:

Replace $x$ in equation $(2)$ with $\frac{1}{x}$ and $-\frac{1}{x}$, respectively:

$\begin{aligned} \ln\left(1+\frac{1}{x}\right) &= \ln\left(\frac{x+1}{x}\right) = \sum_{k=1}^{\infty}\frac{(-1)^k}{k}\frac{1}{x^k} \\ \ln\left(1-\frac{1}{x}\right) &= -\sum_{k=1}^{\infty}\frac{1}{kx^k} \end{aligned}$

Note:

$\begin{aligned} \ln\left(\frac{x-1}{x}\right) &= -\sum_{k=1}^{\infty}\frac{1}{kx^k} \\ \Longrightarrow -\ln\left(\frac{x-1}{x}\right) &= \sum_{k=1}^{\infty}\frac{1}{kx^k} \\ \Longrightarrow \ln\left(\frac{x}{x-1}\right) &= \sum_{k=1}^{\infty}\frac{1}{kx^k} \end{aligned}$

Then:

$\begin{aligned} &\ln\left(\frac{x+1}{x}\right) + \ln\left(\frac{x}{x-1}\right) = \ln\left(\frac{x+1}{x-1}\right) \\ &\sum_{k=1}^n\frac{(-1)^k}{kx^k} + \sum_{k=1}^n\frac{1}{kx^k} = 2\sum_{k=0}^n{\frac{1}{2k+1}\frac{1}{x^{2k+1}}} \end{aligned}$

This yields:

Let $z = \frac{(x+1)}{(x-1)}$, then $x = \frac{(z-1)}{(z+1)}$. Substituting this into the above formula gives equation (4).

The advantage of this formula is that if we expand it to the $n$th term:

The remainder term is $o\left(\frac{(x-1)}{(x+1)^{2k+1}}\right)$, and as $\lim_{x\rightarrow+\infty}\left(\frac{x-1}{x+1}\right)^{2k+1} = 1$, the remainder term ensures that the error remains within the range of $1$.

However, this is still not sufficient, as the error of $1$ is still too large. According to the analysis of the remainder term above, the value of $x$ should be closer to $1$ for higher precision.

Considering the IEEE 754 floating-point standard, we know that any floating-point number is represented in the computer as:

This means that we can directly obtain the exponent ($\text{exp}$) and the fractional part ($1.xxx$) using bitwise operations.

Since the floating-point number is already in this form, we can directly use the fractional part $1.xxx$. If we calculate $\ln(1.xxx)$ using the earlier formula, it is equivalent to calculating $\ln\left(\frac{x}{2^{exp}}\right)$. According to the logarithmic formula:

Therefore:

This completes the calculation of $\ln x$, equivalent to computing the logarithm of any base.

However, this is a general overview, and glibc employs various precision improvement techniques in its implementation.