Third Order Convergence for Reciprocal

09192021, 04:14 PM
(This post was last modified: 09202021 11:18 PM by Albert Chan.)
Post: #1




Third Order Convergence for Reciprocal
Inspired by Namir's thread: Third Order Convergence for Square Roots Using Newton's Method
If we define f(x) = n  1/x, Halleys' method won't work XCAS> f := n  1/x XCAS> factor(f/(f'  (f''/2)*(f/f')) → (n*x1)/n Correction involve division, which we can't do (otherwise, we just evaluate 1/n) Slightly modified f work. XCAS> f /= x XCAS> factor(f/(f'  (f''/2)*(f/f')) → x*(n*x2)*(n*x1) Reciprocal, 3rdorder convergence: x ← x + x*(1n*x)*(2n*x) Compare with 2ndorder Newton's: x ← x + x*(1n*x) We can estimate cost of computation: 3rdorder = 1.5× 2ndorder But, 3rdorder run twice is 9thorder, 2ndorder run 3 times only get 8th. XCAS> N(x) := x + x*(1n*x) XCAS> H(x) := x + x*(1n*x)*(2n*x) XCAS> [ H(H(x)), N(N(N(x))) ]  n=5/8, x=1. [1.59976536036, 1.59937429428] If we try 1 Halley + 1 Newton, regardless of order, we get the same result. XCAS> simplify( N(H(x))  H(N(x)) ) 0 

09192021, 04:47 PM
(This post was last modified: 09192021 07:30 PM by Albert Chan.)
Post: #2




RE: Third Order Convergence for Reciprocal
Newton's or Halley's reciprocal iteration formula required a good guess.
Both diverges if guess is off by 100% or more, due to factor (1n*x) To get a good guess, we decompose number to mantissa, exponents. With a small mantissa range, we can curve fit 1/y with a few points. If we give up accuracy at the edges, we get overall better estimate. Below curve fit 1/y with y = [0.5+ε, 0.75, 1.0ε], where ε = 1/30 1/y ≈ (2.586*y5.819)*y+4.243 // 0.5 ≤ y < 1, max_rel_err ≈ 1.0% Apply 3rd order correction 2 times, max_rel_err ≈ (0.01)^(3*3) = 1E18, below machine epsilon. Code: function rcp(x, verbose) Guess max_rel_err at x = 1/2, 5/8, 7/8, and slightly below 1.0 lua> rcp(0.5, true) 1.9800000000000004 1.999998 2 lua> rcp(5/8, true) 1.6162812500000001 1.6000016858668447 1.6 lua> rcp(7/8, true) 1.1312812500000002 1.142855955231404 1.1428571428571428 lua> rcp(11e16, true) 1.0100000000000002 1.0000010000000001 1.0000000000000002 

09192021, 05:59 PM
Post: #3




RE: Third Order Convergence for Reciprocal
This version use a trick from thread Fun Math Algorithm.
(09082020 09:59 PM)Albert Chan Wrote: We can reduce summing terms, from O(n), to O(ln(n)) Unlike Newton's or Halley's, sum is not selfcorrecting. Instead, we add 1 at the very end. 1/(1ε)  1 = (ε+ε²) + ε²(ε+ε²) + ε^{4}[(ε+ε²) + ε²(ε+ε²)] + ... Code: function rcp2(x, verbose) Convergence is 2ndorder (this use 2 mul + 1 add per loop, slightly less costly than Newton's) To speed up convergence, mantissa [1/2, 1) → [2/3, 4/3), thus ε ≤ 1/3 lua> r = rcp2(1/2, true)  mantissa shifted to 1 2 lua> r = rcp2(5/8, true)  mantissa shifted to 5/4 1.625 1.6015625 1.600006103515625 1.6000000000931323 1.6 lua> r = rcp2(7/8, true) 1.140625 1.142822265625 1.1428571343421936 1.1428571428571423 1.1428571428571428 lua> r = rcp2(11e16, true) 1 

09202021, 04:33 AM
Post: #4




RE: Third Order Convergence for Reciprocal
Just for your information,
"Higherorder convergence algorithm for reciprocal and square root" http://www.finetune.co.jp/~lyuka/technot.../sqrt.html 

09202021, 10:14 AM
Post: #5




RE: Third Order Convergence for Reciprocal
There is an old method for matrix inverse that can be used with ordinary numbers. The idea is that 1/(1+x)=1xx^2x^3.... Then this term is collapsed to (1x)(1x^2)(1x^4)(1x^8)... until x^(2k) is small. There's a similar formula for 1/(1x).
The point is that one computes x^2 (1 multiplication) and gets increasingly accurate approximations with each multiplication. I think it's of exponential order but I don't remember. Let's count: 3 multiplications for order 4, 5 multiplications for order 8, 7 multiplications for order 16, 9 multiplications for order 32... 

09202021, 01:20 PM
Post: #6




RE: Third Order Convergence for Reciprocal
(09202021 04:33 AM)lyuka Wrote: "Higherorder convergence algorithm for reciprocal and square root" Thanks, lyuka If x is a guess of 1/n, n*x = 1h: 1/n = x/(n*x) = x/(1h) = x*(1 + h + h² + ... ), where, h = 1  n*x This version does 1/n with 9th order convergence formula (correction, sum to h^8) Because sum is not selfcorrecting, we need a more accurate calculation of h. Code: function rcp3(x) lua> t = {1/2, 5/8, 7/8, 11e16} lua> for i=1,#t do x=t[i]; print(x, rcp3(x)) end 0.5 2 0.625 1.6 0.875 1.1428571428571428 0.9999999999999999 1.0000000000000002 (09202021 10:14 AM)ttw Wrote: The point is that one computes x^2 (1 multiplication) and gets increasingly accurate approximations with each multiplication. I think it's of exponential order but I don't remember. I just learn there is a name for this ... Estrin's scheme 

09222021, 10:32 AM
(This post was last modified: 09222021 11:30 AM by Albert Chan.)
Post: #7




RE: Third Order Convergence for Reciprocal
(09202021 10:14 AM)ttw Wrote: There is an old method for matrix inverse that can be used with ordinary numbers. The idea is that 1/(1+x)=1xx^2x^3.... Then this term is collapsed to (1x)(1x^2)(1x^4)(1x^8)... until x^(2k) is small. Some signs should be "+" 1/(1x) = 1 + x + x^2 + x^3 + ... = (1+x) * (1 + x^2 + x^4 + ...) = (1+x) * (1+x^2) * (1 + x^4 + x^8 + ...) = (1+x) * (1+x^2) * (1+x^4) * (1+x^8) ... Replace x by (x): 1/(1+x) = 1  x + x^2  x^3 + ... = (1x) * (1+x^2) * (1+x^4) * (1+x^8) ... Using the same idea, we can show rcp(x) and rcp3(x) are equivalent. Let n = 1ε, 1st application of Halley, with guess x = 1, gives x' = (1+ε+ε^2) 2nd application of Halley: ε' = 1  n*x' = 1  (1ε)*(1+ε+ε^2) = ε^3 rcp(1ε) = (1 + ε + ε^2) * (1 + ε^3 + ε^6) = 1 + ε + ε^2 + ... + ε^8 rcp3(1ε) = 1 + (ε+ε^2)*(1+ε^2)*(1+ε^4) = 1 + ε + ε^2 + ... + ε^8 

09252021, 09:39 PM
Post: #8




RE: Third Order Convergence for Reciprocal
Your algorithm for calculatingg the reciprocal is good for (old) machines that were slow in diving numbers. Any idea how multiplication and division compare on (a sample of) calculators and computers?
Namir 

« Next Oldest  Next Newest »

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