Target: Mathematics of Computation. New version of paper in preparation, with 4.5 for reciprocal, 6.5 for quotient, 5.5 for square root, and 8.5 for exponential. See

http://cr.yp.to/fastnewton.html.

REMOVING REDUNDANCY IN HIGH-PRECISION NEWTON ITERATION

DANIEL J. BERNSTEIN

Abstract. This paper speeds up Brent’s algorithms for various high-precision computations in the power series ring C[[t]]. If it takes time 3 to compute a product then it takes time roughly 5.6 to compute a reciprocal; roughly 8.2 to compute a quotient or a logarithm; roughly 6.5 to compute a square root; roughly 9 to compute both a square root and a reciprocal square root; and roughly 10.4 to compute an exponential. The same ideas apply to approximate computations in R, Q_{p}, etc.

1. Introduction

Let f ∈ C[[t]] be a power series with constant coefficient 1. How can one compute f^{−1}? The standard answer is Newton’s method, which shows how to compute f^{−1} mod t^{2n} from f^{−1} mod t^{n} with a few size-n multiplications: f^{−1} mod t^{2n} = g_{0} + (1 − fg_{0})g_{0} mod t^{2n} where g_{0} = f^{−1} mod t^{n}. One can compute g_{0} by the same method recursively. (This algorithm is sometimes credited to Sieveking, who published it in [8]; Kung in [6] pointed out that Sieveking’s method was an example of Newton’s method.) Similar comments apply to computing f^{1/2}, log f, et al.

Using FFTs one can multiply polynomials of degree up to n with 3kn log n+O(n) arithmetic operations in C for some constant k. Brent in [5] showed that one can compute the first n coefficients of f^{−1} with 9kn log n + O(n) operations, log f with 12kn log n + O(n) operations, f^{−1/2} with 13.5kn log n + O(n) operations, both f^{1/2} and f^{−1/2} with 16.5kn log n + O(n) operations, and exp f with 22kn log n + O(n) operations.

The point of this paper is that some obvious redundancies account for a large fraction of the run time of Brent’s algorithms. Sections 2, 3, 4, and 5 present several streamlined examples of Newton’s method, culminating in an algorithm to compute n coefficients of exp f with only 10.4kn log n + O(n) operations. Section 6 gives implementation results.

Generalizations. Newton’s method is not limited to C[[t]]; it is also used for highprecision computations in k[[t]] where k is a finite field, in the p-adic numbers Q_{p}, in R, and so on. See [3] for a survey of relevant multiplication methods. Roundoff error analysis is typically required for Q_{p}, where it is easy, and for R, where it is not so easy; see, e.g., [2, section 8 and section 21]. Some functions, such as log, require completely different methods for R; see, e.g., [5].

Brent’s methods can be streamlined in all of these situations. I have focused on C[[t]] in this paper since it is the simplest case. It is also the case used in [4].

Date: DRAFT 19980627.

1991 Mathematics Subject Classification. Primary 68Q40; Secondary 65Y20.

The author was supported by the National Science Foundation under grant DMS–9600083.

Notation and terminology. Throughout this paper, n is a positive integer. Subscripted variables such as f_{0}, f_{1}, … refer to polynomials of degree below n. Thus every power series can be written uniquely in the form f_{0} + t^{n}f_{1} + t^{2n}f_{2} + · · ·.

A transform means a size-2n FFT or a size-2n inverse FFT. If p is a polynomial of degree below 2n then p^{∗} means the result of applying a size-2n FFT to p. Recall that one can compute a bilinear form such as f_{0}g_{2} + f_{1}g_{1} + f_{2}g_{0}, given f_{0}^{∗}, f_{1}^{∗}, f_{2}^{∗}, g_{0}^{∗}, g_{1}^{∗}, g_{2}^{∗}, with a single transform plus O(n) operations.

If f is a power series then D(f) means t times the derivative of f. If f is a power series with constant coefficient 0 then I(f) means the integral of f/t.

2. Reciprocals

Let f and g be power series with fg = 1. This section considers the problem of computing g given f.

Write f = f_{0} + f_{1}t^{n} + f_{2}t^{2n} + ⋯ and g = g_{0} + g_{1}t^{n} + g_{2}t^{2n} + · · ·. Define q_{1} and r_{1} by 1 + q_{1}t^{n} = f_{0}g_{0} and r_{1} = f_{1}g_{0} mod t^{n}. Then g_{1} = −(q_{1} + r_{1})g_{0} mod t^{n}.

Given f_{0}, f_{1}, g_{0}, Brent suggested computing f_{0}g_{0}, hence q_{1}; f_{1}g_{0}, hence r_{1}; and then (q_{1} +r_{1})g_{0}, hence g_{1}. Each multiplication can be done with 3 transforms, for a total of 9 transforms to compute g mod t^{2n} given g mod t^{n}. The work to compute g mod t^{n} by the same method recursively is comparable to 9/2 + 9/4 + ⋯ = 9 transforms.

