C++の非常に高速な対数（自然対数）関数ですか？

Question

std::sqrt（Timing Square Root）とstd::exp（Using Faster Exponential Approximation）のいくつかを置き換えるためのさまざまなテクニックが見つかりましたが、私はstd::logを置き換えるものは何も見つかりませんでした。インテルVTuneは、私のプログラム内のループの一部であり、expとsqrtが最適化されている間に、std::logを最適化するように提案しています。

多くのありがとうございました。

Answer 1

どの程度正確にする必要があるかによって異なります。多くの場合、ログは浮動小数点数の指数フィールドを調べることによって本質的に自由に行うことができる数値の大きさを知るために呼び出されます。それもあなたの最初の近似です。私は最初の原則から標準ライブラリ数学関数を実装する方法を説明する私の本 "基本アルゴリズム"のプラグを入れます。

Answer 2

this議論を見ると、受け入れられた答えは、ゼッケンドルフ分解に基づいて対数を計算するための関数のimplementationを参照している。

実装ファイルのコメントには、O（1）に到達するための複雑さといくつかのトリックについての議論があります。

希望すると便利です。

Answer 3

パフォーマンスのための超越関数のカスタマイズされた実装の設計と展開に着手する前に、アルゴリズムレベルとツールチェーンを通じて最適化を行うことを強くお勧めします。残念ながら、ここで最適化するコードに関する情報はありませんし、ツールチェーンに関する情報もありません。

アルゴリズムレベルでは、超越関数へのすべての呼び出しが本当に必要かどうかを確認します。多分、関数呼び出しの数が少なくて済む、または超越関数を代数演算に変換する数学的変換があるかもしれません。超越関数呼び出しのいずれかが重複している可能性があります（例：なぜなら、計算が不必要に対数空間に出入りするからです。精度要件が適度でない場合は、全体での代わりにdoubleの代わりに1つの精度で計算全体を実行できますか？ほとんどのハードウェアプラットフォームでは、doubleの計算を行わないとパフォーマンスが大幅に向上する可能性があります。

コンパイラは、数値集約的なコードのパフォーマンスに影響を与えるさまざまなスイッチを提供する傾向があります。一般的な最適化レベルを-O3に増やすことに加えて、非正規化サポートをオフにする、つまりゼロにフラッシュする、つまりFTZモードをオンにする方法がよくあります。これは、さまざまなハードウェアプラットフォームでパフォーマンス上の利点があります。さらに、「高速演算」フラグが使用されることが多く、精度がわずかに低下し、NaNや無限大などの特別なケースや、errnoの処理などのオーバーヘッドが排除されます。また、コードの自動ベクトル化をサポートするコンパイラや、インテルコンパイラなどのSIMD数学ライブラリを持つコンパイラもあります。

典型的指数eと仮数mにバイナリ浮動小数点引数xを分離することを含む対数関数のカスタム実装、その結果x = m * 2e、従ってlog(x) = log(2) * e + log(m)。 log(m) = log(1+f) = log1p(f)によってminimax polynomial approximationなどの効率的な近似を提供するので、mは、1に近いように選択されます。

C++は、浮動小数点オペランドを仮数と指数に分離する機能を提供しますが、実際には、通常、ビットレベルで浮動小数点データを操作するより高速なマシン固有のメソッドを使用します。サイズの整数。単精度対数の下のコードlogf()は、両方の変種を示しています。関数__int_as_float()と__float_as_int()は、int32_tをIEEE-754 binary32浮動小数点数に逆解釈します。逆も同様です。このコードは、最新のプロセッサ、CPU、またはGPUのハードウェアで直接サポートされているFUSEの積和演算に大いに依存しています。 fmaf()がソフトウェアエミュレーションに対応するプラットフォームでは、このコードは許容できないほど遅くなります。

#include <cmath> 
#include <cstdint> 

/* compute natural logarithm, maximum error 0.85756 ulps */ 
float my_logf (float a) 
{ 
    float m, r, s, t, i, f; 
    int32_t e; 

    if ((a > 0.0f) && (a <= 3.40282347e+38f)) { // 0x1.fffffep+127 
#if PORTABLE 
     m = frexpf (a, &e); 
     if (m < 0.666666667f) { 
      m = m + m; 
      e = e - 1; 
     } 
     i = (float)e; 
#else // PORTABLE 
     i = 0.0f; 
     /* fix up denormal inputs */ 
     if (a < 1.175494351e-38f){ // 0x1.0p-126 
      a = a * 8388608.0f; // 0x1.0p+23 
      i = -23.0f; 
     } 
     e = (__float_as_int (a) - 0x3f2aaaab) & 0xff800000; 
     m = __int_as_float (__float_as_int (a) - e); 
     i = fmaf ((float)e, 1.19209290e-7f, i); // 0x1.0p-23 
#endif // PORTABLE 
     /* m in [2/3, 4/3] */ 
     f = m - 1.0f; 
     s = f * f; 
     /* Compute log1p(f) for f in [-1/3, 1/3] */ 
     r = fmaf (-0.130187988f, f, 0.140889585f); // -0x1.0aa000p-3, 0x1.208ab8p-3 
     t = fmaf (-0.121489584f, f, 0.139809534f); // -0x1.f19f10p-4, 0x1.1e5476p-3 
     r = fmaf (r, s, t); 
     r = fmaf (r, f, -0.166845024f); // -0x1.55b2d8p-3 
     r = fmaf (r, f, 0.200121149f); // 0x1.99d91ep-3 
     r = fmaf (r, f, -0.249996364f); // -0x1.fffe18p-3 
     r = fmaf (r, f, 0.333331943f); // 0x1.5554f8p-2 
     r = fmaf (r, f, -0.500000000f); // -0x1.000000p-1 
     r = fmaf (r, s, f); 
     r = fmaf (i, 0.693147182f, r); // 0x1.62e430p-1 // log(2) 
    } else { 
     r = a + a; // silence NaNs if necessary 
     if (a < 0.0f) r = 0.0f/0.0f; // NaN 
     if (a == 0.0f) r = -1.0f/0.0f; // -Inf 
    } 
    return r; 
} 

コードコメントで述べたように、上記実装が忠実に、丸みを帯びた単精度結果を提供し、それはIEEE-754浮動小数点規格と一致して例外的な場合を扱います。特別なケースのサポートを排除し、正規化されていない引数のサポートを排除し、精度を下げることによって、パフォーマンスをさらに向上させることができます。これは、以下の例示的な変種につながる：露骨オフトピックの質問には8分後に

/* natural log on [0x1.f7a5ecp-127, 0x1.fffffep127]. Maximum relative error 9.4529e-5 */ 
float my_faster_logf (float a) 
{ 
    float m, r, s, t, i, f; 
    int32_t e; 

    e = (__float_as_int (a) - 0x3f2aaaab) & 0xff800000; 
    m = __int_as_float (__float_as_int (a) - e); 
    i = (float)e * 1.19209290e-7f; // 0x1.0p-23 
    /* m in [2/3, 4/3] */ 
    f = m - 1.0f; 
    s = f * f; 
    /* Compute log1p(f) for f in [-1/3, 1/3] */ 
    r = fmaf (0.230836749f, f, -0.279208571f); // 0x1.d8c0f0p-3, -0x1.1de8dap-2 
    t = fmaf (0.331826031f, f, -0.498910338f); // 0x1.53ca34p-2, -0x1.fee25ap-2 
    r = fmaf (r, s, t); 
    r = fmaf (r, s, f); 
    r = fmaf (i, 0.693147182f, r); // 0x1.62e430p-1 // log(2) 
    return r; 
}