What is the fastest algorithm for finding the natural logarithm of a big number?

Like for $$\pi$$, we have an algorithm/infinite series that can give us the first 50 decimal places in about 3 terms. So if I wasn't to calculate like $$ln(25551879...)$$ (a really huge integer, most likely a prime), upto 100 decimal places, what will be the algorithm I should use or is used worldwide and how efficient is it? I know that the Taylor series is rather slow in its work, so any other algorithm in which this is computed?

• The discrete logarithm tag relates to a problem of solving equations modulo an integer, not to computing the natural logarithm. – hardmath Oct 5 at 16:30
• The Taylor series is very fast when the argument is close to $1$. Can you use a logarithm identity to make the argument remaining in the log be close to $1$? – Eric Towers Oct 5 at 23:12
• 100 decimal places isn't very big in arbitrary-precision computation. Relevant reference: Fredrik Johansson, "Efficient implementation of elementary functions in the medium-precision range." In 22nd IEEE Symposium on Computer Arithmetic (ARITH 22), pp. 83-89. IEEE, 2015. ArXiv. Abstract We describe a new implementation of the elementary transcendental functions exp, sin, cos, log and atan for variable precision up to approximately 4096 bits. Compared to the MPFR library, we achieve a maximum speedup ranging from a factor 3 for cos to30 for atan. – njuffa Oct 6 at 14:49
• 'big number' is not well-defined, please give us some ballpark bounds. $e^{230} \approx 10^{100}$, do you mean more or less than that? And do you only want for integers, or any arbitrary float? – smci Oct 7 at 2:28

Say you need an absolute tolerance of $$2^{-m}$$ for the answer.

Given a number of the form $$x=a \cdot 2^n$$, $$a \in (1/2,1]$$, write $$\ln(x)=\ln(a)+n\ln(2)$$.

Now compute $$\ln(a)$$ by taking $$m$$ terms of the Maclaurin series of $$\ln(1+x)$$ with $$x=a-1$$, and compute $$\ln(2)$$ as $$-\ln(1/2)$$ by taking $$m \lceil \log_2(n) \rceil$$ terms of the Maclaurin series of $$\ln(1+x)$$ with $$x=-1/2$$.

This way is a little bit fussy in terms of working with decimal numbers vs. binary numbers, but it has the advantage that the $$\ln(a)$$ term converges at worst like $$2^{-m}$$ rather than like $$(9/10)^m$$ like the analogous approach with decimal does. It has the downside that you have to precompute $$\ln(2)$$ to better accuracy since $$n$$ will be larger, but that doesn't matter that much because it's not a "live" problem (provided you enforce some cap on the size of the input and the size of its reciprocal).

This is generally not how people implement library functions in programming languages like C. See, for example, e_log.c at http://www.netlib.org/fdlibm/. This begins with an argument reduction similar to the one I suggested above (where the lower bound and the upper bound for $$a$$ differ by a factor of $$2$$), then converts the problem to $$\ln(1+x)=\ln(1+y)-\ln(1-y)$$ where $$y=\frac{x}{2+x}$$. This conversion leads to some series acceleration, since the series for the difference has only odd powers, and since $$|y|<|x|$$. (You could proceed with a Taylor series approach from here. If you did, it would use around $$m/4$$ terms, due to the aforementioned cancellations and the fact that $$y$$ is in the ballpark of $$x/2$$, taking into account that argument reduction has already been done.)

They then use a minimax polynomial to approximate $$\frac{\ln(1+y)-\ln(1-y)-2y}{y}$$. This kind of approach is what I usually see when I check source code for fast implementations of standard library functions. The coefficients of this minimax polynomial were probably relatively expensive to calculate, but again that's not a "live" problem so its speed doesn't really matter all that much.

• +1 Covering all the bases (not literally) :P – Simply Beautiful Art Oct 6 at 18:44
• To add on a little to how you could construct a slightly faster converging series for $\ln(a)$ and $\ln(2)$, rewrite $$\ln(x)=\int_1^x\frac{\mathrm dt}t$$ and Taylor expand $\frac1t$ at $\frac{x+1}2$. The result in the case of $x=2$ is $$\ln(2)=\frac23\sum_{k=0}^\infty\frac1{(2k+1)9^k}$$which converges much faster than your provided $-\ln(1/2)$ as some intuition on the last paragraph perhaps. – Simply Beautiful Art Oct 6 at 18:55
• @SimplyBeautifulArt The last paragraph does seem to be along these same lines since the idea is to choose an interval centered at $1$ (i.e. centered at the point of expansion) such that integrating $\frac{1}{t}$ there gives the same as integrating over $[1,1+x]$. (Also your answer is a lot better than mine overall; the only advantage my answer has over yours is the comment about netlib and the use of minimax polynomials.) – Ian Oct 7 at 17:45