However, each multiplication uses the same intermediate result g_{0}^{∗}, so only 7 transforms are required: f_{0}^{∗}, f_{1}^{∗}, g_{0}^{∗}, f_{0}g_{0}, f_{1}g_{0}, (q_{1} + r_{1})^{∗}, and (q_{1} + r_{1})g_{0}.

Higher-order iterations. Define q_{2}, q_{3}, r_{2}, r_{3} by (f_{0} + f_{1}t^{n})(g_{0} + g_{1}t^{n}) = 1 + q_{2}t^{2n} + q_{3}t^{3n} and r_{2} + r_{3}t^{n} = (f_{2} + f_{3}t^{n})(g_{0} + g_{1}t^{n}) mod t^{2n}. Then g_{2} + g_{3}t^{n} = −(q_{2} + r_{2} + (q_{3} + r_{3})t^{n})(g_{0} + g_{1}t^{n}) mod t^{2n}.

Brent suggested first computing g_{1} as discussed above, then using size-4n FFTs to multiply f_{0} + f_{1}t^{n}, f_{2} + f_{3}t^{n}, and q_{2} + r_{2} + (q_{3} + r_{3})t^{n} by g_{0} + g_{1}t^{n}.

However, it is wasteful to feed f_{0} and g_{0} to two different sizes of FFTs. One can do better by sticking with size-2n FFTs: compute g_{1}^{∗}, f_{1}g_{0} + f_{0}g_{1}, f_{1}g_{1}, f_{2}^{∗}, f_{3}^{∗}, f_{2}g_{0}, f_{2}g_{1} + f_{3}g_{0}, (q_{2} + r_{2})^{∗}, (q_{3} + r_{3})^{∗}, (q_{2} + r_{2})g_{0}, and (q_{2} + r_{2})g_{1} + (q_{3} + r_{3})g_{0}.

The total is 18 transforms to compute g mod t^{4n} from g mod t^{n}. The work to compute g mod t^{n} by the same method recursively is comparable to 18/4+18/16+ ⋯ = 6 transforms.

The same idea can be applied repeatedly. One can compute g mod t^{8n} from g mod t^{n} with just 40 transforms, for example. The work to compute g mod t^{n} recursively is comparable to 5.5+1.5/(2^{e} −1) transforms with an order-2^{e} iteration. One should select e as a function of n to avoid excessive overhead.

Notes. The idea of reusing transforms is standard. For example, it is well known that squaring takes only 2 transforms while multiplication takes 3. See [3, section 12] for references.

A 6-transform bound for reciprocals, presumably using an order-4 iteration as above, was announced by Schönhage, Grotefeld, and Vetter in [7, page 213].

3. Quotients

Let f, g, h be power series with fg = 1. This section considers the problem of computing hg given h and f.

Brent suggested computing g mod t^{n}, as discussed above, and then multiplying g mod t^{n} by h mod t^{n}. This takes time comparable to 8.6 transforms if g mod t^{n} is computed with an order-16 iteration. However, one can do better by reusing the intermediate results from the computation of g.

Write g = g_{0} + g_{1}t^{n} + g_{2}t^{2n} + · · ·. Section 2 shows how to compute g_{0} in time comparable to 5.6 transforms, and then g_{0}^{∗}, g_{1}^{∗}, g_{2}^{∗}, g_{3}^{∗}, g_{1}, g_{2}, g_{3}, g_{4}, g_{5}, g_{6}, g_{7} in 40 transforms. Then computing g_{4}^{∗}, g_{5}^{∗}, g_{6}^{∗}, g_{7}^{∗} takes just 4 more transforms, and computing hg mod t^{8n} takes 16 more transforms.

The upshot is that computing hg mod t^{n} takes time comparable to 8.2 transforms.

Logarithms. Let f be a power series with constant coefficient 1. Then one can compute log f = I(D(f)/f) in time comparable to 8.2 transforms.

4. Square roots

Let f, g, h be power series with fg = 1 and h = f^{2}. This section considers two problems: computing f and g given h; and computing f given h.

Write f = f_{0} + f_{1}t^{n} + f_{2}t^{2n} + …, g = g_{0} + g_{1}t^{n} + g_{2}t^{2n} + …, and h = h_{0} + h_{1}t^{n} + h_{2}t^{2n} + . . .. Brent suggested (among other techniques) computing f from the standard Newton iteration: f_{1}t^{n} = −(f_{0}^{2}−h)/2f_{0} mod t^{2n}; f_{2}t^{2n}+f_{3}t^{3n} = −((f_{0} + f_{1}t^{n})^{2} − h)/2(f_{0} + f_{1}t^{n}) mod t^{4n}; etc. However, there is some overlap between each step and the next. For example, one can reuse g_{0} = 1/f_{0} mod t^{n} in the computation of g_{0} + g_{1}t^{n} = 1/(f_{0} + f_{1}t^{n}) mod t^{2n}.

So define q_{i}, r_{i} as in section 2. Define s_{1} by s_{1}t^{n} = f_{0}^{2} − h_{0} − h_{1}t^{n}; s_{2} and s_{3} by s_{2}t^{2n} + s_{3}t^{3n} = (f_{0} + f_{1}t^{n})^{2} − h_{0} − h_{1}t^{n} − h_{2}t^{2n} − h_{3}t^{3n}; and so on. Then f_{1} = (−1/2)s_{1}g_{0} mod t^{n}, f_{2} + f_{3}t^{n} = (−1/2)(s_{2} + s_{3}t^{n})(g_{0} + g_{1}t^{n}) mod t^{2n}, etc.

