[Python] How can I approximate 255/sqrt(x) using Newton's method?

I am trying to approximate 255 / sqrt(x) using Newton's method to avoid using division or the sqrt() operation on an architecture where those operations are expensive.

My derivation is:

y = 255 / sqrt(x)
255 / y^2 = x
f(y) = 255 / y^2 - x = 0
y[n+1] = y[n] - f(y) / f’(y)
f’(y) = -510 / y^3
y[n+1] = y[n] - (255 / y[n]^2 - x) / (-510 / y[n] ^ 3)
y[n+1] = 3/2 * y[n] - x * y[n]^2 / 510


I tested this with the following code

# what I want is 255 * isqrt(x) or 255 / sqrt(x)
def inverse_sqrt_newton_255(x, iterations=5):
y = 255. / x
for _ in range(iterations):
y = (1.5 * y) - (x * y * y / 510.)
return y


However, I find that the result is always off by a factor of sqrt(x)

inverse_sqrt_newton_255(25) = 10.2 (should be 51)

I could obviously multiply the result by sqrt(x) but the whole point is to avoid the sqrt() operation.

For reference, the following works correctly for computing 1 / sqrt(x):

def inverse_sqrt_newton(x, iterations=5):
y = 1.0 / x
for _ in range(iterations):
y = y * (1.5 - 0.5 * x * y * y)
return y


What am I doing wrong?

Continue reading...