THE COMPLEXITY OF MULTIPLE-PRECISION

ARITHMETIC ^{1}

Richard P Brent

Computer Centre, Australian National University

In studying the complexity of iterative processes it is usually assumed that the arithmetic operations of addition, multiplication, and division can be performed in certain constant times. This assumption is invalid if the precision required increases as the computation proceeds. We give upper and lower bounds on the number of singleprecision operations required to perform various multiple-precision operations, and deduce some interesting consequences concerning the relative efficiencies of methods for solving nonlinear equations using variable-length multiple-precision arithmetic.

1 Introduction

Traub [28] defines analytic computational complexity to be the optimality theory of analytic or continuous processes. Apart from some work by Schultz [24] on differential equations, most recent results have concerned iterative methods for the solution of nonlinear equations or systems of equations. See, for example, Brent [1,3,6], Brent, Winograd and Wolfe [7], Kung [14, 15], Kung and Traub [16–17], Paterson [21], Rissanen [22], Traub [28–32] and Wozniakowski [35, 36].

The authors just cited make the (usually implicit) assumption that arithmetic is performed with a fixed precision throughout a given computation. This is probably true for most computations programmed in Fortran or Algol 60. Suppose, though, that we are concerned with an iterative process for approximating an irrational number ζ (for example, \sqrt{2}, π or e) to arbitrary accuracy. The iterative process should (theoretically) generate a sequence (x_{i}) of real numbers, such that ζ = lim_{i→∞} x_{i}, provided no rounding errors occur. On a computing machine each x_{i} has to be approximated by a finite-precision machine-representable number x̃_{i}, and ζ = lim_{i→∞} x̃_{i} can only hold if the precision increases indefinitely as i → ∞. In practice, only a finite number of members of the sequence (x̃_{i}) will ever be generated, but if an accurate approximation to ζ is required it may be possible to save a large amount of computational work by using variable precision throughout the computation. This is likely to become easier to program as new languages (and possibly hardware), which allow the precision of floating-point numbers to be varied dynamically, are developed.

1First appeared in The Complexity of Computational Problem Solving (edited by R S Anderssen and R P Brent), Univ. of Queensland Press, 1976, 126–165. Retyped with minor corrections by Frances Page at Oxford University Computing Laboratory, 1999.

Copyright \begin{matrix}c \\ ◯ \end{matrix} 1976, 1999, R. P. Brent and University of Queensland Press. rpb032 typeset using L^{A}TEX.

In Section 7 we discuss the effect of using variable precision when solving nonlinear equations. Before doing so, we consider the complexity of the basic multiple-precision arithmetic operations. We assume that a standard floating-point number representation is used, with a binary fraction of n bits. (Similar results apply for any fixed base, for example, 10.) We are interested in the case where n is much greater than the wordlength of the machine, so the fraction occupies several words. For simplicity, we assume that the exponent field has a fixed length and that numbers remain in the allowable range, so problems of exponent overflow and underflow may be neglected. Note that our assumptions rule out exotic number representations (for example, logarithmic [4] or modular [33, 34] representations) in which it is possible to perform some (but probably not all) of the basic operations faster than with the standard representation. To rule out “table-lookup” methods, we assume that a random-access memory of bounded size and a bounded number of sequential tape units are available. (Formally, our results apply to multitape Turing machines.)

In Sections 2 to 6 we ignore “constant” factors, that is factors which are bounded as n → ∞. Although the constant factors are of practical importance, they depend on the computer and implementation as well as on details of the analysis. Certain machine-independent constants are studied in Sections 7 and 8.

If B is a multiple-precision operation, with operands and result represented as above (that is, “precision n” numbers), then t_{n}(B) denotes the worst-case time required to perform B, obtaining the result with a relative error at most 2^{−n}c, where c is independent of n. We assume that the computation is performed on a serial machine whose single-precision instructions have certain constant execution times. The following definition follows that in Hopcroft [11].

Definition 1.1 B is linearly reducible to C (written B =^{<} C), if there is a positive constant K such that

t_{n}(B) ≤ Kt_{n}(C)  (1.1)

for all sufficiently large n. B is linearly equivalent to C (written B ≡ C) if B =^{<} C and C =^{<} B.

In Section 2 we consider the complexity of multiple-precision addition and some linearly equivalent operations. Then, in Section 3, we show that multiple-precision division, computation of squares or square roots, and a few other operations are linearly equivalent to multiplication. Most of these results are well known [8, 9].

Sections 4 and 5 are concerned with the “operations” of evaluating exponentials, logarithms, and the standard trigonometric and hyperbolic functions (sin, artan, cosh, and so on). It turns out that most of (and probably all) these operations are linearly equivalent so long as certain restrictions are imposed.

Section 6 deals with the relationship between the four equivalence classes established in Sections 2 to 5, and several upper bounds on the complexity of operations in these classes are given. The best known constants relating operations which are linearly equivalent to multiplication are given in Section 7.

Finally, in Section 8, we compare the efficiencies of various methods for solving nonlinear equations using variable-length multiple-precision arithmetic. The relative efficiencies are different from those for the corresponding fixed-precision methods, and some of the conclusions may be rather surprising. The results of Sections 4 to 8 are mainly new.

In the analysis below, c_{1}, c_{2}, … denote certain positive constants which do not need to be specified further. The notation f ∼ g means that lim_{n→∞} f(n)/g(n) = 1, and f = O(g) means that |f(n)| ≤ Kg(n) for some constant K and all sufficiently large n. Finally, the abbreviation “mp” stands for “variable-length multiple-precision”.

2 Addition and linearly equivalent operations

Let A denote the operation of multiple-precision addition. Any reasonable implementation of floating-point addition, using at least one guard digit to avoid the possible occurrence of large relative errors, gives

t_{n}(A) ≤ c_{1}n .  (2.1)

Conversely, from the assumptions stated in Section 1, it is clear that

t_{n}(A) ≥ c_{2}n .  (2.2)

Hence, the complexity of multiple-precision addition is easily established. (For the operations discussed in Sections 3 and 5 the results are less trivial, in fact the conjectured lower bounds corresponding to (2.2) have not been proved rigorously.)

It is easy to see that bounds like (2.1) and (2.2) hold for multiple-precision subtraction, and multiplication or division of a multiple-precision number by a single-precision number (or even by any rational number with bounded numerator and denominator). Hence, all these operations are linearly equivalent to addition.

3 Multiplication and linearly equivalent operations

Let D, I, M, R and S denote the multiple-precision operations of division, taking reciprocals, multiplication, extraction of square roots and forming squares, respectively. In this section, we show that all these operations are linearly equivalent. The proofs are straightforward, but the result is surprising, as it seems intuitively obvious that taking a square root is inherently “more difficult” than forming a square, and similarly for division versus multiplication. (Some bounds on the relative difficulty of these operations are given in Section 7.)

Lemma 3.1

M =^{>} S =^{>} A .  (3.1)

Proof. Clearly

t_{n}(M) ≥ t_{n}(S) ≥ c_{3}n ,  (3.2)

so the result follows from (2.1).

Sharp upper bounds on t_{n}(M) are not needed in this section, so we defer them until Section 6. Lemmas 3.2 and 3.3, although weak, are sufficient for our present purposes.

Lemma 3.2 For all positive n,

t_{2n}(M) ≤ c_{4}t_{n}(M) .  (3.3)

Proof. First assume that n is divisible by 3, and consider operations on the n-bit fractions only. If we can multiply n-bit numbers with relative error 2^{−n}c_{0} then we can multiply n/3-bit numbers exactly (assuming 2^{n/3} > 2c_{0}). Thus, a 2n-bit fraction x may be split up into

x = λa + λ^{2}b + … + λ^{6}f ,  (3.4)

where λ = 2^{−n/3} and a, b, …, f are integers in [0, 2^{n/3}), and the product of two such 2n-bit fractions may be formed exactly with 36 exact multiplications of n/3-bit numbers and some additions. Thus

t_{2n}(M) ≤ 36t_{n}(M) + c_{5}t_{2n}(A) ,  (3.5)

and the result follows from Lemma 3.1. Trivial modifications to the above proof suffice, if n is not divisible by 3.

Lemma 3.3 For some constant c_{6} < 1,

t_{n}(M) ≤ c_{6}t_{8n}(M)  (3.6)

for all sufficiently large n.

Proof. If a, b, c and d are integers in [0, 2^{n}), the identity

(a + λb)(c + λd) = ac + λ(bc + ad) + λ^{2}bd ,  (3.7)

with λ = 2^{3n}, may be used to obtain the products ac and bd from one 8n-bit product. Thus

2t_{n}(M) ≤ t_{8n}(M) + c_{7} .  (3.8)

The result (with c_{6} = 3/4) follows if n is sufficiently large that t_{8n}(M) ≥ 2c_{7}. (We have assumed that the time required for one n-bit multiplication is half the time required for two independent n-bit multiplications, but much weaker assumptions would be sufficient.)

The following lemma will be used to estimate the work required for multiple-precision divisions and square roots.

Lemma 3.4 Given α ∈ (0, 1), there is a constant c_{8} such that, for any integers n_{0}, …, n_{p} satisfying

1 ≤ n_{j} ≤ α^{j}n  (3.9)

for j = 0, 1, …, p, we have

∑_{j=0}^{p} t_{n_{j}}(M) ≤ c_{8}t_{n}(M) .  (3.10)

Proof. Let k be large enough that

α^{k} ≤ 1/8 .  (3.11)

From (3.9) and (3.11),

t_{n_{jk}}(M) ≤ α_{6}^{j}t_{n}(M)  (3.12)

