Third Order Convergence for Square Roots Using Newton's Method

08272019, 06:32 PM
(This post was last modified: 08282019 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 pseudocode for the second algorithm is: Given N and X (initial guess) and tolerance toler: Code:
Enjoy! Namir 

08272019, 07:49 PM
Post: #2




RE: Third Order Convergence for Square Roots Using Newton's Method
(08272019 06:32 PM)Namir Wrote: I stumbled on a third order converging algorithm in an book about ODEs and PDEs. The algorithm is: Quote:Given N and X (initial guess) and tolerance toler: Hi Namir, Could you clarify what variable A is used for? Thanks! — Ian Abbott 

08272019, 09:44 PM
(This post was last modified: 08272019 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.


08272019, 09:55 PM
(This post was last modified: 08282019 03:53 AM by Albert Chan.)
Post: #4




RE: Third Order Convergence for Square Roots Using Newton's Method
(08272019 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 3nbit precise answer, all 3 factors required 3nbits This may be better, correction terms requiring only 2nbits: Halley's formula for \(\Large \sqrt N: x _{new} = x  {2x(x^2N) \over 3x^2 + N}\) Example: N=12345, √N = 111.10 80555 13540 51124 ... Try 5digits decimal guess = 111.11 , √N → 111.10 80555 13540 66013 Try 16bits binary guess = 0xde37p9, √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^3N) \over 2x^3 + N}\) 

08282019, 01:39 AM
Post: #5




RE: Third Order Convergence for Square Roots Using Newton's Method
(08272019 07:49 PM)ijabbott Wrote:(08272019 06:32 PM)Namir Wrote: I stumbled on a third order converging algorithm in an book about ODEs and PDEs. The algorithm is: A is a typo .. should be N. 

08282019, 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.


08282019, 01:52 AM
Post: #7




RE: Third Order Convergence for Square Roots Using Newton's Method
(08282019 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 

08282019, 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 HP67?
Regards, Bob 

09162021, 01:39 PM
Post: #9




RE: Third Order Convergence for Square Roots Using Newton's Method
(08272019 06:32 PM)Namir Wrote: X(n+1) = X(n)*(X(n)^2 + 3*N)(3*X(n)^2 + 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 bitshift (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) + (NX(n)^2)/2 / X(n) Apply this again, we can avoid calculating (NX(n+1)^2) with this equivalent formula: X(n+2) = X(n) + (NX(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) + (NX(n)^2)/2 * (X(n)+X(n+1)) / (N+X(n)^2) Example, guess 7 for √51 ≈ 7.14142842854285 (NX(0)^2)/2 = (5149)/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) 

09162021, 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 higherorder derivatives disappear; on the other hand, this makes memoryless highorder methods tricky.)
To get Sqrt(S), one uses two auxiliary variables, A and B. A=(SX^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. 

09162021, 07:45 PM
(This post was last modified: 09162021 07:48 PM by C.Ret.)
Post: #11




RE: Third Order Convergence for Square Roots Using Newton's Method
(08282019 07:32 PM)bshoring Wrote: Could we have an example of this in RPN, say for HP67? You are welcome. Here is a code with three subroutines that compute square root using the first three methods respectively:
Code: 001 31 25 11 f LBL A 010 31 25 12 f LBL B 024 31 25 13 f LBL C 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, ... 

09162021, 07:58 PM
Post: #12




RE: Third Order Convergence for Square Roots Using Newton's Method
(09162021 05:18 PM)ttw Wrote: To get Sqrt(S), one uses two auxiliary variables, A and B. Very compact formula I believe there is a typo, should be A = (SX^2)/(2X) S = 2XA + X^2 = X*(2A+X) = (BA)*(B+A) Newton's from B: (B + S/B)/2 = (B^2 + (BA)*(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 + (SX*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 overestimate (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=BA)) (A^2+2*B^2)/(2*B) (09162021 01:39 PM)Albert Chan Wrote: X(n+2) = X(n) + (NX(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=BA)) (A^2+2*B^2)/(2*B) 

09172021, 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/thread13750.html Gamo 

09172021, 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({2n1 \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*n1)/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 As expected, for n>1, approx formula overestimated 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*2e4) lua> approx(2, 2e4) * 50/7 7.141428428542851 

09182021, 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 noninteger powers).
Namir 

09202021, 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 256bit 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 biginteger or bigreal 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. 

« Next Oldest  Next Newest »

User(s) browsing this thread: