How Do Computer Calculate the Log-Function

    Take the Glibc's Mindset for Example
    Dec 7, 2023

    Brief Introduction to the Implementation Approach in glibc:

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

    logkx=lnxlnk(1)\log_k{x} = \frac{\ln{x}}{\ln k} \tag{1}

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

    Computation of lnx\ln x

    In advanced calculus, we learned the following Taylor series:

    ln(1+x)=x12x2+13x314x4+=k=1(1)kkxk(2)\ln(1+x) = x - \frac{1}{2}x^2 + \frac{1}{3}x^3 - \frac{1}{4}x^4 + \dots = \sum_{k=1}^{\infty}\frac{(-1)^k}{k}x^k \tag{2}

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

    ln(x)=ln(1+t)=k=1(1)kktk(2)\ln(x) = \ln(1+t) = \sum_{k=1}^{\infty}\frac{(-1)^k}{k}t^k \tag{2}

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

    ln(1+t)=t12t2+13t314t4+=k=1n(1)kktk+o(tn)(3)\ln(1+t) = t - \frac{1}{2}t^2 + \frac{1}{3}t^3 - \frac{1}{4}t^4 + \dots = \sum_{k=1}^{n}\frac{(-1)^k}{k}t^k + o(t^n) \tag{3}

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

    However, another formula is available:

    lnx=2(x1x+1+13(x1x+1)3+15(x1x+1)5+)=2k=012k+1(x1x+1)2k+1(4)\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 + \dots\right) = 2\sum_{k=0}^{\infty}\frac{1}{2k+1}\left(\frac{x-1}{x+1}\right)^{2k+1} \tag{4}

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

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

    ln(1+1x)=ln(x+1x)=k=1(1)kk1xkln(11x)=k=11kxk\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:

    ln(x1x)=k=11kxkln(x1x)=k=11kxkln(xx1)=k=11kxk\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:

    ln(x+1x)+ln(xx1)=ln(x+1x1)k=1n(1)kkxk+k=1n1kxk=2k=0n12k+11x2k+1\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:

    ln(x+1x1)=2k=0n12k+11x2k+1\ln\left(\frac{x+1}{x-1}\right) = 2\sum_{k=0}^n{\frac{1}{2k+1}\frac{1}{x^{2k+1}}}

    Let z=(x+1)(x1)z = \frac{(x+1)}{(x-1)}, then x=(z1)(z+1)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 nnth term:

    lnx=2k=0n12k+1(x1x+1)2k+1+o((x1x+1)2k+1)(5)\ln x = 2\sum_{k=0}^{n}\frac{1}{2k+1}\left(\frac{x-1}{x+1}\right)^{2k+1} + o\left(\left(\frac{x-1}{x+1}\right)^{2k+1}\right) \tag{5}

    The remainder term is o((x1)(x+1)2k+1)o\left(\frac{(x-1)}{(x+1)^{2k+1}}\right), and as limx+(x1x+1)2k+1=1\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 11.

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

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

    x=1.xxx×2expx = 1.xxx \times 2^{exp}

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

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

    ln(x2exp)=lnxexp×ln2\ln\left(\frac{x}{2^{exp}}\right) = \ln x - \text{exp}\times \ln 2

    Therefore:

    lnx=ln(x2exp)+exp×ln2\ln x = \ln\left(\frac{x}{2^{exp}}\right) + \text{exp} \times \ln 2

    This completes the calculation of lnx\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.