This is essentially a more in depth discussion on the efficiency and accuracy of various methods.

Reducing the argument so that it is near $$1$$:

Fundamentally, most of the answers aim for the same goal: reduce the arguments to small values and use the Taylor expansion for $$\ln(x)$$. Thus far we've seen 3 approaches:

1) Factor out a power of 2 and use $$\ln(a\cdot2^n)=\ln(a)+n\ln(2)$$.

2) Factor out a power of 10 and use $$\ln(a\cdot10^n)=\ln(a)+n\ln(10)$$.

3) Reduce by square rooting using $$\ln(x)=2\ln(\sqrt x)$$.

We can note that square rooting the argument repeatedly reduces the argument much faster than the other methods, which divide the argument by a constant repeatedly, since $$\sqrt x\ll x/10 for large $$x$$. Realistically, if your input doesn't have more than, say, 1000 digits, then you only have to square root about 10 times at worst. However, this comes at the cost of having to perform more computations to find the square root itself, which are not easy. On the other hand though, performing the divisions is extremely easy. Due to the nature of how we write/store numbers, all the divisions can be computed at once by simply moving the decimal point. We can then easily truncate to whatever accuracy is desired. With the square rooting, the error is harder to manage, and since the log gets multiplied by 2 every time, the absolute error gets multiplied by 2 as well.

Of course, the choice of writing the argument as a multiple of a power of 2 or a power of 10 depends on whether or not this is being done by a computer or by a human. You will likely prefer to work in base 10.

There is also the additional question of what our range of $$a$$ should be. Since we want this to be as close to 1 as possible, we can do some algebra and see. For powers of 2, we want $$a\in(a_0,2a_0]$$ such that $$2a_0-1=1-a_0$$. Solving this gives $$a\in[\frac23,\frac43]$$. For powers of 10, we want $$a\in[\frac2{11},\frac{20}{11}]$$.

Using series expansions:

From here we could use the standard Taylor expansion for the natural logarithm:

$$\ln(a)=\sum_{k=1}^\infty\frac{(-1)^{k+1}}k(a-1)^k=(a-1)-\frac{(a-1)^2}2+\frac{(a-1)^3}3-\frac{(a-1)^4}4+\mathcal O(a-1)^5$$

however this does not converge as fast as one could manage by performing a Taylor expansion closer to $$a$$. The above is given by using the integral definition of the natural logarithm and Taylor expanding the integrand at $$1$$:

$$\ln(a)=\int_1^a\frac{\mathrm dt}t=\sum_{k=0}^\infty(-1)^k\int_1^a(t-1)^k~\mathrm dt$$

But we can improve on this by expanding in the middle of $$1$$ and $$a$$:

\begin{align}\ln(a)=\int_1^a\frac{\mathrm dt}t&=\sum_{k=0}^\infty(-1)^k\left(\frac2{a+1}\right)^{k+1}\int_1^a\left(t-\frac{a+1}2\right)^k~\mathrm dt\\&=\sum_{k=0}^\infty\frac{(-1)^k}{k+1}\left(\frac{a-1}{a+1}\right)^{k+1}\left(1-(-1)^{k+1}\right)\\&=\sum_{k=0}^\infty\frac2{2k+1}\left(\frac{a-1}{a+1}\right)^{2k+1}\end{align}

For $$a$$ near $$1$$, the above has a rough error of $$\mathcal O((a-1)/2)^{2n+1}$$ when using $$n$$ terms. An algebraic derivation of the above is provided by Wikipedia but doesn't really show off just how fast this one converges. Since we're twice as close to the farthest bound on the integral, we gain an additional binary digit of accuracy per term. But since half of the terms vanish, this means we can basically compute twice as many digits per term! This is the method mentioned by Ian's answer.

Here is a rough Ruby program computing the logarithm using series.

Using root-finding methods:

While the series methods are really nice and converge decently fast, Wikipedia provides two more methods for even higher precision evaluation. The first is provided by Eric Towers and it involves computation of the logarithm via exponential functions. Of course, since the computation is cheap and the accuracy is not effected so much, I would recommend reducing the argument so that it is once again close to $$1$$. This will mean that $$y$$ as defined below will be close to $$0$$, allowing for faster computation of the exponential. This also means we can simply use $$y_0=0$$ as our initial guess.

$$y=\ln(x)\Rightarrow x=\exp(y)\Rightarrow x-\exp(y)=0$$

