Draft.

COMPUTING LOGARITHM INTERVALS

WITH THE ARITHMETIC-GEOMETRIC-MEAN ITERATION

DANIEL J. BERNSTEIN

Abstract. This paper presents a fast algorithm that, given a tight interval around a positive real number x, computes a tight interval around log x. To obtain p bits of precision for typical values of x, the algorithm uses about 2 lg p square roots and about 5 lg p multiplications (or fewer for subsequent logarithms). Here log is the natural logarithm, and lg is the base-2 logarithm. This paper also presents short proofs of all necessary properties of complete elliptic integrals.

1. Introduction

This paper presents an algorithm that, given a positive real number x to high accuracy, computes log x to high accuracy. Here log is the natural logarithm. The algorithm has several useful features:

• It is extremely fast. For p bits of precision, and for x between (e.g.) 2 and 2^{p}, the algorithm time is dominated by about 7 lg p operations on p-bit numbers, specifically 2 lg p square roots and 5 lg p multiplications. Here lg is the base-2 logarithm.

• It computes π to high precision, practically for free, as a side effect of computing log.

• It computes subsequent logarithms at even higher speed: asymptotically, per logarithm, about 2 lg p square roots and about 2 lg p multiplications.

• Given a tight interval containing x, it computes tight intervals containing log x and π. See [\underline{4}] for an application. For large p, the extra cost of handling intervals is negligible.

See Section 5 for details of the algorithm.

The algorithm relies on various properties of complete elliptic integrals and the arithmetic-geometric-mean (AGM) iteration. This paper presents streamlined proofs of those properties. Section 2 introduces the elliptic integrals I and I_{1} used throughout the paper; Section 3 relates I and I_{1} to logarithms; Section 4 introduces the AGM iteration.

Previous work. Salamin pointed out more than thirty years ago that the elliptic integral I could be used to quickly compute high-precision approximations to log x and π; see [\underline{2}, Item 143]. Salamin’s algorithm is dominated by Θ(lg p) high-precision operations, with larger constant factors than the algorithm here.

Permanent document ID: 8f92b1e3ec7918d37b28b9efcee5e97f.

Date: 20030717.

2000 Mathematics Subject Classification. Primary 65D20; Secondary 33E05, 65G30, 11Y60.

The author was supported by the National Science Foundation under grant DMS–0140542, and by the Alfred P. Sloan Foundation.

Subsequent refinements of Salamin’s algorithm included almost all of the ideas necessary for minimizing the constant factors, but the ideas were never combined properly. See Section 6 for a survey of the literature.

Brent in [\underline{9}, Section 6] introduced another method of computing high-precision exponentials and logarithms, without relying on elliptic integrals. See [\underline{3}, Sections 12–16] for further discussion. Brent’s method may be faster than the algorithm here for small p, but for large p it is slower by a factor roughly proportional to lg p.

Logarithm algorithms are occasionally presented with explicit bounds suitable for interval computation. See [\underline{14}] and [\underline{13}].

The results of Sections 2, 3, and 4 (modulo language, and modulo the difference between explicit inequalities and order-of-magnitude estimates) have been known for two centuries, and are a small fraction of the wealth of material collected in [\underline{8}]. However, I have not seen such short proofs of the I_{1} properties, particularly Theorem 2.6 and Theorem 4.3.

2. Elliptic integrals: basic properties

For positive real numbers a, b define I(a, b) = ∫_{0}^{∞}(x^{2} + a^{2})^{−1/2}(x^{2} + b^{2})^{−1/2} dx. Also define I_{1}(a, b) = ∂I(a, b)/∂a = ∫_{0}^{∞} −a(x^{2} + a^{2})^{−3/2}(x^{2} + b^{2})^{−1/2} dx.

Integrability, differentiation under the integral sign, etc. follow from the fact that all the integrands in this section are constant in sign, continuous in x, and differentiable in a, b.

Theorem 2.1. π/2a ≤ I(a, b) ≤ π/2b if b ≤ a.

Proof. x^{2} + b^{2} ≤ x^{2} + a^{2} so ∫_{0}^{∞}(x^{2} + a^{2})^{−1} dx ≤ I(a, b) ≤ ∫_{0}^{∞}(x^{2} + b^{2})^{−1} dx. □

Theorem 2.2. 0 ≤ −aI_{1}(a, b) ≤ I(a, b).

Proof. 0 ≤ a^{2} ≤ x^{2} + a^{2}. □

Theorem 2.3. I(a, b) = I(1, b/a)/a.

Proof. Substitute x = au: aI(a, b) = ∫_{0}^{∞} a^{2}(a^{2}u^{2} + a^{2})^{−1/2}(a^{2}u^{2} + b^{2})^{−1/2} du = ∫_{0}^{∞}(u^{2} + 1)^{−1/2}(u^{2} + b^{2}/a^{2})^{−1/2} du = I(1, b/a). □

Theorem 2.4. I(a, b) + aI_{1}(a, b) + b(∂I(a, b))/∂b = 0.

Proof. Apply a(∂/∂a) + b(∂/∂b) to Theorem 2.3. □

Theorem 2.5. I(a, b) = I(\frac{a+b}{2}, \sqrt{ab}).

Proof. Substitute u = (x − ab/x)/2, using (2u)^{2} + (a + b)^{2} = (x^{2} + a^{2})(x^{2} + b^{2})/x^{2} and x du = (u^{2} + ab)^{1/2} dx: I(a, b) = ∫_{−∞}^{∞}((2u)^{2} + (a + b)^{2})^{−1/2}(u^{2} + ab)^{−1/2} du = 2 ∫_{0}^{∞}(2^{2}(u^{2} + (\frac{a+b}{2})^{2}))^{−1/2}(u^{2} + ab)^{−1/2} du = I(\frac{a+b}{2}, \sqrt{ab}). □

Theorem 2.6. I(a, b) + 2aI_{1}(a, b) = I_{1}(\frac{a+b}{2}, \sqrt{ab})(a − b)/2.

Proof. Apply a(∂/∂a) − b(∂/∂b) to Theorem 2.5: aI_{1}(a, b) − b(∂I(a, b)/∂b) = I(\frac{a+b}{2}, \sqrt{ab})(a − b)/2. By Theorem 2.4, I(a, b) + aI_{1}(a, b) + b(∂I(a, b)/∂b) = 0. □

Theorem 2.7. I(1, b) = 2 ∫_{0}^{√\overline{b}}(x^{2} + 1)^{−1/2}(x^{2} + b^{2})^{−1/2} dx.

Proof. ∫_{√\overline{b}}^{∞}(x^{2} + 1)^{−1/2}(x^{2} + b^{2})^{−1/2} dx = ∫_{0}^{√\overline{b}}((\frac{b}{u})^{2} + 1)^{−1/2}((\frac{b}{u})^{2} + b^{2})^{−1/2} \frac{b}{u^{2}} du =

∫_{0}^{√\overline{b}}(b^{2} + u^{2})^{−1/2}(1 + u^{2})^{−1/2} du. □

Theorem 2.8. (1 + b)^{−1} + I(1, b) + I_{1}(1, b) = 2 ∫_{0}^{√\overline{b}} b^{2}(x^{2} + 1)^{−1/2}(x^{2} + b^{2})^{−3/2} dx.

Proof. By Theorem 2.4, I(1, b) + I_{1}(1, b) + b(∂I(1, b)/∂b) = 0. By Theorem 2.7,

∂I(1, b)/∂b = b^{−1/2}(b+1)^{−1/2}(b+b^{2})^{−1/2}+2 ∫_{0}^{√\overline{b}} −b(x^{2}+1)^{−1/2}(x^{2}+b^{2})^{−3/2} dx. □

Theorem 2.9. a^{2}I(a, b) + (a^{2} − b^{2})aI_{1}(a, b) = ∫_{0}^{∞} a^{2}(x^{2} + a^{2})^{−3/2}(x^{2} + b^{2})^{1/2} dx.

Proof. x^{2} + a^{2} − (a^{2} − b^{2}) = x^{2} + b^{2}. □

