Takayuki HOSODA

Nov. 16, 2009

The Newton-Raphson method is a first-order approximation of a function f(x) around the approximate value of x

Its convergence speed can be increased by higher order approximation.

Japanese edition is here.

XThe significant digit is_{n+1}= 2 * X_{n}- A * X_{n}^{2}

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

hRecurrence formula for 1 / A of_{n}= 1 - A * x_{n}

x_{n+1}= x_{n}* (1 + h_{n}+ h_{n}^{2})

The number of significant digits is tripled with each repetition. h_{n}is equal to the error of x_{n}.

hRecurrence formula for 1 / A of_{n}= 1 - A * x_{n}

x_{n+1}= x_{n}* (1 + h_{n}) * (1 + h_{n}^{2})

h_{n}= 1 - A * x_{n}

x_{n+1}= x_{n}* (1 + (1 + h_{n}^{2}) * (h_{n}+ h_{n}^{2}))

XHowever, this formula has the disadvantage of having to perform each A / X_{n+1}= (X_{n}+ A / X_{n}) / 2

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

XThat is, the above equation (since 3/2 and 1/2 are constants) only requires multiplication._{n+1}= X_{n}* (3 - A * X_{n}^{2}) / 2

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

hRecurrence formula for 1 / √A of_{n}= 1 - A * x_{n}^{2}

x_{n+1}= x_{n}* (1 + h_{n}* (4 + 3 * h_{n}) / 8)

The number of significant digits is tripled with each repetition. h_{n}is equal to the error of x_{n}.

hRecurrence formula for 1 / √A of_{n}= 1 - A * x_{n}^{2}

x_{n+1}= x_{n}* (1 + h_{n}* (8 + h_{n}* (6 + 5 * h_{n})) / 16)

hRecurrence formula for 1 / √A of_{n}= 1 - A * x_{n}^{2}

x_{n+1}= x_{n}* (1 + (64 * h_{n}+ h_{n}^{2}* (48 + 40 * h_{n}+ 35 * h_{n}^{2})) / 128)

h_{n}= 1 - A * x_{n}^{2}

x_{n+1}= x_{n}* (1 + h_{n}* (1/2 + h_{n}* (3/8 + h_{n}* (5/16 + h_{n}* (35/128 + h_{n}* 63/256)))))

Recurrence formula for

h_{n}= 1 - A / x_{n}^{2}

x_{n+1}= x_{n}* (1 - h_{n}* (1/2 + h_{n}* (1/8 + h_{n}* (1/16 + h_{n}* (5/128 + h_{n}* 7/256)))))

Recurrence formula for

h_{n}= 1 - A * x_{n}^{3}

x_{n+1}= x_{n}+ x_{n}* h_{n}* (1/3 + h_{n}* (2/9 + h_{n}* (14/81 + h_{n}* (35/243 + h_{n}* 91/729))))

Recurrence formula for

h_{n}= 1 - A * x_{n}^{4}

x_{n+1}= x_{n}+ x_{n}* h_{n}* (1/4 + h_{n}* (5/32 + h_{n}* (15/128 + h_{n}* (195/2048 + h_{n}* 663/8192))))

Example of a C program with fourth-order convergence to find √x*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

#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)); }^{-16}, three calculations are enough.

SEE ALSO:

- Formula for calculating π by arithmetic mean

- Algorithm for finding the square-root, cubic-root or hypotenuse of an integer
- Higher-order convergence algorithm for hypotenuse