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

Question

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?

Answer 1

per a comment, my math was wrong.

Updated math:

Y = 255 / sqrt(x)
255 / y = sqrt(x)
(255 / y)^2 = x

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

This results in the following code, which works

def inverse_sqrt_newton_255(x, iterations=5):
    y = 255. / x
    for _ in range(iterations):
        # works:
        # y = y - (255**2 / y**2 - x) / (-2 * 255**2 / y**3)

        # works, but dividing by a constant
        y = y - ((x * y**3) - (255**2 * y)) / (2 * 255**2)
    return y

Answer 2

Compare your good code with bad code:

Good code:

    y = y * (1.5 - 0.5 * x * y * y)

Bad code:

    y = (1.5 * y) - (x * y * y / 510.)

Looks like badly placed parentheses.

---

A programmer's way to transform your good code into what you need is to introduce two types of y:

y1 = 1 / sqrt(x)
y2 = 255 / sqrt(x)

Then y2 = y1 * 255; invert that to get y1 = y2 / 255.

The code which works:

y1_new = y1_old * (1.5 - 0.5 * x * y1_old * y1_old)

Replace y1 by its expression:

y2_new / 255 = y2_old / 255 * (1.5 - 0.5 * x * y2_old / 255 * y2_old / 255)

Multiply by 255 and get rid of "new" and "old" subscripts:

y2 = y2 * (1.5 - 0.5 * x * y2 * y2 / 255)