2

平面制限3体問題を解くためのプログラムを書いています。その方程式は以下のとおりです。この関数は、位置と速度の導関数を計算し、それらを配列に書き込みます。

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);
}

したがって、位置ベクトルの2乗を計算してから、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

ファミリ15hAMD64FSQRTプロセッサでは52サイクルかかります。SSE2バリアントは、スカラー値の場合は29、パック操作の場合は38を取ります。のCバージョンsqrt()はおそらくいくつかの追加の指示ですが、それ以上のものではないかと思います。

比較的正確な結果が必要な場合は、他の操作から得られる方がはるかに良いとは思えません。ほとんどの場合、、、などを含む高精度を生成するものはすべてpow()時間exp()がlog()かかります。

ただし、インターネットで質問しても、自分のベンチマークに勝るものはありません。これがコードの重要な部分である場合は、いくつかの異なるバリアントを試して、何が得られるかを確認してください。

于 2013-03-09T09:57:50.520 に答える