Third Order Convergence for Reciprocal

[Albert Chan]

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

[Albert Chan]

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

[Albert Chan]

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+ε) + ε⁴[(1+ε) + ε²(1+ε)] + ...

Unlike Newton's or Halley's, sum is not self-correcting. Instead, we add 1 at the very end.

1/(1-ε) - 1 = (ε+ε²) + ε²(ε+ε²) + ε⁴[(ε+ε²) + ε²(ε+ε²)] + ...

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

[lyuka]

Just for your information,
"Higher-order convergence algorithm for reciprocal and square root"
http://www.finetune.co.jp/~lyuka/technot.../sqrt.html

[ttw]

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...

[Albert Chan]

(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

[Albert Chan]

(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

[Namir]

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