Higher-order convergence algorithm for reciprocal and square root

example of sqrt(2) by Newton-Raphson method
Takayuki HOSODA
Nov. 16, 2009

The Newton-Raphson method is a first-order approximation of a function f(x) around the approximate value of xn .
Its convergence speed can be increased by higher order approximation.

Japanese Japanese edition is here.

Reciprocal

The recurrence formula for reciprocal 1 / A by Newton-Raphson method is Given by
xn+1 = 2 xn - A xn2
The significant digit is Doubles with each iteration.

Since the calculation time of multiplication is dominant in the calculation with a large number of digits, in order to reduce the number of multiplications required to obtain the calculation result of the required number of digits, it is desirable to use a recurrence formula that converges in a higher order.

Recurrence formula for 1 / A of cubic convergence
hn = 1 - A xn
xn+1 = xn (1 + hn + hn2)

The number of significant digits is tripled with each repetition. hn is equal to the error of xn .
Recurrence formula for 1 / A of fourth order convergence
hn = 1 - A xn
xn+1 = xn (1 + hn) (1 + hn2)
Recurrence formula for 1 / A of fifth order convergence
hn = 1 - A xn
xn+1 = xn (1 + (1 + hn2) (hn + hn2))

Square root

A similar method can be applied to the square root calculation. The recurrence formula by the Newton-Raphson method for √A (square root of A) is given by,
xn+1 = (xn + A / xn) / 2
However, this formula has the disadvantage of having to perform each A / xn division, which requires a lot of computation time.
Then, when we consider the Newton-Raphson method of 1 / √A, we find that a better recurrence formulas are obtained.

Recurrence formula for 1 / √A by Newton-Raphson method
xn+1 = xn (3 - A xn2) / 2
That is, the above equation (since 3/2 and 1/2 are constants) only requires multiplication.
To find √A, find 1 / √A and then √A = A (1 / √A) to find √A.
To find 1 / √A, it is better to use the recurrence formula of higher-order convergence, as in the case of the reciprocal calculation above.
Here we use a series expansion of 1 / √(1 - h).

Recurrence formula for 1 / √A of cubic convergence
hn = 1 - A xn2
xn+1 = xn (1 + hn (4 + 3 hn) / 8)

The number of significant digits is tripled with each repetition. hn is equal to the error of xn .
Recurrence formula for 1 / √A of quartic convergence
hn = 1 - A xn2
xn+1 = xn (1 + hn (8 + hn (6 + 5 hn)) / 16)
Recurrence formula for 1 / √A of 5th-order convergence
hn = 1 - A xn2
xn+1 = xn (1 + (64 hn + hn2 (48 + 40 hn + 35 hn2)) / 128)
Recurrence formula for 1 / √A of 6th-order convergence
hn = 1 - A xn2
xn+1 = xn (1 + hn (1/2 + hn (3/8 + hn (5/16 + hn (35/128 + hn 63/256)))))

Appendix ‐ Other recurrence formulas of 6th-order convergence

As in the case of 1 / √A, we use the series expansion of √(1 - h) .

Recurrence formula for √A of 6th-order convergence
hn = 1 - A / xn2
xn+1 = xn (1 - hn (1/2 + hn (1/8 + hn (1/16 + hn (5/128 + hn 7/256)))))

Recurrence formula for 1 / ∛A (= A-1/3 ) of 6th-order convergence
hn = 1 - A xn3
xn+1 = xn + xn hn (1/3 + hn (2/9 + hn (14/81 + hn (35/243 + hn 91/729))))

Recurrence formula for 1 / ∜A (= A-1/4 ) of 6th-order convergence
hn = 1 - A xn4
xn+1 = xn + xn hn (1/4 + hn (5/32 + hn (15/128 + hn (195/2048 + hn 663/8192))))

A program example

Example of a C program with fourth-order convergence to find √x 
#include <stdio.h>
#include <math.h>
#ifndef DBL_EPSILON
#define DBL_EPSILON     2.2204460492503131E-16
#endif

double
_sqrtinv(a)
    double a;
{
    double x, h, g;
    int    e;

    /* enhalf the exponent for a half digit initial accuracy */
    frexp(a, &e);
    x = ldexp(1.0, -e >> 1);

    /* 1/sqrt(a), fourth-order convergence */
    g = 1.0;
    while(fabs( h = 1.0 - a * x * x) < fabs(g)) {
        x += x * (h * (8.0 + h * (6.0 + 5.0 * h)) / 16.0);
        g = h;
    }
    return(x);
}

double
sqrt(a)
    double(a);
{
    if(a < 0){
        fprintf(stderr, "sqrt: domain error\n");
        return(0);
    }
    return(a * _sqrtinv(a));
}
*note1 : It is better to expand the loop as current processors can calculate the contents while determining the jump condition of the loop. If the machine ε is about 10-16, three calculations are enough.
SEE ALSO:
www.finetune.co.jp [Mail]