Post Reply 
Third Order Convergence for Reciprocal
09-19-2021, 04:14 PM (This post was last modified: Yesterday 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
Yesterday, 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
Yesterday, 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
Yesterday, 01:20 PM
Post: #6
RE: Third Order Convergence for Reciprocal
(Yesterday 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

(Yesterday 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
Post Reply 




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