Theorem 2.9 explains why I(a, b), a^{2}I(a, b) + (a^{2} − b^{2})aI_{1}(a, b), etc. are called elliptic integrals, and why related algebraic objects are called elliptic curves: substitute x = a tan θ to see that a^{2}I(a, b) + (a^{2} − b^{2})aI_{1}(a, b) is the arc length of one quadrant of the ellipse θ ↦ (a cos θ, b sin θ). Theorem 2.9 is not used elsewhere in this paper.

3. Elliptic integrals: logarithm bounds

This section presents precise bounds along the following lines: “If b ≈ 0 then I(1, b) ≈ log(4/b) and I(1, b) + I_{1}(1, b) ≈ 1.” Write L(b) = log(\sqrt{b^{−1}} + \sqrt{b^{−1} + 1}).

The simple proof technique used here is not new, and explicit bounds in this context are not new, but I am not aware of previous simple proofs of explicit bounds. For inexplicit bounds using the same proof technique, see [\underline{18}, page 522] and [\underline{15}]; for explicit bounds with longer proofs, see [\underline{7}, pages 355–356] and [\underline{8}, Theorem 7.2].

Theorem 3.1. If 0 < b ≤ 1 then (2 + \frac{1}{2}b^{2})L(b) − (\frac{1}{2}b)(1 + b)^{1/2} ≤ I(1, b) ≤ (2 + \frac{1}{2}b^{2} + \frac{9}{32}b^{4})L(b) − (\frac{1}{2}b − \frac{3}{16}b^{2} + \frac{9}{32}b^{3})(1 + b)^{1/2}.

The difference between lower and upper bounds is on the scale of b^{2}. The bounds can easily be made tighter by the same technique.

Proof. I(1, b) = 2 ∫_{0}^{√\overline{b}}(1 + x^{2})^{−1/2}(x^{2} + b^{2})^{−1/2} dx by Theorem 2.7. Recall that (1+x^{2})^{−1/2} ≥ 1−x^{2}/2 for 0 ≤ x ≤ 1, since 1 ≥ 1−(3−x^{2})x^{4}/4 = (1−x^{2}/2)^{2}(1+x^{2}). Thus

I(1, b) ≥ ∫_{0}^{√\overline{b}}(2 − x^{2})(x^{2} + b^{2})^{−1/2} dx

= (2 + \frac{1}{2}b^{2}) log(x + (x^{2} + b^{2})^{1/2}) − (\frac{1}{2}x)(x^{2} + b^{2})^{1/2}|_{0}^{√\overline{b}}

= (2 + \frac{1}{2}b^{2}) log (\frac{\sqrt{b} + \sqrt{b + b^{2}}}{\sqrt{0} + \sqrt{0 + b^{2}}} ) − (\frac{1}{2} \sqrt{b})(b + b^{2})^{1/2}

= (2 + \frac{1}{2}b^{2})L(b) − (\frac{1}{2}b)(1 + b)^{1/2}.

Similarly, (1 + x^{2})^{−1/2} ≤ 1 − \frac{1}{2}x^{2} + \frac{3}{8}x^{4} for 0 ≤ x, so

I(1, b) ≤ ∫_{0}^{√\overline{b}}(2 − x^{2} + \frac{3}{4}x^{4})(x^{2} + b^{2})^{−1/2} dx

= (2 + \frac{1}{2}b^{2} + \frac{9}{32}b^{4}) log(x + (x^{2} + b^{2})^{1/2}) − (\frac{1}{2} − \frac{3}{16}x^{2} + \frac{9}{32}b^{2})x(x^{2} + b^{2})^{1/2}|_{0}^{√\overline{b}}

= (2 + \frac{1}{2}b^{2} + \frac{9}{32}b^{4})L(b) − (\frac{1}{2}b − \frac{3}{16}b^{2} + \frac{9}{32}b^{3})(1 + b)^{1/2}.
□

Theorem 3.2. If 0 < b ≤ 1 then (2+b^{2})(1+b)^{−1/2} −b^{2}L(b) ≤ (1+b)^{−1} +I(1, b)+ I_{1}(1, b) ≤ (2 + b^{2} + \frac{3}{8}b^{3} + \frac{9}{8}b^{4})(1 + b)^{−1/2} − (b^{2} + \frac{9}{8}b^{4})L(b).

