9

最近我在分析一个程序，其中的热点肯定是这个

double d = somevalue();
double d2=d*d;
double c = 1.0/d2   // HOT SPOT

之后不使用值 d2，因为我只需要值 c。前段时间我读过关于快速反平方根的卡马克方法，显然不是这种情况，但我想知道类似的算法是否可以帮助我计算 1/x^2。

我需要相当准确的精度，我已经检查过我的程序没有使用 gcc -ffast-math 选项给出正确的结果。(g++-4.5)

4

3 回答 3

19

做快速平方根之类的技巧通过牺牲精度来获得它们的性能。（嗯，他们中的大多数。）

  1. 你确定你需要double精度吗？您可以很容易地牺牲精度：

    double d = somevalue();
float c = 1.0f / ((float) d * (float) d);

    在这种1.0f情况下，绝对是强制性的，如果您使用它，1.0您将获得double精确度。

  2. 您是否尝试过在编译器上启用“草率”数学？在 GCC 上您可以使用-ffast-math，其他编译器也有类似的选项。对于你的应用程序来说，草率的数学可能已经足够好了。（编辑：我没有看到生成的程序集有任何区别。）

  3. 如果您使用的是 GCC，您是否考虑过使用 GCC -mrecip？有一个“倒数估计”函数，它只有大约 12 位的精度，但速度要快得多。您可以使用 Newton-Raphson 方法来提高结果的精度。该-mrecip选项将导致编译器自动为您生成倒数估计和 Newton-Raphson 步骤，但如果您想微调性能-精度权衡，您始终可以自己编写程序集。（Newton-Raphson 收敛速度非常快。）（编辑：我无法让 GCC 生成 RCPSS。见下文。）

我发现一篇博客文章（来源）讨论了您正在经历的确切问题，作者的结论是，像 Carmack 方法这样的技术与 RCPSS 指令（-mrecipGCC 上的标志使用）没有竞争力。

除法如此缓慢的原因是处理器通常只有一个除法单元，而且它通常不是流水线的。因此，您可以在管道中同时执行一些乘法运算，但在前一个除法完成之前不能发出除法。

不起作用的技巧

  1. Carmack 的方法：它在现代处理器上已经过时了，它具有互易估计操作码。对于倒数，我见过的最好的版本只给出了一点精度——与RCPSS. 我认为这个技巧对倒数平方根非常有效，这是一个巧合。一个不太可能重演的巧合。

  2. 重新标记变量。1.0/(x*x)就编译器而言，和之间几乎没有区别double x2 = x*x; 1.0/x2。如果你发现一个编译器可以为两个版本生成不同的代码，并且优化到最低级别，我会感到惊讶。

  3. 使用pow. 库pow函数是一个彻头彻尾的怪物。-ffast-math关闭GCC后，库调用相当昂贵。打开GCC后，您-ffast-math将获得与 完全相同的汇编代码，因此没有任何好处。pow(x, -2)1.0/(x*x)

更新

这是双精度浮点值的平方反比的 Newton-Raphson 近似示例。

static double invsq(double x)
{
    double y;
    int i;
    __asm__ (
        "cvtpd2ps %1, %0\n\t"
        "rcpss %0, %0\n\t"
        "cvtps2pd %0, %0"
        : "=x"(y)
        : "x"(x));
    for (i = 0; i < RECIP_ITER; ++i)
        y *= 2 - x * y;
    return y * y;
}

不幸的是，RECIP_ITER=1在我的计算机上进行基准测试时，它比简单版本稍慢（~5%）1.0/(x*x)。零迭代速度更快（快 2 倍），但您只能获得 12 位的精度。我不知道12位对你来说是否足够。

我认为这里的问题之一是微优化太小了。在这个规模上，编译器编写者与汇编黑客几乎处于同等地位。也许如果我们有更大的图景，我们可以找到一种让它更快的方法。

例如，您说这-ffast-math导致了不希望的精度损失；这可能表明您正在使用的算法存在数值稳定性问题。通过正确选择算法，很多问题都可以用float代替来解决double。（当然，您可能只需要超过 24 位。我不知道。）

RCPSS如果您想并行计算其中的几个，我怀疑该方法会发光。

于 2012-03-16T11:57:17.507 回答
5

是的，您当然可以尝试解决问题。让我给你一些大致的想法，你可以填写细节。

首先，让我们看看为什么 Carmack 的 root 有效：

我们用通常的方式写x  =  M  × 2 E。现在回想一下，IEEE 浮点数存储了一个偏差的指数偏移量：如果e表示指数字段，我们有 e = Bias +  E  ≥ 0。重新排列，我们得到E  = e - Bias。

现在对于逆平方根：x -1/2  =  M -1/2  × 2 - E /2。新的指数字段是：

       e'  = 偏差 -  E /2 = 3/2 偏差 - e/2

通过位摆弄，我们可以通过移位从e中得到e /2的值，而 3/2 Bias 只是一个常数。

此外，尾数M存储为 1.0 +  x且x  < 1，我们可以将M -1/2近似为 1 + x/2。同样，只有x以二进制形式存储的事实意味着我们可以通过简单的位移得到除以 2。


现在我们看x -2：这等于M -2  × 2 -2 E，我们正在寻找一个指数场：

       e'  = 偏差 - 2  E  = 3 偏差 - 2  e

同样，3 Bias 只是一个常数，您可以 通过位移从e得到 2 e。至于尾数，您可以将 (1 + x) -2近似为 1 - 2  x，因此问题简化为 从x获得 2 x。


请注意，Carmack 的魔术浮点运算实际上并没有立即计算结果：相反，它产生了非常准确的估计值，用作传统迭代计算的起点。但是因为估计很好，你只需要很少的几轮后续迭代就可以得到一个可以接受的结果。

于 2012-03-16T11:57:09.930 回答
1

对于您当前的程序，您已经确定了热点 - 很好。作为加速 1/d^2 的替代方法，您可以选择更改程序，使其不经常计算 1/d^2。你能把它吊出一个内循环吗？对于多少个不同的 d 值，您计算 1/d^2？你能预先计算出你需要的所有值，然后查找结果吗？这对于 1/d^2 来说有点麻烦，但如果 1/d^2 是更大代码块的一部分，那么将这个技巧应用于此可能是值得的。你说如果你降低精度，你就得不到足够好的答案。有什么方法可以改写代码，这可能会提供更好的行为？数值分析非常微妙，可能值得尝试一些事情并看看发生了什么。

当然，理想情况下，您会发现一些基于多年研究的优化例程——lapack 或 linpack 中有什么可以链接到的吗？

于 2012-03-16T18:21:34.200 回答