At which we may apply standard root-finding methods, such as Newton's method (doubles accurate digits per step) or Halley's method (triples accurate digits per step).

The calculation of the exponential functions can be performed using the Maclaurin expansion:

$$\operatorname{exmp1}(y)=\exp(y)-1=\sum_{n=1}^\infty\frac{y^n}{n!}=y+\frac{y^2}2+\mathcal O(y^3)$$

Since $$y$$ is near $$0$$, there is large floating point error in computing $$\exp(y)$$, which has a dominant term of $$1$$, so we use $$\operatorname{expm1}(y)$$ to circumvent this.

One may also note that since $$\Delta y_n\to0$$, it is easier to compute $$\exp(\Delta y_n)$$ than $$\exp(y_{n+1})$$ directly, and use $$\exp(y_{n+1})=\exp(\Delta y_n)\exp(y_n)=\exp(y_n)+\exp(y_n)\operatorname{expm1}(\Delta y_n)$$. This reduces the main exponentiation down to the first step, which is trivial since $$\exp(0)=1$$.

Let $$y_0=0$$ and $$\operatorname{expy}_0=1$$.

For Newton's method, let $$\displaystyle\Delta y_0=x\operatorname{expy}_0-1$$ and:

$$\begin{cases}\Delta y_n=x\operatorname{expy}_n-1,\\\operatorname{expy}_{n+1}=\operatorname{expy}_n+\operatorname{expy}_n\operatorname{expm1}(-\Delta y_n),\\y_{n+1}=y_n+\Delta y_n\end{cases}$$

For Halley's method, let $$\displaystyle\Delta y_0=2\cdot\frac{x-\operatorname{expy}_0}{x+\operatorname{expy}_0}$$ and:

$$\begin{cases}\displaystyle\Delta y_n=2\cdot\frac{x-\operatorname{expy}_n}{x+\operatorname{expy}_n},\\\operatorname{expy}_{n+1}=\operatorname{expy}_n+\operatorname{expy}_n\operatorname{expm1}(\Delta y_n),\\y_{n+1}=y_n+\Delta y_n\end{cases}$$

Here is a rough Ruby program computing the logarithm with Newton's method and here is a rough Ruby program computing the logarithm with Halley's method.

Using the agm:

The arithmetic-geometric mean is a powerful tool which can be used here to quickly compute the logarithm as well as $$\pi$$ and certain integrals. It is defined as:

$$a_0=a,b_0=b\\a_{n+1}=\frac{a_n+b_n}2,b_{n+1}=\sqrt{a_nb_n}\\M(a,b)=\lim_{n\to\infty}a_n=\lim_{n\to\infty}b_n$$

According to Wikipedia, this is so fast and cheap to compute that this can be used to compute the exponential function using logarithms faster than series approximating the exponential function! To achieve $$p$$ bits of accuracy, take an $$m$$ so that $$s=x2^m$$ is greater than $$2^{p/2}$$. We may then compute the natural logarithm as:

$$\ln(x)=\lim_{m\to\infty}\frac{\pi x2^{m-2}}{2M(x2^{m-2},1)}-m\ln(2)$$

which is essentially a restatement of the convergence rate of $$M(1,t)$$ as $$t\to\infty$$. For this method, reducing the argument isn't even necessary. We may in fact apply this directly to large arguments!

There are some drawbacks to this method, however. The computation requires us to compute some square roots on large floats, but this can be handled with a specially defined float classes and the respective functions.

Alternatively, of course, one could simply reduce the argument down to avoid large floats like before.

Here is a rough Ruby program computing the logarithm using the arithmetic-geometric mean.

I don't know what the fastest way is, but here's one reasonably efficient approach:

• You can divide numbers with Newton-Raphson;
• Once you know how to do that, you can also take square roots with Newton-Raphson;
• You can use $$\ln x=2\ln\sqrt{x}$$ as often as you need to get the logarithm's argument close to $$1$$ before using the Maclaurin series of $$\ln(1+x)$$;
• If you need logarithms in another base, use $$\log_ax=\frac{\ln x}{\ln a}$$.

Halley's method is iterative and its convergence is cubic. Applied here, we would invert to use the exponential (which happily is its own derivative): $$y_{n+1} = y_n+2 \cdot \frac{x- \mathrm{e}^{y_n}}{x + \mathrm{e}^{y_n}} \text{.}$$ For instance, with $$x = 25551879$$ and $$y_0 = 2$$ (i.e., not close), the iterates (all computed with 15-ish digit intermediates) are $$2$$, $$4.$$, $$5.99999$$, $$7.99993$$, $$9.99946$$, $$11.996$$, $$13.9708$$, $$15.7959$$, $$16.9122$$, $$17.056$$, $$17.0562$$. My point is not that 15 digits is enough, but that the method reached the shown precision in only ten steps.