for j = 0, 1, …, ⌊p/k⌋, provided n_{jk} is sufficiently large for Lemma 3.3 to be applicable. Thus,

∑_{j=0}^{p} t_{n_{j}}(M) ≤ kt_{n}(M) (1 + c_{6} + c_{6}^{2} + …) + c_{7} ,  (3.13)

where the term c_{7} allows for those t_{n_{j}}(M) for which Lemma 3.3 is not applicable. If

c_{8} = k/ (1 − c_{6}) + c_{7} ,  (3.14)

the result follows from (3.13).

The following lemma shows that multiple-precision multiplication is linearly equivalent to squaring. This result is essentially due to Floyd [9].

Lemma 3.5

M ≡ S .  (3.15)

Proof. Since squaring is a special case of multiplication,

M =^{>} S .  (3.16)

Conversely, we may use the identity

4λab = (a + λb)^{2} − (a − λb)^{2} ,  (3.17)

where λ is a power of 2 chosen so that

\frac{1}{2} ≤ |λb/a| ≤ 2  (3.18)

(unless a = 0 or b = 0). This scaling is necessary to avoid excessive cancellation in (3.17). (A detailed discussion of a similar situation is given in Brent [5].) From (3.17),

t_{n}(M) ≤ 2t_{n}(S) + 3t_{n}(A) + c_{9},  (3.19)

so M =^{<} S follows from Lemma 3.1.

The next two lemmas show that multiple-precision multiplication is linearly equivalent to taking reciprocals and to division. The idea of the proof of Lemma 3.6 is to use a Newton iteration involving only multiplications and additions to approximate 1/a. Computational work is saved by starting with low precision and approximately doubling the precision at each iteration. The basic idea is well-known and has even been implemented in hardware.

The possibility of saving work by increasing the precision at each iteration is examined more closely in Sections 7 and 8.

Lemma 3.6

I =^{<} D =^{<} M .  (3.20)

Proof. Consider the iteration

x_{j+1} = x_{j} (2 − ax_{j})  (3.21)

obtained by applying Newton’s method to the equation x^{−1} − a = 0. If

x_{j} = (1 − ε_{j}) a^{−1} ,  (3.22)

then substitution in (3.21) shows that

ε_{j+1} = ε_{j}^{2} ,  (3.23)

so the order of convergence is two. A single-precision computation is sufficient to give an initial approximation such that |ε_{0}| ≤ \frac{1}{2}, and it follows from (3.23) that

|ε_{j}| ≤ 2^{−2^{j}}  (3.24)

for all j ≥ 0.

In deriving (3.24) we have assumed that (3.21) is satisfied exactly, but a result like (3.24) holds so long as the right hand side of (3.21) is evaluated using a precision of at least 2^{j+1} bits. Thus, an n-bit approximation to a^{−1} can be obtained by performing ⌈log_{2} n⌉ iterations of (3.22) with precision at least 2, 2^{2}, 2^{3}, …, 2^{⌈log_{2} n⌉−1}, n at each iteration. From Lemma 3.4 (with α = \frac{1}{2}), this gives

t_{n}(I) ≤ c_{10}t_{n}(M) .  (3.25)

Since b/a = b(1/a), it follows that

t_{n}(D) ≤ c_{11}t_{n}(M) ,  (3.26)

so D =^{<} M. Since I =^{<} D is trivial, the proof is complete.

From ab = a/(1/b) it is clear that M =^{<} D. The proof that M =^{<} I is not quite so obvious, and uses the equivalences of multiplication and squaring (Lemma 3.5).

Lemma 3.7

M =^{<} I .  (3.27)

Proof. We may apply the identity

a^{2}(1 − λa)^{−1} = λ^{−2} [(1 − λa)^{−1} − (1 + λa)]  (3.28)

to obtain an approximation to a^{2}, using only the operation of taking reciprocals, addition (or subtraction) and multiplication by powers of two. If a ≠ 0, choose λ to be a power of two such that

2^{−n/3−1} < |λa| < 2^{1−n/3} ,  (3.29)

and evaluate the right hand side of (3.28), using precision n. This gives an approximation to a^{2} with precision ⌈n/3⌉, so

S_{⌈n/3⌉} =^{<} I_{n} ,  (3.30)

where the subscripts denote the precision. Thus, the result follows from Lemmas 3.2 and 3.5.

To conclude this section we consider the complexity of multiple-precision square roots. Results like Lemmas 3.8 and 3.9 actually hold if x^{\frac{1}{2}} is replaced by x^{p} for any fixed rational p ≠ 0 or 1 (we have already shown this for p = −1).

Lemma 3.8

M =^{<} R .  (3.31)

Proof. The proof is similar to that of Lemma 3.7, using the approximation 2λ^{−2} [1 + λa − (1 + 2λa)^{\frac{1}{2}} ] to a^{2}.

Lemma 3.9

R =^{<} M .  (3.32)

Proof. The proof is similar to that of Lemma 3.6, using Newton’s iteration

x_{j+1} = \frac{1}{2} (x_{j} + a/x_{j}) ,  (3.33)

with precision increasing at each iteration, to approximate \sqrt{a}. Alternatively, it is possible to avoid multiple-precision division by using the iteration

x_{j+1} = x_{j} (3 − ax_{j}^{2}) /2  (3.34)

to approximate a^{− \frac{1}{2}} , and then use \sqrt{a} = a.a^{− \frac{1}{2}} to evaluate \sqrt{a}.

The results of Lemmas 3.5 to 3.9 may be summarized in the following:

Theorem 3.1

D ≡ I ≡ M ≡ R ≡ S .  (3.35)

4 Some regularity conditions

Before discussing the complexity of multiple-precision evaluation of exponentials, trigonometric functions, etc., we need some definitions. Throughout this section, let ϕ(x) be a real-valued function which is positive and monotonic increasing for all sufficiently large positive x.

Definition 4.1 ϕ ∈ Φ_{1} iff, for all α ∈ (0, 1), for some positive K, for all sufficiently large x and all x_{0}, …, x_{J} satisfying

1 ≤ x_{j} ≤ α^{j}x  (4.1)

for j = 0, …, J, we have

∑_{j=0}^{J}ϕ(x_{j}) ≤ Kϕ(x) .  (4.2)

ϕ ∈ Φ_{2} iff, for some α, β ∈ (0, 1) and all sufficiently large x,

ϕ(αx) ≤ βϕ(x) .  (4.3)

ϕ ∈ Φ_{3} iff, for some positive K_{1}, K_{2} and p, there is a monotonic increasing function ψ such that

K_{1}x^{p}ψ(x) ≤ ϕ(x) ≤ K_{2}x^{p}ψ(x)  (4.4)

for all sufficiently large x.

Note the similarity between the definition of Φ_{1} and the statement of Lemma 3.4. In Section 5, we need to assume that the time ϕ(n) required to perform certain operations with precision n satisfies (4.2). The following lemmas make this assumption highly plausible. Lemma 4.1 shows that “for all α” in the definition of Φ_{1} may be replaced by “for some α”.

Lemma 4.1 If, for some α ∈ (0, 1) and some positive K, for all sufficiently large x and all x_{0}, …, x_{J} satisfying (4.1), we have (4.2), then ϕ ∈ Φ_{1}.

Proof. Take any α_{1} and α_{2} in (0, 1), and suppose that (4.1) with α replaced by α_{2} implies (4.2) with K replaced by K_{2}. Let m be a positive integer such that α_{1}^{m} ≤ α_{2}. If (4.1) holds with α replaced by α_{1} for a sequence (x_{0}, x_{1}, …, x_{J}), then (4.1) also holds with α replaced by α_{2} for each of the m subsequences

(x_{0}, x_{m}, …), (x_{1}, x_{m+1}, …), …, (x_{m−1}, x_{2m−1}, …) ,

so (4.2) holds with K replaced by K_{1} = mK_{2}.

Lemmas 4.2 and 4.3 show that ϕ ∈ Φ_{2} or ϕ ∈ Φ_{3} is a sufficient condition for ϕ ∈ Φ_{1}. The proof of Lemma 4.2 is similar to that of Lemma 3.4 (using Lemma 4.1), so is omitted.

Lemma 4.2

Φ_{2} ⊆ Φ_{1} .  (4.5)

Lemma 4.3

Φ_{2} = Φ_{3} .  (4.6)

Proof. First suppose that ϕ ∈ Φ_{3}, so (4.4) holds for some function ψ and some positive K_{1}, K_{2} and p. Choose α ∈ (0, 1) such that β = α^{p}K_{2}/K_{1} < 1. For all sufficiently large x, we have

ϕ(αx) ≤ K_{2}α^{p}x^{p}ψ(αx) ≤ K_{2}α^{p}x^{p}ψ(x) ≤ (K_{2}α^{p}/K_{1}) ϕ(x) ≤ βϕ(x)  (4.7)

(using (4.4) and the monotonicity of ψ), so ϕ ∈ Φ_{2}.

Conversely, suppose that ϕ ∈ Φ_{2}, so (4.3) holds for all sufficiently large x (say x ≥ x_{0} > 0) and some α, β ∈ (0, 1). Choose p small enough that β ≤ α^{p}, so

ϕ(αx) ≤ α^{p}ϕ(x)  (4.8)

for x ≥ x_{0}. Since ϕ(x) is positive for sufficiently large x, we may assume that ϕ(x_{0}) > 0. Let K_{1} = α^{p}, K_{2} = 1, and

ψ(x) = sup_{x_{0}≤y≤x} ϕ(y)/y^{p}  (4.9)

for x ≥ x_{0}. Thus, ψ(x) is monotonic increasing and

ψ(x) ≥ ϕ(x)/x^{p}  (4.10)

so