Proof. (1+b)^{−1}+I(1, b)+I_{1}(1, b) = 2 ∫_{0}^{√\overline{b}} b^{2}(x^{2}+1)^{−1/2}(x^{2}+b^{2})^{−3/2} dx by Theorem 2.8. Recall that (1 + x^{2})^{−1/2} ≥ 1 − x^{2}/2 for 0 ≤ x ≤ 1, as in Theorem 3.1, so

(1 + b)^{−1} + I(1, b) + I_{1}(1, b) ≥ ∫_{0}^{√\overline{b}} b^{2}(2 − x^{2})(x^{2} + b^{2})^{−3/2} dx

= (2 + b^{2})x(x^{2} + b^{2})^{−1/2} − b^{2} log(x + (x^{2} + b^{2})^{1/2})|_{0}^{√\overline{b}}

= (2 + b^{2})(1 + b)^{−1/2} − b^{2}L(b).

Similarly, (1 + x^{2})^{−1/2} ≤ 1 − \frac{1}{2}x^{2} + \frac{3}{8}x^{4} for 0 ≤ x, so

(1 + b)^{−1} + I(1, b) + I_{1}(1, b) ≤ ∫_{0}^{√\overline{b}} b^{2}(2 − x^{2} + \frac{3}{4}x^{4})(x^{2} + b^{2})^{−3/2} dx

= (2 + b^{2} + \frac{3}{8}b^{2}x^{2} + \frac{9}{8}b^{4})x(x^{2} + b^{2})^{−1/2} − (b^{2} + \frac{9}{8}b^{4}) log(x + (x^{2} + b^{2})^{1/2})|_{0}^{√\overline{b}}

= (2 + b^{2} + \frac{3}{8}b^{3} + \frac{9}{8}b^{4})(1 + b)^{−1/2} − (b^{2} + \frac{9}{8}b^{4})L(b).
□

4. Elliptic integrals: AGM iteration

Throughout this section, a_{0} and b_{0} are positive real numbers; a_{n} and b_{n} are defined recursively by a_{n+1} = (a_{n} + b_{n})/2 and b_{n+1} = \sqrt{a_{n}b_{n}}. As n increases, both a_{n} and b_{n} converge rapidly to π/2I(a_{0}, b_{0}), the arithmetic-geometric mean of a_{0} and b_{0}; see Theorems 4.2 and 4.5.

Theorem 4.1. If n ≥ 1 then a_{n} ≥ b_{n}.

Proof. The geometric mean cannot exceed the arithmetic mean: if n ≥ 0 then a_{n+1}^{2} − b_{n+1}^{2} = (a_{n} − b_{n})^{2}/4 ≥ 0 so a_{n+1} ≥ b_{n+1}. □

Theorem 4.2. π/2a_{n} ≤ I(a_{0}, b_{0}) ≤ π/2b_{n} for n ≥ 1.

Proof. By Theorem 2.5, I(a_{0}, b_{0}) = I(a_{1}, b_{1}) = I(a_{2}, b_{2}) = ⋯ = I(a_{n}, b_{n}). By Theorem 4.1, a_{n} ≥ b_{n}. Apply Theorem 2.1. □

Theorem 4.3. Define ϵ_{n} = 2^{n}(a_{n}^{2} − b_{n}^{2})(−a_{n}(I_{1}/I)(a_{n}, b_{n})). Then 0 ≤ ϵ_{n} ≤ 2^{n}(a_{n}^{2}−b_{n}^{2}) for n ≥ 1, and (a_{0}^{2}−b_{0}^{2})(−a_{0}(I_{1}/I)(a_{0}, b_{0})) = ϵ_{n}+∑_{0≤i<n} 2^{i−1}(a_{i}^{2}−b_{i}^{2}).

Proof. 0 < b_{n} ≤ a_{n} for n ≥ 1, so 0 ≤ −a_{n}(I_{1}/I)(a_{n}, b_{n}) ≤ 1 by Theorem 2.2; i.e., 0 ≤ ϵ_{n} ≤ 2^{n}(a_{n}^{2} − b_{n}^{2}).

Substitute (a_{n+1}^{2} − b_{n+1}^{2})a_{n+1} = (a_{n}^{2} − b_{n}^{2})(a_{n} − b_{n})/8:

