Third Order Convergence for Reciprocal

[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