Third Order Convergence for Reciprocal

[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