Third Order Convergence for Reciprocal

[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