方程式を解く数値コードがあり、累乗する必要がありf(x) = 0
ます。いろいろ使って解いていくのですが、結局はニュートン法です。解はたまたま に等しく、それが問題の原因です。反復解が例えば に近づくと、計算に必要な時間が途方もなく大きくなり、簡単に 100 倍になり、コードが使用できなくなります。x
p
x = 1
1
x = 1 + 1e-13
std::pow(x, p)
このことを実行しているマシンは、CentOS の AMD64 (Opteron 6172) であり、コマンドは単純y = std::pow(x, p);
です。私のすべてのマシン、すべて x64 で同様の動作が見られます。hereに記載されているように、これは私だけの問題ではなく (つまり、他の誰かも腹を立てています)、x64 でのみ発生しx
、1.0
. 同様のことが起こりexp
ます。
この問題を解決することは、私にとって非常に重要です。この遅さを回避する方法があるかどうか誰かが知っていますか?
編集: ジョンは、これは非正規化によるものだと指摘しました。問題は、これをどのように修正するかです。コードは C++ で、g++
で使用するために でコンパイルされていますGNU Octave
。とCXXFLAGS
を含めるように設定しましたが、それは役に立たず、コードの実行が遅くなるようです。-mtune=native
-ffast-math
今のところ疑似的な解決策: この問題を気にかけているすべての人に、以下に提案されている解決策は個人的にはうまくいきませんでした。の通常の速度が本当に必要ですが、std::pow()
周りのもたつきはありませんx = 1
。個人的な解決策は、次のハックを使用することです。
inline double mpow(double x, double p) __attribute__ ((const));
inline double mpow(double x, double p)
{
double y(x - 1.0);
return (std::abs(y) > 1e-4) ? (std::pow(x, p)) : (1.0 + p * y * (1.0 + (p - 1.0) * y * (0.5 + (1.0 / 6.0) * (p - 2.0) * y)));
}
範囲は変更できますが、-40 < p < 40 の場合、誤差は約 1e-11 よりも小さく、これで十分です。オーバーヘッドは私が見つけたものから最小限であるため、問題は解決します。