         Путь: Математика » Быстрые вычисления » Обратное число, корень m-й степени
Точное вычисление обратного числа, корня m-й степени  High precision inverse and square roots  Статья любезно предоставлена:
web: numbers.computation.free.fr

Inverting a number of size n can be done in time O(n2) with the classical school division. This time cost can be reduced to O(nlog(n)) when using a fast multiplication, thanks to a Newton process. The same type of techniques apply to the n-th root computation of a number. For those operations we also describe some very efficient algorithms with high order of convergence (that is quartic, quintic, sextic, ... algorithms).

1.  High precision inversion  In this section we want to compute the inverse 1/A of the number A to a great accuracy.

### 1.1  Newton's iteration

Starting from an approximation x0 of 1/A, we apply Newton iteration  to the function f(x) = 1/x-A, which gives the approximating sequence

 xn+1 = xn- f(xn)fў(xn) = 2xn-Axn2.     (*)
It's a property of Newton iteration to converge quadratically near the root, that is the number of right digits doubles at each iteration. Thus, a number of log(n) iterations suffices to have n digits of 1/A.

The iteration (*) have a nice feature : the only needed operations are the product by 2, the difference and the product of high precision numbers. The Newton convergence being quadratic, the precision needed to compute iteration (*) is not the full final precision but only with the precision wanted at step n.

As an example, we compute the inverse of p, starting from the approximation x0 = 0.31831 (5 digits). The first three iterations become

 x1
 =
 2x0-px02 = 0.3183098861 (operations done with 10 digits of precision)
 x2
 =
 2x1-px12 = 0.31830988618379067151 (20 digits of precision)
 x3
 =
 2x2-px22 = 0.3183098861837906715377675267450287240665 (40 digits)

The first 38 digits of x3 agree with those of 1/p.

The cost of this process is of the order of the multiplication cost : since each iteration involves two multiplications on high precision numbers, to find 1/A we need to compute 2 multiplications of size n (last step), 2 multiplications of size n/2 (previous step), 2 multiplications of size n/4 , etc. Using for example an FFT multiplication, with cost O(nlogn), the final cost is of order

 2nlog(n)+2 n2 log жз и n2 цч ш +2 n4 log жз и n4 цч ш +ј Ј 2n жз и 1+ 12 + 14 +ј цч ш log(n) = 4nlog(n).

Note : a better iteration is obtained by writing (*) in the form

 hn
 =
 1-Axn
 xn+1
 =
 xn+xnhn.
Since cancellations occurs in hn = 1-Axn, the second multiplication xnhn can be done by taking into account only the non zero part of hn. This has also another advantage : the number of vanishing terms in front of hn gives the approximation precision of xn.

A division B/A is performed by multiplying B by the inverse of A.

### 1.2  A cubic iteration

From the general Householder's iteration 

 xn+1 = xn- f(xn)fў(xn) жз и 1+ f(xn)fўў(xn)2fў(xn)2 цч ш ,

applied on f(x) = 1/x-A, we find

 hn
 =
 1-Axn
 xn+1
 =
 xn+xn(hn+hn2).

Starting with a good initial point x0, the convergence is cubical, that is, the number of good digits is multiplied by 3 at each step of the process. (The second form of this iteration is given in ). Since hn tends to zero, the new term hn2 can be computed at a low cost.

As an example, we compute again, the inverse of p, starting from the approximation x0 = 0.31831 (5 digits). The two first iterations are

 x1
 =
 0.3183098861837906715(523...),(19 good digits)
 x2
 =
 0.3183098861837906715377675267450287240689...,(58 good digits)

and x3has 174 good digits...

### 1.3  High order iterations

Quartic iteration

It's possible to find a quartic iteration, just apply twice Newton's iteration on f(x) = 1/x-A, we find

 hn
 =
 1-Axn
 xn+1
 =
 xn+xn( hn+hn2+hn3) .

If we take a look at our example (inverse of p), starting as usual with x0 = 0.31831, we find that x1 has 26 good digits, x2 is 103 digits exact, x3 has 413 good digits ... The sequence of correct digits is then for this example

 26,103,413,1650,6601,26405,...

Quintic iteration

