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 Newton's iteration. 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). Computing square roots and inverses to a great accuracy is also useful to compute many usual mathematical functions like logarithms and trigonometric functions (see Brent's remarkable paper [2]) or to compute constants like p,log(2),...([2], [6]).
In this section we want to compute the inverse 1/A of the number A to a great accuracy.
Starting from a few digits approximation x0 of 1/A, we apply Newton's iteration to the function f(x) = 1/x-A, which gives the approximating sequence
|
The iteration (*) have a nice feature : the only needed operations are the product by 2, the difference and the product of large 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
|
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 (that is xn×xn and A×xn2), 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
|
Note : a better iteration is obtained by writing (*) in the form
|
A division B/A is performed by multiplying B by the inverse of A.
A practical example of inverse of large numbers from FFT multiplication can be found as an easy C sample code at Easy programs for constants computation.
From the general Householder's iteration ([3] and [1] p. 216)
|
applied on f(x) = 1/x-A, we find
|
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 (A similar form of this iteration is given in [1] p. 216 and [4] p. 279). Since hn tends to zero, the new term hn2 can be computed at a lower cost.
As an example, we compute again, the inverse of p, starting from the approximation x0 = 0.31831 (5 digits). The first two iterations are
|
and x3 has 174 good digits...
Quartic iteration
It's possible to find a quartic iteration, just apply twice Newton's iteration on f(x) = 1/x-A, we find
|
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
|
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 our essay on Newton's methods) is given by
|
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 ... The second factorization of this iteration requires only 4 multiplications (A×xn,hn×hn,(1+hn2)×(hn+hn2),xn×(1+hn2)(hn+hn2)) at each step!
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
|
where Pm(u) is a polynomial of degree m, given by
|
For example, for m = 3,m = 4, m = 5, m = 6 and m = 7 (respectively quartic, quintic, sextic, septic and octic iterations)
|
In practical cases the number hn tends to zero, therefore the evaluation of the powers hnk is more and more efficient when k increases. Observe that some of the given factorizations may reduce the number of required multiplications. A careful implementation of those high order iterations may be very efficient.
The same technique as for the inverses applies for square roots.
Using the function f(x) = x2-A, Newton iteration gives
|
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
|
Note : as for the inverse, the best way to compute the iteration (**) is to write it in the form
|
Again, from the Householder's iteration, applied on the function f(x) = 1/x2-A, we find after some easy manipulations
|
This process converges cubically to the inverse of the square root of the number A (This form of the iteration is also given in [1] p. 217). Just like for the quadratic iteration, a more efficient way is to write it
|
Again like for the inverse, it's possible to find a quartic iteration, just apply twice Newton's iteration on f(x) = 1/x2-A :
|
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)
|
This simpler 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
|
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 iteration of order r is given by
|
with Pm(u) being a polynomial of degree m with rational coefficients. The first polynomials are given with factorization to save some multiplications, the best choice depends on your implementation.
Quintic iteration
|
Sextic iteration
|
Septic iteration
|
Octic iteration
|
For example to evaluate P7(u), you may compute in this order
|
that is, four large multiplications and don't forget that, in our cases, u is small so that those multiplications are even faster to perform.
Expression of the polynomials Pm(u)
The coefficients of the polynomials Pm(u) are in fact given by the series expansion of
|
this allow to define easily any order algorithm to compute square roots.
The previous 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, we use Newton's method on the function f(x) = 1/xm-A and this will produce the process
|
|
Householder's iteration also becomes
|
which converges cubically to A-1/m.
The general pattern for any order iteration is given by
|
and to find P(u) you need to compute the following series expansion
|
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.
For m = 4 the following algorithm converges quartically to A-1/4
|
and A1/4 is given by
|