Third Order Convergence for Square Roots Using Newton's Method

[Namir]

Hi All,

Most of us are familiar with the second order converging algorithm for obtaining the square root of N:

X(n+1) = (X(n) + N/X(n)) / 2
where X(0) is the initial guess for the square root of N.

I stumbled on a third order converging algorithm in an book about ODEs and PDEs. The algorithm is:


X(n+1) = X(n)*(X(n)^2 + 3*N)(3*X(n)^2 + N)
where X(0) is the initial guess for the square root of N.

I compared the two algorithms using Excel and the second one does converge faster than the first one.

The pseudo-code for the second algorithm is:

Given N and X (initial guess) and tolerance toler:

Code:


Repeat
  Y = X*X;
  X = X * (Y + 3*N) / (3*Y + N)
Until Abs(X*X-N) <= toler

Enjoy!

Namir

[ijabbott]

(08-27-2019 06:32 PM)Namir Wrote:  I stumbled on a third order converging algorithm in an book about ODEs and PDEs. The algorithm is:


X(n+1) = X(n)*(X(n)^2 + 3*N)(3*X(n)^2 + A)
where X(0) is the initial guess for the square root of N.

Quote:Given N and X (initial guess) and tolerance toler:

Code:


Repeat
  Y = X*X;
  X = X * (Y + 3*N) / (3*Y + A)
Until Abs(X*X-A) <= toler

Hi Namir,

Could you clarify what variable A is used for? Thanks!

[KeithB]

It looks to be the root you are looking for, since you are comparing it to the guess * the guess.

[Albert Chan]

(08-27-2019 06:32 PM)Namir Wrote:  X(n+1) = X(n)*(X(n)^2 + 3*N)/(3*X(n)^2 + N)

Hi, ijabbott. Typo fixed

However, this setup may not be ideal for arbitrary precision math.
For a 3n-bit precise answer, all 3 factors required 3n-bits

This may be better, correction terms requiring only 2n-bits:

Halley's formula for \(\Large \sqrt N: x _{new} = x - {2x(x^2-N) \over 3x^2 + N}\)

Example: N=12345, √N = 111.10 80555 13540 51124 ...

Try 5-digits decimal guess = 111.11 ,   √N → 111.10 80555 13540 66013
Try 16-bits binary guess = 0xde37p-9, √N → 111.10 80555 13540 50609

You can DIY other functions, see Method of obtaining Third Order Iterative Formulas

Example: Halley's formula for \(\Large \sqrt[3] N: x _{new} = x - {x(x^3-N) \over 2x^3 + N}\)

[Namir]

(08-27-2019 07:49 PM)ijabbott Wrote:  

(08-27-2019 06:32 PM)Namir Wrote:  I stumbled on a third order converging algorithm in an book about ODEs and PDEs. The algorithm is:


X(n+1) = X(n)*(X(n)^2 + 3*N)(3*X(n)^2 + A)
where X(0) is the initial guess for the square root of N.

Quote:Given N and X (initial guess) and tolerance toler:

Code:


Repeat
  Y = X*X;
  X = X * (Y + 3*N) / (3*Y + A)
Until Abs(X*X-A) <= toler

Hi Namir,

Could you clarify what variable A is used for? Thanks!

A is a typo .. should be N.

[Namir]

I corected my original post to replace A with N.

[rprosperi]

(08-28-2019 01:40 AM)Namir Wrote:  I corected my original post to replace A with N.

Albert's typo comments also inserted a division symbol in the middle. Does that belong there as well?

Also, does anyone know what the ":x" symbology means in front of the equations Albert proposed?

[bshoring]

Could we have an example of this in RPN, say for HP-67?

[Albert Chan]

(08-27-2019 06:32 PM)Namir Wrote:  X(n+1) = X(n)*(X(n)^2 + 3*N)(3*X(n)^2 + N)

where X(0) is the initial guess for the square root of N.

I wonder if this is more efficient than applying Newton's method twice.
Assuming X(n)^2 is calculated and saved, above required 2 add, 3 mul, 1 div

Quote:X(n+1) = (X(n) + N/X(n)) / 2

Newton's (applied twice): 2 add, 2 div, 2 bit-shift (or, replace /2 by *0.5, 2 mul)
But, 2 Newton's give quartic convergence !

---

Apply Newton's mentally, division part is not easy.
Rewrite Newton's method as corrections to guess maybe easier.

