
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 x_{0} of 1/A, we apply Newton
iteration [1] to the function f(x) = 1/xA, which gives the
approximating sequence
x_{n+1} = x_{n} 
f(x_{n}) f^{ў}(x_{n})

= 2x_{n}Ax_{n}^{2}. (*) 

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 x_{0} = 0.31831 (5 digits). The first three iterations become



2x_{0}px_{0}^{2} = 0.3183098861 (operations done with 10 digits of precision) 
 

2x_{1}px_{1}^{2} = 0.31830988618379067151 (20 digits of precision) 
 

2x_{2}px_{2}^{2} = 0.3183098861837906715377675267450287240665 (40 digits) 

 

The first 38 digits of x_{3} 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 
n 2

log 
ж з
и


n 2

ц ч
ш

+2 
n 4

log 
ж з
и


n 4

ц ч
ш

+ј Ј 2n 
ж з
и

1+ 
1 2

+ 
1 4

+ј 
ц ч
ш

log(n) = 4nlog(n). 

Note : a better iteration is obtained by writing (*) in the form
Since cancellations occurs in h_{n} = 1Ax_{n}, the second multiplication x_{n}h_{n} can be done by taking into account only the non zero part of h_{n}. This has also another advantage : the number of vanishing terms in
front of h_{n} gives the approximation precision of x_{n}.
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 [2]
x_{n+1} = x_{n} 
f(x_{n}) f^{ў}(x_{n})


ж з
и

1+ 
f(x_{n})f^{ўў}(x_{n}) 2f^{ў}(x_{n})^{2}

ц ч
ш

, 

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



 

x_{n}+x_{n}(h_{n}+h_{n}^{2}). 

 

Starting with a good initial point x_{0}, 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 [3]).
Since h_{n} tends to zero, the new term h_{n}^{2} can be computed at a
low cost.
As an example, we compute again, the inverse of p, starting from the
approximation x_{0} = 0.31831 (5 digits). The two first iterations are



0.3183098861837906715(523...),(19 good digits) 
 

0.3183098861837906715377675267450287240689...,(58 good digits) 

 

and x_{3}has 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/xA, we find



 

x_{n}+x_{n}( h_{n}+h_{n}^{2}+h_{n}^{3}) . 

 

If we take a look at our example (inverse of p), starting as usual with
x_{0} = 0.31831, we find that x_{1} has 26 good digits, x_{2} is 103
digits exact, x_{3} 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



 

x_{n}+x_{n}( h_{n}+h_{n}^{2}+h_{n}^{3}+h_{n}^{4}) , 

 

the rate of convergence is impressive, again on the previous example, x_{1}
has 32 correct digits and x_{2} has 161 good digits, x_{3} 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



 

x_{n}+x_{n}P^{r1}(h_{n}), 

 

where P^{m}(u) is a polynomial of degree m, given by
For example, for m = 3,m = 4 and m = 5 (respectively quartic, quintic, sextic
iterations)



 

 

u+u^{2}+u^{3}+u^{4}+u^{5}, 
 


 

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


In this section we apply the same technique to compute square roots.
2.1 Newton's iteration
Using the function f(x) = x^{2}A, Newton iteration gives
x_{n+1} = x_{n} 
f(x_{n}) f^{ў}(x_{n})

= 
x_{n} 2

+ 
A 2x_{n}

. 

(In other words, the iteration is the mean value of x_{n} and A/x_{n}).
This iteration has one disadvantage : one should compute the division A/x_{n}. When one wants to compute 1/ЦA, a better iteration is
obtained from the function f(x) = 1/x^{2}A, yielding the iteration
x_{n+1} = 
3 2

x_{n} 
1 2

Ax_{n}^{3}. (**) 

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
since cancellations occur in h_{n} = 1Ax_{n}^{2}.
2.2 A cubic iteration
Again, from the Householder's iteration, applied on the function f(x) = 1/x^{2}A, we find after some easy manipulations
x_{n+1} = 
x_{n} 8

(1510Ax_{n}^{2}+3A^{2}x_{n}^{4}). 

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



 

x_{n}+x_{n} 
ж з
и


1 2

h_{n}+ 
3 8

h_{n}^{2} 
ц ч
ш

. 

 

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/x^{2}A :
x_{n+1} = x_{n} 
x_{n} 16

(1Ax_{n}^{2})(20+19Ax_{n}^{2}8A^{2}x_{n}^{4}+A^{3}x_{n}^{6}). 

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)
x_{n+1} = x_{n} 
x_{n} 16

(1Ax_{n}^{2})(19+16Ax_{n}^{2}5A^{2}x_{n}^{4}). 

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



 

x_{n}+x_{n} 
ж з
и


1 2

h_{n}+ 
3 8

h_{n}^{2}+ 
5 16

h_{n}^{3} 
ц ч
ш

. 

 

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
with P(u) being a polynomial with rational coefficients.
Quintic iteration
P(u) = 
ж з
и


1 2

u+ 
3 8

u^{2}+ 
5 16

u^{3}+ 
35 128

u^{4} 
ц ч
ш

. 

Sextic iteration
P(u) = 
ж з
и


1 2

u+ 
3 8

u^{2}+ 
5 16

u^{3}+ 
35 128

u^{4}+ 
63 256

u^{5} 
ц ч
ш

. 

Septic iteration
P(u) = 
ж з
и


1 2

u+ 
3 8

u^{2}+ 
5 16

u^{3}+ 
35 128

u^{4}+ 
63 256

u^{5}+ 
231 1024

u^{6} 
ц ч
ш

. 

Octic iteration
P(u) = 
ж з
и


1 2

u+ 
3 8

u^{2}+ 
5 16

u^{3}+ 
35 128

u^{4}+ 
63 256

u^{5}+ 
231 1024

u^{6}+ 
429 2048

u^{7} 
ц ч
ш

. 

The coefficient of the polynomials P(u) are in fact given by the series
expansion of
this allow to define easily an algorithm at any order to compute square
roots.


The same techniques may be used in order to compute the mth 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/x^{m}A and this will produce the process
which converges quadratically to A^{1/m}. To find A^{1/m} just compute
at the end of the N previous iterations
Householder's iteration also becomes



 

x_{n}+x_{n} 
ж з
и


1 m

h_{n}+ 
1+m 2m^{2}

h_{n}^{2} 
ц ч
ш



 

which converges cubically to A^{1/m}.
3.1 High order iteration
The general pattern for any order iteration is given by
and to find P(u) you need to compute the following series expansion
(1u)^{1/m}1 = 
1 m

u+ 
1+m 2m^{2}

u^{2}+ 
(1+m)(1+2m) 3!m^{3}

u^{3}+... 

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}



 

x_{n}+x_{n} 
ж з
и


1 4

h_{n}+ 
5 32

h_{n}^{2}+ 
15 128

h_{n}^{3} 
ц ч
ш



 

and A^{1/4} is given by at the end of the process by