ϕ(x) ≤ K_{2}x^{p}ψ(x)  (4.11)

for x ≥ x_{0}.

By repeated application of (4.8) we have, for k ≥ 0,

ϕ(x)/x^{p} ≥ ϕ(α^{k}x)/(α^{k}x)^{p}  (4.12)

provided α^{k}x ≥ x_{0}. Thus, from (4.9),

ψ(x) = sup_{αx≤y≤x} ϕ(y)/y^{p}  (4.13)

≤ ϕ(x)/(αx)^{p}  (4.14)

for x ≥ x_{0}/α. Thus,

ϕ(x) ≥ K_{1}x^{p}ψ(x)  (4.15)

and, in view of (4.11), ϕ ∈ Φ_{2}.

5 Linear equivalence of various elementary functions

In this section, we consider the multiple-precision “operations” of evaluating certain elementary functions (log, exp, sin, artan, etc). First we prove three theorems which apply under fairly general conditions. Theorem 5.1 is a generalization of Lemmas 3.7 and 3.8, and gives a simple condition under which the evaluation of f(x) is at least as difficult as a multiplication (in the sense of Definition 1.1).

NOTATION. If f is a real-valued function defined on some finite interval [a, b], the operation of evaluating f(x) to (relative) precision n for x ∈ [a, b] is denoted by E_{[a,b]}^{(n)}(f). If there is no risk of confusion, we write simply E_{[a,b]}(f) or E(f). We sometimes write t_{n}(f) for t_{n} (E(f)). LC^{(m)}[a, b] is the class of functions with Lipschitz continuous m-th derivatives on [a, b]. We always assume that b > a.

Theorem 5.1 If f ∈ LC^{(2)}[a, b] and there is a point x_{0} ∈ (a, b) such that f′′(x_{0}) ≠ 0, then

E(f) ≥ M .  (5.1)

Proof. For all sufficiently small h, we have (from [2, Lemma 3.2])

f(x_{0} + h) + f(x_{0} − h) − 2f(x_{0}) = h^{2}f′′(x_{0}) + R(x) ,  (5.2)

where

|R(x)| ≤ c_{12}|h|^{3} .  (5.3)

Let c = f′′(x_{0}) ≠ 0. Three evaluations of f and some additions may be used to approximate ch^{2}, using (5.2). If h is of order 2^{−n/3}, the resulting approximation to ch^{2} has relative error of order 2^{−n/3}. Proceeding as in the proof of Lemma 3.5, we see that six evaluations of f and some additions may be used to approximate cxy to precision n/3, for any x and y. Applying this result, with x replaced by the stored constant c^{−2}, and y replaced by the computed cxy, shows that 12 evaluations of f give c (c^{−2}) (cxy) = xy to precision n/3. The result now follows from Lemma 3.2.

REMARK. If f′′(x) is not constant on [a, b], the point x_{0} may be chosen so that f′′(x) is rational, so (5.2) may be used to approximate h^{2}, and the result follows more easily (as in the proof of Lemma 3.7).

Theorem 5.2 gives conditions under which the multiple-precision evaluation of the inverse function g = f^{(−1)} of a function f is linearly reducible to the evaluation of f. (The inverse function satisfies g (f(x)) = x.) The condition 0 ∉ [a, b] could be dropped, if we only required the computation of g with an absolute (rather than relative) error of order 2^{−n}.

Theorem 5.2 If 0 ∉ [a, b], f ∈ LC^{(1)}[a, b], f′(x) ≠ 0 on [a, b], E(f) =^{>} M, and

t_{n}(f) ∈ Φ_{1} ,  (5.4)

then

E(g) =^{<} E(f) ,  (5.5)

where g = f^{(−1)} and Φ_{1} is as in Definition 4.1.

Proof. Since f′(x) is continuous and nonzero on [a, b], there is no loss of generality in assuming that

f′(x) ≥ c_{13} > 0  (5.6)

on [a, b]. Thus, g(y) exists on [c, d] = [f(a), f(b)]. Also, since 0 ≠ [a, b], we have

|g(y)| ≥ c_{14} > 0  (5.7)

on [c, d].

To estimate g(y) we may solve ψ(x) = 0 by a discrete version of Newton’s method, where

ψ(x) = f(x) − y .  (5.8)

Consider the iteration

x_{j+1} = x_{j} − ψ(x_{j})/μ_{j} ,  (5.9)

where

μ_{j} = (ψ(x_{j} + h_{j}) − ψ(x_{j})) /h_{j} ,  (5.10)

and the computation of μ_{j} and x_{j+1} is performed with precision n_{j} ≤ n, giving computed values μ̂_{j} and x̂_{j+1} respectively. If h_{j} is of order 2^{−n_{j}/2}, then

|μ̂_{j} − ψ′(x̂_{j})| ≤ 2^{−n_{j}/2}c_{15} ,  (5.11)

and it is easy to show that

|x̂_{j+1} − g(y)| ≤ |x̂_{j} − g(y)|^{2} c_{16} + 2^{−n_{j}/2} |x̂_{j} − g(y)| c_{17} + 2^{−n_{j}} c_{18} .  (5.12)

Since a sufficiently good starting approximation x_{0} may be found using single-precision (or at most bounded-precision) computation, (5.12) ensures that

|x̂_{j+1} − g(y)| ≤ |x̂_{j} − g(y)|^{2} c_{19} ,  (5.13)

provided

|x̂_{j} − g(y)| ≥ 2^{−n_{j}/2} .  (5.14)

Hence, we may approximately double the precision at each iteration, and (5.13) guarantees convergence of order two. A final iteration with h_{j} = 2^{−n/2} will be sufficient to give

|x̂_{j+1} − g(y)| ≤ 2^{−n}c_{20} .  (5.15)

Since E(f) ≥ M, the result follows from (5.4), (5.7), (5.15), and Lemma 3.6.

Theorem 5.3 If 0 ∉ [a, b], f ∈ LC^{(1)}[a, b], f(x)f′(x) ≠ 0 on [a, b], g = f^{(−1)}, E(f) =^{>} M, E(g) =^{>} M, t_{n}(f) ∈ Φ_{1}, and t_{n}(g) ∈ Φ_{1}, then

E(f) ≡ E(g) .  (5.16)

Proof. Since t_{n}(f) ∈ Φ_{1}, Theorem 5.2 applied to f gives E(g) =^{<} E(f). Similarly, applying Theorem 5.2 to f^{(−1)} gives E(f) =^{<} E(g), so the result follows.

We are now ready to deduce the linear equivalence of mp evaluation of various elementary functions f_{i}, assuming that t_{n}(f_{i}) ∈ Φ_{1}. In view of Lemmas 4.2 and 4.3, this assumption is very plausible.

Corollary 5.1 If 0 < a < b, c < d, 1 ∉ [a, b], t_{n} (E_{[a,b]}(log)) ∈ Φ_{1}, and t_{n} (E_{[c,d]}(exp)) ∈ Φ_{1}, then

E_{[a,b]}(log) ≡ E_{[c,d]}(exp) .  (5.17)

Proof. From Theorem 5.1, E_{[a,b]}(log) =^{>} M and E_{[c,d]}(exp) =^{>} M. Also, the identities

exp(−x) = 1/ exp(x)  (5.18)

and

exp(λx) = (exp(x))^{λ}  (5.19)

(for suitable rational λ) may be used to show that E_{[c,d]}(exp) ≡ E_{[c′,d′]}(exp) for any c′ < d′. Hence, the result follows from Theorem 5.3.

REMARK. If 1 ∈ [a, b], then Theorem 5.2 shows that

E_{[c,d]}^{(n)}(exp) =^{<} E_{[a,b]}^{(n)}(log) ,  (5.20)

and a proof like that of Theorem 5.2 shows that

E_{[a,b]}^{(n)}(log) =^{<} E_{[c,d]}^{(2n)}(exp) ,  (5.21)

so the conclusion of Corollary 5.1 follows, if

E_{[c,d]}^{(2n)}(exp) ≡ E_{[c,d]}^{(n)}(exp) .  (5.22)

Although (5.22) is plausible, no proof of it is known. (The corresponding result for multiplication is given in Lemma 3.2.)

Corollary 5.2

E(sinh) ≡ E(cosh) ≡ E(tanh) ≡ E(arsinh)
≡ E(arcosh) ≡ E(artanh) ≡ E(exp) ≡ E(log)
(5.23)

on any nontrivial closed intervals on which the respective functions are bounded and nonzero, assuming t_{n}(sinh) ∈ Φ_{1} etc.

Corollary 5.3

E(sin) ≡ E(cos) ≡ E(tan) ≡ E(arsin) ≡ E(arcos) ≡ E(artan)  (5.24)

on any nontrivial closed intervals on which the respective functions are bounded and nonzero, assuming t_{n}(sin) ∈ Φ_{1} etc.

REMARKS. The proofs of Corollaries 5.2 and 5.3 are similar to that of Corollary 5.1 (using well-known identities), so are omitted. Since exp(ix) = cos(x) + i sin(x), it is plausible that E(exp) ≡ E(sin), but we have not proved this. (It is just conceivable that the evaluation of exp(x) for complex x is not linearly reducible to the evaluation of exp(x) for real x.)

6 Upper and lower bounds

In this section we give some upper and lower bounds on t_{n}(A), t_{n}(M), t_{n}(exp) and t_{n}(sin). Since the multiplicative constants are not specified, the bounds apply equally well to the operations which are linearly equivalent to addition, multiplication, etc. (see Sections 2 to 5). The lower bounds are trivial: t_{n}\binom{exp}{sin} ≥ c_{21}t_{n}(M) ≥ c_{22}t_{n}(A) ≥ c_{23}n (from (2.2), Lemma 3.1 and Theorem 5.1). The upper bounds are more interesting.

