Third Order Convergence for Square Roots Using Newton's Method - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Calculators (and very old HP Computers) (/forum-3.html) +--- Forum: General Forum (/forum-4.html) +--- Thread: Third Order Convergence for Square Roots Using Newton's Method (/thread-13524.html) Third Order Convergence for Square Roots Using Newton's Method - Namir - 08-27-2019 06:32 PM 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 RE: Third Order Convergence for Square Roots Using Newton's Method - ijabbott - 08-27-2019 07:49 PM (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! RE: Third Order Convergence for Square Roots Using Newton's Method - KeithB - 08-27-2019 09:44 PM It looks to be the root you are looking for, since you are comparing it to the guess * the guess. RE: Third Order Convergence for Square Roots Using Newton's Method - Albert Chan - 08-27-2019 09:55 PM (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}$$ RE: Third Order Convergence for Square Roots Using Newton's Method - Namir - 08-28-2019 01:39 AM (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. RE: Third Order Convergence for Square Roots Using Newton's Method - Namir - 08-28-2019 01:40 AM I corected my original post to replace A with N. RE: Third Order Convergence for Square Roots Using Newton's Method - rprosperi - 08-28-2019 01:52 AM (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? RE: Third Order Convergence for Square Roots Using Newton's Method - bshoring - 08-28-2019 07:32 PM Could we have an example of this in RPN, say for HP-67? RE: Third Order Convergence for Square Roots Using Newton's Method - Albert Chan - 09-16-2021 01:39 PM (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) RE: Third Order Convergence for Square Roots Using Newton's Method - ttw - 09-16-2021 05:18 PM 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. RE: Third Order Convergence for Square Roots Using Newton's Method - C.Ret - 09-16-2021 07:45 PM (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: The second order converging algorithm indicated by Namir aka the Héron method: $$x_{k+1} = \frac{x_i + \frac{N}{x_k}}{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}$$ 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, ... RE: Third Order Convergence for Square Roots Using Newton's Method - Albert Chan - 09-16-2021 07:58 PM (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) RE: Third Order Convergence for Square Roots Using Newton's Method - Gamo - 09-17-2021 12:30 AM This one used " Successive Approximation Algorithms " https://www.hpmuseum.org/forum/thread-13750.html Gamo RE: Third Order Convergence for Square Roots Using Newton's Method - Albert Chan - 09-17-2021 01:33 AM 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 RE: Third Order Convergence for Square Roots Using Newton's Method - Namir - 09-18-2021 10:15 PM 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 RE: Third Order Convergence for Square Roots Using Newton's Method - ttw - 09-20-2021 10:04 AM 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.