Starting from f_{0}, g_{0}, h_{0}, h_{1}, one can compute f_{1} with 5 transforms: f_{0}^{∗}, f_{0}^{2}, s_{1}^{∗}, g_{0}^{∗}, s_{1}g_{0}. One can then compute g_{1} with 5 transforms: f_{1}^{∗}, f_{0}g_{0}, f_{1}g_{0}, (q_{1} + r_{1})^{∗}, (q_{1} + r_{1})g_{0}. One can then compute f_{2}, f_{3} with 7 transforms; then g_{2}, g_{3} with 10 transforms; and so on.

The work to compute both f mod t^{n} and g mod t^{n} recursively is comparable to 9 transforms with an order-4 iteration. The work to compute f mod t^{n} is comparable to 6.5 transforms.

Notes. Bailey in [1] reported a square root time comparable to 21 transforms.

5. Exponentials

Let f, g, h be power series with fg = 1 and f = exp h. This section considers the problem of computing f and g given h.

Write f = f_{0} + f_{1}t^{n} + f_{2}t^{2n} + ⋯, g = g_{0} + g_{1}t^{n} + g_{2}t^{2n} + ⋯, h = h_{0} + h_{1}t^{n} + h_{2}t^{2n} + ⋯, D(f) = a_{0} + a_{1}t^{n} + a_{2}t^{2n} + ⋯, and D(h) = b_{0} + b_{1}t^{n} + b_{2}t^{2n} + · · ·.

Define q_{i}, r_{i} as in section 2. Define s_{1} by s_{1}t^{n} = h − I(a_{0}g_{0} − b_{0}q_{1}t^{n}) mod t^{2n}; define s_{2}, s_{3} by

s_{2}t^{2n} + s_{3}t^{3n} = h − I((a_{0} + a_{1}t^{n})(g_{0} + g_{1}t^{n}) − (b_{0} + b_{1}t^{n})(q_{2}t^{2n} + q_{3}t^{3n})) mod t^{4n};

and so on. Then f_{1} = f_{0}s_{1} mod t^{n}; f_{2} + f_{3}t^{n} = (f_{0} + f_{1}t^{n})(s_{2} + s_{3}t^{n}) mod t^{2n}; etc. XXX still have to check these run times Starting from f_{0}, g_{0}, h_{0}, h_{1}, one can compute f_{1} with 10 transforms: f_{0}^{∗}, g_{0}^{∗}, f_{0}g_{0},

a_{0}^{∗}, a_{0}g_{0}, b_{0}^{∗}, q_{1}^{∗}, b_{0}q_{1}, s_{1}^{∗}, f_{0}s_{1}. One can then compute g_{1} with 4 transforms: f_{1}^{∗}, f_{1}g_{0}, (q_{1} + r_{1})^{∗}, (q_{1} + r_{1})g_{0}. One can then compute f_{2}, f_{3} with 15 transforms; then g_{2}, g_{3} with 8 transforms; and so on.

The work to compute both f mod t^{n} and g mod t^{n} recursively is comparable to 37/3 transforms with an order-4 iteration. The work to compute f mod t^{n} is comparable to 31/3 transforms.

6. Implementation results

XXX This is a first draft! But presumably the timings will be in line with the theoretical estimates shown below.

recip: (9, 6) 5.6 quo or log: (12, 9) 8.2 sqrt: (16.5) 6.5 sqrt and isqrt: (16.5) 9 exp: (22) 31/3 I also still have to analyze “natural” order-3 and order-4 iterations.

References

[1] David H. Bailey, The computation of π to 29,360,000 decimal digits using Borweins’ quartically convergent algorithm, Mathematics of Computation 50 (1988), 283–296.

[2] Daniel J. Bernstein, Detecting perfect powers in essentially linear time, to appear; available

   through http://pobox.com/~djb/papers.html.

[3] Daniel J. Bernstein, Multidigit multiplication for mathematicians, submitted for publication;

   available through http://pobox.com/~djb/papers.html.

[4] Daniel J. Bernstein, Bounding smooth integers (extended abstract), available through http://

   pobox.com/~djb/papers.html.

[5] Richard P. Brent, Multiple-precision zero-finding methods and the complexity of elementary function evaluation, in [9], 151–176.

[6] H. T. Kung, On computing reciprocals of power series, Numerische Mathematik 22 (1974), 341–348.

[7] Arnold Schönhage, Andreas F. W. Grotefeld, Ekkehart Vetter, Fast algorithms: a multitape Turing machine implementation, Bibliographisches Institut, Mannheim, 1994.

[8] Malte Sieveking, An algorithm for division of powerseries, Computing 10 (1972), 153–156.

[9] J. F. Traub, Analytic computational complexity, Academic Press, New York, 1976.

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

E-mail address: djb@pobox.com