ϵ_{n+1} = 2^{n}(a_{n}^{2} − b_{n}^{2}) (\frac{−(a_{n} − b_{n})}{4} (\frac{I_{1}}{I} ) (a_{n+1}, b_{n+1})) .

Then apply Theorems 2.6 and 2.5:

ϵ_{n+1} = 2^{n}(a_{n}^{2} − b_{n}^{2}) (\frac{−1}{2} − a_{n} (\frac{I_{1}}{I} ) (a_{n}, b_{n})) = −2^{n−1}(a_{n}^{2} − b_{n}^{2}) + ϵ_{n}.

Thus ϵ_{0} = ϵ_{n} + ∑_{0≤i<n} 2^{i−1}(a_{i}^{2} − b_{i}^{2}) by induction. □

Theorem 4.4. If 1 ≤ a_{0}/b_{0} ≤ 1 + 2^{2^{m}} then 1 ≤ a_{n}/b_{n} ≤ 1 + 2^{2^{m−n}} for n ≥ 0.

Proof. If 1 ≤ a_{n}/b_{n} ≤ 1 + 2^{2^{m−n}} then a_{n+1}/b_{n+1} = (\sqrt{a_{n}/b_{n}} + \sqrt{b_{n}/a_{n}})/2 ≤ \sqrt{a_{n}/b_{n}} ≤ \sqrt{1 + 2^{2^{m−n}}} ≤ 1 + 2^{2^{m−n−1}}. □

Theorem 4.5. If m ≥ 0 and 1 ≤ a_{0}/b_{0} ≤ 1 + 2^{2^{m}} then 1 ≤ a_{n}/b_{n} ≤ 1 + 2^{3−2^{n+1−m}} for n ≥ m.

Proof. The base case n = m follows from Theorem 4.4 since 2^{2^{m−m}} = 2^{3−2^{m+1−m}}.

If a_{n}/b_{n} = 1+ϵ with 0 ≤ ϵ ≤ 2^{3−2^{n+1−m}} then 4+4ϵ+ϵ^{2} ≤ (4+ϵ^{2} +ϵ^{4}/16)(1+ϵ), so (1 + ϵ) + 2 + 1/(1 + ϵ) ≤ (2 + ϵ^{2}/4)^{2}, so \sqrt{a_{n}/b_{n}} + \sqrt{b_{n}/a_{n}} ≤ 2 + ϵ^{2}/4, so a_{n+1}/b_{n+1} ≤ 1 + ϵ^{2}/8 ≤ 1 + 2^{3−2^{n+2−m}}. □

5. Computing logarithm intervals

This section presents several algorithms that, given an interval containing x, compute an interval containing log x. See [\underline{6}] for fast subroutines to compute sums, differences, products, quotients, and square roots of intervals.

These algorithms use a parameter p to decide when to stop. For the “arbitrary x” algorithms, the output interval has approximately p bits of precision if the input interval does. For the “super-size x” algorithms, the output precision is limited to about 4 lg x bits, even if p is much larger. The “super-size x” algorithms also slow down as lg lg x grows past lg p; if lg lg x is much larger than lg p, one should use the “arbitrary x” algorithms instead.

Beware that, as discussed in [\underline{5}], the usual algorithms for arithmetic operations such as square root—and for sequences of arithmetic operations—contain many redundancies that can be eliminated. A sequence of AGM steps can be made almost three times faster than reported in [\underline{11}, Theorem 9.1], for example. Note also that Borwein and Borwein in [\underline{8}, page 222] observed speed improvements from a “quartic AGM” in which one computes \sqrt{a_{n+2}} and \sqrt{b_{n+2}} directly from \sqrt{a_{n}} and \sqrt{b_{n}}; my impression is that these speedups become slowdowns when the square-root algorithms are properly optimized, but I will leave experiments to the reader.

Computing log x for a super-size x. Let x > 4 be a real number. Define a_{0} = 1 and b_{0} = b = (2x/(x^{2}−1))^{2}. Note that \sqrt{b^{−1}}+\sqrt{b^{−1} + 1} = x, so L(b) = log x, where L is defined in Section 3; note also that 0 < b < 1/2. Define a_{n+1} = (a_{n} + b_{n})/2 and b_{n+1} = \sqrt{a_{n}b_{n}} for n ≥ 0, as in Section 4. Define c = 1 + (I_{1}/I)(1, b); note that 0 ≤ c ≤ 1 by Theorem 2.2.