UPPER BOUNDS ON t_{n}(M)

The obvious algorithm for multiplication of multiple-precision numbers gives

t_{n}(M) ≤ c_{24}n^{2} ,  (6.1)

but this is not the best possible upper bound. Karatsuba and Ofman [12] showed that

t_{n}(M) ≤ c_{25}n^{1.58...} ,  (6.2)

where 1.58 … = log_{2} 3. The idea of the proof is that, to compute

(a + λb)(c + λd) = ac + λ(ad + bc) + λ^{2}bd ,  (6.3)

where λ is a suitable power of two, we compute the three products m_{1} = ac, m_{2} = bd, and m_{3} = (a + b)(c + d), and use the identity

ad + bc = m_{3} − (m_{1} + m_{2}) .  (6.4)

Thus, 2n-bit integers can be multiplied with three multiplications of (at most) (n + 1)-bit integers, some multiplications by powers of two, and six additions of (at most) 4n-bit integers. This observation leads to a recurrence relation from which (6.2) follows.

More complicated identities like (6.4) may be used to reduce the exponent in (6.2). Recently Schönhage and Strassen [23] showed that the exponent can be taken arbitrarily close to unity. Their method gives the best known upper bound

t_{n}(M) ≤ c_{26}n log(n) log log(n) ,  (6.5)

and uses an alogrithm related to the fast Fourier transform to compute certain convolutions. For a description of this and earlier methods see Knuth [13 (revised)]. Knuth conjectures that (6.5) is optimal, though the term log log(n) is rather dubious. (It may be omitted if a machine with random-access memory of size O(n^{p}) for some fixed positive p is assumed.) From results of Morgenstern [19] and Cook and Aanderaa [8], it is extremely probable that

lim_{n→∞} t_{n}(M)/n = ∞ ,  (6.6)

which implies that M ≢ A, but more work remains to be done to establish this rigorously.

UPPER BOUNDS ON t_{n}(exp) AND t_{n}(sin)

To evaluate exp(x) to precision n from the power series

exp(±x) = ∑_{j=0}^{∞}(±x)^{j}/j! ,  (6.7)

it is sufficient to take c_{27}n/ log(n) terms, so

t_{n}(exp) ≤ c_{28}t_{n}(M)n/ log(n) .  (6.8)

Theorem 6.1 shows that the bound (6.8) may be reduced by a factor of order \sqrt{n}/ log(n).

Theorem 6.1

t_{n}(exp) ≤ c_{29}\sqrt{n} t_{n}(M)  (6.9)

and

t_{n}(sin) ≤ c_{30}\sqrt{n} t_{n}(M) .  (6.10)

Proof. To establish (6.9), we use the identity

exp(x) = (exp(x/λ))^{λ}  (6.11)

with λ = 2^{q}, where q = ⌊n^{\frac{1}{2}} ⌋. If [a, b] is the domain of x, and c = max(|a|, |b|), then

|(x/λ)^{r}/r!| ≤ 2^{−qr} ,  (6.12)

if r is large enough that

c^{r} ≤ r! .  (6.13)

Hence, it is sufficient to take r = ⌈n/q⌉ terms in the power series for exp(x/λ) to give an absolute error of order 2^{−n} in the approximation to exp(x/λ). Since exp(x/λ) is close to unity, the relative error will also be of the order 2^{−n} for large n. From (6.11), q squarings may be used to compute exp(x) once exp(x/λ) is known.

The method just described gives exp(x) to precision n − n^{\frac{1}{2}} , for the relative error in exp(x/λ) is amplified by the factor λ. This may be avoided by taking r = ⌈n/q⌉ + 1, and either working with precision n + n^{\frac{1}{2}} , or evaluating

exp(|x/λ|) − 1 ≃ ∑_{j=1}^{r} |x/λ|^{j}/j!  (6.14)

and then using the identity

(1 + ε)^{2} − 1 = 2ε + ε^{2}  (6.15)

to evaluate exp(|x|) − 1 without appreciable loss of significant figures. Thus, (6.9) follows (using Lemma 3.2 if necessary).

The proof of (6.10) is similar, using the identity

sin(x) = ±2 sin(x/2)\sqrt{1 − sin^{2}(x/2)} ,  (6.16)

q times to reduce the computation of sin(x) to that of sin(x/λ) (recall Lemma 3.9).

REMARKS. If x is a rational number with small numerator and denominator, the time required to sum r terms in the power series for exp(x/λ) is O(rn), and the time required for q squarings is O (qt_{n}(M)). Thus, choosing r = ⌊\sqrt{t_{n}(M)}⌋ and q = ⌈n/r⌉ gives total time O (n\sqrt{t_{n}(M)}). It is also possible to evaluate exp(x) in this time for general x, by using a form of preconditioning to reduce the number of multiplications required to evaluate the power series for exp(x/λ).

A NUMERICAL EXAMPLE

The following example illustrates the ideas of Theorem 6.1. Suppose we wish to calculate e to 30 decimal places. The obvious method is to use the approximation

e ≃ ∑_{j=0}^{28} 1/j!  (6.17)

(since 29! ≃ 8.8 × 10^{30}). On the other hand

e ≃ (∑_{j=0}^{10} \frac{1}{j! 256^{j}} )^{256}  (6.18)

also gives the desired accuracy (since 11! 256^{10} ≃ 4.8 × 10^{31}). Thus, the computation of 18 inverse factorials may be saved at the expense of 8 squarings.

Similarly, the computation of e to 10^{6} decimal places by the obvious method requires the sum of about 205,030 inverse factorials, but the approximation

e ≃ (∑_{j=0}^{1819} \frac{1}{j! 2^{1820j}} )^{2^{1820}} ,  (6.19)

requiring only 1820 terms and 1820 squarings, is sufficiently accurate.

BASE CONVERSION

Schönhage has shown that conversion from binary to decimal or vice versa may be done in time O (n (log(n))^{2} log (log(n))) (see Knuth [13, ex. 4.4.14 (revised)]). We describe his method here, as a similar idea is used below to improve Theorem 6.1.

Let β > 1 be a fixed base (e.g. β = 10), and suppose we know the base β representation of an

integer x, i.e. we know the digits d_{0}, …, d_{t−1}, where 0 ≤ d_{i} < β and x = ∑_{0}^{t−1} d_{i}β^{i}. Suppose that

n-bit binary numbers can be multiplied exactly in time M(n), where

2M(n) ≤ M(2n)  (6.20)

for all sufficiently large n. (This is certainly true if the Schönhage-Strassen method [13, 23] is used.) We describe how the binary representation of x may be found in time O (M(n) log(n)), where n is sufficiently large for x to be representable as an n-bit number (i.e. 2^{n} ≥ β^{t}).

Without changing the result, we may suppose t = 2^{k} for some positive integer k. Let the time for conversion to binary and computation of β^{2^{k}} be C(k). Thus, we can compute β^{t/2} and

convert the numbers x_{1} = ∑_{0}^{t/2−1} d_{i}β^{i} and x_{2} = ∑_{t/2}^{t−1} d_{i}β^{i−t/2} to binary in time 2C(k − 1), and then

x = x_{1} + β^{t/2}x_{2} and β^{t} = (β^{t/2})^{2} may be computed in time 2M(n/2) + O(n). Thus

C(k) ≤ 2C(k − 1) + 2M(n/2) + O(n) ,  (6.21)

so

C(k) ≤ 2M(n/2) + 4M(n/4) + 8M(n/8) + … + O(n log(n))
≤ O(M(n) log(n))  (6.22)

(using (6.20)).

The proof that conversion from base 3 to base β may be done in time (6.22) is similar, and once we can convert integers it is easy to convert floating-point numbers.

COMPUTATION OF e AND π

We may regard e − 2 = 1/2! + 1/3! + … as given by a mixed-base fraction 0.111 …, where the base is 2, 3, … Hence, it is possible to evaluate e to precision n, using a slight modification of the above base-conversion method, in time O(M(n) log(n)).

Similarly, artan(1/j) may be computed to precision n in time O(M(n) log^{2}(n)), for any small integer j ≥ 2, and then π may be computed from well-known identities such as Machin’s

π = 16 artan(1/5) − 4 artan(1/239) .  (6.23)

The methods just described are asymptotically faster than the O(n^{2}) methods customarily used in multiple-precision calculations of e and π (see, for example, Shanks and Wrench [25, 26]). It would be interesting to know how large n has to be before the asymptotically faster methods are actually faster. A proof that even faster methods are impossible would be of great interest, for it would imply the transcendence of e and π.

IMPROVED UPPER BOUNDS ON t_{n}(exp) AND t_{n}(sin)

The following lemma uses an idea similar to that described above for base conversion and computation of e.

Lemma 6.1 If p and q are positive integers such that p^{2} ≤ q ≤ 2^{n}, then exp(p/q) may be computed to precision n in time O(M(n) log(n)).

Proof. The approximation

exp(p/q) ≃ ∑_{j=0}^{k} \frac{(p/q)^{j}}{j!}  (6.24)

is sufficiently accurate if k is chosen so that

\frac{(p/q)^{k+1}}{(k + 1)!} ≤ 2^{−n} ≤ \frac{(p/q)^{k}}{k!} .  (6.25)

Since p^{2} ≤ q, (6.25) gives k!q^{k/2} ≤ 2^{n}, so certainly

k!q^{k} ≤ 2^{2n} ,  (6.26)

Hence, a method like that described above for the computation of e may be used, and (6.26) ensures that the integers in intermediate computations do not grow too fast.

