Third Order Convergence for Reciprocal
09-19-2021, 04:14 PM (This post was last modified: 09-20-2021 11:18 PM by Albert Chan.)
Post: #1
 Albert Chan Senior Member Posts: 1,659 Joined: Jul 2018
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*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
09-19-2021, 04:47 PM (This post was last modified: 09-19-2021 07:30 PM by Albert Chan.)
Post: #2
 Albert Chan Senior Member Posts: 1,659 Joined: Jul 2018
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 (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 end

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(1-1e-16, true)
1.0100000000000002
1.0000010000000001
1.0000000000000002
09-19-2021, 05:59 PM
Post: #3
 Albert Chan Senior Member Posts: 1,659 Joined: Jul 2018
RE: Third Order Convergence for Reciprocal
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 end

Convergence 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
09-20-2021, 04:33 AM
Post: #4
 lyuka Junior Member Posts: 48 Joined: May 2014
RE: Third Order Convergence for Reciprocal
"Higher-order convergence algorithm for reciprocal and square root"
http://www.finetune.co.jp/~lyuka/technot.../sqrt.html
09-20-2021, 10:14 AM
Post: #5
 ttw Member Posts: 241 Joined: Jun 2014
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)=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...
09-20-2021, 01:20 PM
Post: #6
 Albert Chan Senior Member Posts: 1,659 Joined: Jul 2018
RE: Third Order Convergence for Reciprocal
(09-20-2021 04:33 AM)lyuka Wrote:  "Higher-order convergence algorithm for reciprocal and square root"
http://www.finetune.co.jp/~lyuka/technot.../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) end

lua> 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
09-22-2021, 10:32 AM (This post was last modified: 09-22-2021 11:30 AM by Albert Chan.)
Post: #7
 Albert Chan Senior Member Posts: 1,659 Joined: Jul 2018
RE: Third Order Convergence for Reciprocal
(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
09-25-2021, 09:39 PM
Post: #8
 Namir Senior Member Posts: 810 Joined: Dec 2013
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)