Compute an interval containing b. Successively compute intervals containing a_{1}, b_{1}, a_{2}, b_{2}, …, a_{n}, b_{n}, stopping when a_{n} −b_{n} is no longer clearly larger than 1/2^{p}. The number of steps n is at most about lg lg x + lg p by Theorem 4.5.

Compute an interval containing [0, 2^{n}(a_{n}^{2} − b_{n}^{2})] + ∑_{0≤i<n} 2^{i−1}(a_{i}^{2} − b_{i}^{2}). By Theorem 4.3, this interval contains (1 − b^{2})(1 − c). Divide by an interval containing 1 − b^{2}, and subtract from 1, to obtain an interval containing c. Finally, compute an interval containing

\begin{matrix}[\frac{(\frac{1}{2}b − \frac{3}{16}b^{2} + \frac{9}{32}b^{3})c(1 + b)^{1/2} − (1 + b)^{−1} + (2 + b^{2})(1 + b)^{−1/2}}{(2 + \frac{1}{2}b^{2} + \frac{9}{32}b^{4})c + b^{2}} , \\ \frac{(\frac{1}{2}b)c(1 + b)^{1/2} − (1 + b)^{−1} + (2 + b^{2} + \frac{3}{8}b^{3} + \frac{9}{8}b^{4})(1 + b)^{−1/2}}{(2 + \frac{1}{2}b^{2})c + (b^{2} + \frac{9}{8}b^{4})} ]. \end{matrix}

By Theorems 3.1 and 3.2, this interval contains L(b) = log x. Note that terms such as \frac{9}{8}b^{4} have very little effect on the output and can be replaced by crude bounds.

At this point one can also compute an interval containing π with a few additional operations: the interval

[(2 + \frac{1}{2}b^{2})L(b) − (\frac{1}{2}b)(1 + b)^{1/2}, (2 + \frac{1}{2}b^{2} + \frac{9}{32}b^{4})L(b) − (\frac{1}{2}b − \frac{3}{16}b^{2} + \frac{9}{32}b^{3})(1 + b)^{1/2}]

contains I(1, b) by Theorem 3.1, and the interval [2b_{n}I(1, b), 2a_{n}I(1, b)] contains π by Theorem 4.2. This very fast computation of π will be exploited later.

Numerical stability: The computation of c as 1−(1−c) loses about lg log(4/b) ≈ lg lg x bits of precision, since c ≈ 1/ log(4/b). The other arithmetic operations are stable, each losing only a bounded number of bits of precision. However, the intervals in Theorems 3.1 and 3.2 are inherently limited to about lg(1/b^{2}) ≈ 4 lg x bits of precision.

Computing log x for several super-size x’s. After computing a super-size log as explained above, one can compute another super-size log at somewhat higher speed, by taking advantage of the π interval obtained from the first computation.

Starting from x, define and compute b, a_{1}, b_{1}, …, a_{n}, b_{n} as above, but skip the computation of ∑_{i} 2^{i−1}(a_{i}^{2} − b_{i}^{2}). Compute an interval containing [π/2a_{n}, π/2b_{n}]; by Theorem 4.2, this interval contains I(1, b). Compute an interval containing

[\frac{I(1, b) + (\frac{1}{2}b − \frac{3}{16}b^{2} + \frac{9}{32}b^{3})(1 + b)^{1/2}}{2 + \frac{1}{2}b^{2} + \frac{9}{32}b^{4}} , \frac{I(1, b) + (\frac{1}{2}b)(1 + b)^{1/2}}{2 + \frac{1}{2}b^{2}} ];

by Theorem 3.1, this interval contains L(b) = log x.

A further improvement is available in applications that compute logarithms only to divide them by each other. The above method is to use I_{1}/I to compute log x and π, then use I to compute π/ log y, then divide to obtain (log y)/ log x; it is somewhat faster to use I twice to compute π/ log x and π/ log y, then divide to obtain (log y)/ log x.

