C/C++中的快速反转平方

Question

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

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

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

我需要非常准确的精度,我已经检查过我的程序没有使用gcc -ffast-math选项给出正确的结果.(克++ - 4.5)

Answer 1

做快速平方根等的技巧通过牺牲精度来获得它们的性能.(好吧,大多数人.)

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

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

   该1.0f是在这种情况下绝对强制性的,如果你使用的1.0不是你会得到double精度.

2. 您是否尝试在编译器上启用"草率"数学？在GCC上你可以使用-ffast-math,其他编译器有类似的选项.草率的数学运算可能足以满足您的应用需求.(编辑:我没有看到生成的程序集有任何不同.)

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

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

究其原因为何师可以如此缓慢是因为处理器一般只有一个划分单元,它的往往不是流水线.因此,您可以在管道中进行一些全部同时执行的乘法运算,但在前一个分段完成之前不能发布任何除法.

不起作用的技巧

1. Carmack的方法:它在现代处理器上已经过时,它具有相互估计的操作码.对于倒数,我见过的最好的版本只能提供一点精度 - 与12位相比没什么RCPSS.我认为这个技巧对于互惠的平方根很有效,这是巧合; 一个不太可能重复的巧合.

2. 重新标记变量.就编译器而言,1.0/(x*x)和之间几乎没有什么区别double x2 = x*x; 1.0/x2.如果您发现编译器为两个版本生成不同的代码,并且优化已打开,即使是最低级别,我会感到惊讶.

3. 用pow.该pow库函数是一个总的怪物.-ffast-math关闭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如果你想并行计算其中的几个,这个方法就会闪耀.

Answer 2

是的,你当然可以尝试解决问题.让我给你一些一般性的想法,你可以填写详细信息.

首先,让我们看看为什么Carmack的根工作原理:

我们用通常的方式写x  =  M  ×2 ^E. 现在回想一下,在IEEE浮点存储指数由偏置偏移:如果é表示的指数领域,我们有E =偏差+  Ë  ≥0花事,我们得到ë  = E -偏差.

现在对于平方根: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个偏差只是一个常数,你可以通过位移从e获得2  e.作为尾数,可以近似(1 + ^x)的^-2由1 - 2  X,所以问题简化为获得2  X从X.

请注意,Carmack的魔术浮点小提琴实际上并不能正确计算结果:相反,它会产生非常准确的估计,用作传统迭代计算的起点.但由于估算非常好,您只需要几轮后续迭代即可得到可接受的结果.