You might ask, how do I get those exponentials? Use the power series. That converges quickly for any argument you are likely to see. As a rule of thumb, start with twice as many terms as your argument, so for $$\mathrm{e}^{17.0562}$$, start with $$34$$ terms in this Taylor series. This gives $$2.5549{\dots}\times 10^{7}$$ with error $$2355.61{\dots}$$. Then increase the number of terms in the exponentials by (in this case) $$34$$ as long as your estimate for $$y_{n+1}$$ still changes within your target precision. When that stops happening, take that as your final $$y_{n+1}$$ and repeat the process of extending the exponential series until your $$y_{n+2}$$ stabilizes. Continue until you get two $$y$$s in a row that agree to your target precision (plus enough extra unchanging bits that at least one of them is a zero so that you know which way to round the last bit in your reported answer).

• For an even faster $\exp$, one can typically use $s = \mathrm{sinh}(x); \exp(x) =s+\sqrt{1+s^{2}}$. Cuts the number of terms in the Taylor expansion in half. – njuffa Oct 7 at 2:21
• CORDIC is also a candidate. – smci Oct 7 at 2:24
• @njuffa According to Wikipedia, one would usually use $\exp(x)=1+\frac{2\tanh(x/2)}{1-\tanh(x/2)}$ instead. – Simply Beautiful Art Oct 7 at 22:55
• Thinking about what @smci said, we can note that $\Delta y_n\to0$ very fast, and hence $\exp(y_{n+1})=\exp(\Delta y_n)\exp(y_n)$ can be computed quickly and accurately. However this also means computing each to more than the desired accuracy, since we lose a bit of accuracy each time we multiply. – Simply Beautiful Art Oct 7 at 23:32
• @njuffa Note though that the argument of mine is halved whereas yours is not. Mine also does not involve squaring or square roots. – Simply Beautiful Art Oct 8 at 1:22

Well$$\ln(25551879) = \ln(0.25551879 \times 10^{8})$$ $$= \ln(0.25551879) + \ln(10^8)$$ $$= 8 \times \ln(10) + \ln(0.25551879)$$

Since $$\ln(10)$$ is a constant that can be precomputed to a huge number of decimal places we only need a method that converges quickly for values less than $$1.0$$. I don't know if Taylor is good enough on that restricted range or if there is some other better method.

This addresses the issue you raised about the argument being a large number. As to generating lots of digits, there are a lot of good answers on this previous question.

• Yes but my question still remains, one way or another. You still need to compute $\ln(0.25551879)$ and also $\ln(10)$. And yes log of 10 is already computed, but what with what algorithm? That is my question. Any number can be broken down to smaller terms but say we need to compute the log of that big value at to like a million digits, then what? – Rayreware Oct 5 at 14:28
• Edited to address comment. – JonathanZ Oct 5 at 14:36
• I would prefer $7\ln 10 +1+\ln (1-d)$ where $1-d=2.5551879/e$, because $d$ is fairly close to $0$. – DanielWainfleet Oct 5 at 14:56
• We can tabulate precomputed constants like ln(10) to arbitary precision with whatever algorithm we need, it doesn't really matter. Here we only care about the dynamic part (the logarithm of the mantissa). – smci Oct 7 at 2:22

The Taylor series is "slow" for radius close to $$1$$ (since it's centered on $$x=1$$, radius $$1$$ corresponds to $$x=0$$ and $$x=2$$), and doesn't work for radius greater than or equal to $$1$$. Pretty much any Taylor series will be fast for small radii. The log Taylor series has decreasing coefficients and it is alternating, so its error term can be given an upper bound of $$x^n$$. So if you have a radius less than $$x^k$$, then you'll be getting at least $$k$$ digits of accuracy for every term.

$$25551879$$ is well outside the radius of $$1$$, so a direct application of the Taylor series won't work. You'll have to write it in terms of some $$x$$ that is close to $$1$$. You could write it as $$10^8 * 0.255$$, then find its logarithm as $$\ln(10^8) + \ln(0.255) = 8\ln(10)+\ln(1-0.745)$$, and then use a precomputed value for $$\ln(10)$$ and the Taylor series for $$\ln(1-0.745)$$, but $$0.745$$ is a large radius, so you can get faster convergence by writing it in another form. For instance, if you've precomputed $$\ln(2)$$, then you can write it as $$10^7*2*1.2775$$, and now your radius is only $$0.2775$$.