Computing log x for an arbitrary x. Given an interval containing a positive real number x, find an integer k such that 2^{k}x is larger, but not much larger, than 2^{p/4}. Then use the algorithm above to compute intervals containing the logs of the super-size numbers 2^{k}x and 2^{⌈p/4⌉}. Extract intervals containing log 2 = (log 2^{⌈p/4⌉})/ ⌈p/4⌉ and log x = log 2^{k}x − k log 2.

The time taken by this algorithm is practically independent of x. However, readers trying to formulate a theorem along these lines are cautioned to consider the possibility that k has many digits.

When x is between (for example) 2 and 2^{p}, the following algorithm is faster: select an integer m > lg p − 2 − lg lg x, compute an interval containing the supersize number x^{2^{m}} by repeated squaring, compute an interval containing log x^{2^{m}}, and divide by 2^{m}. Beware that this approach is slow when x is very close to 1.

Computing log x for several arbitrary x’s. The general problem of computing logarithms of ℓ numbers can be reduced to the problem of computing logarithms of ℓ + 1 super-size numbers, one of which is a power of 2, as in the case ℓ = 1 above.

When all the numbers are in reasonable ranges, and when ℓ is small, it is faster to repeatedly square each number, obtaining ℓ super-size numbers. However, for large ℓ, this is slower than the power-of-2 method.

6. Previous logarithm algorithms using elliptic integrals

In the following survey, the word “optimal” means “as fast as possible among all the techniques that I know.” It is not meant to exclude the possibility of future improvements.

Salamin, as reported in [\underline{2}, Item 143], proposed using the AGM iteration to compute π/ log x for super-size x, and thus to compute log x using π. This is the optimal strategy when x is super-size and π is already known.

Salamin proposed computing π as follows: compute exp 1 by a different method, then compute (exp 1)^{2^{m}} = exp 2^{m} by repeated squaring, then use the AGM iteration to compute π/ log exp 2^{m} = π/2^{m}. Schroeppel in [\underline{2}, Item 144] proposed using 2^{2^{m}} and 2^{2^{m}} exp 1 instead of exp 2^{m}. Both methods are considerably slower than optimal.

For x in a reasonable range such as [2, 4], Salamin proposed computing log x by computing log x^{2^{m}}. This is the optimal strategy when only one log is to be computed, although it is not optimal when many logs are to be computed.

Brent in [\underline{10}] proposed a different log algorithm using incomplete elliptic integrals. Brent’s algorithm is somewhat slower than optimal: it saves half the AGM steps by focusing on moderate values of b, but it works with more than two values of b.

Brent, and independently Salamin in [\underline{16}], also proposed using the Legendre-Gauss formula I(1, 2^{−1/2})(I(1, 2^{−1/2}) + I_{1}(1, 2^{−1/2})) = π/2 to compute π. This is faster than Salamin’s previous method of computing π, but in the context of log computation it is not optimal. The same comment applies to several subsequent methods of computing π, not cited here.

Brent in [\underline{11}, Section 9] proposed computing log x for arbitrary x by computing log 2^{k}x where k is chosen so that 2^{k}x is super-size. This is the optimal strategy when many logs are to be computed.

Borwein and Borwein in [\underline{7}, Section 4] proposed another log algorithm in the spirit of [\underline{10}] but relying solely on complete elliptic integrals. The strategy is somewhat slower than optimal.

Newman in [\underline{15}] proposed computing π and log x for super-size x by using the AGM iteration to compute (log x)/π, using it again to compute (log(x + 1))/π ≈ (log x + 1/x)/π, and subtracting. Beware that the subtraction loses about lg x bits of precision. A better (but still suboptimal) strategy is to use x exp 1 in place of x + 1, with exp 1 computed by a different method.

Newman’s goal was not maximum speed, but maximum simplicity. In particular, Newman avoided the Legendre-Gauss formula and any other use of I_{1}. However, I disagree with Newman’s view of “Landen’s transformation law” (Theorem 2.6 here) as “heavy use of elliptic function theory”; I hope I have convinced the reader that the basic properties of I_{1} can be established just as easily as the basic properties of I. (The Legendre-Gauss formula is not much more difficult.)

Borwein and Borwein in [\underline{8}, Algorithm 7.2] proposed the following method of computing log x for, e.g., x between 2 and 4: use I_{1}/I to compute two super-size logarithms, namely log 2^{k} and log 2^{k}x for an appropriate k. The use of I_{1}/I is optimal for computing one super-size logarithm, but it is suboptimal for computing two or more.

References

[\underline{1}] Robert S. Anderssen, Richard P. Brent (editors), The complexity of computational problem solving, University of Queensland Press, Brisbane, 1976. ISBN 0–7022–1213–X. Available

from http://web.comlab.ox.ac.uk/oucl/work/richard.brent/pub/pub031.html.

[\underline{2}] M. Beeler, R. W. Gosper, R. Schroeppel, HAKMEM, Artificial Intelligence Memo No. 239, Massachusetts Institute of Technology, 1972. Available from http://www.inwap.com/pdp10/

hbaker/hakmem/hakmem.html.

[\underline{3}] Daniel J. Bernstein, Fast multiplication and its applications. Available from http://cr.yp.

to/papers.html.

[\underline{4}] Daniel J. Bernstein, Computing logarithm floors in essentially linear time, draft. Available

from http://cr.yp.to/papers.html.

[\underline{5}] Daniel J. Bernstein, Removing redundancy in high-precision Newton iteration, draft. [\underline{6}] Daniel J. Bernstein, Fast, harmless intervals around high-precision floating-point numbers,

draft. [\underline{7}] Jonathan M. Borwein, Peter B. Borwein, The arithmetic-geometric mean and fast computa-

tion of elementary functions, SIAM Review 26 (1984), 351–366. MR 86d:65029. [\underline{8}] Jonathan M. Borwein, Peter B. Borwein, Pi and the AGM, Wiley, New York, 1987. ISBN

0–471–83138–7. MR 89a:11134. [\underline{9}] Richard P. Brent, The complexity of multiple-precision arithmetic, in [\underline{1}] (1976), 126–165.

Available from http://web.comlab.ox.ac.uk/oucl/work/richard.brent/pub/pub032.html.

[\underline{10}] Richard P. Brent, Fast multiple-precision evaluation of elementary functions, Journal of the ACM 23 (1976), 242–251. ISSN 0004–5411. MR 52:16111. Available from http://web.

comlab.ox.ac.uk/oucl/work/richard.brent/pub/pub034.html.

[\underline{11}] Richard P. Brent, Multiple-precision zero-finding methods and the complexity of elementary function evaluation, in [\underline{17}] (1976), 151–176. MR 54:11843. Available from http://web.

comlab.ox.ac.uk/oucl/work/richard.brent/pub/pub028.html.

[\underline{12}] Krzysztof Diks, Wojciech Ritter (editors), Mathematical foundations of computer science 2002: 27th international symposium, MFCS 2002, Warsaw, Poland, 26–30.08.2002: proceedings, Lecture Notes in Computer Science, 2420, Springer, Berlin, 2002.

[\underline{13}] Mika Hirvensalo, Juhani Karhumäki, Computing partial information out of intractable one—

the first digit of 2^{n} at base 3 as an example, in [\underline{12}] (2002), 319–327. [\underline{14}] Wolfram Luther, Werner Otten, Computation of standard interval functions in multiple-

precision interval arithmetic, Interval Computations (1994), 78–99. ISSN 0135–4868. MR 96g:65047.

[\underline{15}] Donald J. Newman, A simplified version of the fast algorithms of Brent and Salamin, Math-

ematics of Computation 44 (1985), 207–210. MR 86e:65030. [\underline{16}] Eugene Salamin, Computation of π using arithmetic-geometric mean, Mathematics of Com-

putation 30 (1976), 565–570. ISSN 0025–5718. MR 53:7928. [\underline{17}] Joseph F. Traub, Analytic computational complexity, Academic Press, New York, 1976. MR

52:15938. [\underline{18}] E. T. Whittaker, G. N. Watson, A course of modern analysis, Cambridge University Press,

Cambridge, 1927. ISBN 0–521–58807–3. MR 97k:01072.

Department of Mathematics, Statistics, and Computer Science (M/C 249), The University of Illinois at Chicago, Chicago, IL 60607–7045

E-mail address: djb@cr.yp.to