From Lemma 6.1 it is easy to deduce Theorem 6.2, which is an improvement of Theorem 6.1 for large n. The methods used in the proof of Theorem 6.1 and the following remarks are, however, faster than that of Theorem 6.2 for small and moderate values of n.

Theorem 6.2 If M(n) satisfies (6.20) then

t_{n}(exp) ≤ c_{32}M(n) log^{2}(n)  (6.27)

and

t_{n}(sin) ≤ c_{33}M(n) log^{2}(n) .  (6.28)

Proof. Without affecting the result (6.27) we may assume that n = 2^{k} for some positive integer k. (This assumption simplifies the proof, but it is not essential.) Given an n-bit fraction x ∈ [0, 1), we write

x = ∑_{i=0}^{k} p_{i}/q_{i} ,  (6.29)

where q_{i} = 2^{2^{i}} and 0 ≤ p_{i} < 2^{2^{i−1}} for i = 0, 1, …, k. By Lemma 6.1, exp(p_{i}/q_{i}) can be computed, to sufficient precision, in time O(M(n) log(n)), so

exp(x) = ∏_{i=0}^{k} exp(p_{i}/q_{i})  (6.30)

can be computed in time O(M(n)(log(n))^{2}). This establishes (6.27), and the proof of (6.28) is similar.

Corollary 6.1

t_{n}(exp) ≤ c_{34}n(log(n))^{3} log log(n)  (6.31)

and

t_{n}(sin) ≤ c_{35}n(log(n))^{3} log log(n) .  (6.32)

Proof. This is immediate from the bound (6.5) and Theorem 6.2.

Corollary 6.2

t_{n}(E_{[a,b]}(f)) ≤ c_{36}n(log(n))^{3} log log(n) ,

where

f(x) = log(x), exp(x), sin(x), cos(x), tan(x), sinh(x),

cosh(x), tanh(x), arsin(x), artan(x), arsinh(x),

etc, and [a, b] is any finite interval on which f(x) is bounded.

Proof. This follows from (6.5), Corollaries 5.1 (and the note following), 5.2, 6.1, and Lemma 3.2.

7 Best constants for operations equivalent to multiplication

In this section, we consider in more detail the relationship between the mp operations D, I, M, R, and S defined in Section 3. It is convenient to consider also the operation Q of forming inverse square roots (i.e., y ← x^{− \frac{1}{2}} ). From Theorem 3.1, if we can perform any one of these operations (say Y) to precision n in time t_{n}(Y), then the time required to perform any of the other operations to precision n is at most a constant multiple of t_{n}(Y).

Definition 7.1 C_{XY} is the minimal constant such that, for all positive ε and all sufficiently large n, the operation X can be performed (to precision n) in time (C_{XY} + ε)t_{n}(Y) if Y can be performed in time t_{n}(Y), where X, Y = D, I, M, Q, R or S.

The following inequalities are immediate consequences of Definition 7.1:

C_{XY} C_{Y Z} ≥ C_{XZ}  (7.1)

and

C_{XY} C_{Y X} ≥ C_{XX} = 1 .  (7.2)

ASSUMPTIONS

To enable us to give moderate upper bounds on the constants C_{XY}, it is necessary to make the following plausible assumption (compare (4.3), (6.20)) throughout this section: for all positive α and ε, and all sufficiently large n,

t_{αn}(Y) ≤ (α + ε)t_{n}(Y)  (7.3)

for Y = D, I, M, Q, R and S. We also assume (6.6).

Table 7.1 gives the best known upper bounds on the constants C_{XY}. Space does not permit a detailed proof of all these upper bounds, but the main ideas of the proof are sketched below.

TABLE 7.1 Upper bounds on C_{XY}

X = D  I  M  Q  R  S
Y = D  1.0  1.0 2.0 3.0 2.0 2.0
I  7.0  1.0 6.0 15.0 14.0 3.0
M  4.0  3.0 1.0 4.5 5.5 1.0
Q 10.0  4.0 6.0 1.0 5.0 3.0
R  7.5  6.0 6.0 3.0 1.0 3.0
S  7.5  5.5 2.0 7.0 9.0 1.0

\underline{C_{IM} ≤ 3}

Use the Newton iteration

x_{i+1} = x_{i} − x_{i}(ax_{i} − 1)  (7.4)

to approximate 1/a using multiplications. At the last iteration it is necessary to compute ax_{i} to precision n, but x_{i}(ax_{i} − 1) only to (relative) precision n/2. Since the order of convergence is 2, the assumptions (7.3) (with α = \frac{1}{2}) and (6.6) give

C_{IM} ≤ (1 + \frac{1}{2})(1 + \frac{1}{2} + \frac{1}{2} + …) = 3 .  (7.5)

\underline{C_{QM} ≤ 4.5}

Use the third-order iteration

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

where

ε_{i} = ax_{i}^{2} − 1  (7.7)

to approximate a^{− \frac{1}{2}} . At the last iteration it is necessary to compute 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

C_{QM} ≤ (2 + \frac{1}{3} + \frac{2}{3})(1 + \frac{1}{3} + \frac{1}{9} + …) = \frac{9}{2} .  (7.8)

Note that this bound is sharper than the bound C_{QM} ≤ 5 which may be obtained from the second-order iteration

x_{i+1} = x_{i} − \frac{1}{2}x_{i}ε_{i} .  (7.9)

\underline{C_{RD} ≤ 2} Use Newton’s iteration

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

to approximate \sqrt{a}. \underline{C_{MS} ≤ 2} This follows from (3.19) and our assumptions. \underline{C_{IS} ≤ 5.5} Use the third-order iteration

x_{i+1} = x_{i} − x_{i} (ε_{i} − ε_{i}^{2})  (7.11)

where

ε_{i} = ax_{i} − 1  (7.12)

to approximate 1/a. \underline{C_{QS} ≤ 7} Use the third-order iteration (7.6). \underline{C_{SI} ≤ 3} From the proof of Lemma 3.7,

t_{n/3}(S) ≤ t_{n}(I) + O(n) .  (7.13)

The result follows from the assumption (7.3) with α = 3. (This is the first time we have used (7.3) with α > 1. The assumption is plausible in view of the Schönhage-Strassen bound (6.5).) Upper bounds on C_{SQ} and C_{SR} follow similarly. \underline{C_{MI} ≤ 6} This follows from (7.1) and our bounds on C_{MS} and C_{SI}. Similarly for the bounds on C_{MQ}, C_{MR} and C_{RI}. \underline{C_{QR} ≤ 3} Use the identity

a^{− \frac{1}{2}} = \frac{1}{λ} (\sqrt{a + λ} − \sqrt{a − λ}) + O (λ^{2}/a^{5/2}) ,  (7.14)

where λ is a power of 2 such that

2^{−n/3−1} ≤ λ/a ≤ 2^{1−n/3} .  (7.15)

Thus

t_{2n/3}(Q) ≤ 2t_{n}(R) + O(n) ,  (7.16)

and the result follows from (7.3).

\underline{C_{DR} ≤ 7.5}

Use the identity

b/a = \frac{1}{λ} (\sqrt{a^{2} + λb} − \sqrt{a^{2} − λb}) + O (λ^{2}b^{3}/a^{5}) ,  (7.17)

where λ is a power of 2 such that (for b ≠ 0)

2^{−n/3−1} ≤ λb/a^{2} ≤ 2^{1−n/3} .  (7.18)

Thus

t_{2n/3}(D) ≤ t_{n}(S) + 2t_{n}(R) + O(n) ,  (7.19)

and the result follows.

\underline{C_{IR} ≤ 6}

a^{−1} = (a^{2})^{− \frac{1}{2}} ,  (7.20)

so

C_{IR} ≤ C_{SR} + C_{QR} ≤ 6 .  (7.21)

The bound on C_{IQ} also follows from (7.20), and then the bound on C_{RQ} follows from a^{\frac{1}{2}} = (a^{−1})^{− \frac{1}{2}} .

8 Comparison of some mp methods for nonlinear equations

In this section, we briefly consider methods for finding multiple-precision solutions of non-linear equations of the form

f(x) = 0,  (8.1)

where f(x) can be evaluated for any x in some domain. Additional results are given in [38].

There are many well-known results on the efficiency of various methods for solving (8.1), e.g., Hindmarsh [10], Ostrowski [20], Traub [27] and the references given in Section 1, but the results are only valid if arithmetic operations (in particular the evaluation of f(x), f′(x) etc.) require certain constant times. The examples given below demonstrate that different considerations are relevant when multiple-precision arithmetic of varying precision is used.

For simplicity, we restrict attention to methods for finding a simple zero ζ of f by evaluating f at various points. We assume that f has sufficiently many continuous derivatives in a neighbourhood of ζ, but the methods considered do not require the evaluation of these derivatives.

Since f(x) is necessarily small near ζ, it is not reasonable to assume that f(x) can be evaluated to within a small relative error near ζ. In this section, an evaluation of f “with precision n” means with an absolute error of order 2^{−n}. We suppose that such an evaluation requires time w(n) = t_{n}(E(f)), where

w(cn) ∼ c^{α}w(n)  (8.2)

for some constant α > 1 and all positive c. Since α > 1, the bound (6.5) and condition (8.2) give

lim_{n→∞} t_{n}(M)/w(n) = 0 ,  (8.3)

so we may ignore the time required for a fixed number of multiplications and divisions per iteration, and merely consider the time required for function evaluations. Our results also apply if α = 1, so long as (8.3) holds. (For example, the evaluation of exp(x) by the method of Corollary 6.1 requires time w(n) ∼ c_{37}n(log(n))^{3} log log(n), which satisfies (8.2) with α = 1, and also satisfies (8.3).)

Definition 8.1 If an mp zero-finding method requires time t(n) ∼ C(α)w(n) to approximate ζ ≠ 0 with precision n, where w(n) and ζ are as above, then C(α) is the asymptotic constant of the method. (Not to be confused with the asymptotic error constant as usually defined for fixed-precision methods [2].)

Given several mp methods with various asymptotic constants, it is clear that the method with minimal asymptotic constant is the fastest (for sufficiently large n). The method which is fastest may depend on α, as the following examples show.

DISCRETE NEWTON mp METHODS

Consider iterative methods of the form

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

where g_{i} is a finite-difference approximation to f′(x_{i}). If ε_{i} = |x_{i} − ζ| is sufficiently small, f(x_{i}) is evaluated with absolute error O (ε_{i}^{2}), and

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

then

|x_{i+1} − ζ| = O (ε_{i}^{2}) ,  (8.6)

so the method has order (at least) 2.

The simplest method of estimating f′(x_{i}) to sufficient accuracy is to use the one-sided difference

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

where h_{i} is of order ε_{i}, and the evaluation of f(x_{i} +h_{i}) and f(x_{i}) are performed with an absolute error O (ε_{i}^{2}). Thus, to obtain ζ to precision n by this method (N_{1}), we need two evaluations of f to precision n (at the last iteration), preceded by two evaluations to precision n/2, etc. (The same idea is used above, in the proof of Theorem 5.2.) The time required is

t(n) ∼ 2w(n) + 2w(n/2) + 2w(n/4) + ….  (8.8)

Thus, from (8.2) and Definition 8.1, the asymptotic constant is

C_{N_{1}}(α) = 2(1 + 2^{−α} + 2^{−2α} + …) = 2/(1 − 2^{−α}) .  (8.9)

Since

2 < C_{N_{1}}(α) ≤ 4 ,  (8.10)

the time required to solve (8.1) to precision n is only a small multiple of the time required to evaluate f to the same precision. The same applies for the methods described below.

Using (8.7) is not necessarily the best way to estimate f′(x_{i}). Let p be a fixed positive integer, and consider estimating f′(x_{i}) by evaluating f at the points

x_{i} − ⌊p/2⌋h_{i}, x_{i} − (⌊p/2⌋ − 1)h_{i}, …, x_{i} + ⌈p/2⌉h_{i} .

(The points need not be equally spaced so long as their minimum and maximum separations are of order h_{i}.) Let g_{i} be the derivative (at x_{i}) of the Lagrange interpolating polynomial agreeing with f at these points. Since estimates f′(x_{i}) with truncation error O (h_{i}^{p}), we need h_{i} of order ε_{i}^{1/p}. Then, to ensure that (8.5) holds, the function evaluations at the above points must be made with absolute error O (ε_{i}^{1+1/p}). Thus to obtain ζ to precision n by this method (N_{p}) we need one evaluation f to precision n and p evaluations to precision n(1 + 1/p)/2, preceded by one evaluation precision n/2 and p to precision n(1 + 1/p)/4, etc. The asymptotic constant is

C_{N}(p, α) = (1 + p (\frac{p + 1}{2p})^{α}) /(1 − 2^{−α}) .  (8.11)

Let

C_{N}(α) = min_{p=1,2,...} C_{N}(p, α) ,  (8.12)

so the “optimal mp discrete Newton method” has asymptotic constant C_{N}(α). From (8.11), the p which minimizes C_{N}(p, α) also minimizes p^{1/α}(1 + 1/p), so the minimum for α > 1 occurs at p = ⌊α − 1⌋ or ⌈α − 1⌉. In fact, p = 1 is optimal if

1 ≤ α < log(2)/ log(4/3) = 2.4094 …,  (8.13)

and p ≥ 2 is optimal if

\frac{log(1 − p^{−1})}{log(1 − p^{−2})} < α < \frac{log(1 + p^{−1})}{log(1 + 1/(p(p + 2)))} .  (8.14)

The result that method N_{2} is more efficient than method N_{1} if α > 2.4094 … is interesting, for N_{2} requires one more function evaluation per iteration than N_{1}, and has the same order of convergence. The reason is that not all the function evaluations need to be as accurate for method N_{2} as for method N_{1}. We give below several more examples where methods with lower order and/or more function evaluations per iteration are more efficient than methods with higher order and/or less function evaluations per iteration.

For future reference, we note that

1 < C_{N}(α) ≤ 4 ,  (8.15)

C_{N}(1) = 4 ,  (8.16)

and

C_{N}(α) − 1 ∼ eα2^{−α}  (8.17)

as α → ∞.

A CLASS OF mp SECANT METHODS

It is well-known that the secant method is more efficient that the discrete Newton method for solving nonlinear equations with fixed-precision arithmetic [2, 20]. For mp methods the comparison depends on the exponent α in (8.2).

Let k be a fixed positive integer and p_{k} the positive real root of

x^{k+1} = 1 + x^{k} .  (8.18)

The iterative method S_{k} is defined by

x_{i+1} = x_{i} − f(x_{i}) (\frac{x_{i} − x_{i−k}}{f(x_{i}) − f(x_{i−k})}) ,  (8.19)

where the function evaluations are performed to sufficient accuracy to ensure that the order of convergence is at least p_{k}. Thus, S_{1} is the usual secant method with order p_{1} = \frac{1+\sqrt{5}}{2} = 1.618 …; S_{2}, S_{3} etc. are methods with lower orders p_{2} = 1.4655 …, p_{3} = 1.3802 …, etc. With fixed-precision S_{1} is always preferable to S_{2}, S_{3} etc., but this is not always true if mp arithmetic is used.

Suppose i and k fixed, δ > 0 small, and write ε = |x_{i−k} − ζ| and p = p_{k} − δ. Since the order of convergence is at least p, we have

|x_{i} − ζ| = O (ε^{p^{k}}) ,  (8.20)

|x_{i+1} − ζ| = O (ε^{p^{k+1}}) ,  (8.21)

|x_{i} − x_{i−k}| = O(ε) ,  (8.22)

and

|f(x_{i})| = O (ε^{p^{k}}) .  (8.23)

For the approximate evaluation of the right side of (8.19) to give order p, the absolute error in the evaluation of f(x_{i}) must be O (ε^{p^{k+1}}), and the relative error in the evaluation of (f(x_{i}) −

f(x_{i−k}))/(x_{i} − x_{i−k}) must be O (ε^{p^{k+1}−p^{k}}), so the absolute error in the evaluation of f(x_{i−k})

must be O (ε^{p^{k+1}−p^{k}+1}). From (7.18), for δ sufficiently small,

p^{k+1} − p^{k} + 1 > p ,  (8.24)

so the evaluation of ζ to precision n by method S_{k} requires evaluations of f to precision n, n/p, n/p^{2}, …, n/p^{k−1}, 2n/p^{k+1}, 2n/p^{k+2}, etc. Thus, the asymptotic constant is

C_{S}(k, α) = 1 + p^{−α} + … + p^{(1−k)α} + (2p^{−(k+1)})^{α}(1 + p^{−α} + …)
= \frac{1 − p^{−kα} + (2p^{−(k+1)})^{α}}{1 − p^{−α}} ,  (8.25)

where (after letting δ → 0) p = p_{k} satisfies (8.18).

We naturally choose k to minimize C_{S}(k, α), giving the “optimal mp secant method” with asymptotic constant

C_{S}(α) = min_{k=1,2,...} C_{S}(k, α) .  (8.26)

The following lemmas show that the optimal secant method is S_{1} if α < 4.5243 …, and S_{2} if α > 4.5243 …

Lemma 8.1

C_{S}(k, 1) = 3 + p_{k}^{k} − p_{k} .  (8.27)

Proof. Easy from (8.18) and (8.25).

Lemma 8.2

