卡马克/威尔什倒数平方根算法是否有偏差？

Question

在实现"Carmack's Inverse Square Root"算法时，我注意到结果似乎存在偏差。下面的代码似乎可以得到更好的结果：

float InvSqrtF(float x)
{
    // Initial approximation by Greg Walsh.
    int i  = * ( int* ) &x;
    i  = 0x5f3759df - ( i >> 1 );
    float y  = * ( float * ) &i;
    // Two iterations of Newton-Raphson's method to refine the initial estimate.
    x *= 0.5f;
    float f = 1.5F;
    y  = y * ( f - ( x * y * y ) );
    y  = y * ( f - ( x * y * y ) );
    * ( int * )(&y) += 0x13; // More magic.
    return y;
}

关键区别在于倒数第二行的“更多魔法”。由于初始结果过低，这个因素是相当恒定的，所以只需一条指令，就可以将19 * 2^(exponent(y)-bias)加到结果中。这似乎给我额外增加了3位，但我是否忽略了什么？

- MSalters

---

你的问题是什么？“19是正确的加数吗？”还是其他什么？ - Mats Petersson

这段代码不是注定要受到严格别名规则的未定义行为吗？ - comocomocomocomo

@comocomocomocomo 你可以安全地使用联合类型进行类型转换。 - Daniel Fischer

@MatsPetersson：不知道，目标是ARM。 - MSalters

@MatsPetersson N1570，6.5.2.3（3），脚注95：“如果用于读取联合体对象内容的成员与上次用于在对象中存储值的成员不同，则该值的对象表示的适当部分将被重新解释为新类型的对象表示，如6.2.6所述（这个过程有时称为“类型切换”）。除了陷阱表示外，根据标准，通过联合体进行类型切换是安全的[是的，脚注不是规范性的，但很明显它的意图是有效的]。 - Daniel Fischer

Answer 1

牛顿法会产生偏差。要找到零点的函数，

f(y) = x - 1/y²

这个方程是凹的，所以除非你从y ≥ √(3/x)开始，否则牛顿迭代法只产生近似值≤ 1/√x（而且精度比较差，除非你从精确结果开始）。

浮点运算有时会产生过大的近似值，但通常前两次迭代不会出现这种情况（因为初始猜测通常不够接近）。

所以，是的，存在偏差，添加一个小量通常可以改进结果。但并不总是如此。例如在1.25或0.85附近的区域，没有进行调整的结果比有调整的结果更好。在其他区域，调整可以获得额外一位的精度，而在另一些区域，则可以获得更多的精度。

在任何情况下，为了获得最佳结果，应该根据取值最频繁的x的区域来调整添加的魔术常数。

- Daniel Fischer

Answer 2

由于这种方法是一种近似，结果有时会高估，有时会低估。您可以在麦肯尼里的论文中找到一些关于不同配置下误差分布以及背后的数学知识的好图表。

因此，除非您有确凿的证据表明在您的应用领域中结果明显存在偏差，否则我建议按照Lomont的文档中建议的调整魔法常数 :-)

- dunadar

---

+1。真的很喜欢Lomont对于0x5f375a86常数的分析。 - Brett Hale

实际上，使用0x5f3759df和0x5f375a86我得到了完全相同的偏置。 - MSalters

好的，我维持我之前说的：如果在你特定的应用领域中你知道存在偏差，你可以进行补偿。Daniel Fisher的回答非常清楚 ^_^ - dunadar