“提高”单位长度向量长度的快速方法答案

【问题标题】：Fast way to "improve" the length of a unit length vector“提高”单位长度向量长度的快速方法
【发布时间】：2016-04-30 03:18:41
【问题描述】：

当知道向量已经几乎是单位长度时，为性能关键代码中的完整向量归一化付费似乎很浪费。

有谁知道一种快速实用的方法可以使双精度 3D 向量的长度更接近 1？我正在想象一种基于 Newton-Raphson 迭代或 1 左右的有限泰勒展开的迭代方法。

Here 是一个实际的现实情况，这样的例程可能很有用。 incoming 向量已经几乎是单位长度，但如果没有明确的规范化，它仍然会触发断言。

使用 SSE 2、SSE 4.2 或 AVX 内部函数是可以的。

【问题讨论】：

我不太确定一个小平方根是否如此昂贵。但无论如何，在 x=1 附近，x^2 的图具有导数 2，即它就像一条线，每一个单位在 x 方向上向前上升 2 个单位。因此，您可以通过减去 1、除以 2 并加 1 来估计该区域中的 x（从欧几里得距离的平方）。我还没有检查过这会如何在性能方面发挥作用，只有您可以知道您的“已经几乎”已经足够接近了。
尝试和基准测试，然后询问详细信息。 SO 不是讨论论坛。
请注意，@Cheersandhth.-Alf 建议的结果是x 附近的平方根函数的一阶泰勒展开，尽管显然是通过不同的途径达到的。您可以通过添加更多项来改进它（下一个是 -((x-1)^2)/8，但我预计您很快就会通过 sqrt() 函数更有效的点。还要注意，您可以使用泰勒级数误差界限项来确定这些近似值是否足够接近您。
@JohnBollinger：谢谢，我不这么认为。但是现在你提到它，我记得麦克劳林级数是通过考虑导数来实现的，而泰勒级数是通过考虑麦克劳林与偏移量来实现的？很久以前，我几乎什么都不记得了，除了当一位讲师在与学生进行某种交流时，我在课堂上开始推导出泰勒级数时，我感到多么困惑。我想知道他到底是怎么做到的，这有多棒。现在，我认为他可能已经记住了这个场合...... :)
我应该说sqrt(x)在1附近的扩展，当然。 “Near x”没有任何意义。

标签： c++ c optimization precision intrinsics

【解决方案1】：

手头的问题归结为找到（近似）reciprocal square root。

SSE 和 AVX 包含一个近似倒数平方根机器指令 rsqrt，它特别适合于此。根据原始AMD64 Architecture Programmer's Manual, volume 1，倒数平方根变体的最大相对误差最多为 1.5×2^-12，或小于 0.0004。

如果你使用GCC，你可以使用__builtin_ia32_rsqrtss()SSE内置函数计算向量平方长度的倒数平方根，并将向量分量乘以结果，得到一个“几乎单位”向量。

请注意，SSE 和 AVX 都提供了加速平方长度计算以及将每个分量相乘的函数。（不过，您需要将比例因子复制到大小相等的向量中。）

没有 SSE/AVX，一般的问题是我们希望向量分量乘以 f(S) ≃ sqrt(1/S) == 1/ sqrt(S)，其中S为向量与自身的内积（点积），即其长度的平方；但是 sqrt() 被认为太慢了，并且已知 S 已经接近 1。

在我们认为“接近 1”的范围内，任何值介于 1 和 sqrt(1/S) 之间的函数 f(S) 都可以工作。我能想到的最简单的函数是 f(S) = (C + 1 - S) / C我>。对于 S = 0.5² 到 2²（即对于长度在 1/2 和 2 之间的向量），C 是 6。

如果我们没有任何硬件支持倒数平方根，我会尝试的第一个近似值将遵循以下几行：

计算向量的平方长度S
计算 M = 0.125 * (9 - S)

请注意，任何常量对 C₁ 和 C₂ = 1 + 1 / C ₁ 应该可以工作，只是收敛的范围和速率会有所不同。我为这个例子选择了 C₁ = 1/8 仅仅是因为它在 IEEE-754 浮点表示中是精确的，并且通常乘法比除法快得多。其他值（例如我上面提到的 0.5 到 2 范围内的 1/6）是不精确的，可能需要手动调整（在两个常数中以一种或另一种方式调整最低有效单位）。
将向量的每个分量乘以 M。

如果这没有产生足够好的结果，我会停止担心，而是使用（硬件）平方根。（在某些架构上，将平方长度转换为单精度以计算比例因子可以显着提高速度。但是，在 x86/AMD64 上却没有。）

【讨论】：

感谢您的详细回答，非常感谢。实际上，我通过简单的 Newton-Raphson 迭代得到了很好的结果：M = (3 - S) * 0.5。这似乎比 M = (9 - S) * 0.125 收敛得更快。
经常让人们认为 sqrt 很慢的一件事是他们没有将-fno-math-errno 与 gcc（或其他编译器上的任何等价物）一起使用，因此编译器会生成对 libc 函数的昂贵调用测试退化情况，而不是直接使用硬件指令。
@MarcGlisse 并不是sqrt 很慢；只是它比几个加法和乘法慢。我在帖子中链接的软件在一次运行过程中计算了数千亿平方根。每一点优化都有帮助。使用一次 Newton-Raphson 迭代调整接近单位长度的 3D 双精度向量的长度似乎比普通归一化快约 2.7 倍。
@MarcGlisse：完全同意。此外，GCC 仍然不能很好地向量化操作，因此如果您可以将向量归一化为 2、4、8 或 16 个组（取决于精度以及是否使用 AVX256 或更早版本），您通常可以在大致相同的时间内完成你可以标准化一个。（如果您使用 GCC 的向量内在函数。）也适用于 rsqrt。较新的 GCC 版本应该能够矢量化 sqrt(v) 与 v 是任何矢量类型，而不依赖于特定于架构的内置插件，但我最近没有检查结果。
接近 1 我相信 M = (3 - S)/2 是最好的近似值。它将 1e-5 (10⁻⁵) 的相对误差降低到大约 -1.5e-10，并且它确实具有确保向量长度的有趣特性（在将所有分量乘以 M 之后) 小于一（如果向量的原始长度在 1-1e-5 和 1+1e-5 之间）。