Post Reply 
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
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
Find all posts by this user
Quote this message in a reply
09-19-2021, 04:47 PM (This post was last modified: 09-19-2021 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 (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
Find all posts by this user
Quote this message in a reply
09-19-2021, 05:59 PM
Post: #3
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
Find all posts by this user
Quote this message in a reply
09-20-2021, 04:33 AM
Post: #4
RE: Third Order Convergence for Reciprocal
Just for your information,
"Higher-order convergence algorithm for reciprocal and square root"
http://www.finetune.co.jp/~lyuka/technot.../sqrt.html
Find all posts by this user
Quote this message in a reply
09-20-2021, 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)=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...
Find all posts by this user
Quote this message in a reply
09-20-2021, 01:20 PM
Post: #6
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
Find all posts by this user
Quote this message in a reply
09-22-2021, 10:32 AM (This post was last modified: 09-22-2021 11:30 AM by Albert Chan.)
Post: #7
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
Find all posts by this user
Quote this message in a reply
09-25-2021, 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
Find all posts by this user
Quote this message in a reply
Post Reply 




User(s) browsing this thread: