MULTIPLE-PRECISION ZERO-FINDING METHODS AND THE COMPLEXITY OF ELEMENTARY FUNCTION EVALUATION^{1}

Richard P Brent

Computer Centre, Australian National University, Canberra, A.C.T. 2600, Australia

ABSTRACT

We consider methods for finding high-precision approximations to simple zeros of smooth functions. As an application, we give fast methods for evaluating the elementary functions log(x), exp(x), sin(x) etc. to high precision. For example, if x is a positive floating-point number with an n-bit fraction, then (under rather weak assumptions) an n-bit approximation to log(x) or exp(x) may be computed in time asymptotically equal to 13M(n) log_{2} n as n → ∞, where M(n) is the time required to multiply floating-point numbers with n-bit fractions. Similar results are given for the other elementary functions, and some analogies with operations on formal power series are mentioned.

1 Introduction

When comparing methods for solving nonlinear equations or evaluating functions, it is customary to assume that the basic arithmetic operations (addition, multiplication, etc.) are performed with some fixed precision. However, an irrational number can only be approximated to arbitrary accuracy if the precision is allowed to increase indefinitely. Thus, we shall consider iterative processes using variable precision. Usually the precision will increase as the computation proceeds, and the final result will be obtained to high precision. Of course, we could use the same (high) precision throughout, but then the computation would take longer than with variable precision, and the final result would be no more accurate.

Assumptions

For simplicity we assume that a standard multiple-precision floating-point number representation is used, with a binary fraction of n bits, where n is large. The exponent length is fixed, or may grow as o(n) if necessary. To avoid table-lookup methods, we assume a machine with a finite random-access memory and a fixed number of sequential tape units. Formally, the results hold for multitape Turing machines.

1First appeared in Analytic Computational Complexity (edited by J F Traub), Academic Press, New York, 1975, 151–176. Retyped with minor corrections and postscript by Frances Page at Oxford University Computing Laboratory, 1999 (urls updated 2005).

Copyright \begin{matrix}c \\ ◯ \end{matrix} 1975, 1999, 2005, R. P. Brent and Academic Press, Inc. rpb028 typeset using L^{A}TEX.

Precision n Operations

An operation is performed with precision n if the operands and result are floating-point numbers as above (i.e., precision n numbers), and the relative error in the result is O(2^{−n}).

Precision n Multiplication

Let M(n) be the time required to perform precision n multiplication. (Time may be regarded as the number of single-precision operations, or the number of bit operations, if desired.) The classical method gives M(n) = O(n^{2}), but methods which are faster for large n are known. Asymptotically the fastest method known is that of Schönhage and Strassen [71], which gives

(1.1)  M(n) = O(n log(n) log log(n)) as n → ∞ .

Our results do not depend on the algorithm used for multiplication, provided M(n) satisfies the following two conditions.

(1.2)  n = o(M(n)) , i.e., lim_{n→∞} n/M(n) = 0 ;

and, for any α > 0,

(1.3)  M(αn) ∼ αM(n) , i.e., lim_{n→∞} \frac{M(αn)}{αM(n)} = 1 .

Condition (1.2) enables us to neglect additions, since the time for an addition is O(n), which is asymptotically negligible compared to the time for a multiplication. Condition (1.3) certainly holds if

M(n) ∼ cn[log(n)]^{β}[log log(n)]^{γ} ,

though it does not hold for some implementations of the Schönhage-Strassen method. We need (1.3) to estimate the constants in the asymptotic “O” results: if the constants are not required then much weaker assumptions suffice, as in Brent [75a,b].

The following lemma follows easily from (1.3).

Lemma 1.1

If 0 < α < 1, M(n) = 0 for n < 1, and c_{1} < \frac{1}{1−α} < c_{2}, then

c_{1}M(n) < ∑_{k=0}^{∞} M(α^{k}n) < c_{2}M(n)

for all sufficiently large n.

2 Basic Multiple-precision Operations

In this section we summarize some results on the time required to perform the multiple-precision operations of divison, extraction of square roots, etc. Additional results are given in Brent [75a].

Reciprocals

Suppose a ≠ 0 is given and we want to evaluate a precision n approximation to 1/a. Applying Newton’s method to the equation

f(x) ≡ a − 1/x = 0

gives the well-known iteration

x_{i+1} = x_{i} − x_{i}ε_{i} ,

where

ε_{i} = ax_{i} − 1 .

Since the order of convergence is two, only k ∼ log_{2} n iterations are required if x_{0} is a reasonable approximation to 1/a, e.g., a single-precision approximation.

If ε_{k−1} = O(2^{−n/2}), then ε_{k} = O(2^{−n}), so at the last iteration it is sufficient to perform the multiplication of x_{k−1} by ε_{k−1} using precision n/2, even though ax_{k−1} must be evaluated with precision n. Thus, the time required for the last iteration is M(n) + M(n/2) + O(n). The time for the next to last iteration is M(n/2)+M(n/4)+O(n/2), since this iteration need only give an approximation accurate to O(2^{−n/2}), and so on. Thus, using Lemma 1.1, the total time required is

I(n) ∼ (1 + \frac{1}{2})(1 + \frac{1}{2} + \frac{1}{4} + ⋯)M(n) ∼ 3M(n)

as n → ∞.

Division

Since b/a = b(1/a), precision n division may be done in time

D(n) ∼ 4M(n)

as n → ∞.

Inverse Square Roots

Asymptotically the fastest known method for evaluating a^{− \frac{1}{2}} to precision n is to use the thirdorder iteration