C_{S}(k, α) − 1 ∼ \{{(3 − \sqrt{5})^{α} if k = 1 , \atop p_{k}^{−α}  if k ≥ 2 ,}  (8.28)

as α → ∞.

Proof. From (8.25),

C_{S}(k, α) − 1 ∼ p_{k}^{−α} − p_{k}^{−kα} + (2p_{k}^{−(k+1)})^{α}  (8.29)

as α → ∞. If k ≥ 2 then, from (8.18),

p_{k}^{k} = p_{k}^{−1} + p_{k}^{k−1} ≥ p_{k}^{−1} + p_{k} > 2 ,  (8.30)

so

p_{k}^{−1} > 2p_{k}^{−(k+1)} .  (8.31)

Thus, the result for k ≥ 2 follows from (8.29). The result for k = 1 also follows from (8.29), for 2p_{1}^{−2} = 3 − \sqrt{5}.

Lemma 8.3

C_{S}(α) = \{{C_{S}(1, α) if 1 ≤ α ≤ α_{0} , \atop C_{S}(2, α) if α ≥ α_{0} ,}  (8.32)

where α_{0} = 4.5243 … is the root of

C_{S}(1, α_{0}) = C_{S}(2, α_{0}) .  (8.33)

Proof. The details of the proof are omitted, but we note that the result follows from Lemmas 8.1 and 8.2 for (respectively) small and large values of α.

From (8.25), C_{S}(k, α) is a monotonic decreasing function of α, so the same is true of C_{S}(α). Thus, from Lemmas 8.1, 8.2 and 8.3,

1 < C_{S}(α) ≤ 3 ,  (8.34)

C_{S}(1) = 3 ,  (8.35)

and

C_{S}(α) − 1 ∼ p_{2}^{−α} = (0.6823 …)^{α}  (8.36)

as α → ∞. Comparing these results with (8.15) to (8.17), we see that the optional mp secant method is more efficient than the optimal mp discrete Newton method for small α, but less efficient for large α. (The changeover occurs at α = 8.7143 …)

AN mp METHOD USING INVERSE QUADRATIC INTERPOLATION

For fixed-precision arithmetic the method of inverse quadratic interpolation [2] is slightly more efficient than the secant method, for it has order P_{Q} = 1.8392 … > 1.6180 …, and requires the same number (one) of function evaluations per iteration. For mp arithmetic, it turns out that inverse quadratic interpolation (Q) is always more efficient than the secant method S_{1}, but it is less efficient than the secant method S_{2} if α > 5.0571 … Since the analysis is similar to that for method S_{1} above, the details are omitted. The order p_{Q} is the positive real root of

x^{3} = 1 + x + x^{2} .  (8.37)

For brevity, we write σ = 1/p_{Q} = 0.5436 …

To evaluate ζ to precision n by method Q requires evaluations of f to precision n, (1 − σ + σ^{2})n, and σ^{j}(1 − σ − σ^{2} + 2σ^{3})n for j = 0, 1, 2, … Hence, the asymptotic constant is

C_{Q}(α) = 1 + (1 − σ + σ^{2})^{α} + (1 − σ − σ^{2} + 2σ^{3})^{α}/(1 − σ^{α})
= 1 + (1 − σ + σ^{2})^{α} + (3σ^{3})^{α}/(1 − σ^{α})  (8.38)

from (8.31). Corresponding to the results (8.15) to (8.17) and (8.34) to (8.36), we have that C_{Q}(α) is monotonic decreasing,

1 < C_{Q}(α) ≤ C_{Q}(1) = \frac{1}{2}(7 − 2σ − σ^{2}) = 2.8085 …,  (8.39)

and

C_{Q}(α) − 1 ∼ (1 − σ + σ^{2})^{α} = (0.7519 …)^{α}  (8.40)

as α → ∞. Method Q is more efficient than the optimal mp secant method if α < 5.0571 …, and more efficient than the optimal mp discrete Newton method if α < 7.1349 … We do not know any mp method which is more efficient than method Q for α close to 1.

OTHER mp METHODS USING INVERSE INTERPOLATION

Since inverse quadratic interpolation is more efficient than linear interpolation (at least for α close to 1), it is natural to ask if inverse cubic or higher degree interpolation is even more efficient. Suppose \frac{1}{2} ≤ μ < 1, and consider an inverse interpolation method I_{μ} with order 1/μ. In particular, consider the method I_{μ} which uses inverse interpolation at x_{i}, x_{i−1}, …, x_{i−k} to generate x_{i+1}, where k is sufficiently large, and the function evaluations at x_{i}, …, x_{i−k} are sufficiently accurate to ensure that the order is at least 1/μ and, in general, no more than 1/μ. (The limiting case I_{1/2} is the method which uses inverse interpolation through all previous points x_{0}, x_{1}, …, x_{i} to generate x_{i+1}.)

By an analysis similar to those above, it may be shown that the asymptotic constant of method I_{μ} is

C_{I}(μ, α) = ∑_{j=0}^{∞} (s_{j}(μ))^{α} ,  (8.41)

where s_{0}(μ) = 1 and

s_{j}(μ) = max [μs_{j−1}(μ), 1 + jμ^{j+1} − μ(1 − μ^{j})/(1 − μ)]  (8.42)

for j = 1, 2, … Space does not allow a proof of (8.41), but related results are given in [20, Appendix H]. We note the easily verified special cases

C_{I} (\frac{\sqrt{5} − 1}{2}, α) = C_{S}(1, α)  (8.43)

and

C_{I}(σ, α) = C_{Q}(α) .  (8.44)

The method with maximal order (see [7]) is I_{1/2}, with asymptotic constant

C_{I}(\frac{1}{2}, α) = ∑_{j=2}^{∞}(j2^{1−j})^{α} .  (8.45)

The “optimal mp inverse interpolatory method” is the method I_{μ} with μ(α) chosen to minimize C_{I}(μ, α), so its symptotic constant is

C_{I}(α) = min_{\frac{1}{2} ≤μ≤1} C_{I}(μ, α) .  (8.46)

The following lemma shows that the optimal choice is μ = σ, corresponding to the inverse quadratic method Q discussed above, if α ≤ 4.6056 …

Lemma 8.4 If C_{I}(α) = C_{I}(μ(α), α) then

μ(α) = σ = 0.5436 … if 1 ≤ α ≤ 4.6056 …,  (8.47)

μ(α) is a monotonic decreasing function of α, and

lim_{α→∞} μ(α) = \frac{1}{2} .  (8.48)

From (8.39),

C_{I}(\frac{1}{2}, α) − 1 ∼ (\frac{3}{4})^{α}  (8.49)

as α → ∞, so Lemma 8.4 shows that the optimal inverse interpolatory is more efficient than methods S_{1} and Q (as expected), but less efficient than method S_{2} or the optimal discrete Newton method, for large α. In fact C_{I}(α) < C_{S}(α) for 1 ≤ α < 5.0608 …

A LOWER BOUND FOR C(α)

The following theorem shows that C(α) ≥ 1 for all useful mp methods. The results above (e.g. (7.17)) show that the constant “1” here is best possible, as methods with C(α) → 1 as α → ∞ are possible. The minimal value of C(α) for any finite α is an open question.

Theorem 8.1 If an mp method is well-defined and converges to a zero of the functions f_{1}(x) = F (x) − y and f_{2}(y) = F^{(−1)}(y) − x, where x and y are restricted to nonempty domains D_{x} and D_{y}, and F is some invertible mapping of D_{x} onto D_{y} such that t_{n}(E(F)) satisfies (8.2), then the asymptotic constant of the method satisfies C(α) ≥ 1.

Proof. If C(α) < 1 then, by solving f_{1}(x) = 0, we can evaluate F^{(−1)}(y) (for y in D_{y}) in time less than t_{n}(E(F)), for all sufficiently large n. Applying the same argument to f_{2}(y), we can evaluate F = (F^{(−1)})^{(−1)} in time less than t_{n}(E(F^{(−1)})). Hence, for large n we have

t_{n}(E(F)) < t_{n}(E(F^{(−1)})) < t_{n}(E(F)) ,  (8.50)

a contradiction. Hence, C(α) ≥ 1.

Conjecture 8.1 For all mp methods (using only function evaluations) which are well-defined and convergent for some reasonable class of functions with simple zeros,

C(α) ≥ 1/(1 − 2^{−α}) .  (8.51)

SUMMARY OF mp ZERO-FINDING METHODS

Of the methods described in this section, the most efficient are:

1. optimal inverse interpolation, if 1 ≤ α ≤ 5.0608 … (equivalent to inverse quadratic interpolation, if 1 ≤ α ≤ 4.6056 …) ;

2. optimal secant method (method S_{2}), if 5.0608 … < α ≤ 8.7143 …;

3. optimal discrete Newton, if 8.7143 … < α.

For practical purposes, the inverse quadratic interpolation method is to be recommended, for it is easy to program, and its asymptotic C_{Q}(α) is always within 3.2% of the least constant for the methods above. Numerical values of the asymptotic constants, for various values of α, are given to 4D in Table 8.1. The smallest constant for each α is italicized.

TABLE 8.1 Aysmptotic constants for various mp methods

\begin{matrix}α & C_{N}(α) & C_{S}(1, α) & C_{S}(2, α) & C_{Q}(α) & C_{I}(α) & C_{I}(\frac{1}{2}, α) \\ 1.0 & 4.0000 & 3.0000 & 3.6823 & 2.8085 & 2.8085 & 3.0000 \\ 1.1 & 3.7489 & 2.8093 & 3.4256 & 2.6484 & 2.6484 & 2.8193 \\ 1.5 & 3.0938 & 2.2987 & 2.7241 & 2.2108 & 2.2108 & 2.3219 \\ 2.0 & 2.6667 & 1.9443 & 2.2209 & 1.8954 & 1.8954 & 1.9630 \\ 3.0 & 2.1071 & 1.5836 & 1.6935 & 1.5586 & 1.5586 & 1.5856 \\ 4.0 & 1.6988 & 1.3988 & 1.4248 & 1.3789 & 1.3789 & 1.3898 \\ 5.0 & 1.4260 & 1.2860 & 1.2694 & 1.2677 & 1.2676 & 1.2718 \\ 6.0 & 1.2529 & 1.2105 & 1.1741 & 1.1936 & 1.1930 & 1.1946 \\ 7.0 & 1.1469 & 1.1573 & 1.1137 & 1.1420 & 1.1410 & 1.1416 \\ 8.0 & 1.0838 & 1.1185 & 1.0748 & 1.1051 & 1.1039 & 1.1041 \\ 9.0 & 1.0471 & 1.0898 & 1.0495 & 1.0782 & 1.0770 & 1.0771 \\ 10.0 & 1.0262 & 1.0682 & 1.0328 & 1.0584 & 1.0573 & 1.0573 \\ 15.0 & 1.0012 & 1.0176 & 1.0043 & 1.0139 & 1.0134 & 1.0134 \\ 20.0 & 1.0001 & 1.0046 & 1.0006 & 1.0033 & 1.0032 & 1.0032 \end{matrix}

NOTE ADDED IN PROOF. Theorem 6.2 and its corollaries may be improved by a factor log(n), as described in [37] and [38].

References

[1] Brent, R.P. The computational complexity of iterative methods for systems of nonlinear equations. In Complexity of Computer Computations (edited by R.E. Miller and J.W. Thatcher). Plenum Press, New York, 1972, 61–71.

[2] Brent, R.P. Algorithms for Minimization without Derivatives. Prentice-Hall, Englewood Cliffs, New Jersey, 1973.

[3] Brent, R.P. Some efficient algorithms for solving systems of nonlinear equations. SIAM Numer. Anal. 10, 327–344, 1973.

[4] Brent, R.P. On the precision attainable with various floating-point number systems. IEEE Trans. Comp. C-22, 601–607, 1973.

[5] Brent, R.P. Error analysis of algorithms for matrix multiplication and triangular decomposition using Winograd’s identity. Numer. Math. 16, 145–156, 1970.

[6] Brent, R.P. Numerical solution of nonlinear equations. Computer Sci. Dept., Stanford University, March 1975, 189 pp.

[7] Brent, R.P., Winograd, S. and Wolfe, P. Optimal iterative processes for rootfinding. Numer. Math. 20, 327–341, 1973.

[8] Cook, S.A. and Aanderaa, S.O. On the minimum complexity of functions. Trans. Amer. Math. Soc. 142, 291–314, 1969.

[9] Floyd, R.W. Unpublished notes.

[10] Hindmarsh, A.C. Optimality in a class of rootfinding algorithms. SIAM J. Numer. Anal. 9, 205–214, 1972.

[11] Hopcroft, J.E. Complexity of computer computations. In Information Processing 74. North-Holland, Amsterdam, 1974, 620–626.

[12] Karatsuba, A. and Ofman, Y. Multiplication of multidigit numbers on automata (Russian). Dokl. Akad. Nauk SSSR 145, 293-294, 1962.

[13] Knuth, D.E. The Art of Compter Programming, Vol. II, Seminumerical Algorithms. Addison Wesley, Reading, Massachusetts, 1969. Errata and addenda: Report CS 194, Computer Sci. Department, Stanford University, 1970.

[14] Kung, H.T. The computational complexity of algebraic numbers. SIAM J. Numer. Anal. (to appear).

[15] Kung, H.T. A bound on the multiplicative efficiency of iteration. J. Computer & System Sciences 7, 334–342, 1973.

[16] Kung, H.T. and Traub, J.F. Optimal order of one-point and multipoint iteration. J. ACM 21, 643–651, 1974.

[17] Kung, H.T. and Traub, J.F. Computational complexity of one-point and multipoint iteration. In Complexity of Real Computations (edited by R. Karp). Amer. Math. Soc., Providence, Rhode Island, 1974, 149–160.

[18] Kung, H.T. and Traub, J.F. Optimal order and efficiency for iterations with two evaluations. Tech. Report, Department of Computer Science, Carnegie-Mellon University, 1973.

[19] Morgenstern, J. The linear complexity of computation. J. ACM 20, 305–306 (1973).

[20] Ostrowski, A.M. Solution of Equations in Euclidean and Banach Spaces. Academic Press, New York, 1973.

[21] Paterson, M.S. Efficient iterations for algebraic numbers. In Complexity of Computer Computations (edited by R.E. Miller and J.W. Thatcher). Plenum Press, New York, 1972, 41–52.

[22] Rissanen, J. On optimum root-finding algorithms. J. Math. Anal. Applics. 36, 220–225, 1971.

[23] Schönhage, A. and Strassen, V. Schnelle Multiplikation grosser Zahlen. Computing 7, 281–292, 1971.

[24] Schultz, M.H. The computational complexity of elliptic partial differential equations. In Complexity of Computer Computations (edited by R.E. Miller and J.W. Thatcher). Plenum Press, New York, 1972, 73–83.

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

[26] Shanks, D. and Wrench, J.W. Calculation of e to 100,000 decimals. Math. Comp. 23 679–680, 1969.

[27] Traub, J.F. Iterative Methods for the Solution of Equations. Prentice-Hall, Englewood Cliffs, New Jersey, 1964.

[28] Traub, J.F. Computational complexity of iterative processes. SIAM J. Computing 1, 167–179, 1972.

[29] Traub, J.F. Numerical mathematics and computer science. Comm. ACM 15, 537–541, 1972.

[30] Traub, J.F. Optimal iterative processes: theorems and conjectures. In Information Processing 71 . North-Holland, Amsterdam, 1972, 1273–1277.

[31] Traub, J.F. Theory of optimal algorithms. In Software for Numerical Mathematics (edited by D.J. Evans). Academic Press, 1974.

[32] Traub, J.F. An introduction to some current research in numerical computational complexity. Tech. Report, Department of Computer Science, Carnegie-Mellon University, 1973.

[33] Winograd, S. On the time required to perform addition. J. ACM 12, 277–285, 1965.

[34] Winograd, S. On the time required to perform multiplication. J. ACM 14, 793–802, 1967.

[35] Wozniakowski, H. Generalized information and maximal order of iteration for operator equations. SIAM J. Numer. Anal. 12, 121–135, 1975.

[36] Wozniakowski, H. Maximal stationary iterative methods for the solution of operator equations. SIAM J. Numer. Anal. 11, 934–949, 1974.

[37] Brent, R.P. Fast multiple-precision evaluation of elementary functions. J. ACM 23, 242–251, 1976.

[38] Brent, R.P. Multiple-precision zero-finding methods and the complexity of elementary function evaluation. In Analytic Computational Complexity (edited by J.F. Traub). Academic Press, New York, 1975, 59–73.

Postscript (September 1999)

Historical Notes

This paper was retyped (with minor corrections) in L^{A}TEX during August 1999. It is available electronically in compressed postscript format from ftp://ftp.comlab.ox.ac.uk/pub/

Documents/techpapers/Richard.Brent/rpb032.ps.gz

The related paper Brent [38] is available electronically in compressed postscript format from

ftp://ftp.comlab.ox.ac.uk/pub/Documents/techpapers/Richard.Brent/rpb028.ps.gz

The paper Kung [14] appeared in SIAM J. Numer. Anal. 12 (1975), 89–96.

Sharper Results

Some of the constants given in Table 7.1 can be improved, e.g. C_{DM}, C_{RM}, C_{DS}, C_{RS}. One source of improvement is given in a report by Karp and Markstein^{2}.

For example, consider C_{DM}. We want to compute an n-bit approximation to b/a. If x_{i} → 1/a as in (7.4) and we define y_{i} = bx_{i}, then y_{i} → b/a. Also, if x_{i} satisfies the recurrence (7.4), then y_{i} satisfies

y_{i+1} = y_{i} − x_{i}(ay_{i} − b) .  (7.4′)

Note that (7.4’) is self-correcting because of the computation of the residual ay_{i} − b. Suppose x_{i} has (relative) precision n/2. If we approximate y_{i} = bx_{i} using an \frac{n}{2}-bit multiplication, compute the residual ay_{i} − b using an n-bit multiplication, then its product with x_{i} using an \frac{n}{2}-bit multiplication, we can apply (7.4’) to obtain y_{i+1} with relative precision n. Assuming x_{i} is obtained in time ∼ 3M(n/2) ∼ \frac{3}{2}M(n) (see (7.5)), the time to obtain y_{i+1} is ∼ \frac{7}{2}M(n), i.e. C_{DM} ≤ 3.5, which is sharper than the bound C_{DM} ≤ 4.0 given in Table 7.1.

Similarly, we can obtain C_{RM} ≤ 4.25, which is sharper than the bound C_{RM} ≤ 5.5 given in Table 7.1. If x_{i} → a^{−1/2} and y_{i} = ax_{i} → \sqrt{a}, we compute a precision n/2 approximation x_{i} in time ∼ \frac{9}{2}M(n/2) as in Section 7, then apply a final second-order iteration for

y_{i+1} = y_{i} − x_{i}(y_{i}^{2} − a)/2  (7.9′)

(derived by multiplying (7.9) by a and using (7.7)) to obtain a precision n approximation y_{i+1} to \sqrt{a}.

As a corollary, the time required for an arithmetic-geometric mean iteration [37,38] is reduced from ∼6.5M(n) to ∼5.25M(n).

The Definition of n-bit Multiplication

Our t_{n}(M) (see Sections 1–3) is essentially the time required to compute the most significant n bits in the product of two n-bit numbers. In Brent [38], t_{n}(M) is written as M(n). A related but subtly different function is M^{∗}(n), defined as the time required to compute the full 2n-bit product of n-bit numbers^{3}. Paul Zimmermann^{4} observed that smaller constants can sometimes be obtained in row Y = M of Table 7.1 if we use M^{∗}(n) instead of M(n). (We denote these constants by C_{XM^{∗}} to avoid confusion with the C_{XM} of Table 7.1.) For example, C_{DM^{∗}} < 3.5 and C_{RM^{∗}} < 4.25.

2Alan H. Karp and Peter Markstein, High Precision Division and Square Root, HP Labs Report 93-93-42 (R.1), June 1993, Revised October 1994. Available electronically from http://www.hpl.hp.com/techreports/ 93/HPL-93-42.html

It is an open question whether

M(n) ∼ M^{∗}(n) as n → ∞ ;

with the best available multiplication algorithms (those based on the FFT) this is true^{5}.

Final Comments

Daniel Bernstein^{6} observed that the time required to compute n-bit square roots can be reduced further if the model of computation is relaxed so that redundant FFTs can be eliminated. Similar remarks apply to division, exponentiation etc (and to operations on power series).

In conclusion, 25 years after the paper was written (in 1974), improvements can still be found, and the last word is yet to be written!

3In Brent [37] we (confusingly) used the notation M(n) for M^{∗}(n).

4Personal communication, 1999.

5Similar remarks apply if we consider computing the product of two polynomials of degree n−1, and ask either for the first n terms in the product or the complete product. Although the first computation is faster (by a factor of about two) if the classical order n^{2} algorithms are used, it is not significantly faster if FFT-based algorithms are used.

6Personal communication, 1999.