X(n+1) = X(n) + (N-X(n)^2)/2 / X(n)

Apply this again, we can avoid calculating (N-X(n+1)^2) with this equivalent formula:

X(n+2) = X(n) + (N-X(n)^2)/2 / harmonic_mean(X(n), X(n+1))

Simplify away harmonic_mean(), and 2*X(n)*X(n+1) = N+X(n)^2, we get:

X(n+2) = X(n) + (N-X(n)^2)/2 * (X(n)+X(n+1)) / (N+X(n)^2)

Example, guess 7 for √51 ≈ 7.14142842854285

(N-X(0)^2)/2 = (51-49)/2 = 1

X(1) = 7 + 1/7 = 7.(142857)
X(2) = 7 + 1*(7+(7+1/7)) / (51+49) = 7 + (14+1/7) / 100 = 7.14(142857)

For comparison, Newton's formula from X(1) to X(2):

X(2)
= (7+1/7) + (51-(7+1/7)^2)/(2*(7+1/7))
= (7+1/7) + (51-(49+2+1/49)) / (100/7)
= (7+1/7) + (-1/49)*7/100
= (7+1/7) − 1/700
= 7.14(142857)

[ttw]

Wikipedia gives a simple quartic method from the Bakhtshali Manuscript. It's nice in that it gives rational approximations and could be used on polynomials. It's actually two steps of Newton's method with algebraic simplifications. (Newton's method is nice for square roots because the higher-order derivatives disappear; on the other hand, this makes memoryless high-order methods tricky.)

To get Sqrt(S), one uses two auxiliary variables, A and B.

A=(S-X^2)/2S

B=X+A

X(new)=B-(A^2)/2B equivalent to (X+A)-(A^2)/2*(X+A)

I suppose one could algebraically compound two of these.

The interval of convergence [/quote]isn't given.

[C.Ret]

(08-28-2019 07:32 PM)bshoring Wrote:  Could we have an example of this in RPN, say for HP-67?

You are welcome.

Here is a code with three sub-routines that compute square root using the first three methods respectively:

1. The second order converging algorithm indicated by Namir aka the Héron method:
   \( x_{k+1} = \frac{x_i + \frac{N}{x_k}}{2} \)
2. The third order converging algorithm Namir stumbled on :
   \( x_{k+1} = x_n \frac{x_k^2+ 3.N}{3.x_k^2+N} \)
3. The Halley's formula proposed by Albert Chan:
   \( x_{k+1}=x_n-\frac{2.x_k.(x_k^2−N)}{3.x_k^2+N} \)
   
   
   

Code:

001   31 25 11  f LBL A      010   31 25 12  f LBL B      024   31 25 13  f LBL C
002      34 00    RCL 0      011      33 01    STO 1      025      33 01    STO 1     
003      35 52  h x:y        012   33 71 01    STO×1      026   33 71 01    STO×1        
004         81     ÷         013      34 01    RCL 1      027      34 01    RCL 1
005      35 82  h LSTx       014         03     3         028         02     2
006         61     +         015   33 71 01    STO×1      029   33 71 01    STO×1
007         02     2         016      34 00    RCL 0      030      35 52  h x:y
008         81     ÷         017   33 61 01    STO+1      031   33 61 01    STO+1
009      35 22  h RTN        018         71     ×         032      34 00    RCL 0
                             019         61     +         033   33 61 01    STO+1
                             020      34 01    RCL 1      034         51     - 
                             021         81     ÷         035         71     × 
                             022         71     ×         036         71     ×
                             023      35 22  h RTN        037      34 01    RCL 1
                                                          038         81     ÷
                                                          039         51     -
                                                          040      35 22  h RTN

These programs only run one step of each method. The value n for the expected square root have to be store in register R0 and the first guest or estimate has to be set in register X: of the stack. Each press on A, B or C key affine the estimation displayed in register X: by one step using the corresponding method.

For example :
51 STO 0
7 A display 7.142857145
Further press on A display 7.141428570, 7.141428430, 7.141428430, ...

where
7 B, B, B, ... successively display 7.141414140, 7.141428423, ...

and
7 C, C, C, ... successively display 7.141414141, 7.141428429, 7.141428428, ...

[Albert Chan]

(09-16-2021 05:18 PM)ttw Wrote:  To get Sqrt(S), one uses two auxiliary variables, A and B.

