2

我正在编写一个解决平面受限三体问题的程序。其方程如下。此函数计算位置和速度的导数并将它们写入数组。

valarray<double> force(double t, valarray<double> r)
{
    valarray<double> f(dim);
    valarray<double>r0(r-rb0);
    valarray<double>r1(r-rb1);      

    f[0]=   2 * r[1] + r[2] - (1 - mu)*r0[2]/norm3(r0) - mu*r1[2]/norm3(r1);
    f[1]= - 2 * r[0] + r[3] - mu*r0[3]/norm3(r0) - mu*r1[3]/norm3(r1);
    f[2] = r[0];
    f[3] = r[1];
    return f;
}

double norm3(valarray<double> x)
{
    return pow(x[2]*x[2]+x[3]*x[3],1.5);
}

所以我必须计算位置向量的平方，然后将其提高到 3/2 的幂。我认为这些操作占用了很大一部分计算时间。

现在我使用 math.h 的 pow 函数。是否有另一种更快的算法来计算这种能力？我尝试使用快速反平方根（稍后再进行立方），但它为我的目的提供了太不精确的值并且工作时间更长（可能是因为立方）。

谢谢！

4

2 回答 2

5

一个简单的方法可能是尝试 x*sqrt(x)，但可以肯定地对其进行基准测试。

double norm3(valarray<double> x)
{
    double result=x[2]*x[2]+x[3]*x[3];
    result=result * sqrt(result);
    return result;
}
于 2013-03-09T09:36:34.143 回答
1

FSQRTFamily 15h AMD64 处理器需要 52 个周期。SSE2 变体将 29 用于标量值，将 38 用于打包操作。C 版本的sqrt()可能是一些额外的指令，但我怀疑它更多。

如果您想要相对精确的结果，我怀疑从其他一些操作中可以获得更好的结果。最有可能的是，任何产生良好精度的涉及pow()和exp()等log()的东西都需要更长的时间。

但是，在互联网上询问并不能超过您自己的基准。如果这是您代码的关键部分，请尝试一些不同的变体，看看您会得到什么。

于 2013-03-09T09:57:50.520 回答