如何使用牛顿法近似 255/sqrt(x)？

Question

我尝试使用牛顿法来近似

255 / sqrt(x)

以避免在架构上使用除法或

sqrt()

操作，因为这些操作的成本很高。

我的推导是：

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

我用以下代码对此进行了测试

# 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

但是，我发现结果总是相差

sqrt(x)

inverse_sqrt_newton_255(25) = 10.2

（应该是51）

我显然可以将结果乘以

sqrt(x)

，但重点是避免

sqrt()

操作。

仅供参考，以下计算正确有效：

1 / sqrt(x)

我做错了什么？

Answer 1

好的代码：

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

错误代码：

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

看起来括号位置不对。

程序员将好的代码转变成你需要的代码的方法是引入两种类型的 y：

---

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

则 y2 = y1 * 255；将其反转即可得到 y1 = y2 / 255。

有效的代码：

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

将 y1 替换为其表达式：

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

乘以255并去掉“新”和“旧”下标：

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