A quintic iteration founded by the application of a general quintic modified iteration (here, we omit the details of the demonstration which is deduced from the essay on Newton's methods) is given by

 hn
 =
 1-Axn
 xn+1
 =
 xn+xn( hn+hn2+hn3+hn4) ,

the rate of convergence is impressive, again on the previous example, x1 has 32 correct digits and x2 has 161 good digits, x3 has 806 good digits ...

Higher order

The inverse of a number can be computed with an iteration at any desired order, it's possible to show that the general algorithm of order r is given by

 hn
 =
 1-Axn
 xn+1
 =
 xn+xnPr-1(hn),

where Pm(u) is a polynomial of degree m, given by

 Pm(u)
 =
 m е k = 1 akuk
 ak
 =
 1

For example, for m = 3,m = 4 and m = 5 (respectively quartic, quintic, sextic iterations)

 P3(u)
 =
 u+u2+u3,
 P4(u)
 =
 u+u2+u3+u4,
 P5(u)
 =
 u+u2+u3+u4+u5,
 ...

In practical cases the number hn tends to zero, therefore the evaluation of the powers hnk is more and more efficient when k increases. A careful implementation of those high order iterations is very efficient.

2.  Square roots  In this section we apply the same technique to compute square roots.

### 2.1  Newton's iteration

Using the function f(x) = x2-A, Newton iteration gives

 xn+1 = xn- f(xn)fў(xn) = xn2 + A2xn .
(In other words, the iteration is the mean value of xn and A/xn).

This iteration has one disadvantage : one should compute the division A/xn. When one wants to compute 1/ЦA, a better iteration is obtained from the function f(x) = 1/x2-A, yielding the iteration

 xn+1 = 32 xn- 12 Axn3.     (**)
This formula only requires multiplications ! An easiest way of computing a square root ЦA is then to compute B = 1/ЦA from this latest process and then to perform the product A×B = ЦA.

Note : as for the inverse, the best way to compute the iteration (**) is to write it in the form

 hn
 =
 1-Axn2
 xn+1
 =
 xn+xn hn2 ,
since cancellations occur in hn = 1-Axn2.

### 2.2  A cubic iteration

Again, from the Householder's iteration, applied on the function f(x) = 1/x2-A, we find after some easy manipulations

 xn+1 = xn8 (15-10Axn2+3A2xn4).

This process converges cubically to the inverse of the square root of the number A (This form of the iteration is also given in ). Just like for the quadratic iteration, a more efficient way is to write it

 hn
 =
 1-Axn2
 xn+1
 =
 xn+xn жз и 12 hn+ 38 hn2 цч ш .

### 2.3  High order iterations

Just like for the inverse, it's possible to find a quartic iteration, just apply twice Newton's iteration on f(x) = 1/x2-A :

 xn+1 = xn- xn16 (1-Axn2)(-20+19Axn2-8A2xn4+A3xn6).

This algorithm converges quartically to 1/ЦA.

It is interesting to compare this result to the direct application of the general quartic modified iteration to our function (see the essay on Newton's methods)

 xn+1 = xn- xn16 (1-Axn2)(-19+16Axn2-5A2xn4).

This easier relation suggests that it may be more efficient to use the quartic iteration than twice the Newton iteration. The best form of this quartic algorithm is

 hn
 =
 1-Axn2
 xn+1
 =
 xn+xn жз и 12 hn+ 38 hn2+ 516 hn3 цч ш .

As an application of the modified iterations, we provide the general expression of some high order iterations to compute square roots. The general pattern of a high order process is given by

 hn
 =
 1-Axn2
 xn+1
 =
 xn+xnP(hn)

with P(u) being a polynomial with rational coefficients.

Quintic iteration

 P(u) = жз и 12 u+ 38 u2+ 516 u3+ 35128 u4 цч ш .

Sextic iteration

 P(u) = жз и 12 u+ 38 u2+ 516 u3+ 35128 u4+ 63256 u5 цч ш .

Septic iteration

 P(u) = жз и 12 u+ 38 u2+ 516 u3+ 35128 u4+ 63256 u5+ 2311024 u6 цч ш .

Octic iteration

 P(u) = жз и 12 u+ 38 u2+ 516 u3+ 35128 u4+ 63256 u5+ 2311024 u6+ 4292048 u7 цч ш .

The coefficient of the polynomials P(u) are in fact given by the series expansion of

1
 ___Ц1-u
-1
 =
 Ґ е k = 1 akuk
 ak
 =
 2(2k-1)!k!(k-1)!4k

this allow to define easily an algorithm at any order to compute square roots.

3.  m-th roots  The same techniques may be used in order to compute the m-th root (m being an integer) of the number A. Therefore, to find with a great accuracy the number A-1/m, use Newton's method on the function f(x) = 1/xm-A and this will produce the process

 hn
 =
 1-Axnm
 xn+1
 =
 xn+xn hnm
which converges quadratically to A-1/m. To find A1/m just compute at the end of the N previous iterations
 yN = AxNm-1.

Householder's iteration also becomes

 hn
 =
 1-Axnm
 xn+1
 =
 xn+xn жз и 1m hn+ 1+m2m2 hn2 цч ш

which converges cubically to A-1/m.

### 3.1  High order iteration

The general pattern for any order iteration is given by

 hn
 =
 1-Axnm
 xn+1
 =
 xn+xnP(hn)

and to find P(u) you need to compute the following series expansion

 (1-u)-1/m-1 = 1m u+ 1+m2m2 u2+ (1+m)(1+2m)3!m3 u3+...

The truncated series at order k will provide an algorithm of order k+1. Observe that for the value m = 1, we retrieve the polynomials used to compute the inverse of a number and for m = 2, it gives back the algorithms to find square roots.

#### 3.1.1  Example

For m = 4 the following algorithm converges quartically to A-1/4

 hn
 =
 1-Axn4
 xn+1
 =
 xn+xn жз и 14 hn+ 532 hn2+ 15128 hn3 цч ш

and A1/4 is given by at the end of the process by

 Axn3.

References  I. Newton, Methodus fluxionum et serierum infinitarum, (1664-1671)


A. S. Householder, The Numerical Treatment of a Single Nonlinear Equation, McGraw-Hill, New York, (1970)

J.M. Borwein and P.B. Borwein, ''Pi and the AGM - A study in Analytic Number Theory and Computational Complexity'', A Wiley-Interscience Publication, New York, (1987)