C/C++中双重快速反平方

您确定需要双重精度吗？您可以轻松地牺牲精度：

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

在这种情况下,1.0f绝对是强制性的,如果使用1.0,那么您将获得双精度.
您是否尝试在编译器上启用“马虎”数学？在GCC上,您可以使用-offast-math,其他编译器也有类似的选项.对于您的应用程序来说,数学可能不够好. (编辑：我没有看到结果汇编中有任何区别.)
>如果您使用GCC,是否考虑使用-mrecip？有一个“相互估计”功能只有大约12位的精度,但是要快得多.您可以使用Newton-Raphson方法来提高结果的精度. -mrecip选项将导致编译器为您自动生成相互估计和Newton-Raphson步骤,但如果要精细调整性能 – 精度折衷,您可以随时自行编写程序集. (牛顿 – 拉夫逊收敛非常快)(编辑：我无法让GCC生成RCPSS,见下文)

我发现一篇博客文章(source)讨论了您正在经历的确切问题,作者的结论是,像Carmack方法这样的技术与RCPSS指令(GCC使用的-mrecip标志)不具竞争力.

分区可能如此缓慢的原因是因为处理器通常只有一个分割单元,而且通常不是流水线的.因此,您可以在管道中同时执行几个乘法运算,但是在上一个分区完成之前不能分配.

不工作的技巧

> Carmack的方法：在现代处理器上已经过时了,它们具有相互的估计操作码.对于倒数,我看到的最好的版本只能提供一点精度 – 与RCPSS的12位相比,没有什么.我认为这是巧合,这个伎俩对于相互平方根有好的效果;巧合不大可能重复.
>重新标记变量.就编译器而言,1.0 /(x * x)和double x2 = x * x之间的差异很小. 1.0 / X2.如果您发现一个编译器为两个版本生成不同的代码,即使达到最低级别的优化,我也会感到惊讶.
>使用pow. pow库函数是一个完整的怪物.随着GCC的关闭数学关闭,图书馆的电话是相当昂贵的.随着GCC的数学打开,您将获得与pow(x,-2)完全相同的汇编代码,就像1.0 /(x * x)一样,所以没有任何好处.

更新

以下是双精度浮点值的反平方的牛顿 – 拉夫逊逼近的例子.

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

不幸的是,在我的电脑上使用RECIP_ITER = 1的基准测试,它比简单版本1.0 /(x * x)略慢(约5％).它的速度更快(2x快),但是只有12位的精度.我不知道12位是否足够你.

我认为这里的一个问题是这个微型优化太小了;在这个规模上,编译器作者与装配黑客几乎相同.也许如果我们有更大的图景,我们可以看到一种方法来加快速度.

例如,你说过 – 数学造成了不合理的精度损失;这可能表示您正在使用的算法中的数值稳定性问题.利用算法的正确选择,可以通过float而不是double来解决许多问题. (当然,你可能只需要24位以上,我不知道.)

如果要并行计算其中的几个,我怀疑RCPSS方法是否发光.