+- 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 Reciprocal (/thread-17485.html)
Third Order Convergence for Reciprocal - Albert Chan - 09-19-2021 04:14 PM
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*x-1)/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*x-2)*(n*x-1)
Reciprocal, 3rd-order convergence: x ← x + x*(1-n*x)*(2-n*x)
Compare with 2nd-order Newton's: x ← x + x*(1-n*x)
We can estimate cost of computation: 3rd-order = 1.5× 2nd-order
But, 3rd-order run twice is 9th-order, 2nd-order run 3 times only get 8th.
XCAS> N(x) := x + x*(1-n*x)
XCAS> H(x) := x + x*(1-n*x)*(2-n*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
RE: Third Order Convergence for Reciprocal - Albert Chan - 09-19-2021 04:47 PM
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 (1-n*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*y-5.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) = 1E-18, below machine epsilon.
Code:
function rcp(x, verbose)
local y, p = frexp(x)
y = ldexp((2.586*y-5.819)*y+4.243, -p) -- 1/x guess
if verbose then print(y) end
p = 1-x*y; y = y+y*(p+p*p)
if verbose then print(y) end
p = 1-x*y; y = y+y*(p+p*p)
return y
endGuess 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(1-1e-16, true)
1.0100000000000002
1.0000010000000001
1.0000000000000002
RE: Third Order Convergence for Reciprocal - Albert Chan - 09-19-2021 05:59 PM
This version use a trick from thread Fun Math Algorithm.
(09-08-2020 09:59 PM)Albert Chan Wrote: We can reduce summing terms, from O(n), to O(ln(n))
1/(1-ε) = 1 + ε + ε² + ... = (1+ε) + ε²(1+ε) + ε4[(1+ε) + ε²(1+ε)] + ...
Unlike Newton's or Halley's, sum is not self-correcting. Instead, we add 1 at the very end.
1/(1-ε) - 1 = (ε+ε²) + ε²(ε+ε²) + ε4[(ε+ε²) + ε²(ε+ε²)] + ...
Code:
function rcp2(x, verbose)
local y, p = frexp(x)
if y < 2/3 then y, p = y+y, p-1 end -- y = [2/3, 4/3)
x, p = 1-y, 2^-p -- 1/(1-x) = 1 + x + x^2 + x^3 + ...
local z = x*x
x = x + z
if verbose then print(p+p*x) end
while z >= 0x1p-54 do
x = x + x*z
z = z * z
if verbose then print(p+p*x) end
end
return p+p*x
endConvergence is 2nd-order (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(1-1e-16, true)
1
RE: Third Order Convergence for Reciprocal - lyuka - 09-20-2021 04:33 AM
Just for your information,
"Higher-order convergence algorithm for reciprocal and square root"
http://www.finetune.co.jp/~lyuka/technote/fract/sqrt.html
RE: Third Order Convergence for Reciprocal - ttw - 09-20-2021 10:14 AM
There is an old method for matrix inverse that can be used with ordinary numbers. The idea is that 1/(1+x)=1-x-x^2-x^3.... Then this term is collapsed to (1-x)(1-x^2)(1-x^4)(1-x^8)... until x^(2k) is small. There's a similar formula for 1/(1-x).
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...
RE: Third Order Convergence for Reciprocal - Albert Chan - 09-20-2021 01:20 PM
(09-20-2021 04:33 AM)lyuka Wrote: "Higher-order convergence algorithm for reciprocal and square root"
http://www.finetune.co.jp/~lyuka/technote/fract/sqrt.html
Thanks, lyuka
If x is a guess of 1/n, n*x = 1-h:
1/n = x/(n*x) = x/(1-h) = 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 self-correcting, we need a more accurate calculation of h.
Code:
function rcp3(x)
local y, p = frexp(x)
x = 1-y
y = (2.586*x+0.647)*x+0.01 -- estimated 1/y-1
local h = x-y + x*y -- = 1-(1-x)*(1+y)
x = h * h
h = h + x -- h + h^2
h = h + h*x -- h + h^2 + h^3 + h^4
h = h + h*x*x -- h + h^2 + ... + h^8
return ldexp(y*h+h+y+1, -p)
endlua> t = {1/2, 5/8, 7/8, 1-1e-16}
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
(09-20-2021 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
RE: Third Order Convergence for Reciprocal - Albert Chan - 09-22-2021 10:32 AM
(09-20-2021 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)=1-x-x^2-x^3.... Then this term is collapsed to (1-x)(1-x^2)(1-x^4)(1-x^8)... until x^(2k) is small.
Some signs should be "+"
1/(1-x)
= 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 + ... = (1-x) * (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
RE: Third Order Convergence for Reciprocal - Namir - 09-25-2021 09:39 PM
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