x_{i+1} = x_{i} − \frac{1}{2}x_{i}(ε_{i} − \frac{3}{4}ε_{i}^{2}) ,

where

ε_{i} = ax_{i}^{2} − 1 .

At the last iteration it is sufficient to evaluate ax_{i}^{2} to precision n, ε_{i}^{2} to precision n/3, and x_{i}(ε_{i} − \frac{3}{4}ε_{i}^{2}) to precision 2n/3. Thus, using Lemma 1.1 as above, the total time required is

Q(n) ∼ (2 + \frac{1}{3} + \frac{2}{3})(1 + \frac{1}{3} + \frac{1}{9} + ⋯)M(n) ∼ \frac{9}{2}M(n)

as n → ∞.

Square Roots

Since

\sqrt{a} = \{ \begin{matrix}a.a^{− \frac{1}{2}} & if & a > 0 , \\ 0 & if & a = 0 , \end{matrix}

we can evaluate \sqrt{a} to precision n in time

R(n) ∼ \frac{11}{2}M(n)

as n → ∞. Note the direct use of Newton’s method in the form

(2.1)  x_{i+1} = \frac{1}{2}(x_{i} + a/x_{i})

or

(2.2)  x_{i+1} = x_{i} + (\frac{a − x_{i}^{2}}{2x_{i}} )

is asymptotically slower, requiring time ∼ 8M(n) or ∼ 6M(n) respectively.

3 Variable-precision Zero-finding Methods

Suppose ζ ≠ 0 is a simple zero of the nonlinear equation

f(x) = 0 .

Here, f(x) is a sufficiently smooth function which can be evaluated near ζ, with absolute error O(2^{−n}), in time w(n). We consider some methods for evaluating ζ to precision n. Since we are interested in results for very large n, the time required to obtain a good starting approximation is neglected.

Assumptions

To obtain sharp results we need the following two assumptions, which are similar to (1.2) and (1.3):

(3.1)  M(n) = o(w(n)) , i.e., lim_{n→∞} M(n)/w(n) = 0 ;

and, for some α ≥ 1 and all β > 0,

(3.2)  w(βn) ∼ β^{α}w(n)

as n → ∞.

From (3.1), the time required for a multiplication is negligible compared to the time for a function evaluation, if n is sufficiently large. (3.2) implies (3.1) if α > 1, and (3.2) certainly holds if, for example,

w(n) ∼ cn^{α}[log(n)]^{γ}[log log(n)]^{δ} .

The next lemma follows from our assumptions in much the same way as Lemma 1.1.

Lemma 3.1

If 0 < β < 1, w(n) = 0 for n < 1, and

c_{1} < 1/(1 − β^{α}) < c_{2} ,

then

c_{1}w(n) < ∑_{k=0}^{∞} w(β^{k}n) < c_{2}w(n)

for all sufficiently large n.

A Discrete Newton Method

To illustrate the ideas of variable-precision zero-finding methods, we describe a simple discrete Newton method. More efficient methods are described in the next three sections, and in Brent [75a].

Consider the iteration

x_{i+1} = x_{i} − f(x_{i})/g_{i} ,

where g_{i} is a one-sided difference approximation to f′(x_{i}), i.e.,

g_{i} = \frac{f(x_{i} + h_{i}) − f(x_{i})}{h_{i}} .

If ε = |x_{i} − ζ| is sufficiently small, f(x_{i}) is evaluated with absolute error O(ε_{i}^{2}), and h_{i} is small enough that

(3.3)  g_{i} = f′(x_{i}) + O(ε_{i}) ,

then the iteration converges to ζ with order at least two. To ensure (3.3), take h_{i} of order ε_{i}, e.g. h_{i} = f(x_{i}).

To obtain ζ to precision n, we need two evaluations of f with absolute error O(2^{−n}), preceded by two evaluations with error, O(2^{−n/2}), etc. Thus the time required is

(3.4)  t(n) ∼ 2(1 + 2^{α} + 2^{−2α} + ⋯)w(n) .

Asymptotic Constants

We say that a zero-finding method has asymptotic constant C(α) if, to find a simple zero ζ ≠ 0 to precision n, the method requires time t(n) ∼ C(α)w(n) as n → ∞. (The asymptotic constant as defined here should not be confused with the “asymptotic error constant” as usually defined for single-precision zero-finding methods.)

For example, from (3.4), the discrete Newton method described above has asymptotic constant

C_{N}(α) = 2/(1 − 2^{−α}) ≤ 4 .

Note that the time required to evaluate ζ to precision n is only a small multiple of the time required to evaluate f(x) with absolute error O(2^{−n}). (If we used fixed precision, the time to evaluate ζ would be of order log(n) times the time to evaluate f(x).)

4 A Variable-precision Secant Method

The secant method is known to be more efficient than the discrete Newton method when fixedprecision arithmetic is used. The same is true with variable-precision arithmetic, although the ratio of efficiencies is no longer constant, but depends on the exponent α in (3.2). Several secantlike methods are described in Brent [75a]; here we consider the simplest such method, which is also the most efficient if α < 4.5243 …

The secant iteration is

x_{i+1} = x_{i} − f_{i} (\frac{x_{i} − x_{i−1}}{f_{i} − f_{i−1}} ) ,

where f_{i} = f(x_{i}), and we assume that the function evaluations are performed with suficient accuracy to ensure that the order of convergence is at least ρ = (1 + \sqrt{5})/2 = 1.6180 …, the larger root of

(4.1)  ρ^{2} = ρ + 1 .

Let ε = |x_{i−1} − ζ|. Since the smaller root of (4.1) lies inside the unit circle, we have

x_{i} − ζ = O(ε^{ρ})

and

x_{i+1} − ζ = O(ε^{ρ^{2}}) .

To give order ρ, f_{i} must be evaluated with absolute error O(ε^{ρ^{2}}). Since f_{i} = O(|x_{i} − ζ|) = O(ε^{ρ}), it is also necessary to evaluate (f_{i} − f_{i−1})/(x_{i} − x_{i−1}) with relative error O(ε^{ρ^{2}−ρ}), but |x_{i} −x_{i−1}| ∼ ε, so it is necessary to evaluate f_{i−1} with absolute error O(ε^{ρ^{2}−ρ+1}). [ Since f_{i} must be evaluated with absolute error O(ε^{ρ^{2}}), f_{i−1} must be evaluated with absolute error O(ε^{ρ}), but ρ^{2} − ρ + 1 = 2 > ρ, so this condition is superfluous.]

The conditions mentioned are sufficient to ensure that the order of convergence is at least ρ. Thus, if we replace ε^{ρ^{2}} by 2^{−n}, we see that ζ may be evaluated to precision n if f is evaluated with absolute errors O(2^{−n}), O(2^{−2nρ^{−2}}), O(2^{−2nρ^{−3}}), O(2^{−2nρ^{−4}}), … It follows that the asymptotic constant of the secant method is

C_{S}(α) = 1 + (2ρ^{−2})^{α}/(1 − ρ^{−α}) ≤ C_{S}(1) = 3 .

The following lemma states that the secant method is asymptotically more efficient than the discrete Newton method when variable precision is used.

Lemma 4.1

C_{S}(α) < C_{N}(α) for all α ≥ 1. In fact C_{S}(α)/C_{N}(α) decreases monotonically from \frac{3}{4} (when α = 1) to \frac{1}{2} (as α → ∞).

5 Other Variable-precision Interpolatory Methods

With fixed precision, inverse quadratic interpolation is more efficient than linear interpolation, and inverse cubic interpolation is even more efficient, if the combinatory cost (i.e., “overhead”) is negligible. With variable precision the situation is different. Inverse quadratic interpolation is slightly more efficient than the secant method, but inverse cubic interpolation is not more efficient than inverse quadratic interpolation if α ≤ 4.6056 … Since the combinatory cost of inverse cubic interpolation is considerably higher than that of inverse quadratic interpolation, the inverse cubic method appears even worse if combinatory costs are significant.

Inverse Quadratic Interpolation

The analysis of variable-precision methods using inverse quadratic interpolation is similar to that for the secant method, so we only state the results. The order ρ = 1.8392 … is the positive root of ρ^{3} = ρ^{2}+ρ+1. It is convenient to define σ = 1/ρ = 0.5436 … To evaluate ζ to precision n requires evaluations of f to (absolute) precision n, (1 − σ + σ^{2})n, and σ^{j}(1 − σ − σ^{2} + 2σ^{3})n for j = 0, 1, 2, … Thus, the asymptotic constant is

\begin{aligned}C_{Q}(α) & = & 1 + (1 − σ + σ^{2})^{α} + (3σ^{3})^{α}/(1 − σ^{α}) \\ & ≤ & C_{Q}(1) = \frac{1}{2}(7 − 2σ − σ^{2}) = 2.8085 …. \end{aligned}

Lemma 5.1

C_{Q}(α) < C_{S}(α) for all α ≥ 1. In fact, C_{Q}(α)/C_{S}(α) increases monotonically from 0.9361… (when α = 1) to 1 (as α → ∞).

Inverse Cubic Interpolation, etc

If μ = 0.5187 … is the positive root of μ^{4} +μ^{3} +μ^{2} +μ = 1, then the variable-precision method of order 1/μ = 1.9275 …, using inverse cubic interpolation, has asymptotic constant

\begin{aligned}C_{C}(α) & = & 1 + (1 − μ + μ^{2})^{α} + (1 − μ − μ^{2} + 2μ^{3})^{α} + (4μ^{4})^{α}/(1 − μ^{α}) \\ & ≤ & C_{C}(1) = (13 − 6μ − 4μ^{2} − 2μ^{3})/3 = 2.8438 …. \end{aligned}

Note that C_{C}(1) > C_{Q}(1). Variable-precision methods using inverse interpolation of arbitrary degree are described in Brent [75a]. Some of these methods are slightly more efficient than the inverse quadratic interpolation method if α is large, but inverse quadratic interpolation is the most efficient method known for α < 4.6056 … In practice α is usually 1, 1\frac{1}{2} or 2.

An Open Question

Is there a method with asymptotic constant C(α) such that C(1) < C_{Q}(1)?

6 Variable-precision Methods using Derivatives

In Sections 3 to 5 we considered methods for solving the nonlinear equation f(x) = 0, using only evaluations of f. Sometimes it is easy to evaluate f′(x), f′′(x), … once f(x) has been evaluated, and the following theorem shows that it is possible to take advantage of this. For an application, see Section 10.

Theorem 6.1

If the time to evaluate f(x) with an absolute error O(2^{−n}) is w(n), where w(n) satisfies conditions (3.1) and (3.2), and (for k = 1, 2, …) the time to evaluate F^{(k)}(x) with absolute error O(2^{−n}) is w_{k}(n), where

w_{k}(n) = o(w(n))

as n → ∞, then the time to evaluate a simple zero ζ ≠ 0 of f(x) to precision n is

t(n) ∼ w(n)

as n → ∞.

Proof

For fixed k ≥ 1, we may use a direct or inverse Taylor series method of order k + 1. The combinatory cost is of order k log(k + 1)M(n) (see Brent and Kung [75]). From (3.1), this is o(w(n)) as n → ∞. Thus,

t(n) ≤ [1 − (k + 1)^{−α}]^{−1}w(n) + o(w(n))
≤ (1 + \frac{1}{k} + o(1))w(n) .

For sufficiently large n, the “o(1)” term is less than 1/k, so

t(n) ≤ (1 + \frac{2}{k})w(n) .

Given ε > 0, choose k ≥ 2/ε. Then, for all sufficiently large n,

w(n) ≤ t(n) ≤ (1 + ε)w(n) ,

so t(n) ∼ w(n) as n → ∞.

Corollary 6.1

If the conditions of Theorem 6.1 hold, f: [a, b] → I, f′(x) ≠ 0 for x ∈ [a, b], and g is the inverse function of f, then the time to evaluate g(y) with absolute error O(2^{−n}), for y ∈ I, is

w_{g}(n) ∼ w(n)

as n → ∞.

Note

Corollary 6.1 is optimal in the sense that, if w_{g}(n) ∼ cw(n) for some constant c < 1, then w(n) ∼ cw_{g}(n) by the same argument, so w(n) ∼ c^{2}w(n), a contradiction. Hence, c = 1 is minimal.

7 The Arithmetic-geometric Mean Iteration

Before considering the multiple-precision evaluation of elementary functions, we recall some properties of the arithmetic-geometric (A–G) mean iteration of Gauss [1876]. Starting from any two positive numbers a_{0} and b_{0}, we may iterate as follows:

a_{i+1} = \frac{a_{i} + b_{i}}{2} arithmetic mean

and

b_{i+1} = \sqrt{a_{i}b_{i}} geometric mean

for i = 0, 1, …

Second-order Convergence

The A–G mean iteration is of computational interest because it converges very fast. If b_{i} ≪ a_{i}, then

b_{i+1}/a_{i+1} = \frac{2\sqrt{b_{i}/a_{i}}}{1 + b_{i}/a_{i}} ≃ 2\sqrt{b_{i}/a_{i}} ,

so only about | log_{2}(a_{0}/b_{0})| iterations are required before a_{i}/b_{i} is of order 1. Once a_{i} and b_{i} are close together the convergence is second order, for if b_{i}/a_{i} = 1 − ε_{i} then

ε_{i+1} = 1 − b_{i+1}/a_{i+1} = 1 − 2(1 − ε_{i})^{\frac{1}{2}} /(2 − ε_{i}) = ε_{i}^{2}/8 + O(ε_{i}^{3}) .

Limit of the A–G Mean Iteration

There is no essential loss of generality in assuming that a_{0} = 1 and b_{0} = cos ϕ for some ϕ. If a = lim_{1→∞}a_{i} = lim_{i→∞}b_{i}, then

(7.1)  a = \frac{π}{2K(ϕ)} ,

where K(ϕ) is the complete elliptic integral of the first kind, i.e.,

K(ϕ) = ∫_{0}^{π/2}\frac{1}{\sqrt{1 − sin^{2} ϕ sin^{2} θ}} dθ .

(A simple proof of (7.1) is given in Melzak [73].)

Also, if c_{0} = sin ϕ, c_{i+1} = a_{i} − a_{i+1} (i = 0, 1, …), then

(7.2)  ∑_{i=0}^{∞} 2^{i−1}c_{i}^{2} = 1 − \frac{E(ϕ)}{K(ϕ)} ,

where E(ϕ) is the complete elliptic integral of the second kind, i.e.,

E(ϕ) = ∫_{0}^{π/2}\sqrt{1 − sin^{2} ϕ sin^{2} θ} dθ .

Both (7.1) and (7.2) were known by Gauss.

Legendre’s Identity

For future use, we note the identity

(7.3)  K(ϕ)E(ϕ′) + K(ϕ′)E(ϕ) − K(ϕ)K(ϕ′) = \frac{1}{2}π ,

where ϕ+ϕ′ = \frac{1}{2}π. (Legendre [11] proved by differentiation that the left side of (7.3) is constant, and the constant may be determined by letting ϕ → 0.)

8 Fast Multiple-precision Evaluation of π

The classical methods for evaluating π to precision n take time O(n^{2}): see, for example, Shanks and Wrench [62]. Several methods which are asymptotically faster than O(n^{2}) are known. For example, in Brent [75a] a method which requires time O(M(n) log^{2}(n)) is described. From the bound (1.1) on M(n), this is faster than O(n^{1+ε}) for any ε > 0.

Asymptotically the fastest known methods require time O(M(n) log(n)). One such method is sketched in Beeler et al [72]. The method given here is faster, and does not require the preliminary computation of e.

The Gauss-Legendre Method

Taking ϕ = ϕ′ = π/4 in (7.3), and dividing both sides by π^{2}, we obtain

(8.1)  [2K(π/4)E(π/4) − K^{2}(π/4)]/π^{2} = \frac{1}{2π} .

However, from the A–G mean iteration with a_{0} = 1 and b_{0} = 2^{− \frac{1}{2}} , and the relations (7.1) and (7.2), we can evaluate K(π/4)/π and E(π/4)/π, and thus the left side of (8.1). A division then gives π. (The idea of using (7.3) in this way occurred independently to Salamin [75] and Brent [75b].) After a little simplification, we obtain the following algorithm (written in pseudo-Algol):

A ← 1; B ← 2^{− \frac{1}{2}} ; T ← 1/4; X ← 1;
while A − B > 2^{−n} do

begin Y ← A; A ← \frac{1}{2}(A + B); B ← \sqrt{BY} ;

T ← T − X(A − Y)^{2} ;
X ← 2X

end;

return A^{2}/T [or, better, (A + B)^{2}/(4T)] .

The rate of convergence is illustrated in Table 8.1.

Table 8.1: Convergence of the Gauss-Legendre Method

Iteration  A^{2}/T − π  π − (A + B)^{2}/(4T)
0  8.6 × 10^{−1}  2.3 × 10^{−1}
1  4.6 × 10^{−2}  1.0 × 10^{−3}
2  8.8 × 10^{−5}  7.4 × 10^{−9}
3  3.1 × 10^{−10}  1.8 × 10^{−19}
4  3.7 × 10^{−21}  5.5 × 10^{−41}
5  5.5 × 10^{−43}  2.4 × 10^{−84}
6  1.2 × 10^{−86}  2.3 × 10^{−171}
7  5.8 × 10^{−174}  1.1 × 10^{−345}
8  1.3 × 10^{−348}  1.1 × 10^{−694}
9  6.9 × 10^{−698}  6.1 × 10^{−1393}

Since the A–G mean iteration converges with order 2, we need ∼ log_{2} n iterations to obtain precision n. Each iteration involves one (precision n) square root, one multiplication, one squaring, one multiplication by a power of two, and some additions. Thus from the results of Section 2, the time required to evaluate π is ∼ \frac{15}{2}M(n) log_{2} n.

Comments

1. Unlike Newton’s iteration, the A–G mean iteration is not self-correcting. Thus, we cannot start with low precision and increase it, as was possible in Section 2.

2. Since there are ∼ log_{2} n iterations, we may lose O(log log(n)) bits of accuracy through accumulation of rounding errors, even though the algorithm is numerically stable. Thus, it may be necessary to work with precision n+O(log log(n)). From (1.3), the time required is still ∼ \frac{15}{2}M(n) log_{2} n.

9 Multiple-precision Evaluation of log(x)

There are several algorithms for evaluating log(x) to precision n in time O(M(n) log(n)). For example, a method based on Landen transformations of incomplete elliptic integrals is described in Brent [75b]. The method described here is essentially due to Salamin (see Beeler et al [72]), though the basic relation (9.1) was known by Gauss.

If cos(ϕ) = ε^{\frac{1}{2}} is small, then

(9.1)  K(ϕ) = (1 + O(ε)) log(4ε^{− \frac{1}{2}} )

Thus, taking a_{0} = 1, b_{0} = 4/y, where y = 4ε^{− \frac{1}{2}} , and applying the A–G mean iteration to compute a = lim_{i→∞}a_{i}, gives

log(y) = \frac{π}{2a}(1 + O(y^{−2}))

for large y. Thus, so long as y ≥ 2^{n/2}, we can evaluate log(y) to precision n. If log(y) = O(n) then ∼ 2 log_{2} n iterations are required, so the time is ∼ 13M(n) log_{2} n, assuming π is precomputed.

For example, to find log(10^{6}) we start the A–G mean iteration with a_{0} = 1 and b_{0} = 4 × 10^{−6}. Results of the first seven iterations are given to 10 significant figures in Table 9.1. We find that π/(2a_{7}) = 13.81551056, which is correct.

Table 9.1: Computation of log(10^{6})

\begin{matrix}i & a_{i} & b_{i} \\ 0 & 1.000000000 & 4.000000000 × 10^{−6} \\ 1 & 5.000020000 × 10^{−1} & 2.000000000 × 10^{−3} \\ 2 & 2.510010000 × 10^{−1} & 3.162283985 × 10^{−2} \\ 3 & 1.413119199 × 10^{−1} & 8.909188753 × 10^{−2} \\ 4 & 1.152019037 × 10^{−1} & 1.122040359 × 10^{−1} \\ 5 & 1.137029698 × 10^{−1} & 1.136930893 × 10^{−1} \\ 6 & 1.136980295 × 10^{−1} & 1.136980294 × 10^{−1} \\ 7 & 1.136980295 × 10^{−1} & 1.136980295 × 10^{−1} \end{matrix}

Since log(2) = \frac{1}{n} log(2^{n}), we can evaluate log(2) to precision n in time ∼ 13M(n) log_{2} n. Suppose x ∈ [b, c], where b > 1. We may set y = 2^{n}x, evaluate log(y) as above, and use the identity

log(x) = log(y) − n log(2)

to evaluate log(x). Since log(y) ≃ n log(2), approximately log_{2} n significant bits will be lost through cancellation, so it is necessary to work with precision n + O(log(n)).

If x is very close to 1, we have to be careful in order to obtain log(x) with a small relative error. Suppose x = 1 + δ. If |δ| < 2^{−n/ log(n)} we may use the power series

log(1 + δ) = δ − δ^{2}/2 + δ^{3}/3 − …,

and it is sufficient to take about log(n) terms. If δ is larger, we may use the above A–G mean method, with working precision n + O(n/ log(n)) to compensate for any cancellation.

Finally, if 0 < x < 1, we may use log(x) = − log(1/x), where log(1/x) is computed as above. To summarize, we have proved:

Theorem 9.1

If x > 0 is a precision n number, then log(x) may be evaluated to precision n in time ∼ 13M(n) log_{2} n as n → ∞ [ assuming π and log(2) precomputed to precision n + O(n/ log(n))].

Note

The time required to compute log(x) by the obvious power series method is O(nM(n)). Since 13 log_{2} n < n for n ≥ 83, the method described here may be useful for moderate n, even if the classical O(n^{2}) multiplication algorithm is used.

10 Multiple-precision Evaluation of exp(x)

Corresponding to Theorem 9.1 we have:

Theorem 10.1

If [a, b] is a fixed interval, and x ∈ [a, b] is a precision n number such that exp(x) does not underflow or overflow, the exp(x) can be evaluated to precision n in time ∼ 13M(n) log_{2} n as n → ∞ (assuming π and log(2) are precomputed).

Proof

To evaluate exp(x) we need to solve the equation f(y) = 0, where f(y) = log(y) − x, and x is regarded as constant. Since

f^{(k)}(y) = (−1)^{k−1} (k − 1) ! y^{−k}

can be evaluated in time O(M(n)) = o(M(n) log(n)) for any fixed k ≥ 1, the result follows from Theorems 6.1 and 9.1. We remark that the (k + 1)-th order method in the proof of Theorem 6.1 may simply be taken as

y_{i+1} = y_{i} ∑_{j=0}^{k}(x − log(y_{i}))^{j}/j!

11 Multiple-precision Operations on Complex Numbers

Before considering the multiple-precision evaluation of trigonometric functions, we need to state some results on multiple-precision operations with complex numbers. We assume that a precision n complex number z = x + iy is represented as a pair (x, y) of precision n real numbers. As before, a precision n operation is one which gives a result with a relative error O(2^{−n}). (Now, of course, the relative error may be complex, but its absolute value must be O(2^{−n}).) Note that the smaller component of a complex result may occasionally have a large relative error, or even the wrong sign!

Complex Multiplication

Since z = (t + iu)(v + iw) = (tv − uw) + i(tw + uv), a complex multiplication may be done with four real multiplications and two additions. However, we may use an idea of Karatsuba and Ofman [62] to reduce the work required to three real multiplications and some additions: evalute tv, uw, and (t + u)(v + w), then use

tw + uv = (t + u)(v + w) − (tv + uw) .

Since |t + u| ≤ \sqrt{2} |t + iu| and |v + w| ≤ \sqrt{2} |v + iw|, we have

|(t + u)(v + w)| ≤ 2|z| .

Thus, all rounding errors are of order 2^{−n}|z| or less, and the computed product has a relative error O(2^{−n}). The time for the six additions is asymptotically negligible compared to that for the three multiplications, so precision n complex multiplication may be performed in time ∼ 3M(n).

Complex Squares

Since (v + iw)^{2} = (v − w)(v + w) + 2ivw, a complex square may be evaluated with two real multiplications and additions, in time ∼ 2M(n).

Complex Division

Using complex multiplication as above, and the same division algorithm as in the real case, we can perform complex division in time ∼ 12M(n). However, it is faster to use the identity

\frac{t + iu}{v + iw} = (v^{2} + w^{2})^{−1}[(t + iu)(v − iw)] ,

reducing the problem to one complex multiplication, four real multiplications, one real reciprocal, and some additions. This gives time ∼ 10M(n). For complex reciprocals we have t = 1, u = 0, and time ∼ 7M(n).

Complex Square Roots

Using (2.2) requires, at the last iteration, one precision n complex squaring and one precision n/2 complex division. Thus, the time required is ∼ 2(2 + 10/2)M(n) = 14M(n).

Complex A–G Mean Iteration

From the above results, a complex square root and multiplication may be performed in time ∼ 17M(n). Each iteration transforms two points in the complex plane into two new points, and has an interesting geometric interpretation.

12 Multiple-precision Evaluation of Trigonometric Functions

Since

(12.1)  log(v + iw) = log |v + iw| + i.artan(w/v)

and

(12.2)  exp(iθ) = cos(θ) + i. sin(θ) ,

we can evaluate artan(x), cos(x) and sin(x) if we can evaluate log(z) and exp(z) for complex arguments z. This may be done just as described above for real z, provide we choose the correct value of \sqrt{a_{j}b_{j}}. Some care is necessary to avoid excessive cancellation; for example, we should use the power series for sin(x) if |x| is very small, as described above for log(1 + δ). Since ∼ 2 log_{2} n A–G mean iterations are required to evaluate log(z), and each iteration requires time ∼ 17M(n), we can evaluate log(z) in time ∼ 34M(n) log_{2} n. From the complex version of Theorem 6.1, exp(z) may also be evaluated in time ∼ 34M(n) log_{2} n.

As an example, consider the evaluation of log(z) for z = 10^{6}(2 + i). The A–G mean iteration is started with a_{0} = 1 and b_{0} = 4/z = 1.6 × 10^{−6} − (8.0 × 10^{−7})i. The results of six iterations are given, to 8 significant figures, in Table 12.1.

Table 12.1: Evaluation of log 10^{6}(2 + i)

j  a_{j}  b_{j}
0  (1.0000000,  0.0000000)  (1.6000000 × 10^{−6}, −8.0000000 × 10^{−7})
1  (5.0000080 × 10^{−1}, −4.0000000 × 10^{−7})  (1.3017017 × 10^{−3}, −3.0729008 × 10^{−4})
2  (2.5065125 × 10^{−1}, −1.5384504 × 10^{−4})  (2.5686505 × 10^{−2}, −2.9907884 × 10^{−3})
3  (1.3816888 × 10^{−1}, −1.5723167 × 10^{−3})  (8.0373334 × 10^{−2}, −4.6881008 × 10^{−3})
4  (1.0927111 × 10^{−1}, −3.1302088 × 10^{−3})  (1.0540970 × 10^{−1}, −3.6719673 × 10^{−3})
5  (1.0734040 × 10^{−1}, −3.4010880 × 10^{−3})  (1.0732355 × 10^{−1}, −3.4064951 × 10^{−3})
6  (1.0733198 × 10^{−1}, −3.4037916 × 10^{−3})  (1.0733198 × 10^{−1}, −3.4037918 × 10^{−3})

We find that

\frac{π}{2a_{7}} = 14.620230 + 0.46364761i
≃ log |z| + i.artan(\frac{1}{2})

as expected.

Another method for evaluating trigonometric functions in time O(M(n) log(n)), without using the identities (12.1) and (12.2) is described in Brent [75b].

13 Operations on Formal Power Series

There is an obvious similarity between a multiple-precision number with base β:

β^{e} ∑_{i=1}^{n} a_{i}β^{−i} (0 ≤ a_{i} < β) ,

and a formal power series:

∑_{i=0}^{∞} a_{i}x^{i} (a_{i} real, x an indeterminate) .

Thus, it is not surprising that algorithms similar to those described in Section 2 may be used to perform operations on power series.

In this section only, M(n) denotes the number of scalar operations required to evaluate the first n coefficients c_{0}, …, c_{n−1} in the formal product

(∑_{i=0}^{∞} a_{i}x^{i}) (∑_{i=0}^{∞} b_{i}x^{i}) = ∑_{i=0}^{∞} c_{i}x^{i} .

Clearly, c_{j} depends only on a_{0}, …, a_{j} and b_{0}, …, b_{j}, in fact

c_{j} = ∑_{i=0}^{j} a_{i}b_{j−i} .

The classical algorithm gives M(n) = O(n^{2}), but it is possible to use the fast Fourier transform (FFT) to obtain

M(n) = O(n log(n))

(see Borodin [73]).

If we assume that M(n) satisfies conditions (1.2) and (1.3), then the time bounds given in Section 2 for division, square roots, etc. of multiple-precision numbers also apply for the corresponding operations on power series (where we want the first n terms in the result). For example,

if P (x) = ∑_{i=0}^{∞}a_{i}x^{i} and a_{0} ≠ 0, then the first n terms in the expansion of 1/P (x) may be found

with ∼ 3M(n) operations as n → ∞. However, some operations, e.g. computing exponentials, are much easier for power series than for multiple-precision numbers!

Evaluation of log(P (x))

If a_{0} > 0 we may want to compute the first n terms in the power series Q(x) = log(P (x)). Since Q(x) = log(a_{0}) + log(P (x)/a_{0}), there is no loss of generality in assuming that a_{0} = 1. Suppose

Q(x) = ∑_{i=0}^{∞}b_{i}x^{i}. From the relation

(13.1)  Q′(x) = P′(x)/P (x) ,

where the prime denotes formal differentiation with respect to x, we have

(13.2)  ∑_{i=1}^{∞} ib_{i}x^{i−1} = (∑_{i=1}^{∞} ia_{i}x^{i−1}) [figure](∑_{i=0}^{∞} a_{i}x^{i}) .

The first n terms in the power series for the right side of (13.2) may be evaluated with ∼ 4M(n) operations, and then we need only compare coefficients to find b, …, b_{n−1}. (Since a_{0} = 1, we know that b_{0} = 0.) Thus, the first n terms in log(P (x)) may be found in ∼ 4M(n) operations. It is interesting to compare this result with Theorem 9.1.

Evaluation of exp(P (x))

If R(x) = exp(P (x)) then R(x) = exp(a_{0}) exp(P (x) − a_{0}), so there is no loss of generality in assuming that a_{0} = 0. Now log(R(x)) − P (x) = 0, and we may regard this as an equation for the unknown power series R(x), and solve it by one of the usual iterative methods. For example, Newton’s method gives the iteration

(13.3)  R_{i+1}(x) = R_{i}(x) − R_{i}(x)(log(R_{i}(x)) − P (x)) .

If we use the starting approximation R_{0}(x) = 1, then the terms in R_{k}(x) agree exactly with those in R(x) up to (but excluding) the term O(x^{2^{k}}). Thus, using (13.3), we can find the first n terms of exp(P (x)) in ∼ 9M(n) operations, and it is possible to reduce this to ∼ \frac{22}{3}M(n) operations by using a fourth-order method instead of (13.3). Compare Theorem 10.1.

Evaluation of P^{m}

Suppose we want to evaluate (P (x))^{m} for some large positive integer m. We can assume that a_{0} ≠ 0, for otherwise some power of x may be factored out. Also, since P^{m} = a_{0}^{m}(P/a_{0})^{m}, we can assume that a_{0} = 1. By forming P^{2}, P^{4}, P^{8}, …, and then the appropriate product given by the binary expansion of m, we can find the first n terms of P^{m} in O(M(n) log_{2} m) operations. Surprisingly, this is not the best possible result, at least for large m. From the identity

(13.4)  P^{m} = exp(m log(P))

and the above results, we can find the first n terms of P^{m} in O(M(n)) operations! (If a_{0} ≠ 1, we also need O(log_{2} m) operations to evaluate a_{0}^{m}.) If the methods described above are used to compute the exponential and logarithm in (13.4), then the number of operations is ∼ \frac{34}{3}M(n) as n → ∞.

Other operations on power series

The method used to evaluate log(P (x)) can easily be generalised to give a method for f(P (x)), where df(t)/dt is a function of t which may be written in terms of square roots, reciprocals etc. For example, with f(t) = artan(t) we have df/dt = 1/(1 + t^{2}), so it is easy to evaluate artan(P (x)). Using Newton’s method we can evaluate the inverse function f^{(−1)}(P (x)) if f(P (x)) can be evaluated. Generalizations and applications are given in Brent and Kung [75].

Some operations on formal power series do not correspond to natural operations on multipleprecision numbers. One example, already mentioned above, is formal differentiation. Other interesting examples are composition and reversion. The classical composition and reversion algorithms, as given in Knuth [69], are O(n^{3}), but much faster algorithms exist: see Brent and Kung [75].

References

Beeler, Gosper and Schroeppel [72] Beeler, M., Gosper, R.W., and Schroeppel, R. “Hakmem”. Memo No. 239, M.I.T. Artificial Intelligence Lab., 1972, 70–71.

Borodin [73] Borodin, A., “On the number of arithmetics required to compute certain functions – circa May 1973”. In Complexity of Sequential and Parallel Numerical Algorithms (ed. by J.F. Traub), Academic Press, New York, 1973, 149–180.

Brent [75a] Brent, R.P., “The complexity of multiple-precision arithmetic”. Proc. Seminar on Complexity of Computational Problem Solving (held at the Australian National University, Dec. 1974), Queensland Univ. Press, Brisbane, 1975.

Brent [75b] Brent, R.P., “Fast multiple-precision evaluation of elementary functions”. Submitted to J. ACM.

Brent and Kung [75] Brent, R.P. and Kung, H.T., “Fast algorithms for reversion and composition of power series”. To appear. (A preliminary paper appears in these Proceedings.)

Gauss [1876] Gauss, C.F., “Carl Friedrich Gauss Werke”, (Ed. 3), Göttingen, 1876, 362–403.

Katatsuba and Ofman [62] Karatsuba, A. and Ofman, Y., “Multiplication of multidigit numbers on automata”, (in Russian). Dokl. Akad. Nauk SSSR 146 (1962), 293–294.

Knuth [69] Knuth, D.E., “The Art of Computer Programming”, (Vol. 2), Addison Wesley, Reading, Mass., 1969, Sec. 4.7.

Legendre [11] Legendre, A.M., “Exercices de Calcul Integral”, (Vol. 1), Paris, 1811, 61.

Melzak [73] Melzak, Z.A., “Companion to Concrete Mathematics”, Wiley, New York, 1973, 68–69.

Salamin [75] Salamin, E., “A fast algorithm for the computation of π”. To appear in Math. Comp.

Schönhage and Strassen [71] Schönhage, A. and Strassen, V., “Schnelle Multiplikation grosser Zahlen”. Computing 7 (1971), 281–292.

Shanks and Wrench [62] Shanks, D. and Wrench, J.W., “Calculation of π to 100,000 decimals”. Math. Comp. 16 (1962), 76–99.

Postscript (September 1999)

Historical Notes and References

This paper was retyped in L^{A}TEX with minor corrections in September 1999. It is available electronically in compressed postscript format from http://wwwmaths.anu.edu.au/~brent/

pub/pub028.html

Brent [75a] is available electronically in compressed postscript format from

http://wwwmaths.anu.edu.au/~brent/pub/pub032.html

Brent [75b] appeared in Journal of the ACM 23 (1976), 242–251. See

http://theory.lcs.mit.edu/~jacm/jacm76.html#Brent1976:242

The “preliminary paper” Brent and Kung [75] appeared as “O((n log n)^{3/2}) algorithms for composition and reversion of power series” in Analytic Computational Complexity (edited by J. F. Traub), Academic Press, New York, 1975, 217–225. The final paper appeared as “Fast algorithms for manipulating formal power series”, Journal of the ACM 25 (1978), 581–595. See

http://theory.lcs.mit.edu/~jacm/jacm78.html#BrentK1978:581

A generalisation to multivariate power series appeared as R. P. Brent and H. T. Kung, “Fast algorithms for composition and reversion of multivariate power series (preliminary version)”, in Proceedings of a Conference on Theoretical Computer Science held at the University of Waterloo, Dept. of Computer Science, University of Waterloo, Waterloo, Ontario (August 1977), 149–158. Abstract available electronically from http://wwwmaths.anu.edu.au/~brent/pub/

pub039.html

For more on (generalized) composition of power series, see R. P. Brent and J. F. Traub, “On the complexity of composition and generalized composition of power series”, SIAM J. Computing 9 (1980), 54–66. Abstract available electronically from http://wwwmaths.anu.edu.au/~brent/

pub/pub050.html

Hakmem by Beeler, Gosper and Schroeppel [72] is available electronically in various formats

from http://www.inwap.com/pdp10/hbaker/hakmem/hakmem.html

Salamin [75] appeared as “Computation of π using arithmetic-geometric mean” in Mathematics of Computation 30 (1976), 565–570.

For much more on the arithmetic-geometric mean, see J. M. Borwein and P. B. Borwein, Pi and the AGM, Wiley-Interscience, 1987.

Sharper Results

Some of the constants can be improved. For example, π can be computed in ∼ 6.25M(n) log_{2} n by the Gauss-Legendre method of Section 8, and the constant 13 in Theorem 9.1 can be replaced by 10.5. For more information see the postscript to Brent [75a], available electronically from

http://wwwmaths.anu.edu.au/~brent/pub/pub032.html