A=(S-X^2)/2S

B=X+A

X(new)=B-(A^2)/2B equivalent to (X+A)-(A^2)/2*(X+A)

Very compact formula
I believe there is a typo, should be A = (S-X^2)/(2X)

S = 2XA + X^2 = X*(2A+X) = (B-A)*(B+A)

Newton's from B: (B + S/B)/2 = (B^2 + (B-A)*(B+A))/(2B) = B - A^2/(2B)       QED

Using the same symbols, this is Halley's method for √S

X(new) = X + A/(1+A/(2X)) = X + A/C

C = 1 + A/(2X) = 1 + (S-X*X)/(2X)^2 = 3/4 + S/(2X)^2 ≥ 3/4

Assuming Halley's method is better than Newton's:

A > 0 → C > 1 → B > √S
A < 0 → C < 1 → B > √S

Newton's method for √S over-estimate (convergence always from the high side)

We can apply small correction to Halley's, for quartic convergence.

XCAS> simplify(subst(X + A/(1+A/(X+B)), X=B-A))

(-A^2+2*B^2)/(2*B)

(09-16-2021 01:39 PM)Albert Chan Wrote:  X(n+2) = X(n) + (N-X(n)^2)/2 * (X(n)+X(n+1)) / (N+X(n)^2)

Might as well proof this too. denominator (N+X(n)^2) = (2*X(n)*X(n+1))

X(new) = X + (X*A)*(X+B)/(2*X*B) = X + (A/B)*(X+B)/2

XCAS> simplify(subst(X + (A/B)*(X+B)/2, X=B-A))

(-A^2+2*B^2)/(2*B)

[Gamo]

This one used " Successive Approximation Algorithms "

https://www.hpmuseum.org/forum/thread-13750.html


Gamo

[Albert Chan]

Using XCas pade((1+n*x)^(1/n), 3, 3), I discovered a pattern:

\(\displaystyle \sqrt[n]{1+n\,x} = 1 + x \left(
\frac {1+\left({n+1 \over 6}\right) x}
{1 + \left({2n-1 \over 3}\right) x}
\right) + O(x^4) \)

For n = 1/2 or 1, O(x^4) can be dropped.

lua> exact = function(n,x) return (1+n*x)^(1/n) end
lua> approx = function(n,x) return 1+x*(1+(n+1)/6*x) / (1+(2*n-1)/3*x) end

lua> y = 0.1 -- test approx formula for (1+y)^(1/n)
lua> for i=2,12 do n=i/4; print(n, exact(n,y/n), approx(n,y/n)) end

Code:

0.5     1.2100000000000002      1.21
0.75    1.135508127002004       1.1355072463768117
1       1.1                     1.1
1.25    1.079230345298891       1.0792307692307692
1.5     1.0656022367666107      1.0656028368794326
1.75    1.0559733623179208      1.055974025974026
2       1.0488088481701516      1.0488095238095239
2.25    1.0432700717216588      1.0432707355242568
2.5     1.0388601182540846      1.038860759493671
2.75    1.0352658433367348      1.0352664576802508
3       1.0322801154563672      1.032280701754386

As expected, for n>1, approx formula over-estimated true root.

51/49 = 1 + 2/49
sqrt(51)/7 = sqrt(1 + 2*1/49)

lua> approx(2, 1/49) * 7
7.141428571428571

51*49 / 50^2 = 1 - 1/50^2
sqrt(51)*7/50 = sqrt(1 - 2*2e-4)

lua> approx(2, -2e-4) * 50/7
7.141428428542851

[Namir]

I am a fan of the Pade polynomials they add more flexibility to curve fitting. I have even integrated them with my Shammas Polynomias (non integer power polynomials) to yield a special version of the Pad polynomials (with non-integer powers).

Namir

[ttw]

I wrote a short 50g program for taking exact integer square roots of 256-bit numbers for the 50g. My tests showed it worked for even larger integers. (I got it from the Wikipedia article on Integer Square Roots.)

I'll post it later as I don't have the calculator handy right now. I was surprised at how quickly it works. The first step is an (approximate) conversion to reals which is good to about 36 bits. After the first iteration, it always converges from above so stopping is easy.

In many of the big-integer or big-real computations, the order of convergence isn't so important. The last step increases the multiplication time so much that most of the CPU time is used for the last multiplication or division. This particularly applies to the estimates of Pi.