Post Reply 
Third Order Convergence for Square Roots Using Newton's Method
08-27-2019, 06:32 PM (This post was last modified: 08-28-2019 01:41 AM by Namir.)
Post: #1
Third Order Convergence for Square Roots Using Newton's Method
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
Find all posts by this user
Quote this message in a reply
08-27-2019, 07:49 PM
Post: #2
RE: Third Order Convergence for Square Roots Using Newton's Method
(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!

— Ian Abbott
Find all posts by this user
Quote this message in a reply
08-27-2019, 09:44 PM (This post was last modified: 08-27-2019 09:45 PM by KeithB.)
Post: #3
RE: Third Order Convergence for Square Roots Using Newton's Method
It looks to be the root you are looking for, since you are comparing it to the guess * the guess.
Find all posts by this user
Quote this message in a reply
08-27-2019, 09:55 PM (This post was last modified: 08-28-2019 03:53 AM by Albert Chan.)
Post: #4
RE: Third Order Convergence for Square Roots Using Newton's Method
(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}\)
Find all posts by this user
Quote this message in a reply
08-28-2019, 01:39 AM
Post: #5
RE: Third Order Convergence for Square Roots Using Newton's Method
(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.
Find all posts by this user
Quote this message in a reply
08-28-2019, 01:40 AM
Post: #6
RE: Third Order Convergence for Square Roots Using Newton's Method
I corected my original post to replace A with N.
Find all posts by this user
Quote this message in a reply
08-28-2019, 01:52 AM
Post: #7
RE: Third Order Convergence for Square Roots Using Newton's Method
(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?

--Bob Prosperi
Find all posts by this user
Quote this message in a reply
08-28-2019, 07:32 PM
Post: #8
RE: Third Order Convergence for Square Roots Using Newton's Method
Could we have an example of this in RPN, say for HP-67?


Regards,
Bob
Find all posts by this user
Quote this message in a reply
09-16-2021, 01:39 PM
Post: #9
RE: Third Order Convergence for Square Roots Using Newton's Method
(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)
Find all posts by this user
Quote this message in a reply
09-16-2021, 05:18 PM
Post: #10
RE: Third Order Convergence for Square Roots Using Newton's Method
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.
Find all posts by this user
Quote this message in a reply
09-16-2021, 07:45 PM (This post was last modified: 09-16-2021 07:48 PM by C.Ret.)
Post: #11
RE: Third Order Convergence for Square Roots Using Newton's Method
(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, ...
Find all posts by this user
Quote this message in a reply
09-16-2021, 07:58 PM
Post: #12
RE: Third Order Convergence for Square Roots Using Newton's Method
(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 Smile
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)
Find all posts by this user
Quote this message in a reply
09-17-2021, 12:30 AM
Post: #13
RE: Third Order Convergence for Square Roots Using Newton's Method
This one used " Successive Approximation Algorithms "

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


Gamo
Find all posts by this user
Quote this message in a reply
09-17-2021, 01:33 AM
Post: #14
RE: Third Order Convergence for Square Roots Using Newton's Method
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
Find all posts by this user
Quote this message in a reply
09-18-2021, 10:15 PM
Post: #15
RE: Third Order Convergence for Square Roots Using Newton's Method
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
Find all posts by this user
Quote this message in a reply
09-20-2021, 10:04 AM
Post: #16
RE: Third Order Convergence for Square Roots Using Newton's Method
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.
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: 1 Guest(s)